PolyFEM
Loading...
Searching...
No Matches
Jacobian.cpp
Go to the documentation of this file.
1#include <numeric>
2#include <element_validity.hpp>
3#include <polyfem/utils/Logger.hpp>
4#include <polyfem/utils/par_for.hpp>
5#include "Jacobian.hpp"
6#include <polyfem/autogen/auto_p_bases.hpp>
7#include <polyfem/io/Evaluator.hpp>
8
9using namespace element_validity;
10using namespace polyfem::assembler;
11
12namespace polyfem::utils
13{
14 namespace {
15 template <int n, int s, int p>
16 bool check_static(
17 const int n_threads,
18 const Eigen::MatrixXd &cp,
19 int &invalid_id,
20 Tree &tree)
21 {
22 std::vector<int> hierarchy;
23 StaticValidator<n,s,p> check(n_threads);
24 const auto flag_ = check.isValid(cp, &hierarchy, &invalid_id);
25 const bool flag = flag_ == Validity::valid;
26 if (!flag) {
27 Tree *dst = &tree;
28 for (const auto i : hierarchy) {
29 dst->add_children(1<<n);
30 dst = &(dst->child(i));
31 }
32 }
33 return flag;
34 }
35
36 template <int n, int s, int p>
37 bool check_transient(
38 const int n_threads,
39 const Eigen::MatrixXd &cp1,
40 const Eigen::MatrixXd &cp2,
41 int &invalid_id)
42 {
43 std::vector<int> hierarchy;
44 ContinuousValidator<n,s,p> check(n_threads);
45 check.setPrecisionTarget(1);
46 const auto flag_ = check.maxTimeStep(cp1, cp2, &hierarchy, &invalid_id);
47 return flag_ == 1.;
48 }
49
50 template <int n, int s, int p>
51 void check_transient(
52 const int n_threads,
53 const Eigen::MatrixXd &cp1,
54 const Eigen::MatrixXd &cp2,
55 double &step,
56 int &invalid_id,
57 double &invalid_step,
58 bool &gaveUp,
59 Tree &tree)
60 {
61 std::vector<int> hierarchy;
62 ContinuousValidator<n,s,p> check(n_threads);
63 typename ContinuousValidator<n,s,p>::Info info;
64 step = check.maxTimeStep(
65 cp1, cp2, &hierarchy, &invalid_id, &invalid_step, &info
66 );
67 gaveUp = !info.success();
68 if (step < 1) {
69 Tree *dst = &tree;
70 for (const auto i : hierarchy) {
71 dst->add_children(1<<n);
72 dst = &(dst->child(i));
73 }
74 }
75 }
76 }
77 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)
78 {
79 if (n_elem < 0)
80 n_elem = bases.size();
81 Eigen::MatrixXd local_pts;
82 if (dim == 3)
83 autogen::p_nodes_3d(order, local_pts);
84 else
85 autogen::p_nodes_2d(order, local_pts);
86 const int n_basis_per_cell = local_pts.rows();
87 Eigen::MatrixXd cp = Eigen::MatrixXd::Zero(n_elem * n_basis_per_cell, dim);
88 for (int e = 0; e < n_elem; ++e)
89 {
90 ElementAssemblyValues vals;
91 vals.compute(e, dim == 3, local_pts, bases[e], gbases[e]);
92
93 for (std::size_t j = 0; j < vals.basis_values.size(); ++j)
94 for (const auto &g : vals.basis_values[j].global)
95 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();
96
97 Eigen::MatrixXd mapped;
98 gbases[e].eval_geom_mapping(local_pts, mapped);
99 cp.middleRows(e * n_basis_per_cell, n_basis_per_cell) += mapped;
100 }
101 return cp;
102 }
103
104 Eigen::MatrixXd extract_nodes(const int dim, const basis::ElementBases &basis, const basis::ElementBases &gbasis, const Eigen::VectorXd &u, int order)
105 {
106 Eigen::MatrixXd local_pts;
107 if (dim == 3)
108 autogen::p_nodes_3d(order, local_pts);
109 else
110 autogen::p_nodes_2d(order, local_pts);
111
112 Eigen::MatrixXd cp;
113 gbasis.eval_geom_mapping(local_pts, cp);
114
115 ElementAssemblyValues vals;
116 vals.compute(0, dim == 3, local_pts, basis, gbasis);
117 for (std::size_t j = 0; j < vals.basis_values.size(); ++j)
118 for (const auto &g : vals.basis_values[j].global)
119 cp += g.val * vals.basis_values[j].val * u.segment(g.index * dim, dim).transpose();
120
121 return cp;
122 }
123
124 Eigen::VectorXd robust_evaluate_jacobian(
125 const int order,
126 const Eigen::MatrixXd &cp,
127 const Eigen::MatrixXd &uv)
128 {
129 if (cp.cols() == 2) {
130 switch (order) {
131 case 1: {
132 JacobianEvaluator<2, 2, 1> evaluator(cp);
133 return evaluator.eval(uv, 0);
134 }
135 case 2: {
136 JacobianEvaluator<2, 2, 2> evaluator(cp);
137 return evaluator.eval(uv, 0);
138 }
139 case 3: {
140 JacobianEvaluator<2, 2, 3> evaluator(cp);
141 return evaluator.eval(uv, 0);
142 }
143 case 4: {
144 JacobianEvaluator<2, 2, 4> evaluator(cp);
145 return evaluator.eval(uv, 0);
146 }
147 default: throw std::invalid_argument("Order not supported");
148 }
149 }
150 else {
151 switch (order) {
152 case 1: {
153 JacobianEvaluator<3, 3, 1> evaluator(cp);
154 return evaluator.eval(uv, 0);
155 }
156 case 2: {
157 JacobianEvaluator<3, 3, 2> evaluator(cp);
158 return evaluator.eval(uv, 0);
159 }
160 case 3: {
161 JacobianEvaluator<3, 3, 3> evaluator(cp);
162 return evaluator.eval(uv, 0);
163 }
164 // case 4: {
165 // JacobianEvaluator<3, 3, 4> evaluator(cp);
166 // return evaluator.eval(uv, 0);
167 // }
168 default: throw std::invalid_argument("Order not supported");
169 }
170 }
171 }
172
173 std::vector<int> count_invalid(
174 const int dim,
175 const std::vector<basis::ElementBases> &bases,
176 const std::vector<basis::ElementBases> &gbases,
177 const Eigen::VectorXd &u)
178 {
179 const int order = std::max(bases[0].bases.front().order(), gbases[0].bases.front().order());
180 const int n_basis_per_cell = std::max(bases[0].bases.size(), gbases[0].bases.size());
181 const Eigen::MatrixXd cp = extract_nodes(dim, bases, gbases, u, order);
182
183 std::vector<int> invalidList;
184 const int n_threads = utils::NThread::get().num_threads();
185
186 if (cp.cols() == 2) {
187 switch (order) {
188 case 1: {
189 StaticValidator<2, 2, 1> check(n_threads);
190 check.isValid(cp, nullptr, nullptr, &invalidList);
191 break;
192 }
193 case 2: {
194 StaticValidator<2, 2, 2> check(n_threads);
195 check.isValid(cp, nullptr, nullptr, &invalidList);
196 break;
197 }
198 case 3: {
199 StaticValidator<2, 2, 3> check(n_threads);
200 check.isValid(cp, nullptr, nullptr, &invalidList);
201 break;
202 }
203 case 4: {
204 StaticValidator<2, 2, 4> check(n_threads);
205 check.isValid(cp, nullptr, nullptr, &invalidList);
206 break;
207 }
208 default: throw std::invalid_argument("Order not supported");
209 }
210 }
211 else {
212 switch (order) {
213 case 1: {
214 StaticValidator<3, 3, 1> check(n_threads);
215 check.isValid(cp, nullptr, nullptr, &invalidList);
216 break;
217 }
218 case 2: {
219 StaticValidator<3, 3, 2> check(n_threads);
220 check.isValid(cp, nullptr, nullptr, &invalidList);
221 break;
222 }
223 case 3: {
224 StaticValidator<3, 3, 3> check(n_threads);
225 check.isValid(cp, nullptr, nullptr, &invalidList);
226 break;
227 }
228 // case 4: {
229 // StaticValidator<3, 3, 4> check(n_threads);
230 // check.isValid(cp, nullptr, nullptr, &invalidList);
231 // break;
232 // }
233 default: throw std::invalid_argument("Order not supported");
234 }
235 }
236
237 return invalidList;
238 }
239
240 std::tuple<bool, int, Tree>
241 is_valid(
242 const int dim,
243 const std::vector<basis::ElementBases> &bases,
244 const std::vector<basis::ElementBases> &gbases,
245 const Eigen::VectorXd &u,
246 const double threshold)
247 {
248 const int order = std::max(bases[0].bases.front().order(), gbases[0].bases.front().order());
249 const int n_basis_per_cell = std::max(bases[0].bases.size(), gbases[0].bases.size());
250 Eigen::MatrixXd cp = extract_nodes(dim, bases, gbases, u, order);
251
252 bool flag = false;
253 int invalid_id = 0;
254 Tree tree;
255
256 const int n_threads = utils::NThread::get().num_threads();
257
258 if (dim == 2) {
259 switch (order) {
260 case 1:
261 flag = check_static<2, 2, 1>(n_threads, cp, invalid_id, tree);
262 break;
263 case 2:
264 flag = check_static<2, 2, 2>(n_threads, cp, invalid_id, tree);
265 break;
266 case 3:
267 flag = check_static<2, 2, 3>(n_threads, cp, invalid_id, tree);
268 break;
269 case 4:
270 flag = check_static<2, 2, 4>(n_threads, cp, invalid_id, tree);
271 break;
272 default:
273 throw std::invalid_argument("Order not supported");
274 }
275 }
276 else {
277 switch (order) {
278 case 1:
279 flag = check_static<3, 3, 1>(n_threads, cp, invalid_id, tree);
280 break;
281 case 2:
282 flag = check_static<3, 3, 2>(n_threads, cp, invalid_id, tree);
283 break;
284 case 3:
285 flag = check_static<3, 3, 3>(n_threads, cp, invalid_id, tree);
286 break;
287 // case 4:
288 // flag = check_static<3, 3, 4>(n_threads, cp, invalid_id, tree);
289 // break;
290 default:
291 throw std::invalid_argument("Order not supported");
292 }
293 }
294
295 return {flag, invalid_id, tree};
296 }
297
298 bool is_valid(
299 const int dim,
300 const std::vector<basis::ElementBases> &bases,
301 const std::vector<basis::ElementBases> &gbases,
302 const Eigen::VectorXd &u1,
303 const Eigen::VectorXd &u2,
304 const double threshold)
305 {
306 const int order = std::max(bases[0].bases.front().order(), gbases[0].bases.front().order());
307 const int n_basis_per_cell = std::max(bases[0].bases.size(), gbases[0].bases.size());
308
309 const Eigen::MatrixXd cp1 = extract_nodes(dim, bases, gbases, u1, order);
310 const Eigen::MatrixXd cp2 = extract_nodes(dim, bases, gbases, u2, order);
311
312 int invalid_id = 0;
313 const int n_threads = utils::NThread::get().num_threads();
314
315 if (dim == 2) {
316 switch (order) {
317 case 1:
318 return check_transient<2, 2, 1>(n_threads, cp1, cp2, invalid_id);
319 case 2:
320 return check_transient<2, 2, 2>(n_threads, cp1, cp2, invalid_id);
321 case 3:
322 return check_transient<2, 2, 3>(n_threads, cp1, cp2, invalid_id);
323 case 4:
324 return check_transient<2, 2, 4>(n_threads, cp1, cp2, invalid_id);
325 default:
326 throw std::invalid_argument("Order not supported");
327 }
328 }
329 else {
330 switch (order) {
331 case 1:
332 return check_transient<3, 3, 1>(n_threads, cp1, cp2, invalid_id);
333 case 2:
334 return check_transient<3, 3, 2>(n_threads, cp1, cp2, invalid_id);
335 case 3:
336 return check_transient<3, 3, 3>(n_threads, cp1, cp2, invalid_id);
337 // case 4:
338 // return check_transient<3, 3, 4>(n_threads, cp1, cp2, invalid_id);
339 default:
340 throw std::invalid_argument("Order not supported");
341 }
342 }
343 }
344
345 std::tuple<double, int, double, Tree> max_time_step(
346 const int dim,
347 const std::vector<basis::ElementBases> &bases,
348 const std::vector<basis::ElementBases> &gbases,
349 const Eigen::VectorXd &u1,
350 const Eigen::VectorXd &u2,
351 double precision)
352 {
353 const int order = std::max(bases[0].bases.front().order(), gbases[0].bases.front().order());
354 const int n_basis_per_cell = std::max(bases[0].bases.size(), gbases[0].bases.size());
355 Eigen::MatrixXd cp1 = extract_nodes(dim, bases, gbases, u1, order);
356 Eigen::MatrixXd cp2 = extract_nodes(dim, bases, gbases, u2, order);
357
358 // logger().debug("Jacobian check order {}, number of nodes per cell {}, number of total nodes {}", order, n_basis_per_cell, cp2.rows());
359
360 int invalid_id = -1;
361 bool gaveUp = false;
362 double step = 1;
363 double invalid_step = 1.;
364 Tree tree;
365
366 const int n_threads = utils::NThread::get().num_threads();
367
368 if (dim == 2) {
369 switch (order) {
370 case 1:
371 check_transient<2, 2, 1>(n_threads, cp1, cp2, step, invalid_id, invalid_step, gaveUp, tree);
372 break;
373 case 2:
374 check_transient<2, 2, 2>(n_threads, cp1, cp2, step, invalid_id, invalid_step, gaveUp, tree);
375 break;
376 case 3:
377 check_transient<2, 2, 3>(n_threads, cp1, cp2, step, invalid_id, invalid_step, gaveUp, tree);
378 break;
379 case 4:
380 check_transient<2, 2, 4>(n_threads, cp1, cp2, step, invalid_id, invalid_step, gaveUp, tree);
381 break;
382 default:
383 throw std::invalid_argument("Order not supported");
384 }
385 }
386 else {
387 switch (order) {
388 case 1:
389 check_transient<3, 3, 1>(n_threads, cp1, cp2, step, invalid_id, invalid_step, gaveUp, tree);
390 break;
391 case 2:
392 check_transient<3, 3, 2>(n_threads, cp1, cp2, step, invalid_id, invalid_step, gaveUp, tree);
393 break;
394 case 3:
395 check_transient<3, 3, 3>(n_threads, cp1, cp2, step, invalid_id, invalid_step, gaveUp, tree);
396 break;
397 // case 4:
398 // check_transient<3, 3, 4>(n_threads, cp1, cp2, step, invalid_id, invalid_step, gaveUp, tree);
399 // break;
400 default:
401 throw std::invalid_argument("Order not supported");
402 }
403 }
404
405 // if (gaveUp)
406 // logger().warn("Jacobian check gave up!");
407
408 return {step, invalid_id, invalid_step, tree};
409 }
410}
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:79
polyfem::utils::Tree::child
Tree & child(int i)
Definition Jacobian.hpp:77
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:835
polyfem::autogen::p_nodes_3d
void p_nodes_3d(const int p, Eigen::MatrixXd &val)
Definition auto_p_bases.cpp:2792
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:241
polyfem::utils::robust_evaluate_jacobian
Eigen::VectorXd robust_evaluate_jacobian(const int order, const Eigen::MatrixXd &cp, const Eigen::MatrixXd &uv)
Definition Jacobian.cpp:124
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:345
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:173
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:77
par_for.hpp