18 const Eigen::SparseMatrix<double> &M,
19 const Eigen::SparseMatrix<double> &A,
20 const Eigen::Ref<const Eigen::MatrixXd> &
y)
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",
"");
29 solver = polysolve::linear::Solver::create(
"Eigen::SimplicialLDLT",
"");
32 solver->analyze_pattern(M, 0);
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));
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));
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)
54 std::vector<int> free_nodes;
55 for (
int i = 0, j = 0; i <
y.rows(); ++i)
57 if (boundary_nodes[j] == i)
60 free_nodes.push_back(i);
63 const Eigen::MatrixXd H = M(free_nodes, free_nodes);
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;
72 std::shared_ptr<polysolve::nonlinear::Solver> nl_solver,
74 const Eigen::SparseMatrix<double> &M,
75 const Eigen::SparseMatrix<double> &A,
76 const Eigen::VectorXd &
y,
78 const Eigen::MatrixXd &rest_positions,
79 const Eigen::MatrixXi &elements,
82 const ipc::CollisionMesh &collision_mesh,
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,
92 const std::vector<int> &boundary_nodes,
93 const size_t obstacle_ndof,
94 const Eigen::VectorXd &target_x,
96 const Eigen::VectorXd &x0)
100 assert(M.rows() == M.cols());
101 assert(A.rows() == M.rows());
102 assert(A.cols() ==
y.size());
104 std::vector<std::shared_ptr<Form>> forms;
106 forms.push_back(std::make_shared<L2ProjectionForm>(M, A,
y));
108 forms.push_back(std::make_shared<InversionBarrierForm>(
109 rest_positions, elements, dim, 1e-12));
110 forms.back()->set_weight(barrier_stiffness);
111 assert(forms.back()->is_step_valid(x0, x0));
113 if (collision_mesh.num_vertices() != 0)
115 forms.push_back(std::make_shared<ContactForm>(
116 collision_mesh, dhat, 1.0, use_area_weighting, use_improved_max_operator,
117 use_physical_barrier,
false,
true,
118 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))));
124 const int ndof = x0.size();
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);
132 std::vector<std::shared_ptr<AugmentedLagrangianForm>> bc_forms;
133 bc_forms.push_back(bc_lagrangian_form);
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;
147 bc_forms, al_initial_weight,
148 al_scaling, al_max_weight, al_eta_tol,
149 [&](
const Eigen::MatrixXd &
x) {});
151 Eigen::MatrixXd sol = x0;
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);
157 nl_solver->stop_criteria().iterations = default_max_iterations;
158 al_solver.solve_reduced(problem, sol, nl_solver);
161 assert(forms[1]->is_step_valid(sol, sol));
162 if (collision_mesh.num_vertices() != 0)
163 assert(!ipc::has_intersections(collision_mesh, collision_mesh.displace_vertices(
utils::unflatten(sol, dim))));
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)