Loading [MathJax]/extensions/tex2jax.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
24
25namespace polyfem::solver
26{
27 using namespace polyfem::time_integrator;
28
29 std::vector<std::shared_ptr<Form>> SolveData::init_forms(
30 // General
31 const Units &units,
32 const int dim,
33 const double t,
34
35 // Elastic form
36 const int n_bases,
37 std::vector<basis::ElementBases> &bases,
38 const std::vector<basis::ElementBases> &geom_bases,
39 const assembler::Assembler &assembler,
40 assembler::AssemblyValsCache &ass_vals_cache,
41 const assembler::AssemblyValsCache &mass_ass_vals_cache,
42 const double jacobian_threshold,
43 const ElementInversionCheck check_inversion,
44
45 // Body form
46 const int n_pressure_bases,
47 const std::vector<int> &boundary_nodes,
48 const std::vector<mesh::LocalBoundary> &local_boundary,
49 const std::vector<mesh::LocalBoundary> &local_neumann_boundary,
50 const int n_boundary_samples,
51 const Eigen::MatrixXd &rhs,
52 const Eigen::MatrixXd &sol,
53 const assembler::Density &density,
54
55 // Pressure form
56 const std::vector<mesh::LocalBoundary> &local_pressure_boundary,
57 const std::unordered_map<int, std::vector<mesh::LocalBoundary>> &local_pressure_cavity,
58 const std::shared_ptr<assembler::PressureAssembler> pressure_assembler,
59
60 // Inertia form
61 const bool ignore_inertia,
62 const StiffnessMatrix &mass,
63 const std::shared_ptr<assembler::ViscousDamping> damping_assembler,
64
65 // Lagged regularization form
66 const double lagged_regularization_weight,
67 const int lagged_regularization_iterations,
68
69 // Augemented lagrangian form
70 // const std::vector<int> &boundary_nodes,
71 // const std::vector<mesh::LocalBoundary> &local_boundary,
72 // const std::vector<mesh::LocalBoundary> &local_neumann_boundary,
73 // const int n_boundary_samples,
74 // const StiffnessMatrix &mass,
75 const size_t obstacle_ndof,
76
77 // Contact form
78 const bool contact_enabled,
79 const ipc::CollisionMesh &collision_mesh,
80 const double dhat,
81 const double avg_mass,
82 const bool use_area_weighting,
83 const bool use_improved_max_operator,
84 const bool use_physical_barrier,
85 const json &barrier_stiffness,
86 const ipc::BroadPhaseMethod broad_phase,
87 const double ccd_tolerance,
88 const long ccd_max_iterations,
89 const bool enable_shape_derivatives,
90
91 // Normal Adhesion Form
92 const bool adhesion_enabled,
93 const double dhat_p,
94 const double dhat_a,
95 const double Y,
96
97 // Tangential Adhesion Form
98 const double tangential_adhesion_coefficient,
99 const double epsa,
100 const int tangential_adhesion_iterations,
101
102 // Homogenization
103 const assembler::MacroStrainValue &macro_strain_constraint,
104
105 // Periodic contact
106 const bool periodic_contact,
107 const Eigen::VectorXi &tiled_to_single,
108 const std::shared_ptr<utils::PeriodicBoundary> &periodic_bc,
109
110 // Friction form
111 const double friction_coefficient,
112 const double epsv,
113 const int friction_iterations,
114
115 // Rayleigh damping form
116 const json &rayleigh_damping)
117 {
118 const bool is_time_dependent = time_integrator != nullptr;
119 assert(!is_time_dependent || time_integrator != nullptr);
120 const double dt = is_time_dependent ? time_integrator->dt() : 0.0;
121 const int ndof = n_bases * dim;
122 // if (is_formulation_mixed) // mixed not supported
123 // ndof_ += n_pressure_bases; // pressure is a scalar
124 const bool is_volume = dim == 3;
125
126 std::vector<std::shared_ptr<Form>> forms;
127 al_form.clear();
128
129 elastic_form = std::make_shared<ElasticForm>(
130 n_bases, bases, geom_bases, assembler, ass_vals_cache,
131 t, dt, is_volume, jacobian_threshold, check_inversion);
132 forms.push_back(elastic_form);
133
134 if (rhs_assembler != nullptr)
135 {
136 body_form = std::make_shared<BodyForm>(
137 ndof, n_pressure_bases, boundary_nodes, local_boundary,
138 local_neumann_boundary, n_boundary_samples, rhs, *rhs_assembler,
139 density, /*is_formulation_mixed=*/false,
140 is_time_dependent);
141 body_form->update_quantities(t, sol);
142 forms.push_back(body_form);
143 }
144
145 if (pressure_assembler != nullptr)
146 {
147 pressure_form = std::make_shared<PressureForm>(
148 ndof,
149 local_pressure_boundary,
150 local_pressure_cavity,
151 boundary_nodes,
152 n_boundary_samples, *pressure_assembler,
153 is_time_dependent);
154 pressure_form->update_quantities(t, sol);
155 forms.push_back(pressure_form);
156 }
157
158 inertia_form = nullptr;
159 damping_form = nullptr;
160 if (is_time_dependent)
161 {
162 if (!ignore_inertia)
163 {
164 assert(time_integrator != nullptr);
165 inertia_form = std::make_shared<InertiaForm>(mass, *time_integrator);
166 forms.push_back(inertia_form);
167 }
168
169 if (damping_assembler != nullptr)
170 {
171 damping_form = std::make_shared<ElasticForm>(
172 n_bases, bases, geom_bases, *damping_assembler, ass_vals_cache, t, dt, is_volume);
173 forms.push_back(damping_form);
174 }
175 }
176 else
177 {
178 if (lagged_regularization_weight > 0)
179 {
180 forms.push_back(std::make_shared<LaggedRegForm>(lagged_regularization_iterations));
181 forms.back()->set_weight(lagged_regularization_weight);
182 }
183 }
184
185 if (rhs_assembler != nullptr)
186 {
187 // assembler::Mass mass_mat_assembler;
188 // mass_mat_assembler.set_size(dim);
189 StiffnessMatrix mass_tmp = mass;
190 // mass_mat_assembler.assemble(dim == 3, n_bases, bases, geom_bases, mass_ass_vals_cache, mass_tmp, true);
191 // assert(mass_tmp.rows() == mass.rows() && mass_tmp.cols() == mass.cols());
192
193 al_form.push_back(std::make_shared<BCLagrangianForm>(
194 ndof, boundary_nodes, local_boundary, local_neumann_boundary,
195 n_boundary_samples, mass_tmp, *rhs_assembler, obstacle_ndof, is_time_dependent, t, periodic_bc));
196 forms.push_back(al_form.back());
197 }
198
199 if (macro_strain_constraint.is_active())
200 {
201 // don't push these two into forms because they take a different input x
202 strain_al_lagr_form = std::make_shared<MacroStrainLagrangianForm>(macro_strain_constraint);
203 }
204
205 contact_form = nullptr;
206 periodic_contact_form = nullptr;
207 friction_form = nullptr;
208 if (contact_enabled)
209 {
210 const bool use_adaptive_barrier_stiffness = !barrier_stiffness.is_number();
211
212 if (periodic_contact)
213 {
214 periodic_contact_form = std::make_shared<PeriodicContactForm>(
215 collision_mesh, tiled_to_single, dhat, avg_mass, use_area_weighting, use_improved_max_operator, use_physical_barrier,
216 use_adaptive_barrier_stiffness, is_time_dependent, enable_shape_derivatives, broad_phase, ccd_tolerance,
217 ccd_max_iterations);
218
219 if (use_adaptive_barrier_stiffness)
220 {
221 periodic_contact_form->set_barrier_stiffness(1);
222 // logger().debug("Using adaptive barrier stiffness");
223 }
224 else
225 {
226 assert(barrier_stiffness.is_number());
227 assert(barrier_stiffness.get<double>() > 0);
228 periodic_contact_form->set_barrier_stiffness(barrier_stiffness);
229 // logger().debug("Using fixed barrier stiffness of {}", contact_form->barrier_stiffness());
230 }
231 }
232 else
233 {
234 contact_form = std::make_shared<ContactForm>(
235 collision_mesh, dhat, avg_mass, use_area_weighting, use_improved_max_operator, use_physical_barrier,
236 use_adaptive_barrier_stiffness, is_time_dependent, enable_shape_derivatives, broad_phase, ccd_tolerance * units.characteristic_length(),
237 ccd_max_iterations);
238
239 if (use_adaptive_barrier_stiffness)
240 {
241 contact_form->set_barrier_stiffness(1);
242 // logger().debug("Using adaptive barrier stiffness");
243 }
244 else
245 {
246 assert(barrier_stiffness.is_number());
247 assert(barrier_stiffness.get<double>() > 0);
248 contact_form->set_barrier_stiffness(barrier_stiffness);
249 // logger().debug("Using fixed barrier stiffness of {}", contact_form->barrier_stiffness());
250 }
251
252 if (contact_form)
253 forms.push_back(contact_form);
254
255 // ----------------------------------------------------------------
256 }
257
258 if (friction_coefficient != 0)
259 {
260 friction_form = std::make_shared<FrictionForm>(
261 collision_mesh, time_integrator, epsv, friction_coefficient,
262 broad_phase, *contact_form, friction_iterations);
263 friction_form->init_lagging(sol);
264 forms.push_back(friction_form);
265 }
266
267 if (adhesion_enabled)
268 {
269 normal_adhesion_form = std::make_shared<NormalAdhesionForm>(
270 collision_mesh, dhat_p, dhat_a, Y, is_time_dependent, enable_shape_derivatives,
271 broad_phase, ccd_tolerance * units.characteristic_length(), ccd_max_iterations
272 );
273 forms.push_back(normal_adhesion_form);
274
275 if (tangential_adhesion_coefficient != 0)
276 {
277 tangential_adhesion_form = std::make_shared<TangentialAdhesionForm>(
278 collision_mesh, time_integrator, epsa, tangential_adhesion_coefficient,
279 broad_phase, *normal_adhesion_form, tangential_adhesion_iterations
280 );
281 forms.push_back(tangential_adhesion_form);
282 }
283 }
284
285
286 }
287
288 const std::vector<json> rayleigh_damping_jsons = utils::json_as_array(rayleigh_damping);
289 if (is_time_dependent)
290 {
291 // Map from form name to form so RayleighDampingForm::create can get the correct form to damp
292 const std::unordered_map<std::string, std::shared_ptr<Form>> possible_forms_to_damp = {
293 {"elasticity", elastic_form},
294 {"contact", contact_form},
295 };
296
297 for (const json &params : rayleigh_damping_jsons)
298 {
299 forms.push_back(RayleighDampingForm::create(
300 params, possible_forms_to_damp,
302 }
303 }
304 else if (rayleigh_damping_jsons.size() > 0)
305 {
306 log_and_throw_adjoint_error("Rayleigh damping is only supported for time-dependent problems");
307 }
308
309 update_dt();
310
311 return forms;
312 }
313
314 void SolveData::update_barrier_stiffness(const Eigen::VectorXd &x)
315 {
316 if (contact_form == nullptr || !contact_form->use_adaptive_barrier_stiffness())
317 return;
318
319 Eigen::VectorXd grad_energy = Eigen::VectorXd::Zero(x.size());
320 const std::array<std::shared_ptr<Form>, 4> energy_forms{
322 for (const std::shared_ptr<Form> &form : energy_forms)
323 {
324 if (form == nullptr || !form->enabled())
325 continue;
326
327 Eigen::VectorXd grad_form;
328 form->first_derivative(x, grad_form);
329 grad_energy += grad_form;
330 }
331
332 contact_form->update_barrier_stiffness(x, grad_energy);
333 }
334
336 {
337 if (time_integrator == nullptr) // if is not time dependent
338 return;
339
340 const std::array<std::shared_ptr<Form>, 6> energy_forms{
342 for (const std::shared_ptr<Form> &form : energy_forms)
343 {
344 if (form == nullptr)
345 continue;
346 form->set_weight(time_integrator->acceleration_scaling());
347 }
348 }
349
350 std::vector<std::pair<std::string, std::shared_ptr<solver::Form>>> SolveData::named_forms() const
351 {
352 std::vector<std::pair<std::string, std::shared_ptr<solver::Form>>> res{
353 {"elastic", elastic_form},
354 {"inertia", inertia_form},
355 {"body", body_form},
356 {"contact", contact_form},
357 {"friction", friction_form},
358 {"damping", damping_form},
359 {"pressure", pressure_form},
360 {"strain_augmented_lagrangian_lagr", strain_al_lagr_form},
361 {"periodic_contact", periodic_contact_form},
362 };
363
364 for (const auto &form : al_form)
365 res.push_back({"augmented_lagrangian", form});
366
367 return res;
368 }
369} // 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::vector< std::shared_ptr< Form > > init_forms(const Units &units, const int dim, const double t, 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 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:29
std::shared_ptr< solver::FrictionForm > friction_form
std::shared_ptr< solver::InertiaForm > inertia_form
std::shared_ptr< solver::PeriodicContactForm > periodic_contact_form
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
std::vector< T > json_as_array(const json &j)
Return the value of a json object as an array.
Definition JSONUtils.hpp:38
nlohmann::json json
Definition Common.hpp:9
void log_and_throw_adjoint_error(const std::string &msg)
Definition Logger.cpp:77
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix
Definition Types.hpp:22