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