8#include <polysolve/linear/FEMSolver.hpp>
27#include <ipc/barrier/barrier.hpp>
28#include <ipc/utils/local_to_global.hpp>
29#include <ipc/potentials/friction_potential.hpp>
32#include <unsupported/Eigen/SparseExtra>
47 reduced_mat.resize(mat.rows(), mat.cols());
49 std::vector<bool> mask(mat.rows(),
false);
53 std::vector<Eigen::Triplet<double>> coeffs;
54 for (
int k = 0; k < mat.outerSize(); ++k)
56 for (StiffnessMatrix::InnerIterator it(mat, k); it; ++it)
60 if (it.row() == it.col())
61 coeffs.emplace_back(it.row(), it.col(), 1.0);
64 coeffs.emplace_back(it.row(), it.col(), it.value());
67 reduced_mat.setFromTriplets(coeffs.begin(), coeffs.end());
73 vertices.setZero(
mesh->n_vertices(),
mesh->dimension());
75 for (
int v = 0; v <
mesh->n_vertices(); v++)
76 vertices.row(v) =
mesh->point(v);
81 assert(
mesh->is_simplicial());
86 int dim =
mesh->dimension();
87 elements.setZero(gbases.size(), dim + 1);
88 for (
int e = 0; e < gbases.size(); e++)
91 for (
const auto &gbs : gbases[e].bases)
92 elements(e, i++) = node_to_primitive_map[gbs.global()[0].index];
98 assert(vertex.size() ==
mesh->dimension());
99 mesh->set_point(v_id, vertex);
105 if (current_step == 0)
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;
115 if (!
problem->is_time_dependent() || current_step > 0)
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();
131 if (
problem->is_time_dependent())
133 if (
args[
"time"][
"quasistatic"].get<bool>())
139 Eigen::MatrixXd vel, acc;
140 if (current_step == 0)
150 vel = euler_integrator->
v_prev();
155 acc.setZero(
ndof(), 1);
174 if (
problem->is_time_dependent())
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);
220 assert(force_step > 0);
221 assert(force_step > sol_step);
234 if (
problem->is_time_dependent())
238 if (sol_step == force_step - 1)
244 const Eigen::MatrixXd surface_velocities = (surface_solution - surface_solution_prev) / dt;
245 const double dv_dut = -1 / dt;
252 surface_solution_prev,
256 ipc::FrictionPotential::DiffWRT::LAGGED_DISPLACEMENTS)
261 surface_solution_prev,
265 ipc::FrictionPotential::DiffWRT::VELOCITIES)
317 if (sol_step == force_step - 1)
325 const Eigen::MatrixXd surface_velocities = (surface_solution - surface_solution_prev) / dt;
326 const double dv_dut = -1 / dt;
328 adhesion_hessian_prev =
333 surface_solution_prev,
337 ipc::TangentialPotential::DiffWRT::LAGGED_DISPLACEMENTS)
342 surface_solution_prev,
346 ipc::TangentialPotential::DiffWRT::VELOCITIES)
349 adhesion_hessian_prev *= -1;
351 adhesion_hessian_prev =
collision_mesh.to_full_dof(adhesion_hessian_prev);
353 hessian_prev += adhesion_hessian_prev;
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);
363 hessian_prev += damping_hessian_prev;
366 if (sol_step == force_step - 1)
370 hessian_prev += body_force_hessian;
383 if (
problem->is_time_dependent())
391 Eigen::MatrixXd b = adjoint_rhs;
393 Eigen::MatrixXd adjoint;
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;
403 b.conservativeResizeLike(Eigen::MatrixXd::Zero(A.rows(), b.cols()));
405 std::vector<int> boundary_nodes_tmp;
415 adjoint.setZero(
ndof(), adjoint_rhs.cols());
416 for (
int i = 0; i < b.cols(); i++)
418 Eigen::VectorXd
x, tmp;
423 adjoint.col(i) =
periodic_bc->periodic_to_full(full_size,
x);
430 auto solver = polysolve::linear::Solver::create(
args[
"solver"][
"adjoint_linear"],
adjoint_logger());
440 adjoint.setZero(
ndof(), adjoint_rhs.cols());
441 for (
int i = 0; i < b.cols(); i++)
443 Eigen::VectorXd tmp = b.col(i);
447 x.setZero(tmp.size());
448 dirichlet_solve(*solver, A, tmp,
boundary_nodes,
x, A.rows(),
"",
false,
false,
false);
456 solver->analyze_pattern(A, A.rows());
457 solver->factorize(A);
459 adjoint.setZero(adjoint_rhs.rows(), adjoint_rhs.cols());
460 for (
int i = 0; i < b.cols(); i++)
462 Eigen::MatrixXd tmp = b.col(i);
465 x.setZero(tmp.size());
466 solver->solve(tmp,
x);
467 x.conservativeResize(adjoint.rows());
479 const double dt =
args[
"time"][
"dt"];
480 const int time_steps =
args[
"time"][
"time_steps"];
483 if (
args[
"time"][
"integrator"].is_string())
485 else if (
args[
"time"][
"integrator"][
"type"] ==
"ImplicitEuler")
487 else if (
args[
"time"][
"integrator"][
"type"] ==
"BDF")
488 bdf_order =
args[
"time"][
"integrator"][
"steps"].get<
int>();
492 assert(adjoint_rhs.cols() == time_steps + 1);
494 const int cols_per_adjoint = time_steps + 1;
495 Eigen::MatrixXd adjoints;
496 adjoints.setZero(
ndof(), cols_per_adjoint * 2);
502 Eigen::MatrixXd sum_alpha_p, sum_alpha_nu;
503 for (
int i = time_steps; i >= 0; --i)
506 sum_alpha_p.setZero(
ndof(), 1);
507 sum_alpha_nu.setZero(
ndof(), 1);
509 const int num = std::min(bdf_order, time_steps - i);
511 Eigen::VectorXd bdf_coeffs = Eigen::VectorXd::Zero(num);
512 for (
int j = 0; j < bdf_order && i + j < time_steps; ++j)
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;
519 Eigen::VectorXd rhs_ = -reduced_mass.transpose() * sum_alpha_nu - adjoint_rhs.col(i);
520 for (
int j = 1; j <= bdf_order; j++)
522 if (i + j > time_steps)
529 rhs_ += -gradu_h_prev.transpose() * tmp;
536 rhs_ += (1. / beta_dt) * (
diff_cached.
gradu_h(i) - reduced_mass).transpose() * sum_alpha_p;
540 Eigen::VectorXd b_ = rhs_;
543 auto solver = polysolve::linear::Solver::create(
args[
"solver"][
"adjoint_linear"],
adjoint_logger());
546 dirichlet_solve(*solver, A, b_,
boundary_nodes,
x, A.rows(),
"",
false,
false,
false);
547 adjoints.col(i + cols_per_adjoint) =
x;
552 if (i + 1 < cols_per_adjoint)
554 if (i + 2 < cols_per_adjoint)
559 adjoints.col(i) = beta_dt * adjoints.col(i + cols_per_adjoint) - sum_alpha_p;
563 adjoints.col(i) = -reduced_mass.transpose() * sum_alpha_p;
564 adjoints.col(i + cols_per_adjoint) = rhs_;
577 const int e = lb.element_id();
578 for (
int i = 0; i < lb.size(); ++i)
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);
584 if (boundary_id == surface_selection)
586 for (
long n = 0; n < nodes.size(); ++n)
588 const int g_id = gbases[e].bases[nodes(n)].global()[0].index;
590 if (std::count(node_ids.begin(), node_ids.end(), g_id) == 0)
591 node_ids.push_back(g_id);
605 const int e = lb.element_id();
606 for (
int i = 0; i < lb.size(); ++i)
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);
611 for (
long n = 0; n < nodes.size(); ++n)
613 const int g_id = gbases[e].bases[nodes(n)].global()[0].index;
615 if (std::count(node_ids.begin(), node_ids.end(), g_id) == 0)
616 node_ids.push_back(g_id);
627 for (
int e = 0; e < gbases.size(); e++)
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);
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
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)
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)
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)
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)
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