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