Loading [MathJax]/extensions/tex2jax.js
PolyFEM
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
NLHomoProblem.cpp
Go to the documentation of this file.
1#include "NLHomoProblem.hpp"
2#include <polyfem/State.hpp>
7
8namespace polyfem::solver
9{
10 NLHomoProblem::NLHomoProblem(const int full_size,
11 const assembler::MacroStrainValue &macro_strain_constraint,
12 const State &state,
13 const double t,
14 const std::vector<std::shared_ptr<Form>> &forms,
15 const std::vector<std::shared_ptr<AugmentedLagrangianForm>> &penalty_forms,
16 const bool solve_symmetric_macro_strain,
17 const std::shared_ptr<polysolve::linear::Solver> &solver)
18 : NLProblem(full_size, state.periodic_bc, t, forms, penalty_forms, solver),
19 state_(state),
20 only_symmetric(solve_symmetric_macro_strain),
21 macro_strain_constraint_(macro_strain_constraint)
22 {
24 }
25
27 {
28 const int dim = state_.mesh->dimension();
30 {
31 macro_mid_to_full_.setZero(dim * dim, (dim * (dim + 1)) / 2);
32 macro_full_to_mid_.setZero((dim * (dim + 1)) / 2, dim * dim);
33 for (int i = 0, idx = 0; i < dim; i++)
34 {
35 for (int j = i; j < dim; j++)
36 {
37 macro_full_to_mid_(idx, i * dim + j) = 1;
38
39 macro_mid_to_full_(j * dim + i, idx) = 1;
40 macro_mid_to_full_(i * dim + j, idx) = 1;
41
42 idx++;
43 }
44 }
45 }
46 else
47 {
48 macro_mid_to_full_.setIdentity(dim * dim, dim * dim);
49 macro_full_to_mid_.setIdentity(dim * dim, dim * dim);
50 }
52 }
53
54 Eigen::VectorXd NLHomoProblem::extended_to_reduced(const Eigen::VectorXd &extended) const
55 {
56 const int dim = state_.mesh->dimension();
57 const int dof2 = macro_reduced_size();
58 const int dof1 = reduced_size();
59
60 Eigen::VectorXd reduced(dof1 + dof2);
61 reduced.head(dof1) = NLProblem::full_to_reduced(extended.head(extended.size() - dim * dim));
62 reduced.tail(dof2) = macro_full_to_reduced(extended.tail(dim * dim));
63
64 return reduced;
65 }
66
67 Eigen::VectorXd NLHomoProblem::reduced_to_extended(const Eigen::VectorXd &reduced, bool homogeneous) const
68 {
69 const int dim = state_.mesh->dimension();
70 const int dof2 = macro_reduced_size();
71 const int dof1 = reduced_size();
72 assert(reduced.size() == dof1 + dof2);
73
74 Eigen::VectorXd full = NLProblem::reduced_to_full(reduced.head(dof1));
75 Eigen::VectorXd disp_grad = macro_reduced_to_full(reduced.tail(dof2), homogeneous);
76 Eigen::VectorXd extended(full.size() + disp_grad.size());
77 extended << full, disp_grad;
78
79 return extended;
80 }
81
82 Eigen::VectorXd NLHomoProblem::extended_to_reduced_grad(const Eigen::VectorXd &extended) const
83 {
84 const int dim = state_.mesh->dimension();
85 const int dof2 = macro_reduced_size();
86 const int dof1 = reduced_size();
87
88 Eigen::VectorXd grad(dof1 + dof2);
89 grad.head(dof1) = NLProblem::full_to_reduced_grad(extended.head(extended.size() - dim * dim));
90 grad.tail(dof2) = macro_full_to_reduced_grad(extended.tail(dim * dim));
91
92 return grad;
93 }
94
95 NLHomoProblem::TVector NLHomoProblem::reduced_to_full_shape_derivative(const Eigen::MatrixXd &disp_grad, const TVector &adjoint_full) const
96 {
97 const int dim = state_.mesh->dimension();
98
99 Eigen::VectorXd term;
100 term.setZero(state_.n_bases * dim);
101 for (int i = 0; i < state_.n_bases; i++)
102 term.segment(i * dim, dim) += disp_grad.transpose() * adjoint_full.segment(i * dim, dim);
103
105 }
106
107 double NLHomoProblem::value(const TVector &x)
108 {
110
112 {
113 val += penalty_problem_->value(x);
114 }
115
116 for (auto &form : homo_forms)
117 if (form->enabled())
118 val += form->value(reduced_to_extended(x));
119
120 return val;
121 }
122 void NLHomoProblem::gradient(const TVector &x, TVector &gradv)
123 {
125
126 if (full_size() != current_size())
127 {
128 gradv = full_to_reduced_grad(gradv);
129 }
130 else if (penalty_problem_)
131 {
132 TVector tmp;
133 penalty_problem_->gradient(x, tmp);
134 gradv += tmp;
135 }
136
137 for (auto &form : homo_forms)
138 if (form->enabled())
139 {
140 Eigen::VectorXd grad_extended;
141 form->first_derivative(reduced_to_extended(x), grad_extended);
142 gradv += extended_to_reduced_grad(grad_extended);
143 }
144 }
145 void NLHomoProblem::extended_hessian_to_reduced_hessian(const THessian &extended, THessian &reduced) const
146 {
147 const int dim = state_.mesh->dimension();
148 const int dof2 = macro_reduced_size();
149 const int dof1 = reduced_size();
150
151 Eigen::MatrixXd A12, A22;
152 {
153 Eigen::MatrixXd tmp = Eigen::MatrixXd(extended.rightCols(dim * dim));
154 A12 = macro_full_to_reduced_grad(tmp.topRows(tmp.rows() - dim * dim).transpose()).transpose();
155 A22 = macro_full_to_reduced_grad(macro_full_to_reduced_grad(tmp.bottomRows(dim * dim)).transpose());
156 }
157
158 std::vector<Eigen::Triplet<double>> entries;
159 entries.reserve(extended.nonZeros() + A12.size() * 2 + A22.size());
160
161 for (int k = 0; k < full_size_; ++k)
162 for (StiffnessMatrix::InnerIterator it(extended, k); it; ++it)
163 if (it.row() < full_size_ && it.col() < full_size_)
164 entries.emplace_back(it.row(), it.col(), it.value());
165
166 for (int i = 0; i < A12.rows(); i++)
167 {
168 for (int j = 0; j < A12.cols(); j++)
169 {
170 entries.emplace_back(i, full_size_ + j, A12(i, j));
171 entries.emplace_back(full_size_ + j, i, A12(i, j));
172 }
173 }
174
175 for (int i = 0; i < A22.rows(); i++)
176 for (int j = 0; j < A22.cols(); j++)
177 entries.emplace_back(i + full_size_, j + full_size_, A22(i, j));
178
179 // Eigen::VectorXd tmp = macro_full_to_reduced_grad(Eigen::VectorXd::Ones(dim*dim));
180 // for (int i = 0; i < tmp.size(); i++)
181 // entries.emplace_back(i + full_size_, i + full_size_, 1 - tmp(i));
182
183 reduced.resize(0, 0);
184 reduced.resize(full_size_ + dof2, full_size_ + dof2);
185 reduced.setFromTriplets(entries.begin(), entries.end());
186
187 {
188 assert(dof1 == Q2_.cols());
189 StiffnessMatrix Q2_extended(Q2_.rows() + dof2, Q2_.cols() + dof2);
190
191 {
192 entries.clear();
193 for (int k = 0; k < Q2_.cols(); ++k)
194 for (StiffnessMatrix::InnerIterator it(Q2_, k); it; ++it)
195 {
196 entries.emplace_back(it.row(), it.col(), it.value());
197 }
198
199 for (int k = 0; k < dof2; k++)
200 {
201 entries.emplace_back(Q2_.rows() + k, dof1 + k, 1);
202 }
203
204 Q2_extended.setFromTriplets(entries.begin(), entries.end());
205 }
206
207 reduced = Q2_extended.transpose() * reduced * Q2_extended;
208 // remove numerical zeros
209 reduced.prune([](const Eigen::Index &row, const Eigen::Index &col, const Scalar &value) {
210 return std::abs(value) > 1e-10;
211 });
212 }
213 }
214 Eigen::MatrixXd NLHomoProblem::reduced_to_disp_grad(const TVector &reduced, bool homogeneous) const
215 {
216 const int dim = state_.mesh->dimension();
217 const int dof2 = macro_reduced_size();
218 const int dof1 = reduced_size();
219
220 return utils::unflatten(macro_reduced_to_full(reduced.tail(dof2), homogeneous), dim);
221 }
222 void NLHomoProblem::hessian(const TVector &x, THessian &hessian)
223 {
225
227
228 for (auto &form : homo_forms)
229 if (form->enabled())
230 {
231 THessian hess_extended;
232 form->second_derivative(reduced_to_extended(x), hess_extended);
233
234 THessian hess;
235 extended_hessian_to_reduced_hessian(hess_extended, hess);
236 hessian += hess;
237 }
238 }
239
240 void NLHomoProblem::set_fixed_entry(const Eigen::VectorXi &fixed_entry)
241 {
242 const int dim = state_.mesh->dimension();
243
244 Eigen::VectorXd fixed_mask;
245 fixed_mask.setZero(dim * dim);
246 fixed_mask(fixed_entry.array()).setOnes();
247 fixed_mask = (macro_full_to_mid_ * fixed_mask).eval();
248
249 fixed_mask_.setZero(fixed_mask.size());
250 for (int i = 0; i < fixed_mask.size(); i++)
251 if (abs(fixed_mask(i)) > 1e-8)
252 fixed_mask_(i) = true;
253
254 const int new_reduced_size = fixed_mask_.size() - fixed_mask_.sum();
255 macro_mid_to_reduced_.setZero(new_reduced_size, fixed_mask_.size());
256 for (int i = 0, j = 0; i < fixed_mask_.size(); i++)
257 if (!fixed_mask_(i))
258 macro_mid_to_reduced_(j++, i) = 1;
259 }
260
262 {
263 const int dim = state_.mesh->dimension();
264 const int dof2 = macro_reduced_size();
265 const int dof1 = reduced_size();
266 const int full_size = hessian.rows();
267
268 Eigen::MatrixXd tmp = constraint_grad();
269 Eigen::MatrixXd A12 = hessian * tmp.transpose();
270 Eigen::MatrixXd A22 = tmp * A12;
271
272 std::vector<Eigen::Triplet<double>> entries;
273 entries.reserve(hessian.nonZeros() + A12.size() * 2 + A22.size());
274
275 for (int k = 0; k < hessian.outerSize(); ++k)
276 for (StiffnessMatrix::InnerIterator it(hessian, k); it; ++it)
277 entries.emplace_back(it.row(), it.col(), it.value());
278
279 for (int i = 0; i < A12.rows(); i++)
280 for (int j = 0; j < A12.cols(); j++)
281 {
282 entries.emplace_back(i, full_size + j, A12(i, j));
283 entries.emplace_back(full_size + j, i, A12(i, j));
284 }
285
286 for (int i = 0; i < A22.rows(); i++)
287 for (int j = 0; j < A22.cols(); j++)
288 entries.emplace_back(i + full_size, j + full_size, A22(i, j));
289
290 hessian.resize(0, 0);
291 hessian.resize(full_size + dof2, full_size + dof2);
292 hessian.setFromTriplets(entries.begin(), entries.end());
293 // NLProblem::full_hessian_to_reduced_hessian(hessian);
294
295 {
296 assert(dof1 == Q2_.cols());
297 StiffnessMatrix Q2_extended(Q2_.rows() + dof2, Q2_.cols() + dof2);
298
299 {
300 entries.clear();
301 for (int k = 0; k < Q2_.cols(); ++k)
302 for (StiffnessMatrix::InnerIterator it(Q2_, k); it; ++it)
303 {
304 entries.emplace_back(it.row(), it.col(), it.value());
305 }
306
307 for (int k = 0; k < dof2; k++)
308 {
309 entries.emplace_back(Q2_.rows() + k, dof1 + k, 1);
310 }
311
312 Q2_extended.setFromTriplets(entries.begin(), entries.end());
313 }
314
315 hessian = Q2_extended.transpose() * hessian * Q2_extended;
316 // remove numerical zeros
317 hessian.prune([](const Eigen::Index &row, const Eigen::Index &col, const Scalar &value) {
318 return std::abs(value) > 1e-10;
319 });
320 }
321 }
322
323 NLHomoProblem::TVector NLHomoProblem::full_to_reduced(const TVector &full) const
324 {
325 log_and_throw_error("Invalid function!");
326 return TVector();
327 }
328
329 NLHomoProblem::TVector NLHomoProblem::full_to_reduced(const TVector &full, const Eigen::MatrixXd &disp_grad) const
330 {
331 const int dim = state_.mesh->dimension();
332 const int dof2 = macro_reduced_size();
333 const int dof1 = reduced_size();
334
335 TVector reduced;
336 reduced.setZero(dof1 + dof2);
337
339 reduced.tail(dof2) = macro_full_to_reduced(utils::flatten(disp_grad));
340
341 return reduced;
342 }
343 NLHomoProblem::TVector NLHomoProblem::full_to_reduced_grad(const TVector &full) const
344 {
345 const int dim = state_.mesh->dimension();
346 const int dof2 = macro_reduced_size();
347 const int dof1 = reduced_size();
348
349 TVector reduced;
350 reduced.setZero(dof1 + dof2);
351 reduced.head(dof1) = NLProblem::full_to_reduced_grad(full);
352 reduced.tail(dof2) = constraint_grad() * full;
353
354 return reduced;
355 }
356 NLHomoProblem::TVector NLHomoProblem::reduced_to_full(const TVector &reduced) const
357 {
358 const int dim = state_.mesh->dimension();
359 const int dof2 = macro_reduced_size();
360 const int dof1 = reduced_size();
361
362 Eigen::MatrixXd disp_grad = utils::unflatten(macro_reduced_to_full(reduced.tail(dof2)), dim);
364 }
365
367 {
368 return macro_mid_to_reduced_.rows();
369 }
370 NLHomoProblem::TVector NLHomoProblem::macro_full_to_reduced(const TVector &full) const
371 {
373 }
374 Eigen::MatrixXd NLHomoProblem::macro_full_to_reduced_grad(const Eigen::MatrixXd &full) const
375 {
376 return macro_mid_to_reduced_ * macro_mid_to_full_.transpose() * full;
377 }
378 NLHomoProblem::TVector NLHomoProblem::macro_reduced_to_full(const TVector &reduced, bool homogeneous) const
379 {
380 TVector mid = macro_mid_to_reduced_.transpose() * reduced;
381 const TVector fixed_values = homogeneous ? TVector::Zero(macro_full_to_mid_.rows()) : TVector(macro_full_to_mid_ * utils::flatten(macro_strain_constraint_.eval(t_)));
382 for (int i = 0; i < fixed_mask_.size(); i++)
383 if (fixed_mask_(i))
384 mid(i) = fixed_values(i);
385
386 return macro_mid_to_full_ * mid;
387 }
388
389 void NLHomoProblem::init(const TVector &x0)
390 {
391 for (auto &form : homo_forms)
392 form->init(reduced_to_extended(x0));
394 }
395
396 bool NLHomoProblem::is_step_valid(const TVector &x0, const TVector &x1)
397 {
398 bool flag = NLProblem::is_step_valid(x0.head(reduced_size()), x1.head(reduced_size()));
399 for (auto &form : homo_forms)
400 if (form->enabled())
401 flag &= form->is_step_valid(reduced_to_extended(x0), reduced_to_extended(x1));
402
403 return flag;
404 }
405 bool NLHomoProblem::is_step_collision_free(const TVector &x0, const TVector &x1)
406 {
407 bool flag = NLProblem::is_step_collision_free(x0.head(reduced_size()), x1.head(reduced_size()));
408 for (auto &form : homo_forms)
409 if (form->enabled())
410 flag &= form->is_step_collision_free(reduced_to_extended(x0), reduced_to_extended(x1));
411
412 return flag;
413 }
414 double NLHomoProblem::max_step_size(const TVector &x0, const TVector &x1)
415 {
416 double size = NLProblem::max_step_size(x0.head(reduced_size()), x1.head(reduced_size()));
417 for (auto &form : homo_forms)
418 if (form->enabled())
419 size = std::min(size, form->max_step_size(reduced_to_extended(x0), reduced_to_extended(x1)));
420 return size;
421 }
422
423 void NLHomoProblem::line_search_begin(const TVector &x0, const TVector &x1)
424 {
426 for (auto &form : homo_forms)
427 form->line_search_begin(reduced_to_extended(x0), reduced_to_extended(x1));
428 }
429 void NLHomoProblem::post_step(const polysolve::nonlinear::PostStepData &data)
430 {
431 NLProblem::post_step(polysolve::nonlinear::PostStepData(data.iter_num, data.solver_info, reduced_to_full(data.x.head(reduced_size())), reduced_to_full(data.grad.head(reduced_size()))));
432 for (auto &form : homo_forms)
433 form->post_step(polysolve::nonlinear::PostStepData(
434 data.iter_num, data.solver_info, reduced_to_extended(data.x), reduced_to_extended(data.grad)));
435 }
436
437 void NLHomoProblem::solution_changed(const TVector &new_x)
438 {
440 for (auto &form : homo_forms)
441 form->solution_changed(reduced_to_extended(new_x));
442 }
443
444 void NLHomoProblem::init_lagging(const TVector &x)
445 {
447 for (auto &form : homo_forms)
448 form->init_lagging(reduced_to_extended(x));
449 }
450 void NLHomoProblem::update_lagging(const TVector &x, const int iter_num)
451 {
452 NLProblem::update_lagging(x.head(reduced_size()), iter_num);
453 for (auto &form : homo_forms)
454 form->update_lagging(reduced_to_extended(x), iter_num);
455 }
456
457 void NLHomoProblem::update_quantities(const double t, const TVector &x)
458 {
460 for (auto &form : homo_forms)
461 form->update_quantities(t, reduced_to_extended(x));
462 }
463
464 Eigen::MatrixXd NLHomoProblem::constraint_grad() const
465 {
466 const int dim = state_.mesh->dimension();
467 Eigen::MatrixXd jac; // (dim*dim) x (dim*n_bases)
468
470
471 jac.setZero(dim * dim, full_size_);
472 for (int i = 0; i < X.rows(); i++)
473 for (int j = 0; j < dim; j++)
474 for (int k = 0; k < dim; k++)
475 jac(j * dim + k, i * dim + j) = X(i, k);
476
477 return macro_full_to_reduced_grad(jac);
478 }
479} // namespace polyfem::solver
double val
Definition Assembler.cpp:86
std::vector< Eigen::Triplet< double > > entries
int x
main class that contains the polyfem solver and all its state
Definition State.hpp:79
int n_bases
number of bases
Definition State.hpp:178
std::unique_ptr< mesh::Mesh > mesh
current mesh, it can be a Mesh2D or Mesh3D
Definition State.hpp:471
std::shared_ptr< polyfem::mesh::MeshNodes > mesh_nodes
Mapping from input nodes to FE nodes.
Definition State.hpp:193
StiffnessMatrix basis_nodes_to_gbasis_nodes
Definition State.hpp:706
Eigen::MatrixXd eval(const double t) const
static Eigen::MatrixXd generate_linear_field(const int n_bases, const std::shared_ptr< mesh::MeshNodes > mesh_nodes, const Eigen::MatrixXd &grad)
static Eigen::MatrixXd get_bases_position(const int n_bases, const std::shared_ptr< mesh::MeshNodes > mesh_nodes)
virtual void hessian(const TVector &x, THessian &hessian) override
virtual double value(const TVector &x) override
virtual void init(const TVector &x0) override
virtual void gradient(const TVector &x, TVector &gradv) override
Eigen::MatrixXd reduced_to_disp_grad(const TVector &reduced, bool homogeneous=false) const
void gradient(const TVector &x, TVector &gradv) override
void init(const TVector &x0) override
void post_step(const polysolve::nonlinear::PostStepData &data) override
void set_fixed_entry(const Eigen::VectorXi &fixed_entry)
TVector macro_full_to_reduced(const TVector &full) const
bool is_step_valid(const TVector &x0, const TVector &x1) override
const assembler::MacroStrainValue & macro_strain_constraint_
void line_search_begin(const TVector &x0, const TVector &x1) override
void update_quantities(const double t, const TVector &x) override
TVector extended_to_reduced_grad(const TVector &extended) const
NLHomoProblem(const int full_size, const assembler::MacroStrainValue &macro_strain_constraint, const State &state, const double t, const std::vector< std::shared_ptr< Form > > &forms, const std::vector< std::shared_ptr< AugmentedLagrangianForm > > &penalty_forms, const bool solve_symmetric_macro_strain, const std::shared_ptr< polysolve::linear::Solver > &solver)
void hessian(const TVector &x, THessian &hessian) override
void init_lagging(const TVector &x) override
TVector reduced_to_extended(const TVector &reduced, bool homogeneous=false) const
TVector reduced_to_full(const TVector &reduced) const
TVector extended_to_reduced(const TVector &extended) const
std::vector< std::shared_ptr< Form > > homo_forms
Eigen::MatrixXd constraint_grad() const
TVector reduced_to_full_shape_derivative(const Eigen::MatrixXd &disp_grad, const TVector &adjoint_full) const
Eigen::MatrixXd macro_full_to_reduced_grad(const Eigen::MatrixXd &full) const
double value(const TVector &x) override
TVector full_to_reduced(const TVector &full, const Eigen::MatrixXd &disp_grad) const
TVector full_to_reduced_grad(const TVector &full) const override
void full_hessian_to_reduced_hessian(THessian &hessian) const
void solution_changed(const TVector &new_x) override
double max_step_size(const TVector &x0, const TVector &x1) override
void extended_hessian_to_reduced_hessian(const THessian &extended, THessian &reduced) const
bool is_step_collision_free(const TVector &x0, const TVector &x1) override
void update_lagging(const TVector &x, const int iter_num) override
TVector macro_reduced_to_full(const TVector &reduced, bool homogeneous=false) const
const int full_size_
Size of the full problem.
Definition NLProblem.hpp:71
void line_search_begin(const TVector &x0, const TVector &x1) override
StiffnessMatrix Q2_
Q2 block of the QR decomposition of the constraints matrix.
Definition NLProblem.hpp:91
virtual bool is_step_valid(const TVector &x0, const TVector &x1) override
virtual void post_step(const polysolve::nonlinear::PostStepData &data) override
virtual TVector full_to_reduced_grad(const TVector &full) const
TVector full_to_reduced(const TVector &full) const
virtual void update_quantities(const double t, const TVector &x)
void init_lagging(const TVector &x) override
virtual bool is_step_collision_free(const TVector &x0, const TVector &x1) override
TVector reduced_to_full(const TVector &reduced) const
void update_lagging(const TVector &x, const int iter_num) override
virtual double max_step_size(const TVector &x0, const TVector &x1) override
void solution_changed(const TVector &new_x) override
std::shared_ptr< FullNLProblem > penalty_problem_
Definition NLProblem.hpp:98
Eigen::MatrixXd unflatten(const Eigen::VectorXd &x, int dim)
Unflatten rowwises, so every dim elements in x become a row.
Eigen::VectorXd flatten(const Eigen::MatrixXd &X)
Flatten rowwises.
void log_and_throw_error(const std::string &msg)
Definition Logger.cpp:73
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix
Definition Types.hpp:22