39 const Eigen::VectorXi &in_node_to_node,
43 std::vector<basis::ElementBases> &bases,
44 const std::vector<basis::ElementBases> &geom_bases,
48 const double jacobian_threshold,
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,
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,
67 const bool ignore_inertia,
69 const std::shared_ptr<assembler::ViscousDamping> damping_assembler,
72 const double lagged_regularization_weight,
73 const int lagged_regularization_iterations,
76 const size_t obstacle_ndof,
77 const std::vector<std::string> &hard_constraint_files,
78 const std::vector<json> &soft_constraint_files,
81 const bool contact_enabled,
82 const ipc::CollisionMesh &collision_mesh,
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,
95 const bool adhesion_enabled,
101 const double tangential_adhesion_coefficient,
103 const int tangential_adhesion_iterations,
109 const bool periodic_contact,
110 const Eigen::VectorXi &tiled_to_single,
111 const std::shared_ptr<utils::PeriodicBoundary> &periodic_bc,
114 const double friction_coefficient,
116 const int friction_iterations,
119 const json &rayleigh_damping)
124 const int ndof = n_bases * dim;
127 const bool is_volume = dim == 3;
129 std::vector<std::shared_ptr<Form>> forms;
133 n_bases, bases, geom_bases, assembler, ass_vals_cache,
134 t, dt, is_volume, jacobian_threshold, check_inversion);
140 ndof, n_pressure_bases, boundary_nodes, local_boundary,
141 local_neumann_boundary, n_boundary_samples, rhs, *
rhs_assembler,
152 local_pressure_boundary,
153 local_pressure_cavity,
163 if (is_time_dependent)
172 if (damping_assembler !=
nullptr)
175 n_bases, bases, geom_bases, *damping_assembler, ass_vals_cache, t, dt, is_volume);
181 if (lagged_regularization_weight > 0)
183 forms.push_back(std::make_shared<LaggedRegForm>(lagged_regularization_iterations));
184 forms.back()->set_weight(lagged_regularization_weight);
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));
203 if (periodic_bc !=
nullptr)
205 al_form.push_back(std::make_shared<PeriodicLagrangianForm>(ndof, periodic_bc));
208 for (
const auto &path : hard_constraint_files)
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");
216 if (local2global.empty())
218 local2global.resize(in_node_to_node.size());
220 for (
int i = 0; i < local2global.size(); ++i)
221 local2global[i] = in_node_to_node[i];
225 for (
auto &v : local2global)
226 v = in_node_to_node[v];
229 Eigen::MatrixXd bin = file.readDataset<Eigen::MatrixXd>(
"b");
232 Eigen::MatrixXd b, b_proj;
234 if (!file.findDatasets(
"A").empty())
236 Eigen::MatrixXd Ain = file.readDataset<Eigen::MatrixXd>(
"A");
239 if (!file.findDatasets(
"A_proj").empty())
241 Eigen::MatrixXd A_proj_in = file.readDataset<Eigen::MatrixXd>(
"A_proj");
242 if (file.findDatasets(
"b_proj").empty())
245 Eigen::MatrixXd b_proj_in = file.readDataset<Eigen::MatrixXd>(
"b_proj");
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");
257 if (!file.findGroups(
"A_proj_triplets").empty())
259 if (file.findDatasets(
"b_proj").empty())
261 if (file.findDatasets(
"rows",
"/A_proj_triplets").empty())
263 if (file.findDatasets(
"cols",
"/A_proj_triplets").empty())
265 if (file.findDatasets(
"values",
"/A_proj_triplets").empty())
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");
274 utils::scatter_matrix_col(ndof, dim, shape_proj, rows_proj, cols_proj, values_proj, b_projin, local2global, A_proj, b_proj);
278 al_form.push_back(std::make_shared<MatrixLagrangianForm>(A, b, A_proj, b_proj));
282 for (
const auto &j : soft_constraint_files)
284 const std::string &path = j[
"data"];
285 double weight = j[
"weight"];
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");
293 if (local2global.empty())
295 local2global.resize(in_node_to_node.size());
297 for (
int i = 0; i < local2global.size(); ++i)
298 local2global[i] = in_node_to_node[i];
302 for (
auto &v : local2global)
303 v = in_node_to_node[v];
306 Eigen::MatrixXd bin = file.readDataset<Eigen::MatrixXd>(
"b");
311 if (!file.findDatasets(
"A").empty())
313 Eigen::MatrixXd Ain = file.readDataset<Eigen::MatrixXd>(
"A");
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");
326 forms.push_back(std::make_shared<QuadraticPenaltyForm>(A, b, weight));
332 strain_al_lagr_form = std::make_shared<MacroStrainLagrangianForm>(macro_strain_constraint);
340 const bool use_adaptive_barrier_stiffness = !barrier_stiffness.is_number();
342 if (periodic_contact)
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,
349 if (use_adaptive_barrier_stiffness)
356 assert(barrier_stiffness.is_number());
357 assert(barrier_stiffness.get<
double>() > 0);
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(),
369 if (use_adaptive_barrier_stiffness)
376 assert(barrier_stiffness.is_number());
377 assert(barrier_stiffness.get<
double>() > 0);
388 if (friction_coefficient != 0)
397 if (adhesion_enabled)
400 collision_mesh, dhat_p, dhat_a, Y, is_time_dependent, enable_shape_derivatives,
405 if (tangential_adhesion_coefficient != 0)
408 collision_mesh,
time_integrator, epsa, tangential_adhesion_coefficient,
419 if (is_time_dependent)
422 const std::unordered_map<std::string, std::shared_ptr<Form>> possible_forms_to_damp = {
427 for (
const json ¶ms : rayleigh_damping_jsons)
430 params, possible_forms_to_damp,
434 else if (rayleigh_damping_jsons.size() > 0)
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 ¯o_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.