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