PolyFEM
Loading...
Searching...
No Matches
StateDiff.cpp
Go to the documentation of this file.
1#include <polyfem/State.hpp>
3
6
8#include <polysolve/linear/FEMSolver.hpp>
12
15
26
27#include <ipc/ipc.hpp>
28#include <ipc/barrier/barrier.hpp>
29#include <ipc/utils/local_to_global.hpp>
30#include <ipc/potentials/friction_potential.hpp>
31
32#include <Eigen/Dense>
33#include <unsupported/Eigen/SparseExtra>
34#include <deque>
35#include <map>
36#include <algorithm>
37
38#include <fstream>
39
40using namespace polyfem::basis;
41
42namespace polyfem
43{
44 namespace
45 {
46 void replace_rows_by_identity(StiffnessMatrix &reduced_mat, const StiffnessMatrix &mat, const std::vector<int> &rows)
47 {
48 reduced_mat.resize(mat.rows(), mat.cols());
49
50 std::vector<bool> mask(mat.rows(), false);
51 for (int i : rows)
52 mask[i] = true;
53
54 std::vector<Eigen::Triplet<double>> coeffs;
55 for (int k = 0; k < mat.outerSize(); ++k)
56 {
57 for (StiffnessMatrix::InnerIterator it(mat, k); it; ++it)
58 {
59 if (mask[it.row()])
60 {
61 if (it.row() == it.col())
62 coeffs.emplace_back(it.row(), it.col(), 1.0);
63 }
64 else
65 coeffs.emplace_back(it.row(), it.col(), it.value());
66 }
67 }
68 reduced_mat.setFromTriplets(coeffs.begin(), coeffs.end());
69 }
70 } // namespace
71
72 void State::get_vertices(Eigen::MatrixXd &vertices) const
73 {
74 vertices.setZero(mesh->n_vertices(), mesh->dimension());
75
76 for (int v = 0; v < mesh->n_vertices(); v++)
77 vertices.row(v) = mesh->point(v);
78 }
79
80 void State::get_elements(Eigen::MatrixXi &elements) const
81 {
82 assert(mesh->is_simplicial());
83
84 auto node_to_primitive_map = node_to_primitive();
85
86 const auto &gbases = geom_bases();
87 int dim = mesh->dimension();
88 elements.setZero(gbases.size(), dim + 1);
89 for (int e = 0; e < gbases.size(); e++)
90 {
91 int i = 0;
92 for (const auto &gbs : gbases[e].bases)
93 elements(e, i++) = node_to_primitive_map[gbs.global()[0].index];
94 }
95 }
96
97 void State::set_mesh_vertex(int v_id, const Eigen::VectorXd &vertex)
98 {
99 assert(vertex.size() == mesh->dimension());
100 mesh->set_point(v_id, vertex);
101 }
102
103 void State::cache_transient_adjoint_quantities(const int current_step, const Eigen::MatrixXd &sol, const Eigen::MatrixXd &disp_grad)
104 {
105 StiffnessMatrix gradu_h(sol.size(), sol.size());
106 if (current_step == 0)
107 diff_cached.init(mesh->dimension(), ndof(), problem->is_time_dependent() ? args["time"]["time_steps"].get<int>() : 0);
108
109 ipc::NormalCollisions cur_collision_set;
110 ipc::SmoothCollisions cur_smooth_collision_set;
111 ipc::TangentialCollisions cur_friction_set;
112 ipc::NormalCollisions cur_normal_adhesion_set;
113 ipc::TangentialCollisions cur_tangential_adhesion_set;
114
116 {
117 if (!problem->is_time_dependent() || current_step > 0)
118 compute_force_jacobian(sol, disp_grad, gradu_h);
119
121 {
122 if (const auto barrier_contact = dynamic_cast<const solver::BarrierContactForm*>(solve_data.contact_form.get()))
123 cur_collision_set = barrier_contact->collision_set();
124 else if (const auto smooth_contact = dynamic_cast<const solver::SmoothContactForm*>(solve_data.contact_form.get()))
125 cur_smooth_collision_set = smooth_contact->collision_set();
126 }
127 cur_friction_set = solve_data.friction_form ? solve_data.friction_form->friction_collision_set() : ipc::TangentialCollisions();
128 cur_normal_adhesion_set = solve_data.normal_adhesion_form ? solve_data.normal_adhesion_form->collision_set() : ipc::NormalCollisions();
129 cur_tangential_adhesion_set = solve_data.tangential_adhesion_form ? solve_data.tangential_adhesion_form->tangential_collision_set() : ipc::TangentialCollisions();
130 }
131
132 if (problem->is_time_dependent())
133 {
134 if (args["time"]["quasistatic"].get<bool>())
135 {
136 diff_cached.cache_quantities_quasistatic(current_step, sol, gradu_h, cur_collision_set, cur_smooth_collision_set, cur_normal_adhesion_set, disp_grad);
137 }
138 else
139 {
140 Eigen::MatrixXd vel, acc;
141 if (current_step == 0)
142 {
143 if (dynamic_cast<time_integrator::BDF *>(solve_data.time_integrator.get()))
144 {
145 const auto bdf_integrator = dynamic_cast<time_integrator::BDF *>(solve_data.time_integrator.get());
146 vel = bdf_integrator->weighted_sum_v_prevs();
147 }
148 else if (dynamic_cast<time_integrator::ImplicitEuler *>(solve_data.time_integrator.get()))
149 {
150 const auto euler_integrator = dynamic_cast<time_integrator::ImplicitEuler *>(solve_data.time_integrator.get());
151 vel = euler_integrator->v_prev();
152 }
153 else
154 log_and_throw_error("Differentiable code doesn't support this time integrator!");
155
156 acc.setZero(ndof(), 1);
157 }
158 else
159 {
160 vel = solve_data.time_integrator->compute_velocity(sol);
161 acc = solve_data.time_integrator->compute_acceleration(vel);
162 }
163
164 diff_cached.cache_quantities_transient(current_step, solve_data.time_integrator->steps(), sol, vel, acc, gradu_h, cur_collision_set, cur_smooth_collision_set, cur_friction_set);
165 }
166 }
167 else
168 {
169 diff_cached.cache_quantities_static(sol, gradu_h, cur_collision_set, cur_smooth_collision_set, cur_friction_set, cur_normal_adhesion_set, cur_tangential_adhesion_set, disp_grad);
170 }
171 }
172
173 void State::compute_force_jacobian(const Eigen::MatrixXd &sol, const Eigen::MatrixXd &disp_grad, StiffnessMatrix &hessian)
174 {
175 if (problem->is_time_dependent())
176 {
177 if (assembler->is_linear() && !is_contact_enabled())
178 log_and_throw_adjoint_error("Differentiable transient linear solve is not supported!");
179
180 StiffnessMatrix tmp_hess;
181 solve_data.nl_problem->set_project_to_psd(false);
182 solve_data.nl_problem->FullNLProblem::solution_changed(sol);
183 solve_data.nl_problem->FullNLProblem::hessian(sol, tmp_hess);
184 hessian.setZero();
185 replace_rows_by_identity(hessian, tmp_hess, boundary_nodes);
186 }
187 else // static formulation
188 {
189 if (assembler->is_linear() && !is_contact_enabled() && !is_homogenization())
190 {
191 hessian.setZero();
192 StiffnessMatrix stiffness;
193 build_stiffness_mat(stiffness);
194 replace_rows_by_identity(hessian, stiffness, boundary_nodes);
195 }
196 else
197 {
198 solve_data.nl_problem->set_project_to_psd(false);
199 if (is_homogenization())
200 {
201 Eigen::VectorXd reduced;
202 std::shared_ptr<solver::NLHomoProblem> homo_problem = std::dynamic_pointer_cast<solver::NLHomoProblem>(solve_data.nl_problem);
203 reduced = homo_problem->full_to_reduced(sol, disp_grad);
204 solve_data.nl_problem->solution_changed(reduced);
205 solve_data.nl_problem->hessian(reduced, hessian);
206 }
207 else
208 {
209 StiffnessMatrix tmp_hess;
210 solve_data.nl_problem->FullNLProblem::solution_changed(sol);
211 solve_data.nl_problem->FullNLProblem::hessian(sol, tmp_hess);
212 hessian.setZero();
213 replace_rows_by_identity(hessian, tmp_hess, boundary_nodes);
214 }
215 }
216 }
217 }
218
219 void State::compute_force_jacobian_prev(const int force_step, const int sol_step, StiffnessMatrix &hessian_prev) const
220 {
221 assert(force_step > 0);
222 assert(force_step > sol_step);
223 if (assembler->is_linear() && !is_contact_enabled())
224 {
225 hessian_prev = StiffnessMatrix(ndof(), ndof());
226 }
227 else
228 {
229 const Eigen::MatrixXd u = diff_cached.u(force_step);
230 const Eigen::MatrixXd u_prev = diff_cached.u(sol_step);
231 const double beta = time_integrator::BDF::betas(diff_cached.bdf_order(force_step) - 1);
232 const double dt = solve_data.time_integrator->dt();
233
234 hessian_prev = StiffnessMatrix(u.size(), u.size());
235 if (problem->is_time_dependent())
236 {
238 {
239 if (sol_step == force_step - 1)
240 {
241 Eigen::MatrixXd surface_solution_prev = collision_mesh.vertices(utils::unflatten(u_prev, mesh->dimension()));
242 Eigen::MatrixXd surface_solution = collision_mesh.vertices(utils::unflatten(u, mesh->dimension()));
243
244 // TODO: use the time integration to compute the velocity
245 const Eigen::MatrixXd surface_velocities = (surface_solution - surface_solution_prev) / dt;
246 const double dv_dut = -1 / dt;
247
248 if (const auto barrier_contact = dynamic_cast<const solver::BarrierContactForm*>(solve_data.contact_form.get()))
249 {
250 hessian_prev =
251 solve_data.friction_form->friction_potential().force_jacobian(
254 collision_mesh.rest_positions(),
255 /*lagged_displacements=*/surface_solution_prev,
256 surface_velocities,
257 barrier_contact->barrier_potential(),
258 barrier_contact->barrier_stiffness(),
259 ipc::FrictionPotential::DiffWRT::LAGGED_DISPLACEMENTS)
260 + solve_data.friction_form->friction_potential().force_jacobian(
263 collision_mesh.rest_positions(),
264 /*lagged_displacements=*/surface_solution_prev,
265 surface_velocities,
266 barrier_contact->barrier_potential(),
267 barrier_contact->barrier_stiffness(),
268 ipc::FrictionPotential::DiffWRT::VELOCITIES)
269 * dv_dut;
270 }
271
272 hessian_prev *= -1;
273
274 // {
275 // Eigen::MatrixXd X = collision_mesh.rest_positions();
276 // Eigen::VectorXd x = utils::flatten(surface_solution_prev);
277 // const double barrier_stiffness = solve_data.contact_form->barrier_stiffness();
278 // const double dhat = solve_data.contact_form->dhat();
279 // const double mu = solve_data.friction_form->mu();
280 // const double epsv = solve_data.friction_form->epsv();
281
282 // Eigen::MatrixXd fgrad;
283 // fd::finite_jacobian(
284 // x, [&](const Eigen::VectorXd &y) -> Eigen::VectorXd
285 // {
286 // Eigen::MatrixXd fd_Ut = utils::unflatten(y, surface_solution_prev.cols());
287
288 // ipc::TangentialCollisions fd_friction_constraints;
289 // ipc::NormalCollisions fd_constraints;
290 // fd_constraints.set_use_convergent_formulation(solve_data.contact_form->use_convergent_formulation());
291 // fd_constraints.set_enable_shape_derivatives(true);
292 // fd_constraints.build(collision_mesh, X + fd_Ut, dhat);
293
294 // fd_friction_constraints.build(
295 // collision_mesh, X + fd_Ut, fd_constraints, dhat, barrier_stiffness,
296 // mu);
297
298 // return fd_friction_constraints.compute_potential_gradient(collision_mesh, (surface_solution - fd_Ut) / dt, epsv);
299
300 // }, fgrad, fd::AccuracyOrder::SECOND, 1e-8);
301
302 // logger().trace("force Ut derivative error {} {}", (fgrad - hessian_prev).norm(), hessian_prev.norm());
303 // }
304
305 hessian_prev = collision_mesh.to_full_dof(hessian_prev); // / (beta * dt) / (beta * dt);
306 }
307 else
308 {
309 // const double alpha = time_integrator::BDF::alphas(std::min(diff_cached.bdf_order(force_step), force_step) - 1)[force_step - sol_step - 1];
310 // Eigen::MatrixXd velocity = collision_mesh.map_displacements(utils::unflatten(diff_cached.v(force_step), collision_mesh.dim()));
311 // hessian_prev = diff_cached.friction_collision_set(force_step).compute_potential_hessian( //
312 // collision_mesh, velocity, solve_data.friction_form->epsv(), false) * (-alpha / beta / dt);
313
314 // hessian_prev = collision_mesh.to_full_dof(hessian_prev);
315 }
316 }
317
319 {
320
321 if (sol_step == force_step - 1)
322 {
323 StiffnessMatrix adhesion_hessian_prev(u.size(), u.size());
324
325 Eigen::MatrixXd surface_solution_prev = collision_mesh.vertices(utils::unflatten(u_prev, mesh->dimension()));
326 Eigen::MatrixXd surface_solution = collision_mesh.vertices(utils::unflatten(u, mesh->dimension()));
327
328 // TODO: use the time integration to compute the velocity
329 const Eigen::MatrixXd surface_velocities = (surface_solution - surface_solution_prev) / dt;
330 const double dv_dut = -1 / dt;
331
332 adhesion_hessian_prev =
333 solve_data.tangential_adhesion_form->tangential_adhesion_potential().force_jacobian(
336 collision_mesh.rest_positions(),
337 /*lagged_displacements=*/surface_solution_prev,
338 surface_velocities,
339 solve_data.normal_adhesion_form->normal_adhesion_potential(),
340 1,
341 ipc::TangentialPotential::DiffWRT::LAGGED_DISPLACEMENTS)
342 + solve_data.tangential_adhesion_form->tangential_adhesion_potential().force_jacobian(
345 collision_mesh.rest_positions(),
346 /*lagged_displacements=*/surface_solution_prev,
347 surface_velocities,
348 solve_data.normal_adhesion_form->normal_adhesion_potential(),
349 1,
350 ipc::TangentialPotential::DiffWRT::VELOCITIES)
351 * dv_dut;
352
353 adhesion_hessian_prev *= -1;
354
355 adhesion_hessian_prev = collision_mesh.to_full_dof(adhesion_hessian_prev); // / (beta * dt) / (beta * dt);
356
357 hessian_prev += adhesion_hessian_prev;
358 }
359 }
360
361 if (damping_assembler->is_valid() && sol_step == force_step - 1) // velocity in damping uses BDF1
362 {
363 utils::SparseMatrixCache mat_cache;
364 StiffnessMatrix damping_hessian_prev(u.size(), u.size());
365 damping_prev_assembler->assemble_hessian(mesh->is_volume(), n_bases, false, bases, geom_bases(), ass_vals_cache, force_step * args["time"]["dt"].get<double>() + args["time"]["t0"].get<double>(), dt, u, u_prev, mat_cache, damping_hessian_prev);
366
367 hessian_prev += damping_hessian_prev;
368 }
369
370 if (sol_step == force_step - 1)
371 {
372 StiffnessMatrix body_force_hessian(u.size(), u.size());
373 solve_data.body_form->hessian_wrt_u_prev(u_prev, force_step * dt, body_force_hessian);
374 hessian_prev += body_force_hessian;
375 }
376 }
377 }
378 }
379
380 void State::solve_adjoint_cached(const Eigen::MatrixXd &rhs)
381 {
383 }
384
385 Eigen::MatrixXd State::solve_adjoint(const Eigen::MatrixXd &rhs) const
386 {
387 if (problem->is_time_dependent())
389 else
391 }
392
393 Eigen::MatrixXd State::solve_static_adjoint(const Eigen::MatrixXd &adjoint_rhs) const
394 {
395 Eigen::MatrixXd b = adjoint_rhs;
396
397 Eigen::MatrixXd adjoint;
399 {
400 b(boundary_nodes, Eigen::all).setZero();
401
403 const int full_size = A.rows();
404 const int problem_dim = problem->is_scalar() ? 1 : mesh->dimension();
405 int precond_num = problem_dim * n_bases;
406
407 b.conservativeResizeLike(Eigen::MatrixXd::Zero(A.rows(), b.cols()));
408
409 std::vector<int> boundary_nodes_tmp;
410 if (has_periodic_bc())
411 {
412 boundary_nodes_tmp = periodic_bc->full_to_periodic(boundary_nodes);
413 precond_num = periodic_bc->full_to_periodic(A);
414 b = periodic_bc->full_to_periodic(b, true);
415 }
416 else
417 boundary_nodes_tmp = boundary_nodes;
418
419 adjoint.setZero(ndof(), adjoint_rhs.cols());
420 for (int i = 0; i < b.cols(); i++)
421 {
422 Eigen::VectorXd x, tmp;
423 tmp = b.col(i);
424 dirichlet_solve_prefactorized(*lin_solver_cached, A, tmp, boundary_nodes_tmp, x);
425
426 if (has_periodic_bc())
427 adjoint.col(i) = periodic_bc->periodic_to_full(full_size, x);
428 else
429 adjoint.col(i) = x;
430 }
431 }
432 else
433 {
434 auto solver = polysolve::linear::Solver::create(args["solver"]["adjoint_linear"], adjoint_logger());
435
436 StiffnessMatrix A = diff_cached.gradu_h(0); // This should be transposed, but A is symmetric in hyper-elastic and diffusion problems
437
438 /*
439 For non-periodic problems, the adjoint solution p's size is the full size in NLProblem
440 For periodic problems, the adjoint solution p's size is the reduced size in NLProblem
441 */
442 if (!is_homogenization())
443 {
444 adjoint.setZero(ndof(), adjoint_rhs.cols());
445 for (int i = 0; i < b.cols(); i++)
446 {
447 Eigen::VectorXd tmp = b.col(i);
448 tmp(boundary_nodes).setZero();
449
450 Eigen::VectorXd x;
451 x.setZero(tmp.size());
452 dirichlet_solve(*solver, A, tmp, boundary_nodes, x, A.rows(), "", false, false, false);
453
454 adjoint.col(i) = x;
455 adjoint(boundary_nodes, i) = -b(boundary_nodes, i);
456 }
457 }
458 else
459 {
460 solver->analyze_pattern(A, A.rows());
461 solver->factorize(A);
462
463 adjoint.setZero(adjoint_rhs.rows(), adjoint_rhs.cols());
464 for (int i = 0; i < b.cols(); i++)
465 {
466 Eigen::MatrixXd tmp = b.col(i);
467
468 Eigen::VectorXd x;
469 x.setZero(tmp.size());
470 solver->solve(tmp, x);
471 x.conservativeResize(adjoint.rows());
472
473 adjoint.col(i) = x;
474 }
475 }
476 }
477
478 return adjoint;
479 }
480
481 Eigen::MatrixXd State::solve_transient_adjoint(const Eigen::MatrixXd &adjoint_rhs) const
482 {
483 const double dt = args["time"]["dt"];
484 const int time_steps = args["time"]["time_steps"];
485
486 int bdf_order = 1;
487 if (args["time"]["integrator"].is_string())
488 bdf_order = 1;
489 else if (args["time"]["integrator"]["type"] == "ImplicitEuler")
490 bdf_order = 1;
491 else if (args["time"]["integrator"]["type"] == "BDF")
492 bdf_order = args["time"]["integrator"]["steps"].get<int>();
493 else
494 log_and_throw_adjoint_error("Integrator type not supported for differentiability.");
495
496 assert(adjoint_rhs.cols() == time_steps + 1);
497
498 const int cols_per_adjoint = time_steps + 1;
499 Eigen::MatrixXd adjoints;
500 adjoints.setZero(ndof(), cols_per_adjoint * 2);
501
502 // set dirichlet rows of mass to identity
503 StiffnessMatrix reduced_mass;
504 replace_rows_by_identity(reduced_mass, mass, boundary_nodes);
505
506 Eigen::MatrixXd sum_alpha_p, sum_alpha_nu;
507 for (int i = time_steps; i >= 0; --i)
508 {
509 {
510 sum_alpha_p.setZero(ndof(), 1);
511 sum_alpha_nu.setZero(ndof(), 1);
512
513 const int num = std::min(bdf_order, time_steps - i);
514
515 Eigen::VectorXd bdf_coeffs = Eigen::VectorXd::Zero(num);
516 for (int j = 0; j < bdf_order && i + j < time_steps; ++j)
517 bdf_coeffs(j) = -time_integrator::BDF::alphas(std::min(bdf_order - 1, i + j))[j];
518
519 sum_alpha_p = adjoints.middleCols(i + 1, num) * bdf_coeffs;
520 sum_alpha_nu = adjoints.middleCols(cols_per_adjoint + i + 1, num) * bdf_coeffs;
521 }
522
523 Eigen::VectorXd rhs_ = -reduced_mass.transpose() * sum_alpha_nu - adjoint_rhs.col(i);
524 for (int j = 1; j <= bdf_order; j++)
525 {
526 if (i + j > time_steps)
527 break;
528
529 StiffnessMatrix gradu_h_prev;
530 compute_force_jacobian_prev(i + j, i, gradu_h_prev);
531 Eigen::VectorXd tmp = adjoints.col(i + j) * (time_integrator::BDF::betas(diff_cached.bdf_order(i + j) - 1) * dt);
532 tmp(boundary_nodes).setZero();
533 rhs_ += -gradu_h_prev.transpose() * tmp;
534 }
535
536 if (i > 0)
537 {
538 double beta_dt = time_integrator::BDF::betas(diff_cached.bdf_order(i) - 1) * dt;
539
540 rhs_ += (1. / beta_dt) * (diff_cached.gradu_h(i) - reduced_mass).transpose() * sum_alpha_p;
541
542 {
543 StiffnessMatrix A = diff_cached.gradu_h(i).transpose();
544 Eigen::VectorXd b_ = rhs_;
545 b_(boundary_nodes).setZero();
546
547 auto solver = polysolve::linear::Solver::create(args["solver"]["adjoint_linear"], adjoint_logger());
548
549 Eigen::VectorXd x;
550 dirichlet_solve(*solver, A, b_, boundary_nodes, x, A.rows(), "", false, false, false);
551 adjoints.col(i + cols_per_adjoint) = x;
552 }
553
554 // TODO: generalize to BDFn
555 Eigen::VectorXd tmp = rhs_(boundary_nodes);
556 if (i + 1 < cols_per_adjoint)
557 tmp += (-2. / beta_dt) * adjoints(boundary_nodes, i + 1);
558 if (i + 2 < cols_per_adjoint)
559 tmp += (1. / beta_dt) * adjoints(boundary_nodes, i + 2);
560
561 tmp -= (diff_cached.gradu_h(i).transpose() * adjoints.col(i + cols_per_adjoint))(boundary_nodes);
562 adjoints(boundary_nodes, i + cols_per_adjoint) = tmp;
563 adjoints.col(i) = beta_dt * adjoints.col(i + cols_per_adjoint) - sum_alpha_p;
564 }
565 else
566 {
567 adjoints.col(i) = -reduced_mass.transpose() * sum_alpha_p;
568 adjoints.col(i + cols_per_adjoint) = rhs_; // adjoint_nu[0] actually stores adjoint_mu[0]
569 }
570 }
571 return adjoints;
572 }
573
574 void State::compute_surface_node_ids(const int surface_selection, std::vector<int> &node_ids) const
575 {
576 node_ids = {};
577
578 const auto &gbases = geom_bases();
579 for (const auto &lb : total_local_boundary)
580 {
581 const int e = lb.element_id();
582 for (int i = 0; i < lb.size(); ++i)
583 {
584 const int primitive_global_id = lb.global_primitive_id(i);
585 const int boundary_id = mesh->get_boundary_id(primitive_global_id);
586 const auto nodes = gbases[e].local_nodes_for_primitive(primitive_global_id, *mesh);
587
588 if (boundary_id == surface_selection)
589 {
590 for (long n = 0; n < nodes.size(); ++n)
591 {
592 const int g_id = gbases[e].bases[nodes(n)].global()[0].index;
593
594 if (std::count(node_ids.begin(), node_ids.end(), g_id) == 0)
595 node_ids.push_back(g_id);
596 }
597 }
598 }
599 }
600 }
601
602 void State::compute_total_surface_node_ids(std::vector<int> &node_ids) const
603 {
604 node_ids = {};
605
606 const auto &gbases = geom_bases();
607 for (const auto &lb : total_local_boundary)
608 {
609 const int e = lb.element_id();
610 for (int i = 0; i < lb.size(); ++i)
611 {
612 const int primitive_global_id = lb.global_primitive_id(i);
613 const auto nodes = gbases[e].local_nodes_for_primitive(primitive_global_id, *mesh);
614
615 for (long n = 0; n < nodes.size(); ++n)
616 {
617 const int g_id = gbases[e].bases[nodes(n)].global()[0].index;
618
619 if (std::count(node_ids.begin(), node_ids.end(), g_id) == 0)
620 node_ids.push_back(g_id);
621 }
622 }
623 }
624 }
625
626 void State::compute_volume_node_ids(const int volume_selection, std::vector<int> &node_ids) const
627 {
628 node_ids = {};
629
630 const auto &gbases = geom_bases();
631 for (int e = 0; e < gbases.size(); e++)
632 {
633 const int body_id = mesh->get_body_id(e);
634 if (body_id == volume_selection)
635 for (const auto &gbs : gbases[e].bases)
636 for (const auto &g : gbs.global())
637 node_ids.push_back(g.index);
638 }
639 }
640
641} // namespace polyfem
int x
void get_vertices(Eigen::MatrixXd &vertices) const
Definition StateDiff.cpp:72
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)
assembler::AssemblyValsCache ass_vals_cache
used to store assembly values for small problems
Definition State.hpp:196
void set_mesh_vertex(int v_id, const Eigen::VectorXd &vertex)
Definition StateDiff.cpp:97
const std::vector< basis::ElementBases > & geom_bases() const
Get a constant reference to the geometry mapping bases.
Definition State.hpp:223
std::shared_ptr< assembler::ViscousDamping > damping_assembler
Definition State.hpp:164
std::shared_ptr< assembler::Assembler > assembler
assemblers
Definition State.hpp:155
ipc::CollisionMesh collision_mesh
IPC collision mesh.
Definition State.hpp:520
solver::CacheLevel optimization_enabled
Definition State.hpp:654
bool is_homogenization() const
Definition State.hpp:718
Eigen::MatrixXd solve_static_adjoint(const Eigen::MatrixXd &adjoint_rhs) const
void compute_force_jacobian(const Eigen::MatrixXd &sol, const Eigen::MatrixXd &disp_grad, StiffnessMatrix &hessian)
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
bool has_periodic_bc() const
Definition State.hpp:387
std::shared_ptr< assembler::Problem > problem
current problem, it contains rhs and bc
Definition State.hpp:168
std::vector< int > node_to_primitive() const
Definition State.cpp:234
json args
main input arguments containing all defaults
Definition State.hpp:101
solver::DiffCache diff_cached
Definition State.hpp:656
std::vector< basis::ElementBases > bases
FE bases, the size is #elements.
Definition State.hpp:171
void solve_adjoint_cached(const Eigen::MatrixXd &rhs)
Eigen::MatrixXd solve_adjoint(const Eigen::MatrixXd &rhs) const
void get_elements(Eigen::MatrixXi &elements) const
Definition StateDiff.cpp:80
void build_stiffness_mat(StiffnessMatrix &stiffness)
utility that builds the stiffness matrix and collects stats, used only for linear problems
void compute_volume_node_ids(const int volume_selection, std::vector< int > &node_ids) const
std::vector< mesh::LocalBoundary > total_local_boundary
mapping from elements to nodes for all mesh
Definition State.hpp:436
int ndof() const
Definition State.hpp:660
std::vector< int > boundary_nodes
list of boundary nodes
Definition State.hpp:432
solver::SolveData solve_data
timedependent stuff cached
Definition State.hpp:321
std::unique_ptr< polysolve::linear::Solver > lin_solver_cached
Definition State.hpp:658
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
Eigen::MatrixXd solve_transient_adjoint(const Eigen::MatrixXd &adjoint_rhs) const
void compute_force_jacobian_prev(const int force_step, const int sol_step, StiffnessMatrix &hessian_prev) const
void compute_total_surface_node_ids(std::vector< int > &node_ids) const
void compute_surface_node_ids(const int surface_selection, std::vector< int > &node_ids) const
Eigen::MatrixXd rhs
System right-hand side.
Definition State.hpp:207
const StiffnessMatrix & gradu_h(int step) const
int bdf_order(int step) const
const ipc::TangentialCollisions & friction_collision_set(int step) const
const ipc::TangentialCollisions & tangential_adhesion_collision_set(int step) const
void cache_quantities_static(const Eigen::MatrixXd &u, const StiffnessMatrix &gradu_h, const ipc::NormalCollisions &collision_set, const ipc::SmoothCollisions &smooth_collision_set, const ipc::TangentialCollisions &friction_constraint_set, const ipc::NormalCollisions &normal_adhesion_set, const ipc::TangentialCollisions &tangential_adhesion_set, const Eigen::MatrixXd &disp_grad)
Definition DiffCache.hpp:43
void init(const int dimension, const int ndof, const int n_time_steps=0)
Definition DiffCache.hpp:21
void cache_quantities_quasistatic(const int cur_step, const Eigen::MatrixXd &u, const StiffnessMatrix &gradu_h, const ipc::NormalCollisions &collision_set, const ipc::SmoothCollisions &smooth_collision_set, const ipc::NormalCollisions &normal_adhesion_set, const Eigen::MatrixXd &disp_grad)
Definition DiffCache.hpp:94
void cache_quantities_transient(const int cur_step, const int cur_bdf_order, const Eigen::MatrixXd &u, const Eigen::MatrixXd &v, const Eigen::MatrixXd &acc, const StiffnessMatrix &gradu_h, const ipc::NormalCollisions &collision_set, const ipc::SmoothCollisions &smooth_collision_set, const ipc::TangentialCollisions &friction_collision_set)
Definition DiffCache.hpp:66
void cache_adjoints(const Eigen::MatrixXd &adjoint_mat)
Eigen::VectorXd u(int step) const
std::shared_ptr< solver::FrictionForm > friction_form
std::shared_ptr< solver::BodyForm > body_form
std::shared_ptr< solver::NLProblem > nl_problem
std::shared_ptr< solver::NormalAdhesionForm > normal_adhesion_form
std::shared_ptr< solver::ContactForm > contact_form
std::shared_ptr< time_integrator::ImplicitTimeIntegrator > time_integrator
std::shared_ptr< solver::TangentialAdhesionForm > tangential_adhesion_form
Backward Differential Formulas.
Definition BDF.hpp:14
static double betas(const int i)
Retrieve the value of beta used for BDF with i steps.
Definition BDF.cpp:35
Eigen::VectorXd weighted_sum_v_prevs() const
Compute the weighted sum of the previous velocities.
Definition BDF.cpp:62
static const std::vector< double > & alphas(const int i)
Retrieve the alphas used for BDF with i steps.
Definition BDF.cpp:21
Implicit Euler time integrator of a second order ODE (equivently a system of coupled first order ODEs...
const Eigen::VectorXd & v_prev() const
Get the most recent previous velocity value.
Eigen::MatrixXd unflatten(const Eigen::VectorXd &x, int dim)
Unflatten rowwises, so every dim elements in x become a row.
spdlog::logger & adjoint_logger()
Retrieves the current logger for adjoint.
Definition Logger.cpp:30
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