38 mesh =
Mesh::create(local_mesh.rest_positions(), local_mesh.elements());
39 assert(mesh->n_vertices() == local_mesh.num_vertices());
41 std::vector<int> boundary_ids(mesh->n_boundary_elements(), -1);
42 if constexpr (std::is_same_v<M, WildTriRemesher>)
44 const auto local_boundary_ids = std::get<Remesher::EdgeMap<int>>(local_mesh.boundary_ids());
45 for (
int i = 0; i < mesh->n_edges(); i++)
47 const size_t e0 = mesh->edge_vertex(i, 0);
48 const size_t e1 = mesh->edge_vertex(i, 1);
49 boundary_ids[i] = local_boundary_ids.at({{e0, e1}});
54 const auto local_boundary_ids = std::get<Remesher::FaceMap<int>>(local_mesh.boundary_ids());
55 for (
int i = 0; i < mesh->n_faces(); i++)
57 std::array<size_t, 3> f = {{
58 (size_t)mesh->face_vertex(i, 0),
59 (size_t)mesh->face_vertex(i, 1),
60 (size_t)mesh->face_vertex(i, 2),
62 boundary_ids[i] = local_boundary_ids.at(f);
66 mesh->set_boundary_ids(boundary_ids);
67 mesh->set_body_ids(local_mesh.body_ids());
75 Eigen::VectorXi vertex_to_basis;
77 *mesh, state.
formulation(), bases, local_boundary, vertex_to_basis);
79 assert(m_n_bases == local_mesh.num_local_vertices());
80 m_n_bases = local_mesh.num_vertices();
81 assert(vertex_to_basis.size() == m_n_bases);
83 const int start_i = local_mesh.num_local_vertices();
84 if (start_i < m_n_bases)
87 std::iota(vertex_to_basis.begin() + start_i, vertex_to_basis.end(), start_i);
90 assert(std::all_of(vertex_to_basis.begin(), vertex_to_basis.end(), [](
const int basis_id) {
95 local_mesh.reorder_vertices(vertex_to_basis);
96 problem->update_nodes(vertex_to_basis);
97 mesh->update_nodes(vertex_to_basis);
105 assert(mesh !=
nullptr);
108 std::vector<int> pressure_boundary_nodes;
111 std::vector<basis::ElementBases>(), local_boundary,
112 boundary_nodes, local_neumann_boundary, pressure_boundary_nodes,
113 dirichlet_nodes, neumann_nodes);
115 auto find_node_position = [&](
const int n_id) {
116 for (
const auto &bs : bases)
117 for (
const auto &b : bs.bases)
118 for (
const auto &lg : b.global())
119 if (lg.index == n_id)
125 dirichlet_nodes_position.resize(dirichlet_nodes.size());
126 for (
int n = 0; n < dirichlet_nodes.size(); ++n)
127 dirichlet_nodes_position[n] = find_node_position(dirichlet_nodes[n]);
129 neumann_nodes_position.resize(neumann_nodes.size());
130 for (
int n = 0; n < neumann_nodes.size(); ++n)
131 neumann_nodes_position[n] = find_node_position(neumann_nodes[n]);
134 for (
const int vi : local_mesh.fixed_vertices())
136 for (
int d = 0; d < dim(); ++d)
138 boundary_nodes.push_back(vi * dim() + d);
142 std::sort(boundary_nodes.begin(), boundary_nodes.end());
143 auto it = std::unique(boundary_nodes.begin(), boundary_nodes.end());
144 boundary_nodes.erase(it, boundary_nodes.end());
169 mass_assembly_vals_cache.init(
170 is_volume(), bases, bases,
true);
171 assert(mass_matrix_assembler !=
nullptr);
172 mass_matrix_assembler->assemble(
173 is_volume(), n_bases(), bases, bases,
174 mass_assembly_vals_cache,
179 const int local_ndof = dim() * local_mesh.num_local_vertices();
180 for (
int i = local_ndof; i < ndof(); ++i)
181 mass.coeffRef(i, i) = state.
avg_mass;
187 const double current_time,
188 const bool contact_enabled)
191 const Eigen::MatrixXd target_x = this->sol();
194 assembly_vals_cache.init(is_volume(), bases, bases,
false);
201 collision_mesh = ipc::CollisionMesh::build_from_full_mesh(
202 local_mesh.rest_positions(), local_mesh.boundary_edges(),
203 local_mesh.boundary_faces());
206 std::vector<bool> is_vertex_fixed(local_mesh.num_vertices(),
false);
207 for (
const int vi : local_mesh.fixed_vertices())
208 is_vertex_fixed[vi] =
true;
209 collision_mesh.can_collide = [is_vertex_fixed,
this](
size_t vi,
size_t vj) {
210 return !is_vertex_fixed[this->collision_mesh.to_full_vertex_id(vi)]
211 || !is_vertex_fixed[this->collision_mesh.to_full_vertex_id(vj)];
216 if (state.
problem->is_time_dependent())
219 solve_data.time_integrator =
221 state.
args[
"time"][
"integrator"]);
222 Eigen::MatrixXd x_prevs, v_prevs, a_prevs;
224 local_mesh.projection_quantities(), dim(), x_prevs, v_prevs, a_prevs);
225 solve_data.time_integrator->init(
226 x_prevs, v_prevs, a_prevs, state.
args[
"time"][
"dt"]);
233 json rhs_solver_params = state.
args[
"solver"][
"linear"];
234 if (!rhs_solver_params.contains(
"Pardiso"))
235 rhs_solver_params[
"Pardiso"] = {};
236 rhs_solver_params[
"Pardiso"][
"mtype"] = -2;
238 const int size = state.
problem->is_scalar() ? 1 : dim();
239 solve_data.rhs_assembler = std::make_shared<assembler::RhsAssembler>(
240 *assembler, *mesh,
Obstacle(), dirichlet_nodes, neumann_nodes,
241 dirichlet_nodes_position, neumann_nodes_position, n_bases(),
242 dim(), bases, bases, mass_assembly_vals_cache,
243 *state.
problem, state.
args[
"space"][
"advanced"][
"bc_method"],
246 solve_data.rhs_assembler->assemble(mass_matrix_assembler->density(), rhs);
250 std::vector<std::shared_ptr<solver::Form>> forms;
253 forms = solve_data.init_forms(
255 state.
units, dim(), current_time,
257 n_bases(), bases, bases, *assembler,
258 assembly_vals_cache, assembly_vals_cache,
260 0, boundary_nodes, local_boundary,
262 target_x, mass_matrix_assembler->density(),
264 state.
args.value(
"/time/quasistatic"_json_pointer,
true), mass,
267 state.
args[
"solver"][
"advanced"][
"lagged_regularization_weight"],
268 state.
args[
"solver"][
"advanced"][
"lagged_regularization_iterations"],
272 contact_enabled, collision_mesh, state.
args[
"contact"][
"dhat"],
273 state.
avg_mass, state.
args[
"contact"][
"use_convergent_formulation"],
275 state.
args[
"solver"][
"contact"][
"CCD"][
"broad_phase"],
276 state.
args[
"solver"][
"contact"][
"CCD"][
"tolerance"],
277 state.
args[
"solver"][
"contact"][
"CCD"][
"max_iterations"],
280 state.
args[
"contact"][
"friction_coefficient"],
281 state.
args[
"contact"][
"epsv"],
282 state.
args[
"solver"][
"contact"][
"friction_iterations"],
284 state.
args[
"solver"][
"rayleigh_damping"]);
287 forms.erase(std::remove(forms.begin(), forms.end(), solve_data.al_lagr_form), forms.end());
288 forms.erase(std::remove(forms.begin(), forms.end(), solve_data.al_pen_form), forms.end());
289 solve_data.al_lagr_form =
nullptr;
290 solve_data.al_pen_form =
nullptr;
293 solve_data.nl_problem = std::make_shared<polyfem::solver::StaticBoundaryNLProblem>(
294 ndof(), boundary_nodes, target_x, forms);
296 assert(solve_data.time_integrator !=
nullptr);
297 solve_data.nl_problem->update_quantities(current_time, solve_data.time_integrator->x_prev());
298 solve_data.nl_problem->init(target_x);
299 solve_data.nl_problem->init_lagging(solve_data.time_integrator->x_prev());
300 solve_data.nl_problem->update_lagging(target_x, 0);