Loading [MathJax]/extensions/tex2jax.js
PolyFEM
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
StateSolveNonlinear.cpp
Go to the documentation of this file.
1#include <polyfem/State.hpp>
2
5
14
25
26#include <ipc/ipc.hpp>
27
28namespace polyfem
29{
30 using namespace mesh;
31 using namespace solver;
32 using namespace time_integrator;
33 using namespace io;
34 using namespace utils;
35
36 std::shared_ptr<polysolve::nonlinear::Solver> State::make_nl_solver(bool for_al) const
37 {
38 return polysolve::nonlinear::Solver::create(for_al ? args["solver"]["augmented_lagrangian"]["nonlinear"] : args["solver"]["nonlinear"], args["solver"]["linear"], units.characteristic_length(), logger());
39 }
40
41 void State::solve_transient_tensor_nonlinear(const int time_steps, const double t0, const double dt, Eigen::MatrixXd &sol)
42 {
43 init_nonlinear_tensor_solve(sol, t0 + dt);
44
45 // Write the total energy to a CSV file
46 int save_i = 0;
47
48 std::unique_ptr<EnergyCSVWriter> energy_csv = nullptr;
49 std::unique_ptr<RuntimeStatsCSVWriter> stats_csv = nullptr;
50
51 if (args["output"]["stats"])
52 {
53 logger().debug("Saving nl stats to {} and {}", resolve_output_path("energy.csv"), resolve_output_path("stats.csv"));
54 energy_csv = std::make_unique<EnergyCSVWriter>(resolve_output_path("energy.csv"), solve_data);
55 stats_csv = std::make_unique<RuntimeStatsCSVWriter>(resolve_output_path("stats.csv"), *this, t0, dt);
56 }
57 const bool remesh_enabled = args["space"]["remesh"]["enabled"];
58 // const double save_dt = remesh_enabled ? (dt / 3) : dt;
59
60 // Save the initial solution
61 if (energy_csv)
62 energy_csv->write(save_i, sol);
63 save_timestep(t0, save_i, t0, dt, sol, Eigen::MatrixXd()); // no pressure
64 save_i++;
65
67 cache_transient_adjoint_quantities(0, sol, Eigen::MatrixXd::Zero(mesh->dimension(), mesh->dimension()));
68
69 for (int t = 1; t <= time_steps; ++t)
70 {
71 double forward_solve_time = 0, remeshing_time = 0, global_relaxation_time = 0;
72
73 {
74 POLYFEM_SCOPED_TIMER(forward_solve_time);
76 }
77
78 if (remesh_enabled)
79 {
80 if (energy_csv)
81 energy_csv->write(save_i, sol);
82 // save_timestep(t0 + dt * t, save_i, t0, save_dt, sol, Eigen::MatrixXd()); // no pressure
83 save_i++;
84
85 bool remesh_success;
86 {
87 POLYFEM_SCOPED_TIMER(remeshing_time);
88 remesh_success = this->remesh(t0 + dt * t, dt, sol);
89 }
90
91 // Save the solution after remeshing
92 if (energy_csv)
93 energy_csv->write(save_i, sol);
94 // save_timestep(t0 + dt * t, save_i, t0, save_dt, sol, Eigen::MatrixXd()); // no pressure
95 save_i++;
96
97 // Only do global relaxation if remeshing was successful
98 if (remesh_success)
99 {
100 POLYFEM_SCOPED_TIMER(global_relaxation_time);
101 solve_tensor_nonlinear(sol, t, false); // solve the scene again after remeshing
102 }
103 }
104
105 // Always save the solution for consistency
106 if (energy_csv)
107 energy_csv->write(save_i, sol);
108 save_timestep(t0 + dt * t, t, t0, dt, sol, Eigen::MatrixXd()); // no pressure
109 save_i++;
110
112 {
113 cache_transient_adjoint_quantities(t, sol, Eigen::MatrixXd::Zero(mesh->dimension(), mesh->dimension()));
114 }
115
116 {
117 POLYFEM_SCOPED_TIMER("Update quantities");
118
119 solve_data.time_integrator->update_quantities(sol);
120
121 solve_data.nl_problem->update_quantities(t0 + (t + 1) * dt, sol);
122
125 }
126
127 logger().info("{}/{} t={}", t, time_steps, t0 + dt * t);
128
129 const std::string rest_mesh_path = args["output"]["data"]["rest_mesh"].get<std::string>();
130 if (!rest_mesh_path.empty())
131 {
132 Eigen::MatrixXd V;
133 Eigen::MatrixXi F;
136 resolve_output_path(fmt::format(args["output"]["data"]["rest_mesh"], t)),
137 V, F, mesh->get_body_ids(), mesh->is_volume(), /*binary=*/true);
138 }
139
140 const std::string &state_path = resolve_output_path(fmt::format(args["output"]["data"]["state"], t));
141 if (!state_path.empty())
142 solve_data.time_integrator->save_state(state_path);
143
144 // save restart file
145 save_restart_json(t0, dt, t);
146 if (remesh_enabled && stats_csv)
147 stats_csv->write(t, forward_solve_time, remeshing_time, global_relaxation_time, sol);
148 }
149 }
150
151 void State::init_nonlinear_tensor_solve(Eigen::MatrixXd &sol, const double t, const bool init_time_integrator)
152 {
153 assert(sol.cols() == 1);
154 assert(!is_problem_linear()); // non-linear
155 assert(!problem->is_scalar()); // tensor
156 assert(mixed_assembler == nullptr);
157
159 {
160 if (initial_sol_update.size() == ndof())
161 sol = initial_sol_update;
162 else
163 initial_sol_update = sol;
164 }
165
166 // --------------------------------------------------------------------
167 // Check for initial intersections
168 if (is_contact_enabled())
169 {
170 POLYFEM_SCOPED_TIMER("Check for initial intersections");
171
172 const Eigen::MatrixXd displaced = collision_mesh.displace_vertices(
173 utils::unflatten(sol, mesh->dimension()));
174
175 if (ipc::has_intersections(collision_mesh, displaced, args["solver"]["contact"]["CCD"]["broad_phase"]))
176 {
178 resolve_output_path("intersection.obj"), displaced,
179 collision_mesh.edges(), collision_mesh.faces());
180 log_and_throw_error("Unable to solve, initial solution has intersections!");
181 }
182 }
183
184 // --------------------------------------------------------------------
185 // Initialize time integrator
186 if (problem->is_time_dependent())
187 {
188 if (init_time_integrator)
189 {
190 POLYFEM_SCOPED_TIMER("Initialize time integrator");
192
193 Eigen::MatrixXd solution, velocity, acceleration;
194 initial_solution(solution); // Reload this because we need all previous solutions
195 solution.col(0) = sol; // Make sure the current solution is the same as `sol`
196 assert(solution.rows() == sol.size());
197 initial_velocity(velocity);
198 assert(velocity.rows() == sol.size());
199 initial_acceleration(acceleration);
200 assert(acceleration.rows() == sol.size());
201
203 {
204 if (initial_vel_update.size() == ndof())
205 velocity = initial_vel_update;
206 else
207 initial_vel_update = velocity;
208 }
209
210 const double dt = args["time"]["dt"];
211 solve_data.time_integrator->init(solution, velocity, acceleration, dt);
212 }
213 assert(solve_data.time_integrator != nullptr);
214 }
215 else
216 {
217 solve_data.time_integrator = nullptr;
218 }
219
220 // --------------------------------------------------------------------
221 // Initialize forms
222
223 damping_assembler = std::make_shared<assembler::ViscousDamping>();
225
227
228 // for backward solve
229 damping_prev_assembler = std::make_shared<assembler::ViscousDampingPrev>();
231
232 const ElementInversionCheck check_inversion = args["solver"]["advanced"]["check_inversion"];
233 const std::vector<std::shared_ptr<Form>> forms = solve_data.init_forms(
234 // General
235 units,
236 mesh->dimension(), t, in_node_to_node,
237 // Elastic form
238 n_bases, bases, geom_bases(), *assembler, ass_vals_cache, mass_ass_vals_cache, args["solver"]["advanced"]["jacobian_threshold"], check_inversion,
239 // Body form
242 n_boundary_samples(), rhs, sol, mass_matrix_assembler->density(),
243 // Pressure form
245 // Inertia form
246 args.value("/time/quasistatic"_json_pointer, true), mass,
247 damping_assembler->is_valid() ? damping_assembler : nullptr,
248 // Lagged regularization form
249 args["solver"]["advanced"]["lagged_regularization_weight"],
250 args["solver"]["advanced"]["lagged_regularization_iterations"],
251 // Augmented lagrangian form
252 obstacle.ndof(), args["constraints"]["hard"], args["constraints"]["soft"],
253 // Contact form
254 args["contact"]["enabled"], args["contact"]["periodic"].get<bool>() ? periodic_collision_mesh : collision_mesh, args["contact"]["dhat"],
255 avg_mass, args["contact"]["use_convergent_formulation"] ? bool(args["contact"]["use_area_weighting"]) : false,
256 args["contact"]["use_convergent_formulation"] ? bool(args["contact"]["use_improved_max_operator"]) : false,
257 args["contact"]["use_convergent_formulation"] ? bool(args["contact"]["use_physical_barrier"]) : false,
258 args["solver"]["contact"]["barrier_stiffness"],
259 args["solver"]["contact"]["CCD"]["broad_phase"],
260 args["solver"]["contact"]["CCD"]["tolerance"],
261 args["solver"]["contact"]["CCD"]["max_iterations"],
263 // Normal Adhesion Form
264 args["contact"]["adhesion"]["adhesion_enabled"],
265 args["contact"]["adhesion"]["dhat_p"],
266 args["contact"]["adhesion"]["dhat_a"],
267 args["contact"]["adhesion"]["adhesion_strength"],
268 // Tangential Adhesion Form
269 args["contact"]["adhesion"]["tangential_adhesion_coefficient"],
270 args["contact"]["adhesion"]["epsa"],
271 args["solver"]["contact"]["tangential_adhesion_iterations"],
272 // Homogenization
274 // Periodic contact
275 args["contact"]["periodic"], periodic_collision_mesh_to_basis, periodic_bc,
276 // Friction form
277 args["contact"]["friction_coefficient"],
278 args["contact"]["epsv"],
279 args["solver"]["contact"]["friction_iterations"],
280 // Rayleigh damping form
281 args["solver"]["rayleigh_damping"]);
282
283 for (const auto &form : forms)
284 form->set_output_dir(output_dir);
285
286 if (solve_data.contact_form != nullptr)
287 solve_data.contact_form->save_ccd_debug_meshes = args["output"]["advanced"]["save_ccd_debug_meshes"];
288
289 // --------------------------------------------------------------------
290 // Initialize nonlinear problems
291
292 const int ndof = n_bases * mesh->dimension();
293 solve_data.nl_problem = std::make_shared<NLProblem>(
295 polysolve::linear::Solver::create(args["solver"]["linear"], logger()));
296 solve_data.nl_problem->init(sol);
297 solve_data.nl_problem->update_quantities(t, sol);
298 // --------------------------------------------------------------------
299
300 stats.solver_info = json::array();
301 }
302
303 void State::solve_tensor_nonlinear(Eigen::MatrixXd &sol, const int t, const bool init_lagging)
304 {
305 assert(solve_data.nl_problem != nullptr);
306 NLProblem &nl_problem = *(solve_data.nl_problem);
307
308 assert(sol.size() == rhs.size());
309
310 if (nl_problem.uses_lagging())
311 {
312 if (init_lagging)
313 {
314 POLYFEM_SCOPED_TIMER("Initializing lagging");
315 nl_problem.init_lagging(sol); // TODO: this should be u_prev projected
316 }
317 logger().info("Lagging iteration 1:");
318 }
319
320 // ---------------------------------------------------------------------
321
322 // Save the subsolve sequence for debugging
323 int subsolve_count = 0;
324 save_subsolve(subsolve_count, t, sol, Eigen::MatrixXd()); // no pressure
325
326 // ---------------------------------------------------------------------
327
328 std::shared_ptr<polysolve::nonlinear::Solver> nl_solver = make_nl_solver(true);
329
330 ALSolver al_solver(
332 args["solver"]["augmented_lagrangian"]["initial_weight"],
333 args["solver"]["augmented_lagrangian"]["scaling"],
334 args["solver"]["augmented_lagrangian"]["max_weight"],
335 args["solver"]["augmented_lagrangian"]["eta"],
336 [&](const Eigen::VectorXd &x) {
337 this->solve_data.update_barrier_stiffness(sol);
338 });
339
340 al_solver.post_subsolve = [&](const double al_weight) {
341 stats.solver_info.push_back(
342 {{"type", al_weight > 0 ? "al" : "rc"},
343 {"t", t}, // TODO: null if static?
344 {"info", nl_solver->info()}});
345 if (al_weight > 0)
346 stats.solver_info.back()["weight"] = al_weight;
347 save_subsolve(++subsolve_count, t, sol, Eigen::MatrixXd()); // no pressure
348 };
349
350 Eigen::MatrixXd prev_sol = sol;
351 al_solver.solve_al(nl_problem, sol,
352 args["solver"]["augmented_lagrangian"]["nonlinear"], args["solver"]["linear"], units.characteristic_length());
353
354 al_solver.solve_reduced(nl_problem, sol,
355 args["solver"]["nonlinear"], args["solver"]["linear"], units.characteristic_length());
356
357 if (args["space"]["advanced"]["count_flipped_els_continuous"])
358 {
359 const auto invalidList = count_invalid(mesh->dimension(), bases, geom_bases(), sol);
360 logger().debug("Flipped elements (cnt {}) : {}", invalidList.size(), invalidList);
361 }
362
363 // ---------------------------------------------------------------------
364
365 // TODO: Make this more general
366 const double lagging_tol = args["solver"]["contact"].value("friction_convergence_tol", 1e-2) * units.characteristic_length();
367
369 {
370 // Lagging loop (start at 1 because we already did an iteration above)
371 bool lagging_converged = !nl_problem.uses_lagging();
372 for (int lag_i = 1; !lagging_converged; lag_i++)
373 {
374 Eigen::VectorXd tmp_sol = nl_problem.full_to_reduced(sol);
375
376 // Update the lagging before checking for convergence
377 nl_problem.update_lagging(tmp_sol, lag_i);
378
379 // Check if lagging converged
380 Eigen::VectorXd grad;
381 nl_problem.gradient(tmp_sol, grad);
382 const double delta_x_norm = (prev_sol - sol).lpNorm<Eigen::Infinity>();
383 logger().debug("Lagging convergence grad_norm={:g} tol={:g} (||Δx||={:g})", grad.norm(), lagging_tol, delta_x_norm);
384 if (grad.norm() <= lagging_tol)
385 {
386 logger().info(
387 "Lagging converged in {:d} iteration(s) (grad_norm={:g} tol={:g})",
388 lag_i, grad.norm(), lagging_tol);
389 lagging_converged = true;
390 break;
391 }
392
393 if (delta_x_norm <= 1e-12)
394 {
395 logger().warn(
396 "Lagging produced tiny update between iterations {:d} and {:d} (grad_norm={:g} grad_tol={:g} ||Δx||={:g} Δx_tol={:g}); stopping early",
397 lag_i - 1, lag_i, grad.norm(), lagging_tol, delta_x_norm, 1e-6);
398 lagging_converged = false;
399 break;
400 }
401
402 // Check for convergence first before checking if we can continue
403 if (lag_i >= nl_problem.max_lagging_iterations())
404 {
405 logger().warn(
406 "Lagging failed to converge with {:d} iteration(s) (grad_norm={:g} tol={:g})",
407 lag_i, grad.norm(), lagging_tol);
408 lagging_converged = false;
409 break;
410 }
411
412 // Solve the problem with the updated lagging
413 logger().info("Lagging iteration {:d}:", lag_i + 1);
414 nl_problem.init(sol);
416 nl_problem.normalize_forms();
417 nl_solver->minimize(nl_problem, tmp_sol);
418 nl_problem.finish();
419 prev_sol = sol;
420 sol = nl_problem.reduced_to_full(tmp_sol);
421
422 // Save the subsolve sequence for debugging and info
423 stats.solver_info.push_back(
424 {{"type", "rc"},
425 {"t", t}, // TODO: null if static?
426 {"lag_i", lag_i},
427 {"info", nl_solver->info()}});
428 save_subsolve(++subsolve_count, t, sol, Eigen::MatrixXd()); // no pressure
429 }
430 }
431 }
432} // namespace polyfem
int V
int x
#define POLYFEM_SCOPED_TIMER(...)
Definition Timer.hpp:10
Eigen::MatrixXd initial_vel_update
Definition State.hpp:704
assembler::MacroStrainValue macro_strain_constraint
Definition State.hpp:712
std::shared_ptr< utils::PeriodicBoundary > periodic_bc
periodic BC and periodic mesh utils
Definition State.hpp:386
int n_bases
number of bases
Definition State.hpp:178
void cache_transient_adjoint_quantities(const int current_step, const Eigen::MatrixXd &sol, const Eigen::MatrixXd &disp_grad)
int n_pressure_bases
number of pressure bases
Definition State.hpp:180
assembler::AssemblyValsCache ass_vals_cache
used to store assembly values for small problems
Definition State.hpp:196
int n_boundary_samples() const
quadrature used for projecting boundary conditions
Definition State.hpp:263
const std::vector< basis::ElementBases > & geom_bases() const
Get a constant reference to the geometry mapping bases.
Definition State.hpp:223
std::vector< mesh::LocalBoundary > local_pressure_boundary
mapping from elements to nodes for pressure boundary conditions
Definition State.hpp:442
std::shared_ptr< assembler::ViscousDamping > damping_assembler
Definition State.hpp:164
Eigen::VectorXi in_node_to_node
Inpute nodes (including high-order) to polyfem nodes, only for isoparametric.
Definition State.hpp:455
mesh::Obstacle obstacle
Obstacles used in collisions.
Definition State.hpp:473
std::shared_ptr< assembler::Assembler > assembler
assemblers
Definition State.hpp:155
ipc::CollisionMesh periodic_collision_mesh
IPC collision mesh under periodic BC.
Definition State.hpp:523
void save_subsolve(const int i, const int t, const Eigen::MatrixXd &sol, const Eigen::MatrixXd &pressure)
saves a subsolve when save_solve_sequence_debug is true
std::string resolve_output_path(const std::string &path) const
Resolve output path relative to output_dir if the path is not absolute.
ipc::CollisionMesh collision_mesh
IPC collision mesh.
Definition State.hpp:520
std::string output_dir
Directory for output files.
Definition State.hpp:586
solver::CacheLevel optimization_enabled
Definition State.hpp:654
void set_materials(std::vector< std::shared_ptr< assembler::Assembler > > &assemblers) const
set the material and the problem dimension
void save_timestep(const double time, const int t, const double t0, const double dt, const Eigen::MatrixXd &sol, const Eigen::MatrixXd &pressure)
saves a timestep
void solve_transient_tensor_nonlinear(const int time_steps, const double t0, const double dt, Eigen::MatrixXd &sol)
solves transient tensor nonlinear problem
void initial_solution(Eigen::MatrixXd &solution) const
Load or compute the initial solution.
std::shared_ptr< polysolve::nonlinear::Solver > make_nl_solver(bool for_al) const
factory to create the nl solver depending on input
std::shared_ptr< assembler::Mass > mass_matrix_assembler
Definition State.hpp:157
bool remesh(const double time, const double dt, Eigen::MatrixXd &sol)
Remesh the FE space and update solution(s).
std::unique_ptr< mesh::Mesh > mesh
current mesh, it can be a Mesh2D or Mesh3D
Definition State.hpp:471
StiffnessMatrix mass
Mass matrix, it is computed only for time dependent problems.
Definition State.hpp:202
std::shared_ptr< assembler::Problem > problem
current problem, it contains rhs and bc
Definition State.hpp:168
json args
main input arguments containing all defaults
Definition State.hpp:101
void save_restart_json(const double t0, const double dt, const int t) const
Save a JSON sim file for restarting the simulation at time t.
void initial_velocity(Eigen::MatrixXd &velocity) const
Load or compute the initial velocity.
io::OutStatsData stats
Other statistics.
Definition State.hpp:592
Eigen::VectorXi periodic_collision_mesh_to_basis
index mapping from periodic 2x2 collision mesh to FE periodic mesh
Definition State.hpp:525
void initial_acceleration(Eigen::MatrixXd &acceleration) const
Load or compute the initial acceleration.
std::vector< basis::ElementBases > bases
FE bases, the size is #elements.
Definition State.hpp:171
std::vector< mesh::LocalBoundary > local_boundary
mapping from elements to nodes for dirichlet boundary conditions
Definition State.hpp:438
double avg_mass
average system mass, used for contact with IPC
Definition State.hpp:204
int ndof() const
Definition State.hpp:660
assembler::AssemblyValsCache mass_ass_vals_cache
Definition State.hpp:197
void solve_tensor_nonlinear(Eigen::MatrixXd &sol, const int t=0, const bool init_lagging=true)
solves nonlinear problems
std::vector< int > boundary_nodes
list of boundary nodes
Definition State.hpp:432
solver::SolveData solve_data
timedependent stuff cached
Definition State.hpp:321
void init_nonlinear_tensor_solve(Eigen::MatrixXd &sol, const double t=1.0, const bool init_time_integrator=true)
initialize the nonlinear solver
std::shared_ptr< assembler::ViscousDampingPrev > damping_prev_assembler
Definition State.hpp:165
bool is_contact_enabled() const
does the simulation have contact
Definition State.hpp:556
std::vector< mesh::LocalBoundary > local_neumann_boundary
mapping from elements to nodes for neumann boundary conditions
Definition State.hpp:440
std::shared_ptr< assembler::PressureAssembler > build_pressure_assembler() const
Definition State.hpp:256
std::shared_ptr< assembler::PressureAssembler > elasticity_pressure_assembler
Definition State.hpp:162
void build_mesh_matrices(Eigen::MatrixXd &V, Eigen::MatrixXi &F)
Build the mesh matrices (vertices and elements) from the mesh using the bases node ordering.
bool is_problem_linear() const
Returns whether the system is linear. Collisions and pressure add nonlinearity to the problem.
Definition State.hpp:407
std::shared_ptr< assembler::MixedAssembler > mixed_assembler
Definition State.hpp:159
Eigen::MatrixXd initial_sol_update
Definition State.hpp:704
Eigen::MatrixXd rhs
System right-hand side.
Definition State.hpp:207
std::unordered_map< int, std::vector< mesh::LocalBoundary > > local_pressure_cavity
mapping from elements to nodes for pressure boundary conditions
Definition State.hpp:444
double characteristic_length() const
Definition Units.hpp:22
static void write(const std::string &path, const mesh::Mesh &mesh, const bool binary)
saves the mesh
Definition MshWriter.cpp:7
static bool write(const std::string &path, const Eigen::MatrixXd &v, const Eigen::MatrixXi &e, const Eigen::MatrixXi &f)
Definition OBJWriter.cpp:18
json solver_info
information of the solver, eg num iteration, time, errors, etc the informations varies depending on t...
Definition OutData.hpp:375
void solve_reduced(NLProblem &nl_problem, Eigen::MatrixXd &sol, std::shared_ptr< polysolve::nonlinear::Solver > nl_solver)
Definition ALSolver.hpp:41
std::function< void(const double)> post_subsolve
Definition ALSolver.hpp:53
void solve_al(NLProblem &nl_problem, Eigen::MatrixXd &sol, std::shared_ptr< polysolve::nonlinear::Solver > nl_solver)
Definition ALSolver.hpp:29
virtual void init(const TVector &x0) override
double normalize_forms() override
TVector full_to_reduced(const TVector &full) const
virtual void gradient(const TVector &x, TVector &gradv) override
void init_lagging(const TVector &x) override
TVector reduced_to_full(const TVector &reduced) const
void update_lagging(const TVector &x, const int iter_num) override
std::vector< std::shared_ptr< Form > > init_forms(const Units &units, const int dim, const double t, const Eigen::VectorXi &in_node_to_node, const int n_bases, std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &geom_bases, const assembler::Assembler &assembler, assembler::AssemblyValsCache &ass_vals_cache, const assembler::AssemblyValsCache &mass_ass_vals_cache, const double jacobian_threshold, const solver::ElementInversionCheck check_inversion, const int n_pressure_bases, const std::vector< int > &boundary_nodes, const std::vector< mesh::LocalBoundary > &local_boundary, const std::vector< mesh::LocalBoundary > &local_neumann_boundary, const int n_boundary_samples, const Eigen::MatrixXd &rhs, const Eigen::MatrixXd &sol, const assembler::Density &density, const std::vector< mesh::LocalBoundary > &local_pressure_boundary, const std::unordered_map< int, std::vector< mesh::LocalBoundary > > &local_pressure_cavity, const std::shared_ptr< assembler::PressureAssembler > pressure_assembler, const bool ignore_inertia, const StiffnessMatrix &mass, const std::shared_ptr< assembler::ViscousDamping > damping_assembler, const double lagged_regularization_weight, const int lagged_regularization_iterations, const size_t obstacle_ndof, const std::vector< std::string > &hard_constraint_files, const std::vector< json > &soft_constraint_files, const bool contact_enabled, const ipc::CollisionMesh &collision_mesh, const double dhat, const double avg_mass, const bool use_area_weighting, const bool use_improved_max_operator, const bool use_physical_barrier, const json &barrier_stiffness, const ipc::BroadPhaseMethod broad_phase, const double ccd_tolerance, const long ccd_max_iterations, const bool enable_shape_derivatives, const bool adhesion_enabled, const double dhat_p, const double dhat_a, const double Y, const double tangential_adhesion_coefficient, const double epsa, const int tangential_adhesion_iterations, const assembler::MacroStrainValue &macro_strain_constraint, const bool periodic_contact, const Eigen::VectorXi &tiled_to_single, const std::shared_ptr< utils::PeriodicBoundary > &periodic_bc, const double friction_coefficient, const double epsv, const int friction_iterations, const json &rayleigh_damping)
Initialize the forms and return a vector of pointers to them.
Definition SolveData.cpp:34
void update_dt()
updates the dt inside the different forms
std::shared_ptr< solver::NLProblem > nl_problem
std::shared_ptr< solver::ContactForm > contact_form
std::vector< std::shared_ptr< solver::AugmentedLagrangianForm > > al_form
std::shared_ptr< time_integrator::ImplicitTimeIntegrator > time_integrator
void update_barrier_stiffness(const Eigen::VectorXd &x)
update the barrier stiffness for the forms
static std::shared_ptr< ImplicitTimeIntegrator > construct_time_integrator(const json &params)
Factory method for constructing implicit time integrators from the name of the integrator.
Eigen::MatrixXd unflatten(const Eigen::VectorXd &x, int dim)
Unflatten rowwises, so every dim elements in x become a row.
std::vector< int > count_invalid(const int dim, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const Eigen::VectorXd &u)
Definition Jacobian.cpp:173
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:44
void log_and_throw_error(const std::string &msg)
Definition Logger.cpp:73