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