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 autogen::prism_nodes_3d(disc_orders(i), disc_ordersq(i), local_pts);
227 }
228 else
229 {
230 // not supported for polys
231 continue;
232 }
233
234 vals.compute(i, actual_dim == 3, bases[i], gbases[i]);
235 const quadrature::Quadrature &quadrature = vals.quadrature;
236 const double area = (vals.det.array() * quadrature.weights.array()).sum();
237
238 assembler.compute_scalar_value(OutputData(t, i, bs, gbs, local_pts, fun), tmp_s);
239 assembler.compute_tensor_value(OutputData(t, i, bs, gbs, local_pts, fun), tmp_t);
240
241 for (size_t j = 0; j < bs.bases.size(); ++j)
242 {
243 const Basis &b = bs.bases[j];
244 if (b.global().size() > 1)
245 continue;
246
247 auto &global = b.global().front();
248 areas(global.index) += area;
249 }
250
251 if (avg_scalar.empty())
252 {
253 avg_scalar.resize(tmp_s.size());
254 for (auto &m : avg_scalar)
255 {
256 m.resize(n_bases, 1);
257 m.setZero();
258 }
259 }
260
261 if (avg_tensor.empty())
262 {
263 avg_tensor.resize(tmp_t.size());
264 for (auto &m : avg_tensor)
265 {
266 m.resize(n_bases, actual_dim * actual_dim);
267 m.setZero();
268 }
269 }
270
271 for (int k = 0; k < tmp_s.size(); ++k)
272 {
273 local_val = tmp_s[k].second;
274
275 for (size_t j = 0; j < bs.bases.size(); ++j)
276 {
277 const Basis &b = bs.bases[j];
278 if (b.global().size() > 1)
279 continue;
280
281 auto &global = b.global().front();
282 avg_scalar[k](global.index) += local_val(j) * area;
283 }
284 }
285
286 for (int k = 0; k < tmp_t.size(); ++k)
287 {
288 local_val = tmp_t[k].second;
289
290 for (size_t j = 0; j < bs.bases.size(); ++j)
291 {
292 const Basis &b = bs.bases[j];
293 if (b.global().size() > 1)
294 continue;
295
296 auto &global = b.global().front();
297 avg_tensor[k].row(global.index) += local_val.row(j) * area;
298 }
299 }
300 }
301
302 for (auto &m : avg_scalar)
303 {
304 m.array() /= areas.array();
305 }
306
307 for (auto &m : avg_tensor)
308 {
309 for (int i = 0; i < m.rows(); ++i)
310 {
311 m.row(i).array() /= areas(i);
312 }
313 }
314
315 result_scalar.resize(tmp_s.size());
316 for (int k = 0; k < tmp_s.size(); ++k)
317 {
318 result_scalar[k].first = tmp_s[k].first;
319 interpolate_function(mesh, 1, bases, disc_orders, disc_ordersq, polys, polys_3d, sampler, n_points,
320 avg_scalar[k], result_scalar[k].second, use_sampler, boundary_only);
321 }
322
323 result_tensor.resize(tmp_t.size());
324 for (int k = 0; k < tmp_t.size(); ++k)
325 {
326 result_tensor[k].first = tmp_t[k].first;
327 interpolate_function(mesh, actual_dim * actual_dim, bases, disc_orders, disc_ordersq, polys, polys_3d, sampler, n_points,
328 utils::flatten(avg_tensor[k]), result_tensor[k].second, use_sampler, boundary_only);
329 }
330 }
331
332 void Evaluator::compute_stress_at_quadrature_points(
333 const mesh::Mesh &mesh,
334 const bool is_problem_scalar,
335 const std::vector<basis::ElementBases> &bases,
336 const std::vector<basis::ElementBases> &gbases,
337 const Eigen::VectorXi &disc_orders,
338 const Eigen::VectorXi &disc_ordersq,
339 const assembler::Assembler &assembler,
340 const Eigen::MatrixXd &fun,
341 const double t,
342 Eigen::MatrixXd &result,
343 Eigen::VectorXd &von_mises)
344 {
345 // if (!mesh)
346 // {
347 // logger().error("Load the mesh first!");
348 // return;
349 // }
350 if (fun.size() <= 0)
351 {
352 logger().error("Solve the problem first!");
353 return;
354 }
355 if (is_problem_scalar)
356 {
357 logger().error("Define a tensor problem!");
358 return;
359 }
360
361 const int actual_dim = mesh.dimension();
362 assert(!is_problem_scalar);
363
364 Eigen::MatrixXd local_val, local_stress, local_mises;
365
366 int num_quadr_pts = 0;
367 result.resize(disc_orders.sum(), actual_dim == 2 ? 3 : 6);
368 result.setZero();
369 von_mises.resize(disc_orders.sum(), 1);
370 von_mises.setZero();
371 for (int e = 0; e < mesh.n_elements(); ++e)
372 {
373 // Compute quadrature points for element
374 quadrature::Quadrature quadr;
375 if (mesh.is_simplex(e))
376 {
377 if (mesh.is_volume())
378 {
379 quadrature::TetQuadrature f;
380 f.get_quadrature(disc_orders(e), quadr);
381 }
382 else
383 {
384 quadrature::TriQuadrature f;
385 f.get_quadrature(disc_orders(e), quadr);
386 }
387 }
388 else if (mesh.is_cube(e))
389 {
390 if (mesh.is_volume())
391 {
392 quadrature::HexQuadrature f;
393 f.get_quadrature(disc_orders(e), quadr);
394 }
395 else
396 {
397 quadrature::QuadQuadrature f;
398 f.get_quadrature(disc_orders(e), quadr);
399 }
400 }
401 else if (mesh.is_prism(e))
402 {
403 assert(mesh.is_volume());
404
405 quadrature::PrismQuadrature f;
406 f.get_quadrature(disc_orders(e), disc_ordersq(e), quadr);
407 }
408 else
409 {
410 continue;
411 }
412
413 std::vector<std::pair<std::string, Eigen::MatrixXd>> tmp_s, tmp_t;
414
415 assembler.compute_scalar_value(OutputData(t, e, bases[e], gbases[e], quadr.points, fun), tmp_s);
416 assembler.compute_tensor_value(OutputData(t, e, bases[e], gbases[e], quadr.points, fun), tmp_t);
417
418 local_mises = tmp_s[0].second;
419 local_val = tmp_t[0].second;
420
421 if (num_quadr_pts + local_val.rows() >= result.rows())
422 {
423 result.conservativeResize(
424 std::max(num_quadr_pts + local_val.rows() + 1, 2 * result.rows()),
425 result.cols());
426 von_mises.conservativeResize(result.rows(), von_mises.cols());
427 }
428 flattened_tensor_coeffs(local_val, local_stress);
429 result.block(num_quadr_pts, 0, local_stress.rows(), local_stress.cols()) = local_stress;
430 von_mises.block(num_quadr_pts, 0, local_mises.rows(), local_mises.cols()) = local_mises;
431 num_quadr_pts += local_val.rows();
432 }
433 result.conservativeResize(num_quadr_pts, result.cols());
434 von_mises.conservativeResize(num_quadr_pts, von_mises.cols());
435 }
436
437 void Evaluator::interpolate_function(
438 const mesh::Mesh &mesh,
439 const bool is_problem_scalar,
440 const std::vector<basis::ElementBases> &bases,
441 const Eigen::VectorXi &disc_orders,
442 const Eigen::VectorXi &disc_ordersq,
443 const std::map<int, Eigen::MatrixXd> &polys,
444 const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d,
445 const utils::RefElementSampler &sampler,
446 const int n_points,
447 const Eigen::MatrixXd &fun,
448 Eigen::MatrixXd &result,
449 const bool use_sampler,
450 const bool boundary_only)
451 {
452 int actual_dim = 1;
453 if (!is_problem_scalar)
454 actual_dim = mesh.dimension();
455 interpolate_function(mesh, actual_dim, bases, disc_orders, disc_ordersq,
456 polys, polys_3d, sampler, n_points,
457 fun, result, use_sampler, boundary_only);
458 }
459
460 void Evaluator::mark_flipped_cells(
461 const mesh::Mesh &mesh,
462 const std::vector<basis::ElementBases> &gbasis,
463 const std::vector<basis::ElementBases> &basis,
464 const Eigen::VectorXi &disc_orders,
465 const std::map<int, Eigen::MatrixXd> &polys,
466 const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d,
467 const utils::RefElementSampler &sampler,
468 const int n_points,
469 const Eigen::MatrixXd &fun,
470 Eigen::Vector<bool, -1> &result,
471 const bool use_sampler,
472 const bool boundary_only)
473 {
474 if (fun.size() <= 0)
475 {
476 logger().error("Solve the problem first!");
477 return;
478 }
479
480 result.setZero(n_points);
481
482 int index = 0;
483
484 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
485
486 const auto invalidList = utils::count_invalid(mesh.dimension(), basis, gbasis, fun);
487
488 for (int i = 0; i < int(basis.size()); ++i)
489 {
490 const ElementBases &bs = basis[i];
491 Eigen::MatrixXd local_pts;
492
493 if (boundary_only && mesh.is_volume() && !mesh.is_boundary_element(i))
494 continue;
495
496 if (use_sampler)
497 {
498 if (mesh.is_simplex(i))
499 local_pts = sampler.simplex_points();
500 else if (mesh.is_cube(i))
501 local_pts = sampler.cube_points();
502 else
503 {
504 if (mesh.is_volume())
505 sampler.sample_polyhedron(polys_3d.at(i).first, polys_3d.at(i).second, local_pts, vis_faces_poly, vis_edges_poly);
506 else
507 sampler.sample_polygon(polys.at(i), local_pts, vis_faces_poly, vis_edges_poly);
508 }
509 }
510 else
511 {
512 if (mesh.is_volume())
513 {
514 if (mesh.is_simplex(i))
515 autogen::p_nodes_3d(disc_orders(i), local_pts);
516 else if (mesh.is_cube(i))
517 autogen::q_nodes_3d(disc_orders(i), local_pts);
518 else
519 continue;
520 }
521 else
522 {
523 if (mesh.is_simplex(i))
524 autogen::p_nodes_2d(disc_orders(i), local_pts);
525 else if (mesh.is_cube(i))
526 autogen::q_nodes_2d(disc_orders(i), local_pts);
527 else
528 continue;
529 }
530 }
531
532 if (std::find(invalidList.begin(), invalidList.end(), i) != invalidList.end())
533 result.segment(index, local_pts.rows()).array() = true;
534 index += local_pts.rows();
535 }
536 }
537
538 void Evaluator::interpolate_function(
539 const mesh::Mesh &mesh,
540 const int actual_dim,
541 const std::vector<basis::ElementBases> &basis,
542 const Eigen::VectorXi &disc_orders,
543 const Eigen::VectorXi &disc_ordersq,
544 const std::map<int, Eigen::MatrixXd> &polys,
545 const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d,
546 const utils::RefElementSampler &sampler,
547 const int n_points,
548 const Eigen::MatrixXd &fun,
549 Eigen::MatrixXd &result,
550 const bool use_sampler,
551 const bool boundary_only)
552 {
553 if (fun.size() <= 0)
554 {
555 logger().error("Solve the problem first!");
556 return;
557 }
558 assert(fun.cols() == 1);
559
560 std::vector<AssemblyValues> tmp;
561
562 result.resize(n_points, actual_dim);
563
564 int index = 0;
565
566 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
567
568 for (int i = 0; i < int(basis.size()); ++i)
569 {
570 const ElementBases &bs = basis[i];
571 Eigen::MatrixXd local_pts;
572
573 if (boundary_only && mesh.is_volume() && !mesh.is_boundary_element(i))
574 continue;
575
576 if (use_sampler)
577 {
578 if (mesh.is_simplex(i))
579 local_pts = sampler.simplex_points();
580 else if (mesh.is_cube(i))
581 local_pts = sampler.cube_points();
582 else if (mesh.is_prism(i))
583 local_pts = sampler.prism_points();
584 else
585 {
586 if (mesh.is_volume())
587 sampler.sample_polyhedron(polys_3d.at(i).first, polys_3d.at(i).second, local_pts, vis_faces_poly, vis_edges_poly);
588 else
589 sampler.sample_polygon(polys.at(i), local_pts, vis_faces_poly, vis_edges_poly);
590 }
591 }
592 else
593 {
594 if (mesh.is_volume())
595 {
596 if (mesh.is_simplex(i))
597 autogen::p_nodes_3d(disc_orders(i), local_pts);
598 else if (mesh.is_cube(i))
599 autogen::q_nodes_3d(disc_orders(i), local_pts);
600 else if (mesh.is_prism(i))
601 autogen::prism_nodes_3d(disc_orders(i), disc_ordersq(i), local_pts);
602 else
603 continue;
604 }
605 else
606 {
607 if (mesh.is_simplex(i))
608 autogen::p_nodes_2d(disc_orders(i), local_pts);
609 else if (mesh.is_cube(i))
610 autogen::q_nodes_2d(disc_orders(i), local_pts);
611 else
612 continue;
613 }
614 }
615
616 Eigen::MatrixXd local_res = Eigen::MatrixXd::Zero(local_pts.rows(), actual_dim);
617 bs.evaluate_bases(local_pts, tmp);
618 for (size_t j = 0; j < bs.bases.size(); ++j)
619 {
620 const Basis &b = bs.bases[j];
621
622 for (int d = 0; d < actual_dim; ++d)
623 {
624 for (size_t ii = 0; ii < b.global().size(); ++ii)
625 local_res.col(d) += b.global()[ii].val * tmp[j].val * fun(b.global()[ii].index * actual_dim + d);
626 }
627 }
628
629 result.block(index, 0, local_res.rows(), actual_dim) = local_res;
630 index += local_res.rows();
631 }
632 }
633
634 void Evaluator::interpolate_at_local_vals(
635 const mesh::Mesh &mesh,
636 const bool is_problem_scalar,
637 const std::vector<basis::ElementBases> &bases,
638 const std::vector<basis::ElementBases> &gbases,
639 const int el_index,
640 const Eigen::MatrixXd &local_pts,
641 const Eigen::MatrixXd &fun,
642 Eigen::MatrixXd &result,
643 Eigen::MatrixXd &result_grad)
644 {
645 int actual_dim = 1;
646 if (!is_problem_scalar)
647 actual_dim = mesh.dimension();
648 interpolate_at_local_vals(mesh, actual_dim, bases, gbases, el_index,
649 local_pts, fun, result, result_grad);
650 }
651
652 void Evaluator::interpolate_at_local_vals(
653 const mesh::Mesh &mesh,
654 const int actual_dim,
655 const std::vector<basis::ElementBases> &bases,
656 const std::vector<basis::ElementBases> &gbases,
657 const int el_index,
658 const Eigen::MatrixXd &local_pts,
659 const Eigen::MatrixXd &fun,
660 Eigen::MatrixXd &result,
661 Eigen::MatrixXd &result_grad)
662 {
663 if (fun.size() <= 0)
664 {
665 logger().error("Solve the problem first!");
666 return;
667 }
668
669 assert(local_pts.cols() == mesh.dimension());
670 assert(fun.cols() == 1);
671
672 const ElementBases &gbs = gbases[el_index];
673 const ElementBases &bs = bases[el_index];
674
675 ElementAssemblyValues vals;
676 vals.compute(el_index, mesh.is_volume(), local_pts, bs, gbs);
677
678 result.resize(vals.val.rows(), actual_dim);
679 result.setZero();
680
681 result_grad.resize(vals.val.rows(), mesh.dimension() * actual_dim);
682 result_grad.setZero();
683
684 const int n_loc_bases = int(vals.basis_values.size());
685
686 for (int i = 0; i < n_loc_bases; ++i)
687 {
688 const auto &val = vals.basis_values[i];
689
690 for (size_t ii = 0; ii < val.global.size(); ++ii)
691 {
692 for (int d = 0; d < actual_dim; ++d)
693 {
694 result.col(d) += val.global[ii].val * fun(val.global[ii].index * actual_dim + d) * val.val;
695 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;
696 }
697 }
698 }
699 }
700
701 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)
702 {
703 if (fun.size() <= 0)
704 {
705 logger().error("Solve the problem first!");
706 return;
707 }
708
709 assert(fun.cols() == 1);
710
711 result.resize(vals.val.rows(), actual_dim);
712 result.setZero();
713
714 result_grad.resize(vals.val.rows(), dim * actual_dim);
715 result_grad.setZero();
716
717 const int n_loc_bases = int(vals.basis_values.size());
718
719 for (int i = 0; i < n_loc_bases; ++i)
720 {
721 const auto &val = vals.basis_values[i];
722
723 for (size_t ii = 0; ii < val.global.size(); ++ii)
724 {
725 for (int d = 0; d < actual_dim; ++d)
726 {
727 result.col(d) += val.global[ii].val * fun(val.global[ii].index * actual_dim + d) * val.val;
728 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;
729 }
730 }
731 }
732 }
733
734 bool Evaluator::check_scalar_value(
735 const mesh::Mesh &mesh,
736 const bool is_problem_scalar,
737 const std::vector<basis::ElementBases> &bases,
738 const std::vector<basis::ElementBases> &gbases,
739 const Eigen::VectorXi &disc_orders,
740 const Eigen::VectorXi &disc_ordersq,
741 const std::map<int, Eigen::MatrixXd> &polys,
742 const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d,
743 const assembler::Assembler &assembler,
744 const utils::RefElementSampler &sampler,
745 const Eigen::MatrixXd &fun,
746 const double t,
747 const bool use_sampler,
748 const bool boundary_only)
749 {
750 if (fun.size() <= 0)
751 {
752 logger().error("Solve the problem first!");
753 return true;
754 }
755
756 assert(!is_problem_scalar);
757
758 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
759
760 std::vector<std::pair<std::string, Eigen::MatrixXd>> tmp_s;
761
762 for (int i = 0; i < int(bases.size()); ++i)
763 {
764 if (boundary_only && mesh.is_volume() && !mesh.is_boundary_element(i))
765 continue;
766
767 const ElementBases &bs = bases[i];
768 const ElementBases &gbs = gbases[i];
769 Eigen::MatrixXd local_pts;
770
771 if (use_sampler)
772 {
773 if (mesh.is_simplex(i))
774 local_pts = sampler.simplex_points();
775 else if (mesh.is_cube(i))
776 local_pts = sampler.cube_points();
777 else if (mesh.is_prism(i))
778 local_pts = sampler.prism_points();
779 else
780 {
781 if (mesh.is_volume())
782 sampler.sample_polyhedron(polys_3d.at(i).first, polys_3d.at(i).second, local_pts, vis_faces_poly, vis_edges_poly);
783 else
784 sampler.sample_polygon(polys.at(i), local_pts, vis_faces_poly, vis_edges_poly);
785 }
786 }
787 else
788 {
789 if (mesh.is_volume())
790 {
791 if (mesh.is_simplex(i))
792 autogen::p_nodes_3d(disc_orders(i), local_pts);
793 else if (mesh.is_cube(i))
794 autogen::q_nodes_3d(disc_orders(i), local_pts);
795 else if (mesh.is_prism(i))
796 autogen::prism_nodes_3d(disc_orders(i), disc_ordersq(i), local_pts);
797 else
798 continue;
799 }
800 else
801 {
802 if (mesh.is_simplex(i))
803 autogen::p_nodes_2d(disc_orders(i), local_pts);
804 else if (mesh.is_cube(i))
805 autogen::q_nodes_2d(disc_orders(i), local_pts);
806 else
807 continue;
808 }
809 }
810
811 assembler.compute_scalar_value(OutputData(t, i, bs, gbs, local_pts, fun), tmp_s);
812
813 for (const auto &s : tmp_s)
814 if (std::isnan(s.second.norm()))
815 return false;
816 }
817
818 return true;
819 }
820
821 void Evaluator::compute_scalar_value(
822 const mesh::Mesh &mesh,
823 const bool is_problem_scalar,
824 const std::vector<basis::ElementBases> &bases,
825 const std::vector<basis::ElementBases> &gbases,
826 const Eigen::VectorXi &disc_orders,
827 const Eigen::VectorXi &disc_ordersq,
828 const std::map<int, Eigen::MatrixXd> &polys,
829 const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d,
830 const assembler::Assembler &assembler,
831 const utils::RefElementSampler &sampler,
832 const int n_points,
833 const Eigen::MatrixXd &fun,
834 const double t,
835 std::vector<assembler::Assembler::NamedMatrix> &result,
836 const bool use_sampler,
837 const bool boundary_only)
838 {
839 if (fun.size() <= 0)
840 {
841 logger().error("Solve the problem first!");
842 return;
843 }
844
845 result.clear();
846
847 assert(!is_problem_scalar);
848
849 int index = 0;
850
851 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
852 std::vector<std::pair<std::string, Eigen::MatrixXd>> tmp_s;
853
854 for (int i = 0; i < int(bases.size()); ++i)
855 {
856 if (boundary_only && mesh.is_volume() && !mesh.is_boundary_element(i))
857 continue;
858
859 const ElementBases &bs = bases[i];
860 const ElementBases &gbs = gbases[i];
861 Eigen::MatrixXd local_pts;
862
863 if (use_sampler)
864 {
865 if (mesh.is_simplex(i))
866 local_pts = sampler.simplex_points();
867 else if (mesh.is_cube(i))
868 local_pts = sampler.cube_points();
869 else if (mesh.is_prism(i))
870 local_pts = sampler.prism_points();
871 else
872 {
873 if (mesh.is_volume())
874 sampler.sample_polyhedron(polys_3d.at(i).first, polys_3d.at(i).second, local_pts, vis_faces_poly, vis_edges_poly);
875 else
876 sampler.sample_polygon(polys.at(i), local_pts, vis_faces_poly, vis_edges_poly);
877 }
878 }
879 else
880 {
881 if (mesh.is_volume())
882 {
883 if (mesh.is_simplex(i))
884 autogen::p_nodes_3d(disc_orders(i), local_pts);
885 else if (mesh.is_cube(i))
886 autogen::q_nodes_3d(disc_orders(i), local_pts);
887 else if (mesh.is_prism(i))
888 autogen::prism_nodes_3d(disc_orders(i), disc_ordersq(i), local_pts);
889 else
890 continue;
891 }
892 else
893 {
894 if (mesh.is_simplex(i))
895 autogen::p_nodes_2d(disc_orders(i), local_pts);
896 else if (mesh.is_cube(i))
897 autogen::q_nodes_2d(disc_orders(i), local_pts);
898 else
899 continue;
900 }
901 }
902
903 assembler.compute_scalar_value(OutputData(t, i, bs, gbs, local_pts, fun), tmp_s);
904
905 if (result.empty())
906 {
907 result.resize(tmp_s.size());
908 for (int k = 0; k < tmp_s.size(); ++k)
909 {
910 result[k].first = tmp_s[k].first;
911 result[k].second.resize(n_points, 1);
912 }
913 }
914
915 for (int k = 0; k < tmp_s.size(); ++k)
916 {
917 assert(local_pts.rows() == tmp_s[k].second.rows());
918 result[k].second.block(index, 0, tmp_s[k].second.rows(), 1) = tmp_s[k].second;
919 }
920 index += local_pts.rows();
921 }
922 }
923
924 void Evaluator::compute_tensor_value(
925 const mesh::Mesh &mesh,
926 const bool is_problem_scalar,
927 const std::vector<basis::ElementBases> &bases,
928 const std::vector<basis::ElementBases> &gbases,
929 const Eigen::VectorXi &disc_orders,
930 const Eigen::VectorXi &disc_ordersq,
931 const std::map<int, Eigen::MatrixXd> &polys,
932 const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d,
933 const assembler::Assembler &assembler,
934 const utils::RefElementSampler &sampler,
935 const int n_points,
936 const Eigen::MatrixXd &fun,
937 const double t,
938 std::vector<assembler::Assembler::NamedMatrix> &result,
939 const bool use_sampler,
940 const bool boundary_only)
941 {
942 if (fun.size() <= 0)
943 {
944 logger().error("Solve the problem first!");
945 return;
946 }
947
948 result.clear();
949
950 const int actual_dim = mesh.dimension();
951 assert(!is_problem_scalar);
952
953 int index = 0;
954
955 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
956 std::vector<std::pair<std::string, Eigen::MatrixXd>> tmp_t;
957
958 for (int i = 0; i < int(bases.size()); ++i)
959 {
960 if (boundary_only && mesh.is_volume() && !mesh.is_boundary_element(i))
961 continue;
962
963 const ElementBases &bs = bases[i];
964 const ElementBases &gbs = gbases[i];
965 Eigen::MatrixXd local_pts;
966
967 if (use_sampler)
968 {
969 if (mesh.is_simplex(i))
970 local_pts = sampler.simplex_points();
971 else if (mesh.is_cube(i))
972 local_pts = sampler.cube_points();
973 else if (mesh.is_prism(i))
974 local_pts = sampler.prism_points();
975 else
976 {
977 if (mesh.is_volume())
978 sampler.sample_polyhedron(polys_3d.at(i).first, polys_3d.at(i).second, local_pts, vis_faces_poly, vis_edges_poly);
979 else
980 sampler.sample_polygon(polys.at(i), local_pts, vis_faces_poly, vis_edges_poly);
981 }
982 }
983 else
984 {
985 if (mesh.is_volume())
986 {
987 if (mesh.is_simplex(i))
988 autogen::p_nodes_3d(disc_orders(i), local_pts);
989 else if (mesh.is_cube(i))
990 autogen::q_nodes_3d(disc_orders(i), local_pts);
991 else if (mesh.is_prism(i))
992 autogen::prism_nodes_3d(disc_orders(i), disc_ordersq(i), local_pts);
993 else
994 continue;
995 }
996 else
997 {
998 if (mesh.is_simplex(i))
999 autogen::p_nodes_2d(disc_orders(i), local_pts);
1000 else if (mesh.is_cube(i))
1001 autogen::q_nodes_2d(disc_orders(i), local_pts);
1002 else
1003 continue;
1004 }
1005 }
1006
1007 assembler.compute_tensor_value(OutputData(t, i, bs, gbs, local_pts, fun), tmp_t);
1008
1009 if (result.empty())
1010 {
1011 result.resize(tmp_t.size());
1012 for (int k = 0; k < tmp_t.size(); ++k)
1013 {
1014 result[k].first = tmp_t[k].first;
1015 result[k].second.resize(n_points, actual_dim * actual_dim);
1016 }
1017 }
1018
1019 for (int k = 0; k < tmp_t.size(); ++k)
1020 {
1021 assert(local_pts.rows() == tmp_t[k].second.rows());
1022 result[k].second.block(index, 0, tmp_t[k].second.rows(), tmp_t[k].second.cols()) = tmp_t[k].second;
1023 }
1024 index += local_pts.rows();
1025 }
1026 }
1027
1028 Eigen::MatrixXd Evaluator::get_bases_position(
1029 const int n_bases,
1030 const std::shared_ptr<mesh::MeshNodes> mesh_nodes)
1031 {
1032 Eigen::MatrixXd func;
1033 func.setZero(n_bases, mesh_nodes->node_position(0).size());
1034
1035 for (int i = 0; i < n_bases; i++)
1036 func.row(i) = mesh_nodes->node_position(i);
1037
1038 return func;
1039 }
1040
1041 Eigen::MatrixXd Evaluator::generate_linear_field(
1042 const int n_bases,
1043 const std::shared_ptr<mesh::MeshNodes> mesh_nodes,
1044 const Eigen::MatrixXd &grad)
1045 {
1046 return utils::flatten(get_bases_position(n_bases, mesh_nodes) * grad.transpose());
1047 }
1048
1049 Eigen::VectorXd Evaluator::integrate_function(
1050 const std::vector<basis::ElementBases> &bases,
1051 const std::vector<basis::ElementBases> &gbases,
1052 const Eigen::MatrixXd &fun,
1053 const int dim,
1054 const int actual_dim)
1055 {
1056 Eigen::VectorXd result;
1057 result.setZero(actual_dim);
1058 for (int e = 0; e < bases.size(); ++e)
1059 {
1060 ElementAssemblyValues vals;
1061 vals.compute(e, dim == 3, bases[e], gbases[e]);
1062
1063 Eigen::MatrixXd u, grad_u;
1064 io::Evaluator::interpolate_at_local_vals(e, dim, actual_dim, vals, fun, u, grad_u);
1065 const quadrature::Quadrature &quadrature = vals.quadrature;
1066 Eigen::VectorXd da = vals.det.array() * quadrature.weights.array();
1067 result += u.transpose() * da;
1068 }
1069
1070 return result;
1071 }
1072} // 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:1041
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:634
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:1028
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:821
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:924
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:734
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:460
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:437
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:332
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:1049
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