8#include <polysolve/linear/FEMSolver.hpp>
28#include <ipc/barrier/barrier.hpp>
29#include <ipc/utils/local_to_global.hpp>
30#include <ipc/potentials/friction_potential.hpp>
33#include <unsupported/Eigen/SparseExtra>
48 reduced_mat.resize(mat.rows(), mat.cols());
50 std::vector<bool> mask(mat.rows(),
false);
54 std::vector<Eigen::Triplet<double>> coeffs;
55 for (
int k = 0; k < mat.outerSize(); ++k)
57 for (StiffnessMatrix::InnerIterator it(mat, k); it; ++it)
61 if (it.row() == it.col())
62 coeffs.emplace_back(it.row(), it.col(), 1.0);
65 coeffs.emplace_back(it.row(), it.col(), it.value());
68 reduced_mat.setFromTriplets(coeffs.begin(), coeffs.end());
74 vertices.setZero(
mesh->n_vertices(),
mesh->dimension());
76 for (
int v = 0; v <
mesh->n_vertices(); v++)
77 vertices.row(v) =
mesh->point(v);
82 assert(
mesh->is_simplicial());
87 int dim =
mesh->dimension();
88 elements.setZero(gbases.size(), dim + 1);
89 for (
int e = 0; e < gbases.size(); e++)
92 for (
const auto &gbs : gbases[e].bases)
93 elements(e, i++) = node_to_primitive_map[gbs.global()[0].index];
99 assert(vertex.size() ==
mesh->dimension());
100 mesh->set_point(v_id, vertex);
106 if (current_step == 0)
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;
117 if (!
problem->is_time_dependent() || current_step > 0)
123 cur_collision_set = barrier_contact->collision_set();
125 cur_smooth_collision_set = smooth_contact->collision_set();
132 if (
problem->is_time_dependent())
134 if (
args[
"time"][
"quasistatic"].get<bool>())
140 Eigen::MatrixXd vel, acc;
141 if (current_step == 0)
151 vel = euler_integrator->
v_prev();
156 acc.setZero(
ndof(), 1);
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);
175 if (
problem->is_time_dependent())
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);
221 assert(force_step > 0);
222 assert(force_step > sol_step);
235 if (
problem->is_time_dependent())
239 if (sol_step == force_step - 1)
245 const Eigen::MatrixXd surface_velocities = (surface_solution - surface_solution_prev) / dt;
246 const double dv_dut = -1 / dt;
255 surface_solution_prev,
257 barrier_contact->barrier_potential(),
258 barrier_contact->barrier_stiffness(),
259 ipc::FrictionPotential::DiffWRT::LAGGED_DISPLACEMENTS)
264 surface_solution_prev,
266 barrier_contact->barrier_potential(),
267 barrier_contact->barrier_stiffness(),
268 ipc::FrictionPotential::DiffWRT::VELOCITIES)
321 if (sol_step == force_step - 1)
329 const Eigen::MatrixXd surface_velocities = (surface_solution - surface_solution_prev) / dt;
330 const double dv_dut = -1 / dt;
332 adhesion_hessian_prev =
337 surface_solution_prev,
341 ipc::TangentialPotential::DiffWRT::LAGGED_DISPLACEMENTS)
346 surface_solution_prev,
350 ipc::TangentialPotential::DiffWRT::VELOCITIES)
353 adhesion_hessian_prev *= -1;
355 adhesion_hessian_prev =
collision_mesh.to_full_dof(adhesion_hessian_prev);
357 hessian_prev += adhesion_hessian_prev;
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);
367 hessian_prev += damping_hessian_prev;
370 if (sol_step == force_step - 1)
374 hessian_prev += body_force_hessian;
387 if (
problem->is_time_dependent())
395 Eigen::MatrixXd b = adjoint_rhs;
397 Eigen::MatrixXd adjoint;
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;
407 b.conservativeResizeLike(Eigen::MatrixXd::Zero(A.rows(), b.cols()));
409 std::vector<int> boundary_nodes_tmp;
419 adjoint.setZero(
ndof(), adjoint_rhs.cols());
420 for (
int i = 0; i < b.cols(); i++)
422 Eigen::VectorXd
x, tmp;
427 adjoint.col(i) =
periodic_bc->periodic_to_full(full_size,
x);
434 auto solver = polysolve::linear::Solver::create(
args[
"solver"][
"adjoint_linear"],
adjoint_logger());
444 adjoint.setZero(
ndof(), adjoint_rhs.cols());
445 for (
int i = 0; i < b.cols(); i++)
447 Eigen::VectorXd tmp = b.col(i);
451 x.setZero(tmp.size());
452 dirichlet_solve(*solver, A, tmp,
boundary_nodes,
x, A.rows(),
"",
false,
false,
false);
460 solver->analyze_pattern(A, A.rows());
461 solver->factorize(A);
463 adjoint.setZero(adjoint_rhs.rows(), adjoint_rhs.cols());
464 for (
int i = 0; i < b.cols(); i++)
466 Eigen::MatrixXd tmp = b.col(i);
469 x.setZero(tmp.size());
470 solver->solve(tmp,
x);
471 x.conservativeResize(adjoint.rows());
483 const double dt =
args[
"time"][
"dt"];
484 const int time_steps =
args[
"time"][
"time_steps"];
487 if (
args[
"time"][
"integrator"].is_string())
489 else if (
args[
"time"][
"integrator"][
"type"] ==
"ImplicitEuler")
491 else if (
args[
"time"][
"integrator"][
"type"] ==
"BDF")
492 bdf_order =
args[
"time"][
"integrator"][
"steps"].get<
int>();
496 assert(adjoint_rhs.cols() == time_steps + 1);
498 const int cols_per_adjoint = time_steps + 1;
499 Eigen::MatrixXd adjoints;
500 adjoints.setZero(
ndof(), cols_per_adjoint * 2);
506 Eigen::MatrixXd sum_alpha_p, sum_alpha_nu;
507 for (
int i = time_steps; i >= 0; --i)
510 sum_alpha_p.setZero(
ndof(), 1);
511 sum_alpha_nu.setZero(
ndof(), 1);
513 const int num = std::min(bdf_order, time_steps - i);
515 Eigen::VectorXd bdf_coeffs = Eigen::VectorXd::Zero(num);
516 for (
int j = 0; j < bdf_order && i + j < time_steps; ++j)
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;
523 Eigen::VectorXd rhs_ = -reduced_mass.transpose() * sum_alpha_nu - adjoint_rhs.col(i);
524 for (
int j = 1; j <= bdf_order; j++)
526 if (i + j > time_steps)
533 rhs_ += -gradu_h_prev.transpose() * tmp;
540 rhs_ += (1. / beta_dt) * (
diff_cached.
gradu_h(i) - reduced_mass).transpose() * sum_alpha_p;
544 Eigen::VectorXd b_ = rhs_;
547 auto solver = polysolve::linear::Solver::create(
args[
"solver"][
"adjoint_linear"],
adjoint_logger());
550 dirichlet_solve(*solver, A, b_,
boundary_nodes,
x, A.rows(),
"",
false,
false,
false);
551 adjoints.col(i + cols_per_adjoint) =
x;
556 if (i + 1 < cols_per_adjoint)
558 if (i + 2 < cols_per_adjoint)
563 adjoints.col(i) = beta_dt * adjoints.col(i + cols_per_adjoint) - sum_alpha_p;
567 adjoints.col(i) = -reduced_mass.transpose() * sum_alpha_p;
568 adjoints.col(i + cols_per_adjoint) = rhs_;
581 const int e = lb.element_id();
582 for (
int i = 0; i < lb.size(); ++i)
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);
588 if (boundary_id == surface_selection)
590 for (
long n = 0; n < nodes.size(); ++n)
592 const int g_id = gbases[e].bases[nodes(n)].global()[0].index;
594 if (std::count(node_ids.begin(), node_ids.end(), g_id) == 0)
595 node_ids.push_back(g_id);
609 const int e = lb.element_id();
610 for (
int i = 0; i < lb.size(); ++i)
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);
615 for (
long n = 0; n < nodes.size(); ++n)
617 const int g_id = gbases[e].bases[nodes(n)].global()[0].index;
619 if (std::count(node_ids.begin(), node_ids.end(), g_id) == 0)
620 node_ids.push_back(g_id);
631 for (
int e = 0; e < gbases.size(); e++)
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);
void get_vertices(Eigen::MatrixXd &vertices) const
std::shared_ptr< utils::PeriodicBoundary > periodic_bc
periodic BC and periodic mesh utils
int n_bases
number of bases
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
void set_mesh_vertex(int v_id, const Eigen::VectorXd &vertex)
const std::vector< basis::ElementBases > & geom_bases() const
Get a constant reference to the geometry mapping bases.
std::shared_ptr< assembler::ViscousDamping > damping_assembler
std::shared_ptr< assembler::Assembler > assembler
assemblers
ipc::CollisionMesh collision_mesh
IPC collision mesh.
solver::CacheLevel optimization_enabled
bool is_homogenization() const
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
StiffnessMatrix mass
Mass matrix, it is computed only for time dependent problems.
bool has_periodic_bc() const
std::shared_ptr< assembler::Problem > problem
current problem, it contains rhs and bc
std::vector< int > node_to_primitive() const
json args
main input arguments containing all defaults
solver::DiffCache diff_cached
std::vector< basis::ElementBases > bases
FE bases, the size is #elements.
void solve_adjoint_cached(const Eigen::MatrixXd &rhs)
Eigen::MatrixXd solve_adjoint(const Eigen::MatrixXd &rhs) const
void get_elements(Eigen::MatrixXi &elements) const
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
std::vector< int > boundary_nodes
list of boundary nodes
solver::SolveData solve_data
timedependent stuff cached
std::unique_ptr< polysolve::linear::Solver > lin_solver_cached
std::shared_ptr< assembler::ViscousDampingPrev > damping_prev_assembler
bool is_contact_enabled() const
does the simulation have contact
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.
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)
void init(const int dimension, const int ndof, const int n_time_steps=0)
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)
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)
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.
static double betas(const int i)
Retrieve the value of beta used for BDF with i steps.
Eigen::VectorXd weighted_sum_v_prevs() const
Compute the weighted sum of the previous velocities.
static const std::vector< double > & alphas(const int i)
Retrieve the alphas used for BDF with i steps.
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.
void log_and_throw_adjoint_error(const std::string &msg)
void log_and_throw_error(const std::string &msg)
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix