Loading [MathJax]/extensions/tex2jax.js
PolyFEM
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
L2Projection.cpp
Go to the documentation of this file.
1#include "L2Projection.hpp"
2
10
11#include <ipc/ipc.hpp>
12
13#include <polysolve/linear/Solver.hpp>
14
15namespace polyfem::mesh
16{
18 const Eigen::SparseMatrix<double> &M,
19 const Eigen::SparseMatrix<double> &A,
20 const Eigen::Ref<const Eigen::MatrixXd> &y)
21 {
22 // Construct a linear solver for M
23 std::unique_ptr<polysolve::linear::Solver> solver;
24#ifdef POLYSOLVE_WITH_MKL
25 solver = polysolve::linear::Solver::create("Eigen::PardisoLDLT", "");
26#elif defined(POLYSOLVE_WITH_CHOLMOD)
27 solver = polysolve::linear::Solver::create("Eigen::CholmodSimplicialLDLT", "");
28#else
29 solver = polysolve::linear::Solver::create("Eigen::SimplicialLDLT", "");
30#endif
31
32 solver->analyze_pattern(M, 0);
33 solver->factorize(M);
34
35 const Eigen::MatrixXd rhs = A * y;
36 Eigen::MatrixXd x(rhs.rows(), rhs.cols());
37 for (int i = 0; i < x.cols(); ++i)
38 solver->solve(rhs.col(i), x.col(i));
39
40 double residual_error = (M * x - rhs).norm();
41 logger().debug("residual error in L2 projection: {}", residual_error);
42 assert(residual_error < std::max(1e-12 * rhs.norm(), 1e-12));
43
44 return x;
45 }
46
48 const Eigen::MatrixXd &M,
49 const Eigen::MatrixXd &A,
50 const Eigen::Ref<const Eigen::MatrixXd> &y,
51 const std::vector<int> &boundary_nodes,
52 Eigen::Ref<Eigen::MatrixXd> x)
53 {
54 std::vector<int> free_nodes;
55 for (int i = 0, j = 0; i < y.rows(); ++i)
56 {
57 if (boundary_nodes[j] == i)
58 ++j;
59 else
60 free_nodes.push_back(i);
61 }
62
63 const Eigen::MatrixXd H = M(free_nodes, free_nodes);
64
65 const Eigen::MatrixXd g = -((M * x - A * y)(free_nodes, Eigen::all));
66 const Eigen::MatrixXd sol = H.llt().solve(g);
67 x(free_nodes, Eigen::all) += sol;
68 }
69
71 // Nonlinear solver
72 std::shared_ptr<polysolve::nonlinear::Solver> nl_solver,
73 // L2 projection form
74 const Eigen::SparseMatrix<double> &M,
75 const Eigen::SparseMatrix<double> &A,
76 const Eigen::VectorXd &y,
77 // Inversion-free form
78 const Eigen::MatrixXd &rest_positions,
79 const Eigen::MatrixXi &elements,
80 const int dim,
81 // Contact form
82 const ipc::CollisionMesh &collision_mesh,
83 const double dhat,
84 const double barrier_stiffness,
85 const bool use_area_weighting,
86 const bool use_improved_max_operator,
87 const bool use_physical_barrier,
88 const ipc::BroadPhaseMethod broad_phase_method,
89 const double ccd_tolerance,
90 const int ccd_max_iterations,
91 // Augmented lagrangian form
92 const std::vector<int> &boundary_nodes,
93 const size_t obstacle_ndof,
94 const Eigen::VectorXd &target_x,
95 // Initial guess
96 const Eigen::VectorXd &x0)
97 {
98 using namespace polyfem::solver;
99
100 assert(M.rows() == M.cols());
101 assert(A.rows() == M.rows());
102 assert(A.cols() == y.size());
103
104 std::vector<std::shared_ptr<Form>> forms;
105
106 forms.push_back(std::make_shared<L2ProjectionForm>(M, A, y));
107
108 forms.push_back(std::make_shared<InversionBarrierForm>(
109 rest_positions, elements, dim, /*vhat=*/1e-12));
110 forms.back()->set_weight(barrier_stiffness); // use same weight as barrier stiffness
111 assert(forms.back()->is_step_valid(x0, x0));
112
113 if (collision_mesh.num_vertices() != 0)
114 {
115 forms.push_back(std::make_shared<ContactForm>(
116 collision_mesh, dhat, /*avg_mass=*/1.0, use_area_weighting, use_improved_max_operator,
117 use_physical_barrier, /*use_adaptive_barrier_stiffness=*/false, /*is_time_dependent=*/true,
118 /*enable_shape_derivatives=*/false, broad_phase_method, ccd_tolerance,
119 ccd_max_iterations));
120 forms.back()->set_weight(barrier_stiffness);
121 assert(!ipc::has_intersections(collision_mesh, collision_mesh.displace_vertices(utils::unflatten(x0, dim))));
122 }
123
124 const int ndof = x0.size();
125
126 std::shared_ptr<BCLagrangianForm> bc_lagrangian_form = std::make_shared<BCLagrangianForm>(
127 ndof, boundary_nodes, M, obstacle_ndof, target_x);
128 forms.push_back(bc_lagrangian_form);
129
130 // --------------------------------------------------------------------
131
132 std::vector<std::shared_ptr<AugmentedLagrangianForm>> bc_forms;
133 bc_forms.push_back(bc_lagrangian_form);
134
135 StaticBoundaryNLProblem problem(ndof, target_x, forms, bc_forms);
136
137 // --------------------------------------------------------------------
138
139 // Create augmented Lagrangian solver
140 // AL parameters
141 constexpr double al_initial_weight = 1e6;
142 constexpr double al_scaling = 2.0;
143 constexpr int al_max_weight = 100 * al_initial_weight;
144 constexpr double al_eta_tol = 0.99;
145 constexpr size_t al_max_solver_iter = 1000;
146 ALSolver al_solver(
147 bc_forms, al_initial_weight,
148 al_scaling, al_max_weight, al_eta_tol,
149 /*update_barrier_stiffness=*/[&](const Eigen::MatrixXd &x) {});
150
151 Eigen::MatrixXd sol = x0;
152
153 const size_t default_max_iterations = nl_solver->stop_criteria().iterations;
154 nl_solver->stop_criteria().iterations = al_max_solver_iter;
155 al_solver.solve_al(problem, sol, nl_solver);
156
157 nl_solver->stop_criteria().iterations = default_max_iterations;
158 al_solver.solve_reduced(problem, sol, nl_solver);
159
160#ifndef NDEBUG
161 assert(forms[1]->is_step_valid(sol, sol)); // inversion-free
162 if (collision_mesh.num_vertices() != 0)
163 assert(!ipc::has_intersections(collision_mesh, collision_mesh.displace_vertices(utils::unflatten(sol, dim))));
164#endif
165
166 return sol;
167 }
168} // namespace polyfem::mesh
int y
int x
void solve_al(NLProblem &nl_problem, Eigen::MatrixXd &sol, std::shared_ptr< polysolve::nonlinear::Solver > nl_solver)
Definition ALSolver.hpp:29
Eigen::MatrixXd unconstrained_L2_projection(const Eigen::SparseMatrix< double > &M, const Eigen::SparseMatrix< double > &A, const Eigen::Ref< const Eigen::MatrixXd > &y)
void reduced_L2_projection(const Eigen::MatrixXd &M, const Eigen::MatrixXd &A, const Eigen::Ref< const Eigen::MatrixXd > &y, const std::vector< int > &boundary_nodes, Eigen::Ref< Eigen::MatrixXd > x)
Eigen::VectorXd constrained_L2_projection(std::shared_ptr< polysolve::nonlinear::Solver > nl_solver, const Eigen::SparseMatrix< double > &M, const Eigen::SparseMatrix< double > &A, const Eigen::VectorXd &y, const Eigen::MatrixXd &rest_positions, const Eigen::MatrixXi &elements, const int dim, const ipc::CollisionMesh &collision_mesh, const double dhat, const double barrier_stiffness, const bool use_area_weighting, const bool use_improved_max_operator, const bool use_physical_barrier, const ipc::BroadPhaseMethod broad_phase_method, const double ccd_tolerance, const int ccd_max_iterations, const std::vector< int > &boundary_nodes, const size_t obstacle_ndof, const Eigen::VectorXd &target_x, const Eigen::VectorXd &x0)
Eigen::MatrixXd unflatten(const Eigen::VectorXd &x, int dim)
Unflatten rowwises, so every dim elements in x become a row.
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:44