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