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
11#include <polyfem/utils/BoundarySampler.hpp>
12
13#include <polyfem/autogen/auto_p_bases.hpp>
14#include <polyfem/autogen/auto_q_bases.hpp>
15
16#include <polyfem/utils/Logger.hpp>
17
18#include <igl/AABB.h>
19#include <igl/per_face_normals.h>
20
21namespace polyfem::io
22{
23 using namespace mesh;
24 using namespace assembler;
25 using namespace basis;
26
27 namespace
28 {
29 void flattened_tensor_coeffs(const Eigen::MatrixXd &S, Eigen::MatrixXd &X)
30 {
31 if (S.cols() == 4)
32 {
33 X.resize(S.rows(), 3);
34 X.col(0) = S.col(0);
35 X.col(1) = S.col(3);
36 X.col(2) = S.col(1);
37 }
38 else if (S.cols() == 9)
39 {
40 // [S11, S22, S33, S12, S13, S23]
41 X.resize(S.rows(), 6);
42 X.col(0) = S.col(0);
43 X.col(1) = S.col(4);
44 X.col(2) = S.col(8);
45 X.col(3) = S.col(1);
46 X.col(4) = S.col(2);
47 X.col(5) = S.col(5);
48 }
49 else
50 {
51 logger().error("Invalid tensor dimensions.");
52 }
53 }
54 } // namespace
55
56 void Evaluator::get_sidesets(
57 const mesh::Mesh &mesh,
58 Eigen::MatrixXd &pts,
59 Eigen::MatrixXi &faces,
60 Eigen::MatrixXd &sidesets)
61 {
62 // if (!mesh)
63 // {
64 // logger().error("Load the mesh first!");
65 // return;
66 // }
67
68 if (mesh.is_volume())
69 {
70 const Mesh3D &tmp_mesh = dynamic_cast<const Mesh3D &>(mesh);
71 int n_pts = 0;
72 int n_faces = 0;
73 for (int f = 0; f < tmp_mesh.n_faces(); ++f)
74 {
75 if (tmp_mesh.get_boundary_id(f) > 0)
76 {
77 n_pts += tmp_mesh.n_face_vertices(f) + 1;
78 n_faces += tmp_mesh.n_face_vertices(f);
79 }
80 }
81
82 pts.resize(n_pts, 3);
83 faces.resize(n_faces, 3);
84 sidesets.resize(n_pts, 1);
85
86 n_pts = 0;
87 n_faces = 0;
88 for (int f = 0; f < tmp_mesh.n_faces(); ++f)
89 {
90 const int sideset = tmp_mesh.get_boundary_id(f);
91 if (sideset > 0)
92 {
93 const int n_face_vertices = tmp_mesh.n_face_vertices(f);
94
95 for (int i = 0; i < n_face_vertices; ++i)
96 {
97 if (n_face_vertices == 3)
98 faces.row(n_faces) << ((i + 1) % n_face_vertices + n_pts), (i + n_pts), (n_pts + n_face_vertices);
99 else
100 faces.row(n_faces) << (i + n_pts), ((i + 1) % n_face_vertices + n_pts), (n_pts + n_face_vertices);
101 ++n_faces;
102 }
103
104 for (int i = 0; i < n_face_vertices; ++i)
105 {
106 pts.row(n_pts) = tmp_mesh.point(tmp_mesh.face_vertex(f, i));
107 sidesets(n_pts) = sideset;
108
109 ++n_pts;
110 }
111
112 pts.row(n_pts) = tmp_mesh.face_barycenter(f);
113 sidesets(n_pts) = sideset;
114 ++n_pts;
115 }
116 }
117 }
118 else
119 {
120 const Mesh2D &tmp_mesh = dynamic_cast<const Mesh2D &>(mesh);
121 int n_siteset = 0;
122 for (int e = 0; e < tmp_mesh.n_edges(); ++e)
123 {
124 if (tmp_mesh.get_boundary_id(e) > 0)
125 ++n_siteset;
126 }
127
128 pts.resize(n_siteset * 2, 2);
129 faces.resize(n_siteset, 2);
130 sidesets.resize(n_siteset, 1);
131
132 n_siteset = 0;
133 for (int e = 0; e < tmp_mesh.n_edges(); ++e)
134 {
135 const int sideset = tmp_mesh.get_boundary_id(e);
136 if (sideset > 0)
137 {
138 pts.row(2 * n_siteset) = tmp_mesh.point(tmp_mesh.edge_vertex(e, 0));
139 pts.row(2 * n_siteset + 1) = tmp_mesh.point(tmp_mesh.edge_vertex(e, 1));
140 faces.row(n_siteset) << 2 * n_siteset, 2 * n_siteset + 1;
141 sidesets(n_siteset) = sideset;
142 ++n_siteset;
143 }
144 }
145
146 pts.conservativeResize(n_siteset * 2, 3);
147 pts.col(2).setZero();
148 }
149 }
150
151 void Evaluator::interpolate_boundary_function(
152 const mesh::Mesh &mesh,
153 const bool is_problem_scalar,
154 const std::vector<basis::ElementBases> &bases,
155 const std::vector<basis::ElementBases> &gbases,
156 const Eigen::MatrixXd &pts,
157 const Eigen::MatrixXi &faces,
158 const Eigen::MatrixXd &fun,
159 const bool compute_avg,
160 Eigen::MatrixXd &result)
161 {
162 if (fun.size() <= 0)
163 {
164 logger().error("Solve the problem first!");
165 return;
166 }
167 assert(mesh.is_volume());
168
169 const Mesh3D &mesh3d = dynamic_cast<const Mesh3D &>(mesh);
170
171 Eigen::MatrixXd points, uv;
172 Eigen::VectorXd weights;
173
174 int actual_dim = 1;
175 if (!is_problem_scalar)
176 actual_dim = 3;
177
178 igl::AABB<Eigen::MatrixXd, 3> tree;
179 tree.init(pts, faces);
180
181 result.resize(faces.rows(), actual_dim);
182 result.setConstant(std::numeric_limits<double>::quiet_NaN());
183
184 int counter = 0;
185
186 for (int e = 0; e < mesh3d.n_elements(); ++e)
187 {
188 const basis::ElementBases &gbs = gbases[e];
189 const basis::ElementBases &bs = bases[e];
190
191 for (int lf = 0; lf < mesh3d.n_cell_faces(e); ++lf)
192 {
193 const int face_id = mesh3d.cell_face(e, lf);
194 if (!mesh3d.is_boundary_face(face_id))
195 continue;
196
197 if (mesh3d.is_simplex(e))
198 utils::BoundarySampler::quadrature_for_tri_face(lf, 4, face_id, mesh3d, uv, points, weights);
199 else if (mesh3d.is_cube(e))
200 utils::BoundarySampler::quadrature_for_quad_face(lf, 4, face_id, mesh3d, uv, points, weights);
201 else
202 assert(false);
203
204 ElementAssemblyValues vals;
205 vals.compute(e, true, points, bs, gbs);
206 RowVectorNd loc_val(actual_dim);
207 loc_val.setZero();
208
209 // UIEvaluator::ui_state().debug_data().add_points(vals.val, Eigen::RowVector3d(1,0,0));
210
211 // const auto nodes = bs.local_nodes_for_primitive(face_id, mesh3d);
212
213 // for(long n = 0; n < nodes.size(); ++n)
214 for (size_t j = 0; j < bs.bases.size(); ++j)
215 {
216 // const auto &b = bs.bases[nodes(n)];
217 // const AssemblyValues &v = vals.basis_values[nodes(n)];
218 const AssemblyValues &v = vals.basis_values[j];
219 for (int d = 0; d < actual_dim; ++d)
220 {
221 for (size_t g = 0; g < v.global.size(); ++g)
222 {
223 loc_val(d) += (v.global[g].val * v.val.array() * fun(v.global[g].index * actual_dim + d) * weights.array()).sum();
224 }
225 }
226 }
227
228 int I;
229 Eigen::RowVector3d C;
230 const Eigen::RowVector3d bary = mesh3d.face_barycenter(face_id);
231
232 const double dist = tree.squared_distance(pts, faces, bary, I, C);
233 assert(dist < 1e-16);
234
235 assert(std::isnan(result(I, 0)));
236 if (compute_avg)
237 result.row(I) = loc_val / weights.sum();
238 else
239 result.row(I) = loc_val;
240 ++counter;
241 }
242 }
243
244 assert(counter == result.rows());
245 }
246
247 void Evaluator::interpolate_boundary_function_at_vertices(
248 const mesh::Mesh &mesh,
249 const bool is_problem_scalar,
250 const std::vector<basis::ElementBases> &bases,
251 const std::vector<basis::ElementBases> &gbases,
252 const Eigen::MatrixXd &pts,
253 const Eigen::MatrixXi &faces,
254 const Eigen::MatrixXd &fun,
255 Eigen::MatrixXd &result)
256 {
257 if (fun.size() <= 0)
258 {
259 logger().error("Solve the problem first!");
260 return;
261 }
262 if (!mesh.is_volume())
263 {
264 logger().error("This function works only on volumetric meshes!");
265 return;
266 }
267
268 assert(mesh.is_volume());
269
270 const Mesh3D &mesh3d = dynamic_cast<const Mesh3D &>(mesh);
271
272 Eigen::MatrixXd points;
273
274 int actual_dim = 1;
275 if (!is_problem_scalar)
276 actual_dim = 3;
277
278 igl::AABB<Eigen::MatrixXd, 3> tree;
279 tree.init(pts, faces);
280
281 result.resize(pts.rows(), actual_dim);
282 result.setZero();
283
284 for (int e = 0; e < mesh3d.n_elements(); ++e)
285 {
286 const ElementBases &gbs = gbases[e];
287 const ElementBases &bs = bases[e];
288
289 for (int lf = 0; lf < mesh3d.n_cell_faces(e); ++lf)
290 {
291 const int face_id = mesh3d.cell_face(e, lf);
292 // if (!mesh3d.is_boundary_face(face_id))
293 // continue;
294 int I = -1;
295 Eigen::RowVector3d C;
296 const Eigen::RowVector3d bary = mesh3d.face_barycenter(face_id);
297
298 const double dist = tree.squared_distance(pts, faces, bary, I, C);
299 if (dist > 1e-15)
300 continue;
301
302 if (mesh3d.is_simplex(e))
303 autogen::p_nodes_3d(1, points);
304 else if (mesh3d.is_cube(e))
305 autogen::q_nodes_3d(1, points);
306 else
307 assert(false);
308
309 ElementAssemblyValues vals;
310 vals.compute(e, true, points, bs, gbs);
311 Eigen::MatrixXd loc_val(points.rows(), actual_dim);
312 loc_val.setZero();
313
314 // UIEvaluator::ui_state().debug_data().add_points(vals.val, Eigen::RowVector3d(1,0,0));
315
316 for (size_t j = 0; j < bs.bases.size(); ++j)
317 {
318 const Basis &b = bs.bases[j];
319 const AssemblyValues &v = vals.basis_values[j];
320
321 for (int d = 0; d < actual_dim; ++d)
322 {
323 for (size_t ii = 0; ii < b.global().size(); ++ii)
324 loc_val.col(d) += b.global()[ii].val * v.val * fun(b.global()[ii].index * actual_dim + d);
325 }
326 }
327
328 for (int lv_id = 0; lv_id < faces.cols(); ++lv_id)
329 {
330 const int v_id = faces(I, lv_id);
331 const auto p = pts.row(v_id);
332 const auto &mapped = vals.val;
333
334 bool found = false;
335
336 for (int n = 0; n < mapped.rows(); ++n)
337 {
338 if ((p - mapped.row(n)).norm() < 1e-10)
339 {
340 result.row(v_id) = loc_val.row(n);
341 found = true;
342 break;
343 }
344 }
345
346 assert(found);
347 }
348 }
349 }
350 }
351
352 void Evaluator::interpolate_boundary_tensor_function(
353 const mesh::Mesh &mesh,
354 const bool is_problem_scalar,
355 const std::vector<basis::ElementBases> &bases,
356 const std::vector<basis::ElementBases> &gbases,
357 const assembler::Assembler &assembler,
358 const Eigen::MatrixXd &pts,
359 const Eigen::MatrixXi &faces,
360 const Eigen::MatrixXd &fun,
361 const double t,
362 const bool compute_avg,
363 Eigen::MatrixXd &result,
364 Eigen::MatrixXd &stresses,
365 Eigen::MatrixXd &mises,
366 const bool skip_orientation)
367 {
368 interpolate_boundary_tensor_function(
369 mesh, is_problem_scalar, bases, gbases,
370 assembler,
371 pts, faces, fun, Eigen::MatrixXd::Zero(pts.rows(), pts.cols()), t,
372 compute_avg, result, stresses, mises, skip_orientation);
373 }
374
375 void Evaluator::interpolate_boundary_tensor_function(
376 const mesh::Mesh &mesh,
377 const bool is_problem_scalar,
378 const std::vector<basis::ElementBases> &bases,
379 const std::vector<basis::ElementBases> &gbases,
380 const assembler::Assembler &assembler,
381 const Eigen::MatrixXd &pts,
382 const Eigen::MatrixXi &faces,
383 const Eigen::MatrixXd &fun,
384 const Eigen::MatrixXd &disp,
385 const double t,
386 const bool compute_avg,
387 Eigen::MatrixXd &result,
388 Eigen::MatrixXd &stresses,
389 Eigen::MatrixXd &mises,
390 bool skip_orientation)
391 {
392 if (fun.size() <= 0)
393 {
394 logger().error("Solve the problem first!");
395 return;
396 }
397 if (disp.size() <= 0)
398 {
399 logger().error("Solve the problem first!");
400 return;
401 }
402 if (!mesh.is_volume())
403 {
404 logger().error("This function works only on volumetric meshes!");
405 return;
406 }
407 if (is_problem_scalar)
408 {
409 logger().error("Define a tensor problem!");
410 return;
411 }
412
413 std::vector<std::pair<std::string, Eigen::MatrixXd>> tmp_t, tmp_s;
414
415 assert(mesh.is_volume());
416 assert(!is_problem_scalar);
417
418 const Mesh3D &mesh3d = dynamic_cast<const Mesh3D &>(mesh);
419
420 Eigen::MatrixXd normals;
421 igl::per_face_normals((pts + disp).eval(), faces, normals);
422
423 Eigen::MatrixXd points, uv, tmp_n, loc_v;
424 Eigen::VectorXd weights;
425
426 const int actual_dim = 3;
427
428 igl::AABB<Eigen::MatrixXd, 3> tree;
429 tree.init(pts, faces);
430
431 result.resize(faces.rows(), actual_dim);
432 result.setConstant(std::numeric_limits<double>::quiet_NaN());
433
434 stresses.resize(faces.rows(), actual_dim * actual_dim);
435 stresses.setConstant(std::numeric_limits<double>::quiet_NaN());
436
437 mises.resize(faces.rows(), 1);
438 mises.setConstant(std::numeric_limits<double>::quiet_NaN());
439
440 int counter = 0;
441
442 for (int e = 0; e < mesh3d.n_elements(); ++e)
443 {
444 const ElementBases &gbs = gbases[e];
445 const ElementBases &bs = bases[e];
446
447 for (int lf = 0; lf < mesh3d.n_cell_faces(e); ++lf)
448 {
449 const int face_id = mesh3d.cell_face(e, lf);
450 // if (!mesh3d.is_boundary_face(face_id))
451 // continue;
452
453 int I;
454 Eigen::RowVector3d C;
455 const Eigen::RowVector3d bary = mesh3d.face_barycenter(face_id);
456
457 const double dist = tree.squared_distance(pts, faces, bary, I, C);
458 if (dist > 1e-15)
459 continue;
460
461 int lfid = 0;
462 for (; lfid < mesh3d.n_cell_faces(e); ++lfid)
463 {
464 if (mesh.is_simplex(e))
465 loc_v = utils::BoundarySampler::tet_local_node_coordinates_from_face(lfid);
466 else if (mesh.is_cube(e))
467 loc_v = utils::BoundarySampler::hex_local_node_coordinates_from_face(lfid);
468 else
469 assert(false);
470
471 ElementAssemblyValues vals;
472 vals.compute(e, true, loc_v, bs, gbs);
473
474 int count = 0;
475
476 for (int lv_id = 0; lv_id < faces.cols(); ++lv_id)
477 {
478 const int v_id = faces(I, lv_id);
479 const auto p = pts.row(v_id);
480 const auto &mapped = vals.val;
481 assert(mapped.rows() == faces.cols());
482
483 for (int n = 0; n < mapped.rows(); ++n)
484 {
485 if ((p - mapped.row(n)).norm() < 1e-10)
486 {
487 count++;
488 break;
489 }
490 }
491 }
492
493 if (count == faces.cols())
494 break;
495 }
496 assert(lfid < mesh3d.n_cell_faces(e));
497
498 if (mesh.is_simplex(e))
499 {
500 utils::BoundarySampler::quadrature_for_tri_face(lfid, 4, face_id, mesh3d, uv, points, weights);
501 utils::BoundarySampler::normal_for_tri_face(lfid, tmp_n);
502 }
503 else if (mesh.is_cube(e))
504 {
505 utils::BoundarySampler::quadrature_for_quad_face(lfid, 4, face_id, mesh3d, uv, points, weights);
506 utils::BoundarySampler::normal_for_quad_face(lfid, tmp_n);
507 }
508 else
509 assert(false);
510
511 Eigen::RowVector3d tet_n;
512 tet_n.setZero();
513 ElementAssemblyValues vals;
514 vals.compute(e, true, points, bs, gbs);
515 for (int n = 0; n < vals.jac_it.size(); ++n)
516 {
517 Eigen::RowVector3d tmp = tmp_n * vals.jac_it[n];
518 tmp.normalize();
519 tet_n += tmp;
520 }
521
522 assembler.compute_scalar_value(OutputData(t, e, bs, gbs, points, fun), tmp_s);
523 assembler.compute_tensor_value(OutputData(t, e, bs, gbs, points, fun), tmp_t);
524
525 Eigen::MatrixXd loc_val = tmp_t[0].second, local_mises = tmp_s[0].second;
526 Eigen::VectorXd tmp(loc_val.cols());
527 const double tmp_mises = (local_mises.array() * weights.array()).sum();
528
529 for (int d = 0; d < loc_val.cols(); ++d)
530 tmp(d) = (loc_val.col(d).array() * weights.array()).sum();
531 const Eigen::MatrixXd tensor = Eigen::Map<Eigen::MatrixXd>(tmp.data(), 3, 3);
532
533 const Eigen::RowVector3d tmpn = normals.row(I);
534 const Eigen::RowVector3d tmptf = tmpn * tensor;
535 if (skip_orientation || tmpn.dot(tet_n) > 0)
536 {
537 assert(std::isnan(result(I, 0)));
538 assert(std::isnan(stresses(I, 0)));
539 assert(std::isnan(mises(I)));
540
541 result.row(I) = tmptf;
542 stresses.row(I) = tmp;
543 mises(I) = tmp_mises;
544
545 if (compute_avg)
546 {
547 result.row(I) /= weights.sum();
548 stresses.row(I) /= weights.sum();
549 mises(I) /= weights.sum();
550 }
551 ++counter;
552 }
553 }
554 }
555
556 assert(counter == result.rows());
557 }
558
559 void Evaluator::average_grad_based_function(
560 const mesh::Mesh &mesh,
561 const bool is_problem_scalar,
562 const int n_bases,
563 const std::vector<basis::ElementBases> &bases,
564 const std::vector<basis::ElementBases> &gbases,
565 const Eigen::VectorXi &disc_orders,
566 const std::map<int, Eigen::MatrixXd> &polys,
567 const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d,
568 const assembler::Assembler &assembler,
569 const utils::RefElementSampler &sampler,
570 const double t,
571 const int n_points,
572 const Eigen::MatrixXd &fun,
573 std::vector<assembler::Assembler::NamedMatrix> &result_scalar,
574 std::vector<assembler::Assembler::NamedMatrix> &result_tensor,
575 const bool use_sampler,
576 const bool boundary_only)
577 {
578 result_scalar.clear();
579 result_tensor.clear();
580
581 if (fun.size() <= 0)
582 {
583 logger().error("Solve the problem first!");
584 return;
585 }
586 if (is_problem_scalar)
587 {
588 logger().error("Define a tensor problem!");
589 return;
590 }
591
592 assert(!is_problem_scalar);
593 const int actual_dim = mesh.dimension();
594
595 std::vector<Eigen::MatrixXd> avg_scalar;
596
597 Eigen::MatrixXd areas(n_bases, 1);
598 areas.setZero();
599
600 std::vector<std::pair<std::string, Eigen::MatrixXd>> tmp_s;
601 Eigen::MatrixXd local_val;
602
603 ElementAssemblyValues vals;
604 for (int i = 0; i < int(bases.size()); ++i)
605 {
606 const ElementBases &bs = bases[i];
607 const ElementBases &gbs = gbases[i];
608 Eigen::MatrixXd local_pts;
609
610 if (mesh.is_simplex(i))
611 {
612 if (mesh.dimension() == 3)
613 autogen::p_nodes_3d(disc_orders(i), local_pts);
614 else
615 autogen::p_nodes_2d(disc_orders(i), local_pts);
616 }
617 else if (mesh.is_cube(i))
618 {
619 if (mesh.dimension() == 3)
620 autogen::q_nodes_3d(disc_orders(i), local_pts);
621 else
622 autogen::q_nodes_2d(disc_orders(i), local_pts);
623 }
624 else
625 {
626 // not supported for polys
627 continue;
628 }
629
630 vals.compute(i, actual_dim == 3, bases[i], gbases[i]);
631 const quadrature::Quadrature &quadrature = vals.quadrature;
632 const double area = (vals.det.array() * quadrature.weights.array()).sum();
633
634 assembler.compute_scalar_value(OutputData(t, i, bs, gbs, local_pts, fun), tmp_s);
635
636 // assembler.compute_tensor_value(i, bs, gbs, local_pts, fun, local_val);
637 // MatrixXd avg_tensor(n_points * actual_dim*actual_dim, 1);
638
639 for (size_t j = 0; j < bs.bases.size(); ++j)
640 {
641 const Basis &b = bs.bases[j];
642 if (b.global().size() > 1)
643 continue;
644
645 auto &global = b.global().front();
646 areas(global.index) += area;
647 }
648
649 if (avg_scalar.empty())
650 {
651 avg_scalar.resize(tmp_s.size());
652 for (auto &m : avg_scalar)
653 {
654 m.resize(n_bases, 1);
655 m.setZero();
656 }
657 }
658
659 for (int k = 0; k < tmp_s.size(); ++k)
660 {
661 local_val = tmp_s[k].second;
662
663 for (size_t j = 0; j < bs.bases.size(); ++j)
664 {
665 const Basis &b = bs.bases[j];
666 if (b.global().size() > 1)
667 continue;
668
669 auto &global = b.global().front();
670 avg_scalar[k](global.index) += local_val(j) * area;
671 }
672 }
673 }
674
675 for (auto &m : avg_scalar)
676 {
677 m.array() /= areas.array();
678 }
679
680 result_scalar.resize(tmp_s.size());
681 for (int k = 0; k < tmp_s.size(); ++k)
682 {
683 result_scalar[k].first = tmp_s[k].first;
684 interpolate_function(mesh, 1, bases, disc_orders, polys, polys_3d, sampler, n_points,
685 avg_scalar[k], result_scalar[k].second, use_sampler, boundary_only);
686 }
687 // interpolate_function(n_points, actual_dim*actual_dim, bases, avg_tensor, result_tensor, boundary_only);
688 }
689
690 void Evaluator::compute_vertex_values(
691 const mesh::Mesh &mesh,
692 int actual_dim,
693 const std::vector<basis::ElementBases> &basis,
694 const utils::RefElementSampler &sampler,
695 const Eigen::MatrixXd &fun,
696 Eigen::MatrixXd &result)
697 {
698 // if (!mesh)
699 // {
700 // logger().error("Load the mesh first!");
701 // return;
702 // }
703 if (fun.size() <= 0)
704 {
705 logger().error("Solve the problem first!");
706 return;
707 }
708 if (!mesh.is_volume())
709 {
710 logger().error("This function works only on volumetric meshes!");
711 return;
712 }
713
714 const Mesh3D &mesh3d = dynamic_cast<const Mesh3D &>(mesh);
715
716 result.resize(mesh3d.n_vertices(), actual_dim);
717 result.setZero();
718
719 // std::array<int, 8> get_ordered_vertices_from_hex(const int element_index) const;
720 // std::array<int, 4> get_ordered_vertices_from_tet(const int element_index) const;
721
722 std::vector<AssemblyValues> tmp;
723 std::vector<bool> marked(mesh3d.n_vertices(), false);
724 for (int i = 0; i < int(basis.size()); ++i)
725 {
726 const ElementBases &bs = basis[i];
727 Eigen::MatrixXd local_pts;
728 std::vector<int> vertices;
729
730 if (mesh.is_simplex(i))
731 {
732 local_pts = sampler.simplex_corners();
733 auto vtx = mesh3d.get_ordered_vertices_from_tet(i);
734 vertices.assign(vtx.begin(), vtx.end());
735 }
736 else if (mesh.is_cube(i))
737 {
738 local_pts = sampler.cube_corners();
739 auto vtx = mesh3d.get_ordered_vertices_from_hex(i);
740 vertices.assign(vtx.begin(), vtx.end());
741 }
742 // TODO poly?
743 assert((int)vertices.size() == (int)local_pts.rows());
744
745 Eigen::MatrixXd local_res = Eigen::MatrixXd::Zero(local_pts.rows(), actual_dim);
746 bs.evaluate_bases(local_pts, tmp);
747 for (size_t j = 0; j < bs.bases.size(); ++j)
748 {
749 const Basis &b = bs.bases[j];
750
751 for (int d = 0; d < actual_dim; ++d)
752 {
753 for (size_t ii = 0; ii < b.global().size(); ++ii)
754 local_res.col(d) += b.global()[ii].val * tmp[j].val * fun(b.global()[ii].index * actual_dim + d);
755 }
756 }
757
758 for (size_t lv = 0; lv < vertices.size(); ++lv)
759 {
760 int v = vertices[lv];
761 if (marked[v])
762 {
763 assert((result.row(v) - local_res.row(lv)).norm() < 1e-6);
764 }
765 else
766 {
767 result.row(v) = local_res.row(lv);
768 marked[v] = true;
769 }
770 }
771 }
772 }
773
774 void Evaluator::compute_stress_at_quadrature_points(
775 const mesh::Mesh &mesh,
776 const bool is_problem_scalar,
777 const std::vector<basis::ElementBases> &bases,
778 const std::vector<basis::ElementBases> &gbases,
779 const Eigen::VectorXi &disc_orders,
780 const assembler::Assembler &assembler,
781 const Eigen::MatrixXd &fun,
782 const double t,
783 Eigen::MatrixXd &result,
784 Eigen::VectorXd &von_mises)
785 {
786 // if (!mesh)
787 // {
788 // logger().error("Load the mesh first!");
789 // return;
790 // }
791 if (fun.size() <= 0)
792 {
793 logger().error("Solve the problem first!");
794 return;
795 }
796 if (is_problem_scalar)
797 {
798 logger().error("Define a tensor problem!");
799 return;
800 }
801
802 const int actual_dim = mesh.dimension();
803 assert(!is_problem_scalar);
804
805 Eigen::MatrixXd local_val, local_stress, local_mises;
806
807 int num_quadr_pts = 0;
808 result.resize(disc_orders.sum(), actual_dim == 2 ? 3 : 6);
809 result.setZero();
810 von_mises.resize(disc_orders.sum(), 1);
811 von_mises.setZero();
812 for (int e = 0; e < mesh.n_elements(); ++e)
813 {
814 // Compute quadrature points for element
815 quadrature::Quadrature quadr;
816 if (mesh.is_simplex(e))
817 {
818 if (mesh.is_volume())
819 {
820 quadrature::TetQuadrature f;
821 f.get_quadrature(disc_orders(e), quadr);
822 }
823 else
824 {
825 quadrature::TriQuadrature f;
826 f.get_quadrature(disc_orders(e), quadr);
827 }
828 }
829 else if (mesh.is_cube(e))
830 {
831 if (mesh.is_volume())
832 {
833 quadrature::HexQuadrature f;
834 f.get_quadrature(disc_orders(e), quadr);
835 }
836 else
837 {
838 quadrature::QuadQuadrature f;
839 f.get_quadrature(disc_orders(e), quadr);
840 }
841 }
842 else
843 {
844 continue;
845 }
846
847 std::vector<std::pair<std::string, Eigen::MatrixXd>> tmp_s, tmp_t;
848
849 assembler.compute_scalar_value(OutputData(t, e, bases[e], gbases[e], quadr.points, fun), tmp_s);
850 assembler.compute_tensor_value(OutputData(t, e, bases[e], gbases[e], quadr.points, fun), tmp_t);
851
852 local_mises = tmp_s[0].second;
853 local_val = tmp_t[0].second;
854
855 if (num_quadr_pts + local_val.rows() >= result.rows())
856 {
857 result.conservativeResize(
858 std::max(num_quadr_pts + local_val.rows() + 1, 2 * result.rows()),
859 result.cols());
860 von_mises.conservativeResize(result.rows(), von_mises.cols());
861 }
862 flattened_tensor_coeffs(local_val, local_stress);
863 result.block(num_quadr_pts, 0, local_stress.rows(), local_stress.cols()) = local_stress;
864 von_mises.block(num_quadr_pts, 0, local_mises.rows(), local_mises.cols()) = local_mises;
865 num_quadr_pts += local_val.rows();
866 }
867 result.conservativeResize(num_quadr_pts, result.cols());
868 von_mises.conservativeResize(num_quadr_pts, von_mises.cols());
869 }
870
871 void Evaluator::interpolate_function(
872 const mesh::Mesh &mesh,
873 const bool is_problem_scalar,
874 const std::vector<basis::ElementBases> &bases,
875 const Eigen::VectorXi &disc_orders,
876 const std::map<int, Eigen::MatrixXd> &polys,
877 const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d,
878 const utils::RefElementSampler &sampler,
879 const int n_points,
880 const Eigen::MatrixXd &fun,
881 Eigen::MatrixXd &result,
882 const bool use_sampler,
883 const bool boundary_only)
884 {
885 int actual_dim = 1;
886 if (!is_problem_scalar)
887 actual_dim = mesh.dimension();
888 interpolate_function(mesh, actual_dim, bases, disc_orders,
889 polys, polys_3d, sampler, n_points,
890 fun, result, use_sampler, boundary_only);
891 }
892
893 void Evaluator::interpolate_function(
894 const mesh::Mesh &mesh,
895 const int actual_dim,
896 const std::vector<basis::ElementBases> &basis,
897 const Eigen::VectorXi &disc_orders,
898 const std::map<int, Eigen::MatrixXd> &polys,
899 const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d,
900 const utils::RefElementSampler &sampler,
901 const int n_points,
902 const Eigen::MatrixXd &fun,
903 Eigen::MatrixXd &result,
904 const bool use_sampler,
905 const bool boundary_only)
906 {
907 if (fun.size() <= 0)
908 {
909 logger().error("Solve the problem first!");
910 return;
911 }
912
913 std::vector<AssemblyValues> tmp;
914
915 result.resize(n_points, actual_dim);
916
917 int index = 0;
918
919 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
920
921 for (int i = 0; i < int(basis.size()); ++i)
922 {
923 const ElementBases &bs = basis[i];
924 Eigen::MatrixXd local_pts;
925
926 if (boundary_only && mesh.is_volume() && !mesh.is_boundary_element(i))
927 continue;
928
929 if (use_sampler)
930 {
931 if (mesh.is_simplex(i))
932 local_pts = sampler.simplex_points();
933 else if (mesh.is_cube(i))
934 local_pts = sampler.cube_points();
935 else
936 {
937 if (mesh.is_volume())
938 sampler.sample_polyhedron(polys_3d.at(i).first, polys_3d.at(i).second, local_pts, vis_faces_poly, vis_edges_poly);
939 else
940 sampler.sample_polygon(polys.at(i), local_pts, vis_faces_poly, vis_edges_poly);
941 }
942 }
943 else
944 {
945 if (mesh.is_volume())
946 {
947 if (mesh.is_simplex(i))
948 autogen::p_nodes_3d(disc_orders(i), local_pts);
949 else if (mesh.is_cube(i))
950 autogen::q_nodes_3d(disc_orders(i), local_pts);
951 else
952 continue;
953 }
954 else
955 {
956 if (mesh.is_simplex(i))
957 autogen::p_nodes_2d(disc_orders(i), local_pts);
958 else if (mesh.is_cube(i))
959 autogen::q_nodes_2d(disc_orders(i), local_pts);
960 else
961 continue;
962 }
963 }
964
965 Eigen::MatrixXd local_res = Eigen::MatrixXd::Zero(local_pts.rows(), actual_dim);
966 bs.evaluate_bases(local_pts, tmp);
967 for (size_t j = 0; j < bs.bases.size(); ++j)
968 {
969 const Basis &b = bs.bases[j];
970
971 for (int d = 0; d < actual_dim; ++d)
972 {
973 for (size_t ii = 0; ii < b.global().size(); ++ii)
974 local_res.col(d) += b.global()[ii].val * tmp[j].val * fun(b.global()[ii].index * actual_dim + d);
975 }
976 }
977
978 result.block(index, 0, local_res.rows(), actual_dim) = local_res;
979 index += local_res.rows();
980 }
981 }
982
983 void Evaluator::interpolate_at_local_vals(
984 const mesh::Mesh &mesh,
985 const bool is_problem_scalar,
986 const std::vector<basis::ElementBases> &bases,
987 const std::vector<basis::ElementBases> &gbases,
988 const int el_index,
989 const Eigen::MatrixXd &local_pts,
990 const Eigen::MatrixXd &fun,
991 Eigen::MatrixXd &result,
992 Eigen::MatrixXd &result_grad)
993 {
994 int actual_dim = 1;
995 if (!is_problem_scalar)
996 actual_dim = mesh.dimension();
997 interpolate_at_local_vals(mesh, actual_dim, bases, gbases, el_index,
998 local_pts, fun, result, result_grad);
999 }
1000
1001 void Evaluator::interpolate_at_local_vals(
1002 const mesh::Mesh &mesh,
1003 const int actual_dim,
1004 const std::vector<basis::ElementBases> &bases,
1005 const std::vector<basis::ElementBases> &gbases,
1006 const int el_index,
1007 const Eigen::MatrixXd &local_pts,
1008 const Eigen::MatrixXd &fun,
1009 Eigen::MatrixXd &result,
1010 Eigen::MatrixXd &result_grad)
1011 {
1012 if (fun.size() <= 0)
1013 {
1014 logger().error("Solve the problem first!");
1015 return;
1016 }
1017
1018 assert(local_pts.cols() == mesh.dimension());
1019 assert(fun.cols() == 1);
1020
1021 const ElementBases &gbs = gbases[el_index];
1022 const ElementBases &bs = bases[el_index];
1023
1024 ElementAssemblyValues vals;
1025 vals.compute(el_index, mesh.is_volume(), local_pts, bs, gbs);
1026
1027 result.resize(vals.val.rows(), actual_dim);
1028 result.setZero();
1029
1030 result_grad.resize(vals.val.rows(), mesh.dimension() * actual_dim);
1031 result_grad.setZero();
1032
1033 const int n_loc_bases = int(vals.basis_values.size());
1034
1035 for (int i = 0; i < n_loc_bases; ++i)
1036 {
1037 const auto &val = vals.basis_values[i];
1038
1039 for (size_t ii = 0; ii < val.global.size(); ++ii)
1040 {
1041 for (int d = 0; d < actual_dim; ++d)
1042 {
1043 result.col(d) += val.global[ii].val * fun(val.global[ii].index * actual_dim + d) * val.val;
1044 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;
1045 }
1046 }
1047 }
1048 }
1049
1050 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)
1051 {
1052 if (fun.size() <= 0)
1053 {
1054 logger().error("Solve the problem first!");
1055 return;
1056 }
1057
1058 assert(fun.cols() == 1);
1059
1060 result.resize(vals.val.rows(), actual_dim);
1061 result.setZero();
1062
1063 result_grad.resize(vals.val.rows(), dim * actual_dim);
1064 result_grad.setZero();
1065
1066 const int n_loc_bases = int(vals.basis_values.size());
1067
1068 for (int i = 0; i < n_loc_bases; ++i)
1069 {
1070 const auto &val = vals.basis_values[i];
1071
1072 for (size_t ii = 0; ii < val.global.size(); ++ii)
1073 {
1074 for (int d = 0; d < actual_dim; ++d)
1075 {
1076 result.col(d) += val.global[ii].val * fun(val.global[ii].index * actual_dim + d) * val.val;
1077 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;
1078 }
1079 }
1080 }
1081 }
1082
1083 bool Evaluator::check_scalar_value(
1084 const mesh::Mesh &mesh,
1085 const bool is_problem_scalar,
1086 const std::vector<basis::ElementBases> &bases,
1087 const std::vector<basis::ElementBases> &gbases,
1088 const Eigen::VectorXi &disc_orders,
1089 const std::map<int, Eigen::MatrixXd> &polys,
1090 const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d,
1091 const assembler::Assembler &assembler,
1092 const utils::RefElementSampler &sampler,
1093 const Eigen::MatrixXd &fun,
1094 const double t,
1095 const bool use_sampler,
1096 const bool boundary_only)
1097 {
1098 if (fun.size() <= 0)
1099 {
1100 logger().error("Solve the problem first!");
1101 return true;
1102 }
1103
1104 assert(!is_problem_scalar);
1105
1106 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
1107
1108 std::vector<std::pair<std::string, Eigen::MatrixXd>> tmp_s;
1109
1110 for (int i = 0; i < int(bases.size()); ++i)
1111 {
1112 if (boundary_only && mesh.is_volume() && !mesh.is_boundary_element(i))
1113 continue;
1114
1115 const ElementBases &bs = bases[i];
1116 const ElementBases &gbs = gbases[i];
1117 Eigen::MatrixXd local_pts;
1118
1119 if (use_sampler)
1120 {
1121 if (mesh.is_simplex(i))
1122 local_pts = sampler.simplex_points();
1123 else if (mesh.is_cube(i))
1124 local_pts = sampler.cube_points();
1125 else
1126 {
1127 if (mesh.is_volume())
1128 sampler.sample_polyhedron(polys_3d.at(i).first, polys_3d.at(i).second, local_pts, vis_faces_poly, vis_edges_poly);
1129 else
1130 sampler.sample_polygon(polys.at(i), local_pts, vis_faces_poly, vis_edges_poly);
1131 }
1132 }
1133 else
1134 {
1135 if (mesh.is_volume())
1136 {
1137 if (mesh.is_simplex(i))
1138 autogen::p_nodes_3d(disc_orders(i), local_pts);
1139 else if (mesh.is_cube(i))
1140 autogen::q_nodes_3d(disc_orders(i), local_pts);
1141 else
1142 continue;
1143 }
1144 else
1145 {
1146 if (mesh.is_simplex(i))
1147 autogen::p_nodes_2d(disc_orders(i), local_pts);
1148 else if (mesh.is_cube(i))
1149 autogen::q_nodes_2d(disc_orders(i), local_pts);
1150 else
1151 continue;
1152 }
1153 }
1154
1155 assembler.compute_scalar_value(OutputData(t, i, bs, gbs, local_pts, fun), tmp_s);
1156
1157 for (const auto &s : tmp_s)
1158 if (std::isnan(s.second.norm()))
1159 return false;
1160 }
1161
1162 return true;
1163 }
1164
1165 void Evaluator::compute_scalar_value(
1166 const mesh::Mesh &mesh,
1167 const bool is_problem_scalar,
1168 const std::vector<basis::ElementBases> &bases,
1169 const std::vector<basis::ElementBases> &gbases,
1170 const Eigen::VectorXi &disc_orders,
1171 const std::map<int, Eigen::MatrixXd> &polys,
1172 const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d,
1173 const assembler::Assembler &assembler,
1174 const utils::RefElementSampler &sampler,
1175 const int n_points,
1176 const Eigen::MatrixXd &fun,
1177 const double t,
1178 std::vector<assembler::Assembler::NamedMatrix> &result,
1179 const bool use_sampler,
1180 const bool boundary_only)
1181 {
1182 if (fun.size() <= 0)
1183 {
1184 logger().error("Solve the problem first!");
1185 return;
1186 }
1187
1188 result.clear();
1189
1190 assert(!is_problem_scalar);
1191
1192 int index = 0;
1193
1194 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
1195 std::vector<std::pair<std::string, Eigen::MatrixXd>> tmp_s;
1196
1197 for (int i = 0; i < int(bases.size()); ++i)
1198 {
1199 if (boundary_only && mesh.is_volume() && !mesh.is_boundary_element(i))
1200 continue;
1201
1202 const ElementBases &bs = bases[i];
1203 const ElementBases &gbs = gbases[i];
1204 Eigen::MatrixXd local_pts;
1205
1206 if (use_sampler)
1207 {
1208 if (mesh.is_simplex(i))
1209 local_pts = sampler.simplex_points();
1210 else if (mesh.is_cube(i))
1211 local_pts = sampler.cube_points();
1212 else
1213 {
1214 if (mesh.is_volume())
1215 sampler.sample_polyhedron(polys_3d.at(i).first, polys_3d.at(i).second, local_pts, vis_faces_poly, vis_edges_poly);
1216 else
1217 sampler.sample_polygon(polys.at(i), local_pts, vis_faces_poly, vis_edges_poly);
1218 }
1219 }
1220 else
1221 {
1222 if (mesh.is_volume())
1223 {
1224 if (mesh.is_simplex(i))
1225 autogen::p_nodes_3d(disc_orders(i), local_pts);
1226 else if (mesh.is_cube(i))
1227 autogen::q_nodes_3d(disc_orders(i), local_pts);
1228 else
1229 continue;
1230 }
1231 else
1232 {
1233 if (mesh.is_simplex(i))
1234 autogen::p_nodes_2d(disc_orders(i), local_pts);
1235 else if (mesh.is_cube(i))
1236 autogen::q_nodes_2d(disc_orders(i), local_pts);
1237 else
1238 continue;
1239 }
1240 }
1241
1242 assembler.compute_scalar_value(OutputData(t, i, bs, gbs, local_pts, fun), tmp_s);
1243
1244 if (result.empty())
1245 {
1246 result.resize(tmp_s.size());
1247 for (int k = 0; k < tmp_s.size(); ++k)
1248 {
1249 result[k].first = tmp_s[k].first;
1250 result[k].second.resize(n_points, 1);
1251 }
1252 }
1253
1254 for (int k = 0; k < tmp_s.size(); ++k)
1255 {
1256 assert(local_pts.rows() == tmp_s[k].second.rows());
1257 result[k].second.block(index, 0, tmp_s[k].second.rows(), 1) = tmp_s[k].second;
1258 }
1259 index += local_pts.rows();
1260 }
1261 }
1262
1263 void Evaluator::compute_tensor_value(
1264 const mesh::Mesh &mesh,
1265 const bool is_problem_scalar,
1266 const std::vector<basis::ElementBases> &bases,
1267 const std::vector<basis::ElementBases> &gbases,
1268 const Eigen::VectorXi &disc_orders,
1269 const std::map<int, Eigen::MatrixXd> &polys,
1270 const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d,
1271 const assembler::Assembler &assembler,
1272 const utils::RefElementSampler &sampler,
1273 const int n_points,
1274 const Eigen::MatrixXd &fun,
1275 const double t,
1276 std::vector<assembler::Assembler::NamedMatrix> &result,
1277 const bool use_sampler,
1278 const bool boundary_only)
1279 {
1280 if (fun.size() <= 0)
1281 {
1282 logger().error("Solve the problem first!");
1283 return;
1284 }
1285
1286 result.clear();
1287
1288 const int actual_dim = mesh.dimension();
1289 assert(!is_problem_scalar);
1290
1291 int index = 0;
1292
1293 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
1294 std::vector<std::pair<std::string, Eigen::MatrixXd>> tmp_t;
1295
1296 for (int i = 0; i < int(bases.size()); ++i)
1297 {
1298 if (boundary_only && mesh.is_volume() && !mesh.is_boundary_element(i))
1299 continue;
1300
1301 const ElementBases &bs = bases[i];
1302 const ElementBases &gbs = gbases[i];
1303 Eigen::MatrixXd local_pts;
1304
1305 if (use_sampler)
1306 {
1307 if (mesh.is_simplex(i))
1308 local_pts = sampler.simplex_points();
1309 else if (mesh.is_cube(i))
1310 local_pts = sampler.cube_points();
1311 else
1312 {
1313 if (mesh.is_volume())
1314 sampler.sample_polyhedron(polys_3d.at(i).first, polys_3d.at(i).second, local_pts, vis_faces_poly, vis_edges_poly);
1315 else
1316 sampler.sample_polygon(polys.at(i), local_pts, vis_faces_poly, vis_edges_poly);
1317 }
1318 }
1319 else
1320 {
1321 if (mesh.is_volume())
1322 {
1323 if (mesh.is_simplex(i))
1324 autogen::p_nodes_3d(disc_orders(i), local_pts);
1325 else if (mesh.is_cube(i))
1326 autogen::q_nodes_3d(disc_orders(i), local_pts);
1327 else
1328 continue;
1329 }
1330 else
1331 {
1332 if (mesh.is_simplex(i))
1333 autogen::p_nodes_2d(disc_orders(i), local_pts);
1334 else if (mesh.is_cube(i))
1335 autogen::q_nodes_2d(disc_orders(i), local_pts);
1336 else
1337 continue;
1338 }
1339 }
1340
1341 assembler.compute_tensor_value(OutputData(t, i, bs, gbs, local_pts, fun), tmp_t);
1342
1343 if (result.empty())
1344 {
1345 result.resize(tmp_t.size());
1346 for (int k = 0; k < tmp_t.size(); ++k)
1347 {
1348 result[k].first = tmp_t[k].first;
1349 result[k].second.resize(n_points, actual_dim * actual_dim);
1350 }
1351 }
1352
1353 for (int k = 0; k < tmp_t.size(); ++k)
1354 {
1355 assert(local_pts.rows() == tmp_t[k].second.rows());
1356 result[k].second.block(index, 0, tmp_t[k].second.rows(), tmp_t[k].second.cols()) = tmp_t[k].second;
1357 }
1358 index += local_pts.rows();
1359 }
1360 }
1361
1362 Eigen::MatrixXd Evaluator::get_bases_position(
1363 const int n_bases,
1364 const std::shared_ptr<mesh::MeshNodes> mesh_nodes)
1365 {
1366 Eigen::MatrixXd func;
1367 func.setZero(n_bases, mesh_nodes->node_position(0).size());
1368
1369 for (int i = 0; i < n_bases; i++)
1370 func.row(i) = mesh_nodes->node_position(i);
1371
1372 return func;
1373 }
1374
1375 Eigen::MatrixXd Evaluator::generate_linear_field(
1376 const int n_bases,
1377 const std::shared_ptr<mesh::MeshNodes> mesh_nodes,
1378 const Eigen::MatrixXd &grad)
1379 {
1380 return utils::flatten(get_bases_position(n_bases, mesh_nodes) * grad.transpose());
1381 }
1382
1383 Eigen::VectorXd Evaluator::integrate_function(
1384 const std::vector<basis::ElementBases> &bases,
1385 const std::vector<basis::ElementBases> &gbases,
1386 const Eigen::MatrixXd &fun,
1387 const int dim,
1388 const int actual_dim)
1389 {
1390 Eigen::VectorXd result;
1391 result.setZero(actual_dim);
1392 for (int e = 0; e < bases.size(); ++e)
1393 {
1394 ElementAssemblyValues vals;
1395 vals.compute(e, dim == 3, bases[e], gbases[e]);
1396
1397 Eigen::MatrixXd u, grad_u;
1398 io::Evaluator::interpolate_at_local_vals(e, dim, actual_dim, vals, fun, u, grad_u);
1399 const quadrature::Quadrature &quadrature = vals.quadrature;
1400 Eigen::VectorXd da = vals.det.array() * quadrature.weights.array();
1401 result += u.transpose() * da;
1402 }
1403
1404 return result;
1405 }
1406} // 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:151
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::get_sidesets
static void get_sidesets(const mesh::Mesh &mesh, Eigen::MatrixXd &pts, Eigen::MatrixXi &faces, Eigen::MatrixXd &sidesets)
returns a triangulated representation of the sideset.
Definition Evaluator.cpp:56
polyfem::io::Evaluator::interpolate_boundary_function_at_vertices
static void interpolate_boundary_function_at_vertices(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, Eigen::MatrixXd &result)
computes integrated solution (fun) per surface face vertex.
Definition Evaluator.cpp:247
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:1375
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 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:1083
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:983
polyfem::io::Evaluator::interpolate_boundary_tensor_function
static void interpolate_boundary_tensor_function(const mesh::Mesh &mesh, const bool is_problem_scalar, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const assembler::Assembler &assembler, const Eigen::MatrixXd &pts, const Eigen::MatrixXi &faces, const Eigen::MatrixXd &fun, const Eigen::MatrixXd &disp, const double t, const bool compute_avg, Eigen::MatrixXd &result, Eigen::MatrixXd &stresses, Eigen::MatrixXd &mises, bool skip_orientation=false)
computes traction forces for fun (tensor * surface normal) result, stress tensor, and von mises,...
Definition Evaluator.cpp:375
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 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:774
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:1362
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:151
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 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:871
polyfem::io::Evaluator::compute_vertex_values
static void compute_vertex_values(const mesh::Mesh &mesh, int actual_dim, const std::vector< basis::ElementBases > &bases, const utils::RefElementSampler &sampler, const Eigen::MatrixXd &fun, Eigen::MatrixXd &result)
evaluates the function fun at the vertices on the mesh
Definition Evaluator.cpp:690
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 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:559
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:1383
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 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:1165
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 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:1263
polyfem::mesh::Mesh3D::get_ordered_vertices_from_hex
std::array< int, 8 > get_ordered_vertices_from_hex(const int element_index) const
Definition Mesh3D.cpp:322
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::Mesh3D::get_ordered_vertices_from_tet
virtual std::array< int, 4 > get_ordered_vertices_from_tet(const int element_index) const
Definition Mesh3D.cpp:346
polyfem::mesh::Mesh
Abstract mesh class to capture 2d/3d conforming and non-conforming meshes.
Definition Mesh.hpp:39
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:158
polyfem::mesh::Mesh::n_vertices
virtual int n_vertices() const =0
number of vertices
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:348
polyfem::mesh::Mesh::point
virtual RowVectorNd point(const int global_index) const =0
point coordinates
polyfem::mesh::Mesh::is_boundary_face
virtual bool is_boundary_face(const int face_global_id) const =0
is face boundary
polyfem::mesh::Mesh::get_boundary_id
virtual int get_boundary_id(const int primitive) const
Get the boundary selection of an element (face in 3d, edge in 2d)
Definition Mesh.hpp:478
polyfem::mesh::Mesh::is_simplex
bool is_simplex(const int el_id) const
checks if element is simples compatible
Definition Mesh.cpp:418
polyfem::mesh::Mesh::is_volume
virtual bool is_volume() const =0
checks if mesh is volume
polyfem::mesh::Mesh::edge_vertex
virtual int edge_vertex(const int e_id, const int lv_id) const =0
id of the edge vertex
polyfem::mesh::Mesh::dimension
int dimension() const
utily for dimension
Definition Mesh.hpp:148
polyfem::mesh::Mesh::n_faces
virtual int n_faces() const =0
number of faces
polyfem::mesh::Mesh::n_edges
virtual int n_edges() const =0
number of edges
polyfem::mesh::Mesh::n_face_vertices
virtual int n_face_vertices(const int f_id) const =0
number of vertices of a face
polyfem::mesh::Mesh::is_boundary_element
virtual bool is_boundary_element(const int element_global_id) const =0
is cell boundary
polyfem::mesh::Mesh::face_vertex
virtual int face_vertex(const int f_id, const int lv_id) const =0
id of the face vertex
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:150
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:185
polyfem::utils::BoundarySampler::normal_for_quad_face
static void normal_for_quad_face(int index, Eigen::MatrixXd &normal)
Definition BoundarySampler.cpp:408
polyfem::utils::BoundarySampler::hex_local_node_coordinates_from_face
static Eigen::MatrixXd hex_local_node_coordinates_from_face(int lf)
Definition BoundarySampler.cpp:87
polyfem::utils::BoundarySampler::normal_for_tri_face
static void normal_for_tri_face(int index, Eigen::MatrixXd &normal)
Definition BoundarySampler.cpp:420
polyfem::utils::BoundarySampler::tet_local_node_coordinates_from_face
static Eigen::MatrixXd tet_local_node_coordinates_from_face(int lf)
Definition BoundarySampler.cpp:72
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:368
polyfem::utils::RefElementSampler::simplex_corners
const Eigen::MatrixXd & simplex_corners() const
Definition RefElementSampler.hpp:42
polyfem::utils::RefElementSampler::simplex_points
const Eigen::MatrixXd & simplex_points() const
Definition RefElementSampler.hpp:43
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:409
polyfem::utils::RefElementSampler::cube_points
const Eigen::MatrixXd & cube_points() const
Definition RefElementSampler.hpp:37
polyfem::utils::RefElementSampler::cube_corners
const Eigen::MatrixXd & cube_corners() const
Definition RefElementSampler.hpp:36
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::p_nodes_2d
void p_nodes_2d(const int p, Eigen::MatrixXd &val)
Definition auto_p_bases.cpp:837
polyfem::autogen::p_nodes_3d
void p_nodes_3d(const int p, Eigen::MatrixXd &val)
Definition auto_p_bases.cpp:2805
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::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:42
polyfem::RowVectorNd
Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor, 1, 3 > RowVectorNd
Definition Types.hpp:11