Loading [MathJax]/jax/input/TeX/config.js
PolyFEM
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
SolveData.cpp
Go to the documentation of this file.
1#include "SolveData.hpp"
2
27
28#include <h5pp/h5pp.h>
29
30namespace polyfem::solver
31{
32 using namespace polyfem::time_integrator;
33
34 std::vector<std::shared_ptr<Form>> SolveData::init_forms(
35 // General
36 const Units &units,
37 const int dim,
38 const double t,
39 const Eigen::VectorXi &in_node_to_node,
40
41 // Elastic form
42 const int n_bases,
43 std::vector<basis::ElementBases> &bases,
44 const std::vector<basis::ElementBases> &geom_bases,
45 const assembler::Assembler &assembler,
46 assembler::AssemblyValsCache &ass_vals_cache,
47 const assembler::AssemblyValsCache &mass_ass_vals_cache,
48 const double jacobian_threshold,
49 const ElementInversionCheck check_inversion,
50
51 // Body form
52 const int n_pressure_bases,
53 const std::vector<int> &boundary_nodes,
54 const std::vector<mesh::LocalBoundary> &local_boundary,
55 const std::vector<mesh::LocalBoundary> &local_neumann_boundary,
56 const int n_boundary_samples,
57 const Eigen::MatrixXd &rhs,
58 const Eigen::MatrixXd &sol,
59 const assembler::Density &density,
60
61 // Pressure form
62 const std::vector<mesh::LocalBoundary> &local_pressure_boundary,
63 const std::unordered_map<int, std::vector<mesh::LocalBoundary>> &local_pressure_cavity,
64 const std::shared_ptr<assembler::PressureAssembler> pressure_assembler,
65
66 // Inertia form
67 const bool ignore_inertia,
68 const StiffnessMatrix &mass,
69 const std::shared_ptr<assembler::ViscousDamping> damping_assembler,
70
71 // Lagged regularization form
72 const double lagged_regularization_weight,
73 const int lagged_regularization_iterations,
74
75 // Augemented lagrangian form
76 const size_t obstacle_ndof,
77 const std::vector<std::string> &hard_constraint_files,
78 const std::vector<json> &soft_constraint_files,
79
80 // Contact form
81 const bool contact_enabled,
82 const ipc::CollisionMesh &collision_mesh,
83 const double dhat,
84 const double avg_mass,
85 const bool use_area_weighting,
86 const bool use_improved_max_operator,
87 const bool use_physical_barrier,
88 const json &barrier_stiffness,
89 const ipc::BroadPhaseMethod broad_phase,
90 const double ccd_tolerance,
91 const long ccd_max_iterations,
92 const bool enable_shape_derivatives,
93
94 // Normal Adhesion Form
95 const bool adhesion_enabled,
96 const double dhat_p,
97 const double dhat_a,
98 const double Y,
99
100 // Tangential Adhesion Form
101 const double tangential_adhesion_coefficient,
102 const double epsa,
103 const int tangential_adhesion_iterations,
104
105 // Homogenization
106 const assembler::MacroStrainValue &macro_strain_constraint,
107
108 // Periodic contact
109 const bool periodic_contact,
110 const Eigen::VectorXi &tiled_to_single,
111 const std::shared_ptr<utils::PeriodicBoundary> &periodic_bc,
112
113 // Friction form
114 const double friction_coefficient,
115 const double epsv,
116 const int friction_iterations,
117
118 // Rayleigh damping form
119 const json &rayleigh_damping)
120 {
121 const bool is_time_dependent = time_integrator != nullptr;
122 assert(!is_time_dependent || time_integrator != nullptr);
123 const double dt = is_time_dependent ? time_integrator->dt() : 0.0;
124 const int ndof = n_bases * dim;
125 // if (is_formulation_mixed) // mixed not supported
126 // ndof_ += n_pressure_bases; // pressure is a scalar
127 const bool is_volume = dim == 3;
128
129 std::vector<std::shared_ptr<Form>> forms;
130 al_form.clear();
131
132 elastic_form = std::make_shared<ElasticForm>(
133 n_bases, bases, geom_bases, assembler, ass_vals_cache,
134 t, dt, is_volume, jacobian_threshold, check_inversion);
135 forms.push_back(elastic_form);
136
137 if (rhs_assembler != nullptr)
138 {
139 body_form = std::make_shared<BodyForm>(
140 ndof, n_pressure_bases, boundary_nodes, local_boundary,
141 local_neumann_boundary, n_boundary_samples, rhs, *rhs_assembler,
142 density, /*is_formulation_mixed=*/false,
143 is_time_dependent);
144 body_form->update_quantities(t, sol);
145 forms.push_back(body_form);
146 }
147
148 if (pressure_assembler != nullptr)
149 {
150 pressure_form = std::make_shared<PressureForm>(
151 ndof,
152 local_pressure_boundary,
153 local_pressure_cavity,
154 boundary_nodes,
155 n_boundary_samples, *pressure_assembler,
156 is_time_dependent);
157 pressure_form->update_quantities(t, sol);
158 forms.push_back(pressure_form);
159 }
160
161 inertia_form = nullptr;
162 damping_form = nullptr;
163 if (is_time_dependent)
164 {
165 if (!ignore_inertia)
166 {
167 assert(time_integrator != nullptr);
168 inertia_form = std::make_shared<InertiaForm>(mass, *time_integrator);
169 forms.push_back(inertia_form);
170 }
171
172 if (damping_assembler != nullptr)
173 {
174 damping_form = std::make_shared<ElasticForm>(
175 n_bases, bases, geom_bases, *damping_assembler, ass_vals_cache, t, dt, is_volume);
176 forms.push_back(damping_form);
177 }
178 }
179 else
180 {
181 if (lagged_regularization_weight > 0)
182 {
183 forms.push_back(std::make_shared<LaggedRegForm>(lagged_regularization_iterations));
184 forms.back()->set_weight(lagged_regularization_weight);
185 }
186 }
187
188 if (rhs_assembler != nullptr)
189 {
190 // assembler::Mass mass_mat_assembler;
191 // mass_mat_assembler.set_size(dim);
192 StiffnessMatrix mass_tmp = mass;
193 // mass_mat_assembler.assemble(dim == 3, n_bases, bases, geom_bases, mass_ass_vals_cache, mass_tmp, true);
194 // assert(mass_tmp.rows() == mass.rows() && mass_tmp.cols() == mass.cols());
195
196 if (!boundary_nodes.empty())
197 al_form.push_back(std::make_shared<BCLagrangianForm>(
198 ndof, boundary_nodes, local_boundary, local_neumann_boundary,
199 n_boundary_samples, mass_tmp, *rhs_assembler, obstacle_ndof, is_time_dependent, t));
200 // forms.push_back(al_form.back());
201 }
202
203 if (periodic_bc != nullptr)
204 {
205 al_form.push_back(std::make_shared<PeriodicLagrangianForm>(ndof, periodic_bc));
206 }
207
208 for (const auto &path : hard_constraint_files)
209 {
210 logger().debug("Setting up hard constraints for {}", path);
211 h5pp::File file(path, h5pp::FileAccess::READONLY);
212 std::vector<int> local2global;
213 if (!file.findDatasets("local2global").empty())
214 local2global = file.readDataset<std::vector<int>>("local2global");
215
216 if (local2global.empty())
217 {
218 local2global.resize(in_node_to_node.size());
219
220 for (int i = 0; i < local2global.size(); ++i)
221 local2global[i] = in_node_to_node[i];
222 }
223 else
224 {
225 for (auto &v : local2global)
226 v = in_node_to_node[v];
227 }
228
229 Eigen::MatrixXd bin = file.readDataset<Eigen::MatrixXd>("b");
230
231 StiffnessMatrix A, A_proj;
232 Eigen::MatrixXd b, b_proj;
233
234 if (!file.findDatasets("A").empty())
235 {
236 Eigen::MatrixXd Ain = file.readDataset<Eigen::MatrixXd>("A");
237 utils::scatter_matrix(ndof, dim, Ain, bin, local2global, A, b);
238
239 if (!file.findDatasets("A_proj").empty())
240 {
241 Eigen::MatrixXd A_proj_in = file.readDataset<Eigen::MatrixXd>("A_proj");
242 if (file.findDatasets("b_proj").empty())
243 log_and_throw_error("Missing b_proj in hard constraint file");
244
245 Eigen::MatrixXd b_proj_in = file.readDataset<Eigen::MatrixXd>("b_proj");
246 utils::scatter_matrix_col(ndof, dim, A_proj_in, b_proj_in, local2global, A_proj, b_proj);
247 }
248 }
249 else
250 {
251 std::vector<double> values = file.readDataset<std::vector<double>>("A_triplets/values");
252 std::vector<int> rows = file.readDataset<std::vector<int>>("A_triplets/rows");
253 std::vector<int> cols = file.readDataset<std::vector<int>>("A_triplets/cols");
254 std::vector<long> shape = file.readDataset<std::vector<long>>("A_triplets/shape");
255 utils::scatter_matrix(ndof, dim, shape, rows, cols, values, bin, local2global, A, b);
256
257 if (!file.findGroups("A_proj_triplets").empty())
258 {
259 if (file.findDatasets("b_proj").empty())
260 log_and_throw_error("Missing b_proj in hard constraint file");
261 if (file.findDatasets("rows", "/A_proj_triplets").empty())
262 log_and_throw_error("Missing A_proj_triplets/rows in hard constraint file");
263 if (file.findDatasets("cols", "/A_proj_triplets").empty())
264 log_and_throw_error("Missing A_proj_triplets/cols in hard constraint file");
265 if (file.findDatasets("values", "/A_proj_triplets").empty())
266 log_and_throw_error("Missing A_proj_triplets/values in hard constraint file");
267
268 std::vector<double> values_proj = file.readDataset<std::vector<double>>("A_proj_triplets/values");
269 std::vector<int> rows_proj = file.readDataset<std::vector<int>>("A_proj_triplets/rows");
270 std::vector<int> cols_proj = file.readDataset<std::vector<int>>("A_proj_triplets/cols");
271 Eigen::MatrixXd b_projin = file.readDataset<Eigen::MatrixXd>("b_proj");
272 std::vector<long> shape_proj = file.readDataset<std::vector<long>>("A_proj_triplets/shape");
273
274 utils::scatter_matrix_col(ndof, dim, shape_proj, rows_proj, cols_proj, values_proj, b_projin, local2global, A_proj, b_proj);
275 }
276 }
277
278 al_form.push_back(std::make_shared<MatrixLagrangianForm>(A, b, A_proj, b_proj));
279 // forms.push_back(al_form.back());
280 }
281
282 for (const auto &j : soft_constraint_files)
283 {
284 const std::string &path = j["data"];
285 double weight = j["weight"];
286
287 logger().debug("Setting up soft constraints for {}", path);
288 h5pp::File file(path, h5pp::FileAccess::READONLY);
289 std::vector<int> local2global;
290 if (!file.findDatasets("local2global").empty())
291 local2global = file.readDataset<std::vector<int>>("local2global");
292
293 if (local2global.empty())
294 {
295 local2global.resize(in_node_to_node.size());
296
297 for (int i = 0; i < local2global.size(); ++i)
298 local2global[i] = in_node_to_node[i];
299 }
300 else
301 {
302 for (auto &v : local2global)
303 v = in_node_to_node[v];
304 }
305
306 Eigen::MatrixXd bin = file.readDataset<Eigen::MatrixXd>("b");
307
309 Eigen::MatrixXd b;
310
311 if (!file.findDatasets("A").empty())
312 {
313 Eigen::MatrixXd Ain = file.readDataset<Eigen::MatrixXd>("A");
314 utils::scatter_matrix(ndof, dim, Ain, bin, local2global, A, b);
315 }
316 else
317 {
318 std::vector<double> values = file.readDataset<std::vector<double>>("A_triplets/values");
319 std::vector<int> rows = file.readDataset<std::vector<int>>("A_triplets/rows");
320 std::vector<int> cols = file.readDataset<std::vector<int>>("A_triplets/cols");
321 std::vector<long> shape = file.readDataset<std::vector<long>>("A_triplets/shape");
322
323 utils::scatter_matrix(ndof, dim, shape, rows, cols, values, bin, local2global, A, b);
324 }
325
326 forms.push_back(std::make_shared<QuadraticPenaltyForm>(A, b, weight));
327 }
328
329 if (macro_strain_constraint.is_active())
330 {
331 // don't push these two into forms because they take a different input x
332 strain_al_lagr_form = std::make_shared<MacroStrainLagrangianForm>(macro_strain_constraint);
333 }
334
335 contact_form = nullptr;
336 periodic_contact_form = nullptr;
337 friction_form = nullptr;
338 if (contact_enabled)
339 {
340 const bool use_adaptive_barrier_stiffness = !barrier_stiffness.is_number();
341
342 if (periodic_contact)
343 {
344 periodic_contact_form = std::make_shared<PeriodicContactForm>(
345 collision_mesh, tiled_to_single, dhat, avg_mass, use_area_weighting, use_improved_max_operator, use_physical_barrier,
346 use_adaptive_barrier_stiffness, is_time_dependent, enable_shape_derivatives, broad_phase, ccd_tolerance,
347 ccd_max_iterations);
348
349 if (use_adaptive_barrier_stiffness)
350 {
351 periodic_contact_form->set_barrier_stiffness(1);
352 // logger().debug("Using adaptive barrier stiffness");
353 }
354 else
355 {
356 assert(barrier_stiffness.is_number());
357 assert(barrier_stiffness.get<double>() > 0);
358 periodic_contact_form->set_barrier_stiffness(barrier_stiffness);
359 // logger().debug("Using fixed barrier stiffness of {}", contact_form->barrier_stiffness());
360 }
361 }
362 else
363 {
364 contact_form = std::make_shared<ContactForm>(
365 collision_mesh, dhat, avg_mass, use_area_weighting, use_improved_max_operator, use_physical_barrier,
366 use_adaptive_barrier_stiffness, is_time_dependent, enable_shape_derivatives, broad_phase, ccd_tolerance * units.characteristic_length(),
367 ccd_max_iterations);
368
369 if (use_adaptive_barrier_stiffness)
370 {
371 contact_form->set_barrier_stiffness(1);
372 // logger().debug("Using adaptive barrier stiffness");
373 }
374 else
375 {
376 assert(barrier_stiffness.is_number());
377 assert(barrier_stiffness.get<double>() > 0);
378 contact_form->set_barrier_stiffness(barrier_stiffness);
379 // logger().debug("Using fixed barrier stiffness of {}", contact_form->barrier_stiffness());
380 }
381
382 if (contact_form)
383 forms.push_back(contact_form);
384
385 // ----------------------------------------------------------------
386 }
387
388 if (friction_coefficient != 0)
389 {
390 friction_form = std::make_shared<FrictionForm>(
391 collision_mesh, time_integrator, epsv, friction_coefficient,
392 broad_phase, *contact_form, friction_iterations);
393 friction_form->init_lagging(sol);
394 forms.push_back(friction_form);
395 }
396
397 if (adhesion_enabled)
398 {
399 normal_adhesion_form = std::make_shared<NormalAdhesionForm>(
400 collision_mesh, dhat_p, dhat_a, Y, is_time_dependent, enable_shape_derivatives,
401 broad_phase, ccd_tolerance * units.characteristic_length(), ccd_max_iterations
402 );
403 forms.push_back(normal_adhesion_form);
404
405 if (tangential_adhesion_coefficient != 0)
406 {
407 tangential_adhesion_form = std::make_shared<TangentialAdhesionForm>(
408 collision_mesh, time_integrator, epsa, tangential_adhesion_coefficient,
409 broad_phase, *normal_adhesion_form, tangential_adhesion_iterations
410 );
411 forms.push_back(tangential_adhesion_form);
412 }
413 }
414
415
416 }
417
418 const std::vector<json> rayleigh_damping_jsons = utils::json_as_array(rayleigh_damping);
419 if (is_time_dependent)
420 {
421 // Map from form name to form so RayleighDampingForm::create can get the correct form to damp
422 const std::unordered_map<std::string, std::shared_ptr<Form>> possible_forms_to_damp = {
423 {"elasticity", elastic_form},
424 {"contact", contact_form},
425 };
426
427 for (const json &params : rayleigh_damping_jsons)
428 {
429 forms.push_back(RayleighDampingForm::create(
430 params, possible_forms_to_damp,
432 }
433 }
434 else if (rayleigh_damping_jsons.size() > 0)
435 {
436 log_and_throw_adjoint_error("Rayleigh damping is only supported for time-dependent problems");
437 }
438
439 update_dt();
440
441 return forms;
442 }
443
444 void SolveData::update_barrier_stiffness(const Eigen::VectorXd &x)
445 {
446 if (contact_form == nullptr || !contact_form->use_adaptive_barrier_stiffness())
447 return;
448
449 Eigen::VectorXd grad_energy = Eigen::VectorXd::Zero(x.size());
450 const std::array<std::shared_ptr<Form>, 4> energy_forms{
452 for (const std::shared_ptr<Form> &form : energy_forms)
453 {
454 if (form == nullptr || !form->enabled())
455 continue;
456
457 Eigen::VectorXd grad_form;
458 form->first_derivative(x, grad_form);
459 grad_energy += grad_form;
460 }
461
462 contact_form->update_barrier_stiffness(x, grad_energy);
463 }
464
466 {
467 if (time_integrator == nullptr) // if is not time dependent
468 return;
469
470 const std::array<std::shared_ptr<Form>, 6> energy_forms{
472 for (const std::shared_ptr<Form> &form : energy_forms)
473 {
474 if (form == nullptr)
475 continue;
476 form->set_weight(time_integrator->acceleration_scaling());
477 }
478 }
479
480 std::vector<std::pair<std::string, std::shared_ptr<solver::Form>>> SolveData::named_forms() const
481 {
482 std::vector<std::pair<std::string, std::shared_ptr<solver::Form>>> res{
483 {"elastic", elastic_form},
484 {"inertia", inertia_form},
485 {"body", body_form},
486 {"contact", contact_form},
487 {"friction", friction_form},
488 {"damping", damping_form},
489 {"pressure", pressure_form},
490 {"strain_augmented_lagrangian_lagr", strain_al_lagr_form},
491 {"periodic_contact", periodic_contact_form},
492 };
493
494 for (const auto &form : al_form)
495 res.push_back({"augmented_lagrangian", form});
496
497 return res;
498 }
499} // namespace polyfem::solver
int x
double characteristic_length() const
Definition Units.hpp:22
Caches basis evaluation and geometric mapping at every element.
static std::shared_ptr< RayleighDampingForm > create(const json &params, const std::unordered_map< std::string, std::shared_ptr< Form > > &forms, const time_integrator::ImplicitTimeIntegrator &time_integrator)
std::shared_ptr< solver::FrictionForm > friction_form
std::shared_ptr< solver::InertiaForm > inertia_form
std::shared_ptr< solver::PeriodicContactForm > periodic_contact_form
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::PressureForm > pressure_form
std::shared_ptr< solver::BodyForm > body_form
std::shared_ptr< solver::MacroStrainLagrangianForm > strain_al_lagr_form
std::shared_ptr< solver::NormalAdhesionForm > normal_adhesion_form
std::shared_ptr< solver::ContactForm > contact_form
std::vector< std::pair< std::string, std::shared_ptr< solver::Form > > > named_forms() const
std::shared_ptr< solver::ElasticForm > damping_form
std::shared_ptr< solver::ElasticForm > elastic_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
std::shared_ptr< solver::TangentialAdhesionForm > tangential_adhesion_form
std::shared_ptr< assembler::PressureAssembler > pressure_assembler
std::shared_ptr< assembler::RhsAssembler > rhs_assembler
void scatter_matrix(const int n_dofs, const int dim, const Eigen::MatrixXd &A, const Eigen::MatrixXd &b, const std::vector< int > &local_to_global, StiffnessMatrix &Aout, Eigen::MatrixXd &bout)
std::vector< T > json_as_array(const json &j)
Return the value of a json object as an array.
Definition JSONUtils.hpp:38
void scatter_matrix_col(const int n_dofs, const int dim, const Eigen::MatrixXd &A, const Eigen::MatrixXd &b, const std::vector< int > &local_to_global, StiffnessMatrix &Aout, Eigen::MatrixXd &bout)
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:44
nlohmann::json json
Definition Common.hpp:9
void log_and_throw_adjoint_error(const std::string &msg)
Definition Logger.cpp:79
void log_and_throw_error(const std::string &msg)
Definition Logger.cpp:73
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix
Definition Types.hpp:22