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