Loading [MathJax]/extensions/tex2jax.js
PolyFEM
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Jacobian.cpp
Go to the documentation of this file.
1#include <numeric>
2#ifdef POLYFEM_WITH_BEZIER
3#include <element_validity.hpp>
4using namespace element_validity;
5#endif
6#include <polyfem/utils/Logger.hpp>
7#include <polyfem/utils/par_for.hpp>
8#include "Jacobian.hpp"
9#include <polyfem/autogen/auto_p_bases.hpp>
10#include <polyfem/io/Evaluator.hpp>
11
12using namespace polyfem::assembler;
13
14namespace polyfem::utils
15{
16 namespace
17 {
18 template <int n, int s, int p>
19 bool check_static(
20 const int n_threads,
21 const Eigen::MatrixXd &cp,
22 int &invalid_id,
23 Tree &tree)
24 {
25#ifdef POLYFEM_WITH_BEZIER
26 std::vector<int> hierarchy;
27 StaticValidator<n, s, p> check(n_threads);
28 const auto flag_ = check.isValid(cp, &hierarchy, &invalid_id);
29 const bool flag = flag_ == Validity::valid;
30 if (!flag)
31 {
32 Tree *dst = &tree;
33 for (const auto i : hierarchy)
34 {
35 dst->add_children(1 << n);
36 dst = &(dst->child(i));
37 }
38 }
39 return flag;
40#else
41 log_and_throw_error("Enable Bezier library to allow robust Jacobian check!");
42 return false;
43#endif
44 }
45
46 template <int n, int s, int p>
47 bool check_transient(
48 const int n_threads,
49 const Eigen::MatrixXd &cp1,
50 const Eigen::MatrixXd &cp2,
51 int &invalid_id)
52 {
53#ifdef POLYFEM_WITH_BEZIER
54 std::vector<int> hierarchy;
55 ContinuousValidator<n, s, p> check(n_threads);
56 check.setPrecisionTarget(1);
57 const auto flag_ = check.maxTimeStep(cp1, cp2, &hierarchy, &invalid_id);
58 return flag_ == 1.;
59#else
60 log_and_throw_error("Enable Bezier library to allow robust Jacobian check!");
61 return false;
62#endif
63 }
64
65 template <int n, int s, int p>
66 void check_transient(
67 const int n_threads,
68 const Eigen::MatrixXd &cp1,
69 const Eigen::MatrixXd &cp2,
70 double &step,
71 int &invalid_id,
72 double &invalid_step,
73 bool &gaveUp,
74 Tree &tree)
75 {
76#ifdef POLYFEM_WITH_BEZIER
77 std::vector<int> hierarchy;
78 ContinuousValidator<n, s, p> check(n_threads);
79 typename ContinuousValidator<n, s, p>::Info info;
80 step = check.maxTimeStep(
81 cp1, cp2, &hierarchy, &invalid_id, &invalid_step, &info);
82 gaveUp = !info.success();
83 if (step < 1)
84 {
85 Tree *dst = &tree;
86 for (const auto i : hierarchy)
87 {
88 dst->add_children(1 << n);
89 dst = &(dst->child(i));
90 }
91 }
92#else
93 log_and_throw_error("Enable Bezier library to allow robust Jacobian check!");
94 return;
95#endif
96 }
97 } // namespace
98 Eigen::MatrixXd extract_nodes(const int dim, const std::vector<basis::ElementBases> &bases, const std::vector<basis::ElementBases> &gbases, const Eigen::VectorXd &u, int order, int n_elem)
99 {
100 if (n_elem < 0)
101 n_elem = bases.size();
102 Eigen::MatrixXd local_pts;
103 if (dim == 3)
104 autogen::p_nodes_3d(order, local_pts);
105 else
106 autogen::p_nodes_2d(order, local_pts);
107 const int n_basis_per_cell = local_pts.rows();
108 Eigen::MatrixXd cp = Eigen::MatrixXd::Zero(n_elem * n_basis_per_cell, dim);
109 for (int e = 0; e < n_elem; ++e)
110 {
111 ElementAssemblyValues vals;
112 vals.compute(e, dim == 3, local_pts, bases[e], gbases[e]);
113
114 for (std::size_t j = 0; j < vals.basis_values.size(); ++j)
115 for (const auto &g : vals.basis_values[j].global)
116 cp.middleRows(e * n_basis_per_cell, n_basis_per_cell) += g.val * vals.basis_values[j].val * u.segment(g.index * dim, dim).transpose();
117
118 Eigen::MatrixXd mapped;
119 gbases[e].eval_geom_mapping(local_pts, mapped);
120 cp.middleRows(e * n_basis_per_cell, n_basis_per_cell) += mapped;
121 }
122 return cp;
123 }
124
125 Eigen::MatrixXd extract_nodes(const int dim, const basis::ElementBases &basis, const basis::ElementBases &gbasis, const Eigen::VectorXd &u, int order)
126 {
127 Eigen::MatrixXd local_pts;
128 if (dim == 3)
129 autogen::p_nodes_3d(order, local_pts);
130 else
131 autogen::p_nodes_2d(order, local_pts);
132
133 Eigen::MatrixXd cp;
134 gbasis.eval_geom_mapping(local_pts, cp);
135
136 ElementAssemblyValues vals;
137 vals.compute(0, dim == 3, local_pts, basis, gbasis);
138 for (std::size_t j = 0; j < vals.basis_values.size(); ++j)
139 for (const auto &g : vals.basis_values[j].global)
140 cp += g.val * vals.basis_values[j].val * u.segment(g.index * dim, dim).transpose();
141
142 return cp;
143 }
144
145 Eigen::VectorXd robust_evaluate_jacobian(
146 const int order,
147 const Eigen::MatrixXd &cp,
148 const Eigen::MatrixXd &uv)
149 {
150#ifdef POLYFEM_WITH_BEZIER
151 if (cp.cols() == 2)
152 {
153 switch (order)
154 {
155 case 1:
156 {
157 JacobianEvaluator<2, 2, 1> evaluator(cp);
158 return evaluator.eval(uv, 0);
159 }
160 case 2:
161 {
162 JacobianEvaluator<2, 2, 2> evaluator(cp);
163 return evaluator.eval(uv, 0);
164 }
165 case 3:
166 {
167 JacobianEvaluator<2, 2, 3> evaluator(cp);
168 return evaluator.eval(uv, 0);
169 }
170 case 4:
171 {
172 JacobianEvaluator<2, 2, 4> evaluator(cp);
173 return evaluator.eval(uv, 0);
174 }
175 default:
176 throw std::invalid_argument("Order not supported");
177 }
178 }
179 else
180 {
181 switch (order)
182 {
183 case 1:
184 {
185 JacobianEvaluator<3, 3, 1> evaluator(cp);
186 return evaluator.eval(uv, 0);
187 }
188 case 2:
189 {
190 JacobianEvaluator<3, 3, 2> evaluator(cp);
191 return evaluator.eval(uv, 0);
192 }
193 case 3:
194 {
195 JacobianEvaluator<3, 3, 3> evaluator(cp);
196 return evaluator.eval(uv, 0);
197 }
198 // case 4: {
199 // JacobianEvaluator<3, 3, 4> evaluator(cp);
200 // return evaluator.eval(uv, 0);
201 // }
202 default:
203 throw std::invalid_argument("Order not supported");
204 }
205 }
206#else
207 log_and_throw_error("Enable Bezier library to allow robust Jacobian evaluation!");
208 return Eigen::VectorXd::Zero(uv.rows());
209#endif
210 }
211
212 std::vector<int> count_invalid(
213 const int dim,
214 const std::vector<basis::ElementBases> &bases,
215 const std::vector<basis::ElementBases> &gbases,
216 const Eigen::VectorXd &u)
217 {
218 const int order = std::max(bases[0].bases.front().order(), gbases[0].bases.front().order());
219 const int n_basis_per_cell = std::max(bases[0].bases.size(), gbases[0].bases.size());
220 const Eigen::MatrixXd cp = extract_nodes(dim, bases, gbases, u, order);
221
222 std::vector<int> invalidList;
223 const int n_threads = utils::NThread::get().num_threads();
224
225#ifdef POLYFEM_WITH_BEZIER
226 if (cp.cols() == 2)
227 {
228 switch (order)
229 {
230 case 1:
231 {
232 StaticValidator<2, 2, 1> check(n_threads);
233 check.isValid(cp, nullptr, nullptr, &invalidList);
234 break;
235 }
236 case 2:
237 {
238 StaticValidator<2, 2, 2> check(n_threads);
239 check.isValid(cp, nullptr, nullptr, &invalidList);
240 break;
241 }
242 case 3:
243 {
244 StaticValidator<2, 2, 3> check(n_threads);
245 check.isValid(cp, nullptr, nullptr, &invalidList);
246 break;
247 }
248 case 4:
249 {
250 StaticValidator<2, 2, 4> check(n_threads);
251 check.isValid(cp, nullptr, nullptr, &invalidList);
252 break;
253 }
254 default:
255 throw std::invalid_argument("Order not supported");
256 }
257 }
258 else
259 {
260 switch (order)
261 {
262 case 1:
263 {
264 StaticValidator<3, 3, 1> check(n_threads);
265 check.isValid(cp, nullptr, nullptr, &invalidList);
266 break;
267 }
268 case 2:
269 {
270 StaticValidator<3, 3, 2> check(n_threads);
271 check.isValid(cp, nullptr, nullptr, &invalidList);
272 break;
273 }
274 case 3:
275 {
276 StaticValidator<3, 3, 3> check(n_threads);
277 check.isValid(cp, nullptr, nullptr, &invalidList);
278 break;
279 }
280 // case 4: {
281 // StaticValidator<3, 3, 4> check(n_threads);
282 // check.isValid(cp, nullptr, nullptr, &invalidList);
283 // break;
284 // }
285 default:
286 throw std::invalid_argument("Order not supported");
287 }
288 }
289#else
290 log_and_throw_error("Enable Bezier library to allow robust Jacobian check!");
291#endif
292
293 return invalidList;
294 }
295
296 std::tuple<bool, int, Tree>
297 is_valid(
298 const int dim,
299 const std::vector<basis::ElementBases> &bases,
300 const std::vector<basis::ElementBases> &gbases,
301 const Eigen::VectorXd &u,
302 const double threshold)
303 {
304 const int order = std::max(bases[0].bases.front().order(), gbases[0].bases.front().order());
305 const int n_basis_per_cell = std::max(bases[0].bases.size(), gbases[0].bases.size());
306 Eigen::MatrixXd cp = extract_nodes(dim, bases, gbases, u, order);
307
308 bool flag = false;
309 int invalid_id = 0;
310 Tree tree;
311
312 const int n_threads = utils::NThread::get().num_threads();
313
314 if (dim == 2)
315 {
316 switch (order)
317 {
318 case 1:
319 flag = check_static<2, 2, 1>(n_threads, cp, invalid_id, tree);
320 break;
321 case 2:
322 flag = check_static<2, 2, 2>(n_threads, cp, invalid_id, tree);
323 break;
324 case 3:
325 flag = check_static<2, 2, 3>(n_threads, cp, invalid_id, tree);
326 break;
327 case 4:
328 flag = check_static<2, 2, 4>(n_threads, cp, invalid_id, tree);
329 break;
330 default:
331 throw std::invalid_argument("Order not supported");
332 }
333 }
334 else
335 {
336 switch (order)
337 {
338 case 1:
339 flag = check_static<3, 3, 1>(n_threads, cp, invalid_id, tree);
340 break;
341 case 2:
342 flag = check_static<3, 3, 2>(n_threads, cp, invalid_id, tree);
343 break;
344 case 3:
345 flag = check_static<3, 3, 3>(n_threads, cp, invalid_id, tree);
346 break;
347 // case 4:
348 // flag = check_static<3, 3, 4>(n_threads, cp, invalid_id, tree);
349 // break;
350 default:
351 throw std::invalid_argument("Order not supported");
352 }
353 }
354
355 return {flag, invalid_id, tree};
356 }
357
358 bool is_valid(
359 const int dim,
360 const std::vector<basis::ElementBases> &bases,
361 const std::vector<basis::ElementBases> &gbases,
362 const Eigen::VectorXd &u1,
363 const Eigen::VectorXd &u2,
364 const double threshold)
365 {
366 const int order = std::max(bases[0].bases.front().order(), gbases[0].bases.front().order());
367 const int n_basis_per_cell = std::max(bases[0].bases.size(), gbases[0].bases.size());
368
369 const Eigen::MatrixXd cp1 = extract_nodes(dim, bases, gbases, u1, order);
370 const Eigen::MatrixXd cp2 = extract_nodes(dim, bases, gbases, u2, order);
371
372 int invalid_id = 0;
373 const int n_threads = utils::NThread::get().num_threads();
374
375 if (dim == 2)
376 {
377 switch (order)
378 {
379 case 1:
380 return check_transient<2, 2, 1>(n_threads, cp1, cp2, invalid_id);
381 case 2:
382 return check_transient<2, 2, 2>(n_threads, cp1, cp2, invalid_id);
383 case 3:
384 return check_transient<2, 2, 3>(n_threads, cp1, cp2, invalid_id);
385 case 4:
386 return check_transient<2, 2, 4>(n_threads, cp1, cp2, invalid_id);
387 default:
388 throw std::invalid_argument("Order not supported");
389 }
390 }
391 else
392 {
393 switch (order)
394 {
395 case 1:
396 return check_transient<3, 3, 1>(n_threads, cp1, cp2, invalid_id);
397 case 2:
398 return check_transient<3, 3, 2>(n_threads, cp1, cp2, invalid_id);
399 case 3:
400 return check_transient<3, 3, 3>(n_threads, cp1, cp2, invalid_id);
401 // case 4:
402 // return check_transient<3, 3, 4>(n_threads, cp1, cp2, invalid_id);
403 default:
404 throw std::invalid_argument("Order not supported");
405 }
406 }
407 }
408
409 std::tuple<double, int, double, Tree> max_time_step(
410 const int dim,
411 const std::vector<basis::ElementBases> &bases,
412 const std::vector<basis::ElementBases> &gbases,
413 const Eigen::VectorXd &u1,
414 const Eigen::VectorXd &u2,
415 double precision)
416 {
417 const int order = std::max(bases[0].bases.front().order(), gbases[0].bases.front().order());
418 const int n_basis_per_cell = std::max(bases[0].bases.size(), gbases[0].bases.size());
419 Eigen::MatrixXd cp1 = extract_nodes(dim, bases, gbases, u1, order);
420 Eigen::MatrixXd cp2 = extract_nodes(dim, bases, gbases, u2, order);
421
422 // logger().debug("Jacobian check order {}, number of nodes per cell {}, number of total nodes {}", order, n_basis_per_cell, cp2.rows());
423
424 int invalid_id = -1;
425 bool gaveUp = false;
426 double step = 1;
427 double invalid_step = 1.;
428 Tree tree;
429
430 const int n_threads = utils::NThread::get().num_threads();
431
432 if (dim == 2)
433 {
434 switch (order)
435 {
436 case 1:
437 check_transient<2, 2, 1>(n_threads, cp1, cp2, step, invalid_id, invalid_step, gaveUp, tree);
438 break;
439 case 2:
440 check_transient<2, 2, 2>(n_threads, cp1, cp2, step, invalid_id, invalid_step, gaveUp, tree);
441 break;
442 case 3:
443 check_transient<2, 2, 3>(n_threads, cp1, cp2, step, invalid_id, invalid_step, gaveUp, tree);
444 break;
445 case 4:
446 check_transient<2, 2, 4>(n_threads, cp1, cp2, step, invalid_id, invalid_step, gaveUp, tree);
447 break;
448 default:
449 throw std::invalid_argument("Order not supported");
450 }
451 }
452 else
453 {
454 switch (order)
455 {
456 case 1:
457 check_transient<3, 3, 1>(n_threads, cp1, cp2, step, invalid_id, invalid_step, gaveUp, tree);
458 break;
459 case 2:
460 check_transient<3, 3, 2>(n_threads, cp1, cp2, step, invalid_id, invalid_step, gaveUp, tree);
461 break;
462 case 3:
463 check_transient<3, 3, 3>(n_threads, cp1, cp2, step, invalid_id, invalid_step, gaveUp, tree);
464 break;
465 // case 4:
466 // check_transient<3, 3, 4>(n_threads, cp1, cp2, step, invalid_id, invalid_step, gaveUp, tree);
467 // break;
468 default:
469 throw std::invalid_argument("Order not supported");
470 }
471 }
472
473 // if (gaveUp)
474 // logger().warn("Jacobian check gave up!");
475
476 return {step, invalid_id, invalid_step, tree};
477 }
478} // namespace polyfem::utils
vals
ElementAssemblyValues vals
Definition Assembler.cpp:22
auto_p_bases.hpp
polyfem::assembler::ElementAssemblyValues
stores per element basis values at given quadrature points and geometric mapping
Definition ElementAssemblyValues.hpp:14
polyfem::assembler::ElementAssemblyValues::compute
void compute(const int el_index, const bool is_volume, const Eigen::MatrixXd &pts, const basis::ElementBases &basis, const basis::ElementBases &gbasis)
computes the per element values at the local (ref el) points (pts) sets basis_values,...
Definition ElementAssemblyValues.cpp:164
polyfem::basis::ElementBases
Stores the basis functions for a given element in a mesh (facet in 2d, cell in 3d).
Definition ElementBases.hpp:17
polyfem::basis::ElementBases::eval_geom_mapping
void eval_geom_mapping(const Eigen::MatrixXd &samples, Eigen::MatrixXd &mapped) const
Map the sample positions in the parametric domain to the object domain (if the element has no paramet...
Definition ElementBases.cpp:19
polyfem::utils::NThread::get
static NThread & get()
Definition par_for.hpp:19
polyfem::utils::NThread::num_threads
size_t num_threads() const
Definition par_for.hpp:25
polyfem::utils::Tree::add_children
void add_children(int n)
Definition Jacobian.hpp:84
polyfem::utils::Tree::child
Tree & child(int i)
Definition Jacobian.hpp:82
polyfem::assembler
Used for test only.
Definition AMIPSEnergy.cpp:6
polyfem::autogen::p_nodes_2d
void p_nodes_2d(const int p, Eigen::MatrixXd &val)
Definition auto_p_bases.cpp:275
polyfem::autogen::p_nodes_3d
void p_nodes_3d(const int p, Eigen::MatrixXd &val)
Definition auto_p_bases.cpp:830
polyfem::utils::is_valid
std::tuple< bool, int, Tree > is_valid(const int dim, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const Eigen::VectorXd &u, const double threshold)
Definition Jacobian.cpp:297
polyfem::utils::robust_evaluate_jacobian
Eigen::VectorXd robust_evaluate_jacobian(const int order, const Eigen::MatrixXd &cp, const Eigen::MatrixXd &uv)
Definition Jacobian.cpp:145
polyfem::utils::max_time_step
std::tuple< double, int, double, Tree > max_time_step(const int dim, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const Eigen::VectorXd &u1, const Eigen::VectorXd &u2, double precision)
Definition Jacobian.cpp:409
polyfem::utils::count_invalid
std::vector< int > count_invalid(const int dim, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const Eigen::VectorXd &u)
Definition Jacobian.cpp:212
polyfem::utils::extract_nodes
Eigen::MatrixXd extract_nodes(const int dim, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const Eigen::VectorXd &u, int order, int n_elem)
Definition Jacobian.cpp:98
polyfem::log_and_throw_error
void log_and_throw_error(const std::string &msg)
Definition Logger.cpp:73
par_for.hpp