PolyFEM
Loading...
Searching...
No Matches
Evaluator.cpp
Go to the documentation of this file.
1#include "Evaluator.hpp"
2
3#include <polyfem/mesh/mesh2D/Mesh2D.hpp>
4#include <polyfem/mesh/mesh3D/Mesh3D.hpp>
5
6#include <polyfem/quadrature/HexQuadrature.hpp>
7#include <polyfem/quadrature/QuadQuadrature.hpp>
8#include <polyfem/quadrature/TetQuadrature.hpp>
9#include <polyfem/quadrature/TriQuadrature.hpp>
10
11#include <polyfem/utils/BoundarySampler.hpp>
12#include <polyfem/utils/Jacobian.hpp>
13
14#include <polyfem/autogen/auto_p_bases.hpp>
15#include <polyfem/autogen/auto_q_bases.hpp>
16
17#include <polyfem/utils/Logger.hpp>
18
19#include <igl/AABB.h>
20#include <igl/per_face_normals.h>
21
22namespace polyfem::io
23{
24 using namespace mesh;
25 using namespace assembler;
26 using namespace basis;
27
28 namespace
29 {
30 void flattened_tensor_coeffs(const Eigen::MatrixXd &S, Eigen::MatrixXd &X)
31 {
32 if (S.cols() == 4)
33 {
34 X.resize(S.rows(), 3);
35 X.col(0) = S.col(0);
36 X.col(1) = S.col(3);
37 X.col(2) = S.col(1);
38 }
39 else if (S.cols() == 9)
40 {
41 // [S11, S22, S33, S12, S13, S23]
42 X.resize(S.rows(), 6);
43 X.col(0) = S.col(0);
44 X.col(1) = S.col(4);
45 X.col(2) = S.col(8);
46 X.col(3) = S.col(1);
47 X.col(4) = S.col(2);
48 X.col(5) = S.col(5);
49 }
50 else
51 {
52 logger().error("Invalid tensor dimensions.");
53 }
54 }
55 } // namespace
56
57 void Evaluator::interpolate_boundary_function(
58 const mesh::Mesh &mesh,
59 const bool is_problem_scalar,
60 const std::vector<basis::ElementBases> &bases,
61 const std::vector<basis::ElementBases> &gbases,
62 const Eigen::MatrixXd &pts,
63 const Eigen::MatrixXi &faces,
64 const Eigen::MatrixXd &fun,
65 const bool compute_avg,
66 Eigen::MatrixXd &result)
67 {
68 if (fun.size() <= 0)
69 {
70 logger().error("Solve the problem first!");
71 return;
72 }
73 assert(mesh.is_volume());
74
75 const Mesh3D &mesh3d = dynamic_cast<const Mesh3D &>(mesh);
76
77 Eigen::MatrixXd points, uv;
78 Eigen::VectorXd weights;
79
80 int actual_dim = 1;
81 if (!is_problem_scalar)
82 actual_dim = 3;
83
84 igl::AABB<Eigen::MatrixXd, 3> tree;
85 tree.init(pts, faces);
86
87 result.resize(faces.rows(), actual_dim);
88 result.setConstant(std::numeric_limits<double>::quiet_NaN());
89
90 int counter = 0;
91
92 for (int e = 0; e < mesh3d.n_elements(); ++e)
93 {
94 const basis::ElementBases &gbs = gbases[e];
95 const basis::ElementBases &bs = bases[e];
96
97 for (int lf = 0; lf < mesh3d.n_cell_faces(e); ++lf)
98 {
99 const int face_id = mesh3d.cell_face(e, lf);
100 if (!mesh3d.is_boundary_face(face_id))
101 continue;
102
103 if (mesh3d.is_simplex(e))
104 utils::BoundarySampler::quadrature_for_tri_face(lf, 4, face_id, mesh3d, uv, points, weights);
105 else if (mesh3d.is_cube(e))
106 utils::BoundarySampler::quadrature_for_quad_face(lf, 4, face_id, mesh3d, uv, points, weights);
107 else
108 assert(false);
109
110 ElementAssemblyValues vals;
111 vals.compute(e, true, points, bs, gbs);
112 RowVectorNd loc_val(actual_dim);
113 loc_val.setZero();
114
115 // UIEvaluator::ui_state().debug_data().add_points(vals.val, Eigen::RowVector3d(1,0,0));
116
117 // const auto nodes = bs.local_nodes_for_primitive(face_id, mesh3d);
118
119 // for(long n = 0; n < nodes.size(); ++n)
120 for (size_t j = 0; j < bs.bases.size(); ++j)
121 {
122 // const auto &b = bs.bases[nodes(n)];
123 // const AssemblyValues &v = vals.basis_values[nodes(n)];
124 const AssemblyValues &v = vals.basis_values[j];
125 for (int d = 0; d < actual_dim; ++d)
126 {
127 for (size_t g = 0; g < v.global.size(); ++g)
128 {
129 loc_val(d) += (v.global[g].val * v.val.array() * fun(v.global[g].index * actual_dim + d) * weights.array()).sum();
130 }
131 }
132 }
133
134 int I;
135 Eigen::RowVector3d C;
136 const Eigen::RowVector3d bary = mesh3d.face_barycenter(face_id);
137
138 const double dist = tree.squared_distance(pts, faces, bary, I, C);
139 assert(dist < 1e-16);
140
141 assert(std::isnan(result(I, 0)));
142 if (compute_avg)
143 result.row(I) = loc_val / weights.sum();
144 else
145 result.row(I) = loc_val;
146 ++counter;
147 }
148 }
149
150 assert(counter == result.rows());
151 }
152
153 void Evaluator::average_grad_based_function(
154 const mesh::Mesh &mesh,
155 const bool is_problem_scalar,
156 const int n_bases,
157 const std::vector<basis::ElementBases> &bases,
158 const std::vector<basis::ElementBases> &gbases,
159 const Eigen::VectorXi &disc_orders,
160 const std::map<int, Eigen::MatrixXd> &polys,
161 const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d,
162 const assembler::Assembler &assembler,
163 const utils::RefElementSampler &sampler,
164 const double t,
165 const int n_points,
166 const Eigen::MatrixXd &fun,
167 std::vector<assembler::Assembler::NamedMatrix> &result_scalar,
168 std::vector<assembler::Assembler::NamedMatrix> &result_tensor,
169 const bool use_sampler,
170 const bool boundary_only)
171 {
172 result_scalar.clear();
173 result_tensor.clear();
174
175 if (fun.size() <= 0)
176 {
177 logger().error("Solve the problem first!");
178 return;
179 }
180 if (is_problem_scalar)
181 {
182 logger().error("Define a tensor problem!");
183 return;
184 }
185
186 assert(!is_problem_scalar);
187 const int actual_dim = mesh.dimension();
188
189 std::vector<Eigen::MatrixXd> avg_scalar;
190
191 Eigen::MatrixXd areas(n_bases, 1);
192 areas.setZero();
193
194 std::vector<std::pair<std::string, Eigen::MatrixXd>> tmp_s;
195 Eigen::MatrixXd local_val;
196
197 ElementAssemblyValues vals;
198 for (int i = 0; i < int(bases.size()); ++i)
199 {
200 const ElementBases &bs = bases[i];
201 const ElementBases &gbs = gbases[i];
202 Eigen::MatrixXd local_pts;
203
204 if (mesh.is_simplex(i))
205 {
206 if (mesh.dimension() == 3)
207 autogen::p_nodes_3d(disc_orders(i), local_pts);
208 else
209 autogen::p_nodes_2d(disc_orders(i), local_pts);
210 }
211 else if (mesh.is_cube(i))
212 {
213 if (mesh.dimension() == 3)
214 autogen::q_nodes_3d(disc_orders(i), local_pts);
215 else
216 autogen::q_nodes_2d(disc_orders(i), local_pts);
217 }
218 else
219 {
220 // not supported for polys
221 continue;
222 }
223
224 vals.compute(i, actual_dim == 3, bases[i], gbases[i]);
225 const quadrature::Quadrature &quadrature = vals.quadrature;
226 const double area = (vals.det.array() * quadrature.weights.array()).sum();
227
228 assembler.compute_scalar_value(OutputData(t, i, bs, gbs, local_pts, fun), tmp_s);
229
230 // assembler.compute_tensor_value(i, bs, gbs, local_pts, fun, local_val);
231 // MatrixXd avg_tensor(n_points * actual_dim*actual_dim, 1);
232
233 for (size_t j = 0; j < bs.bases.size(); ++j)
234 {
235 const Basis &b = bs.bases[j];
236 if (b.global().size() > 1)
237 continue;
238
239 auto &global = b.global().front();
240 areas(global.index) += area;
241 }
242
243 if (avg_scalar.empty())
244 {
245 avg_scalar.resize(tmp_s.size());
246 for (auto &m : avg_scalar)
247 {
248 m.resize(n_bases, 1);
249 m.setZero();
250 }
251 }
252
253 for (int k = 0; k < tmp_s.size(); ++k)
254 {
255 local_val = tmp_s[k].second;
256
257 for (size_t j = 0; j < bs.bases.size(); ++j)
258 {
259 const Basis &b = bs.bases[j];
260 if (b.global().size() > 1)
261 continue;
262
263 auto &global = b.global().front();
264 avg_scalar[k](global.index) += local_val(j) * area;
265 }
266 }
267 }
268
269 for (auto &m : avg_scalar)
270 {
271 m.array() /= areas.array();
272 }
273
274 result_scalar.resize(tmp_s.size());
275 for (int k = 0; k < tmp_s.size(); ++k)
276 {
277 result_scalar[k].first = tmp_s[k].first;
278 interpolate_function(mesh, 1, bases, disc_orders, polys, polys_3d, sampler, n_points,
279 avg_scalar[k], result_scalar[k].second, use_sampler, boundary_only);
280 }
281 // interpolate_function(n_points, actual_dim*actual_dim, bases, avg_tensor, result_tensor, boundary_only);
282 }
283
284 void Evaluator::compute_stress_at_quadrature_points(
285 const mesh::Mesh &mesh,
286 const bool is_problem_scalar,
287 const std::vector<basis::ElementBases> &bases,
288 const std::vector<basis::ElementBases> &gbases,
289 const Eigen::VectorXi &disc_orders,
290 const assembler::Assembler &assembler,
291 const Eigen::MatrixXd &fun,
292 const double t,
293 Eigen::MatrixXd &result,
294 Eigen::VectorXd &von_mises)
295 {
296 // if (!mesh)
297 // {
298 // logger().error("Load the mesh first!");
299 // return;
300 // }
301 if (fun.size() <= 0)
302 {
303 logger().error("Solve the problem first!");
304 return;
305 }
306 if (is_problem_scalar)
307 {
308 logger().error("Define a tensor problem!");
309 return;
310 }
311
312 const int actual_dim = mesh.dimension();
313 assert(!is_problem_scalar);
314
315 Eigen::MatrixXd local_val, local_stress, local_mises;
316
317 int num_quadr_pts = 0;
318 result.resize(disc_orders.sum(), actual_dim == 2 ? 3 : 6);
319 result.setZero();
320 von_mises.resize(disc_orders.sum(), 1);
321 von_mises.setZero();
322 for (int e = 0; e < mesh.n_elements(); ++e)
323 {
324 // Compute quadrature points for element
325 quadrature::Quadrature quadr;
326 if (mesh.is_simplex(e))
327 {
328 if (mesh.is_volume())
329 {
330 quadrature::TetQuadrature f;
331 f.get_quadrature(disc_orders(e), quadr);
332 }
333 else
334 {
335 quadrature::TriQuadrature f;
336 f.get_quadrature(disc_orders(e), quadr);
337 }
338 }
339 else if (mesh.is_cube(e))
340 {
341 if (mesh.is_volume())
342 {
343 quadrature::HexQuadrature f;
344 f.get_quadrature(disc_orders(e), quadr);
345 }
346 else
347 {
348 quadrature::QuadQuadrature f;
349 f.get_quadrature(disc_orders(e), quadr);
350 }
351 }
352 else
353 {
354 continue;
355 }
356
357 std::vector<std::pair<std::string, Eigen::MatrixXd>> tmp_s, tmp_t;
358
359 assembler.compute_scalar_value(OutputData(t, e, bases[e], gbases[e], quadr.points, fun), tmp_s);
360 assembler.compute_tensor_value(OutputData(t, e, bases[e], gbases[e], quadr.points, fun), tmp_t);
361
362 local_mises = tmp_s[0].second;
363 local_val = tmp_t[0].second;
364
365 if (num_quadr_pts + local_val.rows() >= result.rows())
366 {
367 result.conservativeResize(
368 std::max(num_quadr_pts + local_val.rows() + 1, 2 * result.rows()),
369 result.cols());
370 von_mises.conservativeResize(result.rows(), von_mises.cols());
371 }
372 flattened_tensor_coeffs(local_val, local_stress);
373 result.block(num_quadr_pts, 0, local_stress.rows(), local_stress.cols()) = local_stress;
374 von_mises.block(num_quadr_pts, 0, local_mises.rows(), local_mises.cols()) = local_mises;
375 num_quadr_pts += local_val.rows();
376 }
377 result.conservativeResize(num_quadr_pts, result.cols());
378 von_mises.conservativeResize(num_quadr_pts, von_mises.cols());
379 }
380
381 void Evaluator::interpolate_function(
382 const mesh::Mesh &mesh,
383 const bool is_problem_scalar,
384 const std::vector<basis::ElementBases> &bases,
385 const Eigen::VectorXi &disc_orders,
386 const std::map<int, Eigen::MatrixXd> &polys,
387 const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d,
388 const utils::RefElementSampler &sampler,
389 const int n_points,
390 const Eigen::MatrixXd &fun,
391 Eigen::MatrixXd &result,
392 const bool use_sampler,
393 const bool boundary_only)
394 {
395 int actual_dim = 1;
396 if (!is_problem_scalar)
397 actual_dim = mesh.dimension();
398 interpolate_function(mesh, actual_dim, bases, disc_orders,
399 polys, polys_3d, sampler, n_points,
400 fun, result, use_sampler, boundary_only);
401 }
402
403 void Evaluator::mark_flipped_cells(
404 const mesh::Mesh &mesh,
405 const std::vector<basis::ElementBases> &gbasis,
406 const std::vector<basis::ElementBases> &basis,
407 const Eigen::VectorXi &disc_orders,
408 const std::map<int, Eigen::MatrixXd> &polys,
409 const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d,
410 const utils::RefElementSampler &sampler,
411 const int n_points,
412 const Eigen::MatrixXd &fun,
413 Eigen::Vector<bool, -1> &result,
414 const bool use_sampler,
415 const bool boundary_only)
416 {
417 if (fun.size() <= 0)
418 {
419 logger().error("Solve the problem first!");
420 return;
421 }
422
423 result.setZero(n_points);
424
425 int index = 0;
426
427 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
428
429 const auto invalidList = utils::count_invalid(mesh.dimension(), basis, gbasis, fun);
430
431 for (int i = 0; i < int(basis.size()); ++i)
432 {
433 const ElementBases &bs = basis[i];
434 Eigen::MatrixXd local_pts;
435
436 if (boundary_only && mesh.is_volume() && !mesh.is_boundary_element(i))
437 continue;
438
439 if (use_sampler)
440 {
441 if (mesh.is_simplex(i))
442 local_pts = sampler.simplex_points();
443 else if (mesh.is_cube(i))
444 local_pts = sampler.cube_points();
445 else
446 {
447 if (mesh.is_volume())
448 sampler.sample_polyhedron(polys_3d.at(i).first, polys_3d.at(i).second, local_pts, vis_faces_poly, vis_edges_poly);
449 else
450 sampler.sample_polygon(polys.at(i), local_pts, vis_faces_poly, vis_edges_poly);
451 }
452 }
453 else
454 {
455 if (mesh.is_volume())
456 {
457 if (mesh.is_simplex(i))
458 autogen::p_nodes_3d(disc_orders(i), local_pts);
459 else if (mesh.is_cube(i))
460 autogen::q_nodes_3d(disc_orders(i), local_pts);
461 else
462 continue;
463 }
464 else
465 {
466 if (mesh.is_simplex(i))
467 autogen::p_nodes_2d(disc_orders(i), local_pts);
468 else if (mesh.is_cube(i))
469 autogen::q_nodes_2d(disc_orders(i), local_pts);
470 else
471 continue;
472 }
473 }
474
475 if (std::find(invalidList.begin(), invalidList.end(), i) != invalidList.end())
476 result.segment(index, local_pts.rows()).array() = true;
477 index += local_pts.rows();
478 }
479 }
480
481 void Evaluator::interpolate_function(
482 const mesh::Mesh &mesh,
483 const int actual_dim,
484 const std::vector<basis::ElementBases> &basis,
485 const Eigen::VectorXi &disc_orders,
486 const std::map<int, Eigen::MatrixXd> &polys,
487 const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d,
488 const utils::RefElementSampler &sampler,
489 const int n_points,
490 const Eigen::MatrixXd &fun,
491 Eigen::MatrixXd &result,
492 const bool use_sampler,
493 const bool boundary_only)
494 {
495 if (fun.size() <= 0)
496 {
497 logger().error("Solve the problem first!");
498 return;
499 }
500
501 std::vector<AssemblyValues> tmp;
502
503 result.resize(n_points, actual_dim);
504
505 int index = 0;
506
507 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
508
509 for (int i = 0; i < int(basis.size()); ++i)
510 {
511 const ElementBases &bs = basis[i];
512 Eigen::MatrixXd local_pts;
513
514 if (boundary_only && mesh.is_volume() && !mesh.is_boundary_element(i))
515 continue;
516
517 if (use_sampler)
518 {
519 if (mesh.is_simplex(i))
520 local_pts = sampler.simplex_points();
521 else if (mesh.is_cube(i))
522 local_pts = sampler.cube_points();
523 else
524 {
525 if (mesh.is_volume())
526 sampler.sample_polyhedron(polys_3d.at(i).first, polys_3d.at(i).second, local_pts, vis_faces_poly, vis_edges_poly);
527 else
528 sampler.sample_polygon(polys.at(i), local_pts, vis_faces_poly, vis_edges_poly);
529 }
530 }
531 else
532 {
533 if (mesh.is_volume())
534 {
535 if (mesh.is_simplex(i))
536 autogen::p_nodes_3d(disc_orders(i), local_pts);
537 else if (mesh.is_cube(i))
538 autogen::q_nodes_3d(disc_orders(i), local_pts);
539 else
540 continue;
541 }
542 else
543 {
544 if (mesh.is_simplex(i))
545 autogen::p_nodes_2d(disc_orders(i), local_pts);
546 else if (mesh.is_cube(i))
547 autogen::q_nodes_2d(disc_orders(i), local_pts);
548 else
549 continue;
550 }
551 }
552
553 Eigen::MatrixXd local_res = Eigen::MatrixXd::Zero(local_pts.rows(), actual_dim);
554 bs.evaluate_bases(local_pts, tmp);
555 for (size_t j = 0; j < bs.bases.size(); ++j)
556 {
557 const Basis &b = bs.bases[j];
558
559 for (int d = 0; d < actual_dim; ++d)
560 {
561 for (size_t ii = 0; ii < b.global().size(); ++ii)
562 local_res.col(d) += b.global()[ii].val * tmp[j].val * fun(b.global()[ii].index * actual_dim + d);
563 }
564 }
565
566 result.block(index, 0, local_res.rows(), actual_dim) = local_res;
567 index += local_res.rows();
568 }
569 }
570
571 void Evaluator::interpolate_at_local_vals(
572 const mesh::Mesh &mesh,
573 const bool is_problem_scalar,
574 const std::vector<basis::ElementBases> &bases,
575 const std::vector<basis::ElementBases> &gbases,
576 const int el_index,
577 const Eigen::MatrixXd &local_pts,
578 const Eigen::MatrixXd &fun,
579 Eigen::MatrixXd &result,
580 Eigen::MatrixXd &result_grad)
581 {
582 int actual_dim = 1;
583 if (!is_problem_scalar)
584 actual_dim = mesh.dimension();
585 interpolate_at_local_vals(mesh, actual_dim, bases, gbases, el_index,
586 local_pts, fun, result, result_grad);
587 }
588
589 void Evaluator::interpolate_at_local_vals(
590 const mesh::Mesh &mesh,
591 const int actual_dim,
592 const std::vector<basis::ElementBases> &bases,
593 const std::vector<basis::ElementBases> &gbases,
594 const int el_index,
595 const Eigen::MatrixXd &local_pts,
596 const Eigen::MatrixXd &fun,
597 Eigen::MatrixXd &result,
598 Eigen::MatrixXd &result_grad)
599 {
600 if (fun.size() <= 0)
601 {
602 logger().error("Solve the problem first!");
603 return;
604 }
605
606 assert(local_pts.cols() == mesh.dimension());
607 assert(fun.cols() == 1);
608
609 const ElementBases &gbs = gbases[el_index];
610 const ElementBases &bs = bases[el_index];
611
612 ElementAssemblyValues vals;
613 vals.compute(el_index, mesh.is_volume(), local_pts, bs, gbs);
614
615 result.resize(vals.val.rows(), actual_dim);
616 result.setZero();
617
618 result_grad.resize(vals.val.rows(), mesh.dimension() * actual_dim);
619 result_grad.setZero();
620
621 const int n_loc_bases = int(vals.basis_values.size());
622
623 for (int i = 0; i < n_loc_bases; ++i)
624 {
625 const auto &val = vals.basis_values[i];
626
627 for (size_t ii = 0; ii < val.global.size(); ++ii)
628 {
629 for (int d = 0; d < actual_dim; ++d)
630 {
631 result.col(d) += val.global[ii].val * fun(val.global[ii].index * actual_dim + d) * val.val;
632 result_grad.block(0, d * val.grad_t_m.cols(), result_grad.rows(), val.grad_t_m.cols()) += val.global[ii].val * fun(val.global[ii].index * actual_dim + d) * val.grad_t_m;
633 }
634 }
635 }
636 }
637
638 void Evaluator::interpolate_at_local_vals(const int el_index, const int dim, const int actual_dim, const assembler::ElementAssemblyValues &vals, const Eigen::MatrixXd &fun, Eigen::MatrixXd &result, Eigen::MatrixXd &result_grad)
639 {
640 if (fun.size() <= 0)
641 {
642 logger().error("Solve the problem first!");
643 return;
644 }
645
646 assert(fun.cols() == 1);
647
648 result.resize(vals.val.rows(), actual_dim);
649 result.setZero();
650
651 result_grad.resize(vals.val.rows(), dim * actual_dim);
652 result_grad.setZero();
653
654 const int n_loc_bases = int(vals.basis_values.size());
655
656 for (int i = 0; i < n_loc_bases; ++i)
657 {
658 const auto &val = vals.basis_values[i];
659
660 for (size_t ii = 0; ii < val.global.size(); ++ii)
661 {
662 for (int d = 0; d < actual_dim; ++d)
663 {
664 result.col(d) += val.global[ii].val * fun(val.global[ii].index * actual_dim + d) * val.val;
665 result_grad.block(0, d * val.grad_t_m.cols(), result_grad.rows(), val.grad_t_m.cols()) += val.global[ii].val * fun(val.global[ii].index * actual_dim + d) * val.grad_t_m;
666 }
667 }
668 }
669 }
670
671 bool Evaluator::check_scalar_value(
672 const mesh::Mesh &mesh,
673 const bool is_problem_scalar,
674 const std::vector<basis::ElementBases> &bases,
675 const std::vector<basis::ElementBases> &gbases,
676 const Eigen::VectorXi &disc_orders,
677 const std::map<int, Eigen::MatrixXd> &polys,
678 const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d,
679 const assembler::Assembler &assembler,
680 const utils::RefElementSampler &sampler,
681 const Eigen::MatrixXd &fun,
682 const double t,
683 const bool use_sampler,
684 const bool boundary_only)
685 {
686 if (fun.size() <= 0)
687 {
688 logger().error("Solve the problem first!");
689 return true;
690 }
691
692 assert(!is_problem_scalar);
693
694 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
695
696 std::vector<std::pair<std::string, Eigen::MatrixXd>> tmp_s;
697
698 for (int i = 0; i < int(bases.size()); ++i)
699 {
700 if (boundary_only && mesh.is_volume() && !mesh.is_boundary_element(i))
701 continue;
702
703 const ElementBases &bs = bases[i];
704 const ElementBases &gbs = gbases[i];
705 Eigen::MatrixXd local_pts;
706
707 if (use_sampler)
708 {
709 if (mesh.is_simplex(i))
710 local_pts = sampler.simplex_points();
711 else if (mesh.is_cube(i))
712 local_pts = sampler.cube_points();
713 else
714 {
715 if (mesh.is_volume())
716 sampler.sample_polyhedron(polys_3d.at(i).first, polys_3d.at(i).second, local_pts, vis_faces_poly, vis_edges_poly);
717 else
718 sampler.sample_polygon(polys.at(i), local_pts, vis_faces_poly, vis_edges_poly);
719 }
720 }
721 else
722 {
723 if (mesh.is_volume())
724 {
725 if (mesh.is_simplex(i))
726 autogen::p_nodes_3d(disc_orders(i), local_pts);
727 else if (mesh.is_cube(i))
728 autogen::q_nodes_3d(disc_orders(i), local_pts);
729 else
730 continue;
731 }
732 else
733 {
734 if (mesh.is_simplex(i))
735 autogen::p_nodes_2d(disc_orders(i), local_pts);
736 else if (mesh.is_cube(i))
737 autogen::q_nodes_2d(disc_orders(i), local_pts);
738 else
739 continue;
740 }
741 }
742
743 assembler.compute_scalar_value(OutputData(t, i, bs, gbs, local_pts, fun), tmp_s);
744
745 for (const auto &s : tmp_s)
746 if (std::isnan(s.second.norm()))
747 return false;
748 }
749
750 return true;
751 }
752
753 void Evaluator::compute_scalar_value(
754 const mesh::Mesh &mesh,
755 const bool is_problem_scalar,
756 const std::vector<basis::ElementBases> &bases,
757 const std::vector<basis::ElementBases> &gbases,
758 const Eigen::VectorXi &disc_orders,
759 const std::map<int, Eigen::MatrixXd> &polys,
760 const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d,
761 const assembler::Assembler &assembler,
762 const utils::RefElementSampler &sampler,
763 const int n_points,
764 const Eigen::MatrixXd &fun,
765 const double t,
766 std::vector<assembler::Assembler::NamedMatrix> &result,
767 const bool use_sampler,
768 const bool boundary_only)
769 {
770 if (fun.size() <= 0)
771 {
772 logger().error("Solve the problem first!");
773 return;
774 }
775
776 result.clear();
777
778 assert(!is_problem_scalar);
779
780 int index = 0;
781
782 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
783 std::vector<std::pair<std::string, Eigen::MatrixXd>> tmp_s;
784
785 for (int i = 0; i < int(bases.size()); ++i)
786 {
787 if (boundary_only && mesh.is_volume() && !mesh.is_boundary_element(i))
788 continue;
789
790 const ElementBases &bs = bases[i];
791 const ElementBases &gbs = gbases[i];
792 Eigen::MatrixXd local_pts;
793
794 if (use_sampler)
795 {
796 if (mesh.is_simplex(i))
797 local_pts = sampler.simplex_points();
798 else if (mesh.is_cube(i))
799 local_pts = sampler.cube_points();
800 else
801 {
802 if (mesh.is_volume())
803 sampler.sample_polyhedron(polys_3d.at(i).first, polys_3d.at(i).second, local_pts, vis_faces_poly, vis_edges_poly);
804 else
805 sampler.sample_polygon(polys.at(i), local_pts, vis_faces_poly, vis_edges_poly);
806 }
807 }
808 else
809 {
810 if (mesh.is_volume())
811 {
812 if (mesh.is_simplex(i))
813 autogen::p_nodes_3d(disc_orders(i), local_pts);
814 else if (mesh.is_cube(i))
815 autogen::q_nodes_3d(disc_orders(i), local_pts);
816 else
817 continue;
818 }
819 else
820 {
821 if (mesh.is_simplex(i))
822 autogen::p_nodes_2d(disc_orders(i), local_pts);
823 else if (mesh.is_cube(i))
824 autogen::q_nodes_2d(disc_orders(i), local_pts);
825 else
826 continue;
827 }
828 }
829
830 assembler.compute_scalar_value(OutputData(t, i, bs, gbs, local_pts, fun), tmp_s);
831
832 if (result.empty())
833 {
834 result.resize(tmp_s.size());
835 for (int k = 0; k < tmp_s.size(); ++k)
836 {
837 result[k].first = tmp_s[k].first;
838 result[k].second.resize(n_points, 1);
839 }
840 }
841
842 for (int k = 0; k < tmp_s.size(); ++k)
843 {
844 assert(local_pts.rows() == tmp_s[k].second.rows());
845 result[k].second.block(index, 0, tmp_s[k].second.rows(), 1) = tmp_s[k].second;
846 }
847 index += local_pts.rows();
848 }
849 }
850
851 void Evaluator::compute_tensor_value(
852 const mesh::Mesh &mesh,
853 const bool is_problem_scalar,
854 const std::vector<basis::ElementBases> &bases,
855 const std::vector<basis::ElementBases> &gbases,
856 const Eigen::VectorXi &disc_orders,
857 const std::map<int, Eigen::MatrixXd> &polys,
858 const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d,
859 const assembler::Assembler &assembler,
860 const utils::RefElementSampler &sampler,
861 const int n_points,
862 const Eigen::MatrixXd &fun,
863 const double t,
864 std::vector<assembler::Assembler::NamedMatrix> &result,
865 const bool use_sampler,
866 const bool boundary_only)
867 {
868 if (fun.size() <= 0)
869 {
870 logger().error("Solve the problem first!");
871 return;
872 }
873
874 result.clear();
875
876 const int actual_dim = mesh.dimension();
877 assert(!is_problem_scalar);
878
879 int index = 0;
880
881 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
882 std::vector<std::pair<std::string, Eigen::MatrixXd>> tmp_t;
883
884 for (int i = 0; i < int(bases.size()); ++i)
885 {
886 if (boundary_only && mesh.is_volume() && !mesh.is_boundary_element(i))
887 continue;
888
889 const ElementBases &bs = bases[i];
890 const ElementBases &gbs = gbases[i];
891 Eigen::MatrixXd local_pts;
892
893 if (use_sampler)
894 {
895 if (mesh.is_simplex(i))
896 local_pts = sampler.simplex_points();
897 else if (mesh.is_cube(i))
898 local_pts = sampler.cube_points();
899 else
900 {
901 if (mesh.is_volume())
902 sampler.sample_polyhedron(polys_3d.at(i).first, polys_3d.at(i).second, local_pts, vis_faces_poly, vis_edges_poly);
903 else
904 sampler.sample_polygon(polys.at(i), local_pts, vis_faces_poly, vis_edges_poly);
905 }
906 }
907 else
908 {
909 if (mesh.is_volume())
910 {
911 if (mesh.is_simplex(i))
912 autogen::p_nodes_3d(disc_orders(i), local_pts);
913 else if (mesh.is_cube(i))
914 autogen::q_nodes_3d(disc_orders(i), local_pts);
915 else
916 continue;
917 }
918 else
919 {
920 if (mesh.is_simplex(i))
921 autogen::p_nodes_2d(disc_orders(i), local_pts);
922 else if (mesh.is_cube(i))
923 autogen::q_nodes_2d(disc_orders(i), local_pts);
924 else
925 continue;
926 }
927 }
928
929 assembler.compute_tensor_value(OutputData(t, i, bs, gbs, local_pts, fun), tmp_t);
930
931 if (result.empty())
932 {
933 result.resize(tmp_t.size());
934 for (int k = 0; k < tmp_t.size(); ++k)
935 {
936 result[k].first = tmp_t[k].first;
937 result[k].second.resize(n_points, actual_dim * actual_dim);
938 }
939 }
940
941 for (int k = 0; k < tmp_t.size(); ++k)
942 {
943 assert(local_pts.rows() == tmp_t[k].second.rows());
944 result[k].second.block(index, 0, tmp_t[k].second.rows(), tmp_t[k].second.cols()) = tmp_t[k].second;
945 }
946 index += local_pts.rows();
947 }
948 }
949
950 Eigen::MatrixXd Evaluator::get_bases_position(
951 const int n_bases,
952 const std::shared_ptr<mesh::MeshNodes> mesh_nodes)
953 {
954 Eigen::MatrixXd func;
955 func.setZero(n_bases, mesh_nodes->node_position(0).size());
956
957 for (int i = 0; i < n_bases; i++)
958 func.row(i) = mesh_nodes->node_position(i);
959
960 return func;
961 }
962
963 Eigen::MatrixXd Evaluator::generate_linear_field(
964 const int n_bases,
965 const std::shared_ptr<mesh::MeshNodes> mesh_nodes,
966 const Eigen::MatrixXd &grad)
967 {
968 return utils::flatten(get_bases_position(n_bases, mesh_nodes) * grad.transpose());
969 }
970
971 Eigen::VectorXd Evaluator::integrate_function(
972 const std::vector<basis::ElementBases> &bases,
973 const std::vector<basis::ElementBases> &gbases,
974 const Eigen::MatrixXd &fun,
975 const int dim,
976 const int actual_dim)
977 {
978 Eigen::VectorXd result;
979 result.setZero(actual_dim);
980 for (int e = 0; e < bases.size(); ++e)
981 {
982 ElementAssemblyValues vals;
983 vals.compute(e, dim == 3, bases[e], gbases[e]);
984
985 Eigen::MatrixXd u, grad_u;
986 io::Evaluator::interpolate_at_local_vals(e, dim, actual_dim, vals, fun, u, grad_u);
987 const quadrature::Quadrature &quadrature = vals.quadrature;
988 Eigen::VectorXd da = vals.det.array() * quadrature.weights.array();
989 result += u.transpose() * da;
990 }
991
992 return result;
993 }
994} // namespace polyfem::io
val
double val
Definition Assembler.cpp:86
da
QuadratureVector da
Definition Assembler.cpp:23
vals
ElementAssemblyValues vals
Definition Assembler.cpp:22
quadrature
Quadrature quadrature
Definition MassMatrixAssembler.cpp:30
faces
std::vector< Eigen::VectorXi > faces
Definition PressureAssembler.cpp:52
auto_p_bases.hpp
auto_q_bases.hpp
polyfem::assembler::Assembler::compute_scalar_value
virtual void compute_scalar_value(const OutputData &data, std::vector< NamedMatrix > &result) const
Definition Assembler.hpp:128
polyfem::assembler::Assembler::compute_tensor_value
virtual void compute_tensor_value(const OutputData &data, std::vector< NamedMatrix > &result) const
Definition Assembler.hpp:133
polyfem::assembler::AssemblyValues
stores per local bases evaluations
Definition AssemblyValues.hpp:14
polyfem::assembler::AssemblyValues::global
std::vector< basis::Local2Global > global
Definition AssemblyValues.hpp:20
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::Basis
Represents one basis function and its gradient.
Definition Basis.hpp:44
polyfem::basis::Basis::global
const std::vector< Local2Global > & global() const
Definition Basis.hpp:104
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::evaluate_bases
void evaluate_bases(const Eigen::MatrixXd &uv, std::vector< assembler::AssemblyValues > &basis_values) const
evaluate stored bases at given points on the reference element saves results to basis_values
Definition ElementBases.hpp:71
polyfem::basis::ElementBases::bases
std::vector< Basis > bases
one basis function per node in the element
Definition ElementBases.hpp:27
polyfem::io::Evaluator::generate_linear_field
static Eigen::MatrixXd generate_linear_field(const int n_bases, const std::shared_ptr< mesh::MeshNodes > mesh_nodes, const Eigen::MatrixXd &grad)
Definition Evaluator.cpp:963
polyfem::io::Evaluator::check_scalar_value
bool check_scalar_value(const mesh::Mesh &mesh, const bool is_problem_scalar, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const Eigen::VectorXi &disc_orders, const std::map< int, Eigen::MatrixXd > &polys, const std::map< int, std::pair< Eigen::MatrixXd, Eigen::MatrixXi > > &polys_3d, const assembler::Assembler &assembler, const utils::RefElementSampler &sampler, const Eigen::MatrixXd &fun, const double t, const bool use_sampler, const bool boundary_only)
checks if mises are not nan
Definition Evaluator.cpp:671
polyfem::io::Evaluator::interpolate_at_local_vals
static void interpolate_at_local_vals(const mesh::Mesh &mesh, const bool is_problem_scalar, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const int el_index, const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &fun, Eigen::MatrixXd &result, Eigen::MatrixXd &result_grad)
interpolate solution and gradient at element (calls interpolate_at_local_vals with sol)
Definition Evaluator.cpp:571
polyfem::io::Evaluator::compute_stress_at_quadrature_points
static void compute_stress_at_quadrature_points(const mesh::Mesh &mesh, const bool is_problem_scalar, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const Eigen::VectorXi &disc_orders, const assembler::Assembler &assembler, const Eigen::MatrixXd &fun, const double t, Eigen::MatrixXd &result, Eigen::VectorXd &von_mises)
compute von mises stress at quadrature points for the function fun, also compute the interpolated fun...
Definition Evaluator.cpp:284
polyfem::io::Evaluator::get_bases_position
static Eigen::MatrixXd get_bases_position(const int n_bases, const std::shared_ptr< mesh::MeshNodes > mesh_nodes)
Definition Evaluator.cpp:950
polyfem::io::Evaluator::interpolate_boundary_function
static void interpolate_boundary_function(const mesh::Mesh &mesh, const bool is_problem_scalar, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const Eigen::MatrixXd &pts, const Eigen::MatrixXi &faces, const Eigen::MatrixXd &fun, const bool compute_avg, Eigen::MatrixXd &result)
computes integrated solution (fun) per surface face.
Definition Evaluator.cpp:57
polyfem::io::Evaluator::mark_flipped_cells
static void mark_flipped_cells(const mesh::Mesh &mesh, const std::vector< basis::ElementBases > &gbasis, const std::vector< basis::ElementBases > &basis, const Eigen::VectorXi &disc_orders, const std::map< int, Eigen::MatrixXd > &polys, const std::map< int, std::pair< Eigen::MatrixXd, Eigen::MatrixXi > > &polys_3d, const utils::RefElementSampler &sampler, const int n_points, const Eigen::MatrixXd &fun, Eigen::Vector< bool, -1 > &result, const bool use_sampler, const bool boundary_only)
Definition Evaluator.cpp:403
polyfem::io::Evaluator::interpolate_function
static void interpolate_function(const mesh::Mesh &mesh, const bool is_problem_scalar, const std::vector< basis::ElementBases > &bases, const Eigen::VectorXi &disc_orders, const std::map< int, Eigen::MatrixXd > &polys, const std::map< int, std::pair< Eigen::MatrixXd, Eigen::MatrixXi > > &polys_3d, const utils::RefElementSampler &sampler, const int n_points, const Eigen::MatrixXd &fun, Eigen::MatrixXd &result, const bool use_sampler, const bool boundary_only)
interpolate the function fun.
Definition Evaluator.cpp:381
polyfem::io::Evaluator::average_grad_based_function
static void average_grad_based_function(const mesh::Mesh &mesh, const bool is_problem_scalar, const int n_bases, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const Eigen::VectorXi &disc_orders, const std::map< int, Eigen::MatrixXd > &polys, const std::map< int, std::pair< Eigen::MatrixXd, Eigen::MatrixXi > > &polys_3d, const assembler::Assembler &assembler, const utils::RefElementSampler &sampler, const double t, const int n_points, const Eigen::MatrixXd &fun, std::vector< assembler::Assembler::NamedMatrix > &result_scalar, std::vector< assembler::Assembler::NamedMatrix > &result_tensor, const bool use_sampler, const bool boundary_only)
computes scalar quantity of funtion (ie von mises for elasticity and norm of velocity for fluid) the ...
Definition Evaluator.cpp:153
polyfem::io::Evaluator::integrate_function
static Eigen::VectorXd integrate_function(const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const Eigen::MatrixXd &fun, const int dim, const int actual_dim)
Definition Evaluator.cpp:971
polyfem::io::Evaluator::compute_scalar_value
static void compute_scalar_value(const mesh::Mesh &mesh, const bool is_problem_scalar, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const Eigen::VectorXi &disc_orders, const std::map< int, Eigen::MatrixXd > &polys, const std::map< int, std::pair< Eigen::MatrixXd, Eigen::MatrixXi > > &polys_3d, const assembler::Assembler &assembler, const utils::RefElementSampler &sampler, const int n_points, const Eigen::MatrixXd &fun, const double t, std::vector< assembler::Assembler::NamedMatrix > &result, const bool use_sampler, const bool boundary_only)
computes scalar quantity of funtion (ie von mises for elasticity and norm of velocity for fluid)
Definition Evaluator.cpp:753
polyfem::io::Evaluator::compute_tensor_value
static void compute_tensor_value(const mesh::Mesh &mesh, const bool is_problem_scalar, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const Eigen::VectorXi &disc_orders, const std::map< int, Eigen::MatrixXd > &polys, const std::map< int, std::pair< Eigen::MatrixXd, Eigen::MatrixXi > > &polys_3d, const assembler::Assembler &assembler, const utils::RefElementSampler &sampler, const int n_points, const Eigen::MatrixXd &fun, const double t, std::vector< assembler::Assembler::NamedMatrix > &result, const bool use_sampler, const bool boundary_only)
compute tensor quantity (ie stress tensor or velocy)
Definition Evaluator.cpp:851
polyfem::mesh::Mesh3D::n_cell_faces
virtual int n_cell_faces(const int c_id) const =0
polyfem::mesh::Mesh3D::cell_face
virtual int cell_face(const int c_id, const int lf_id) const =0
polyfem::mesh::Mesh
Abstract mesh class to capture 2d/3d conforming and non-conforming meshes.
Definition Mesh.hpp:39
polyfem::mesh::Mesh::n_elements
int n_elements() const
utitlity to return the number of elements, cells or faces in 3d and 2d
Definition Mesh.hpp:161
polyfem::mesh::Mesh::face_barycenter
virtual RowVectorNd face_barycenter(const int f) const =0
face barycenter
polyfem::mesh::Mesh::is_cube
bool is_cube(const int el_id) const
checks if element is cube compatible
Definition Mesh.cpp:352
polyfem::mesh::Mesh::is_boundary_face
virtual bool is_boundary_face(const int face_global_id) const =0
is face boundary
polyfem::mesh::Mesh::is_simplex
bool is_simplex(const int el_id) const
checks if element is simples compatible
Definition Mesh.cpp:422
polyfem::mesh::Mesh::is_volume
virtual bool is_volume() const =0
checks if mesh is volume
polyfem::mesh::Mesh::dimension
int dimension() const
utily for dimension
Definition Mesh.hpp:151
polyfem::mesh::Mesh::is_boundary_element
virtual bool is_boundary_element(const int element_global_id) const =0
is cell boundary
polyfem::utils::BoundarySampler::quadrature_for_quad_face
static void quadrature_for_quad_face(int index, int order, const int gid, const mesh::Mesh &mesh, Eigen::MatrixXd &uv, Eigen::MatrixXd &points, Eigen::VectorXd &weights)
Definition BoundarySampler.cpp:150
polyfem::utils::BoundarySampler::quadrature_for_tri_face
static void quadrature_for_tri_face(int index, int order, const int gid, const mesh::Mesh &mesh, Eigen::MatrixXd &uv, Eigen::MatrixXd &points, Eigen::VectorXd &weights)
Definition BoundarySampler.cpp:185
polyfem::utils::RefElementSampler::sample_polygon
void sample_polygon(const Eigen::MatrixXd &poly, Eigen::MatrixXd &pts, Eigen::MatrixXi &faces, Eigen::MatrixXi &edges) const
Definition RefElementSampler.cpp:368
polyfem::utils::RefElementSampler::simplex_points
const Eigen::MatrixXd & simplex_points() const
Definition RefElementSampler.hpp:43
polyfem::utils::RefElementSampler::sample_polyhedron
void sample_polyhedron(const Eigen::MatrixXd &vertices, const Eigen::MatrixXi &f, Eigen::MatrixXd &pts, Eigen::MatrixXi &faces, Eigen::MatrixXi &edges) const
Definition RefElementSampler.cpp:409
polyfem::utils::RefElementSampler::cube_points
const Eigen::MatrixXd & cube_points() const
Definition RefElementSampler.hpp:37
polyfem::autogen::q_nodes_2d
void q_nodes_2d(const int q, Eigen::MatrixXd &val)
Definition auto_q_bases_2d_nodes.cpp:67
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::autogen::q_nodes_3d
void q_nodes_3d(const int q, Eigen::MatrixXd &val)
Definition auto_q_bases_3d_nodes.cpp:149
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::flatten
Eigen::VectorXd flatten(const Eigen::MatrixXd &X)
Flatten rowwises.
Definition MatrixUtils.cpp:34
polyfem::logger
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:42
polyfem::RowVectorNd
Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor, 1, 3 > RowVectorNd
Definition Types.hpp:13