PolyFEM
Loading...
Searching...
No Matches
LagrangeBasis2d.cpp
Go to the documentation of this file.
1
2#include "LagrangeBasis2d.hpp"
3
4#include <polyfem/quadrature/TriQuadrature.hpp>
5#include <polyfem/quadrature/QuadQuadrature.hpp>
6#include <polyfem/autogen/auto_p_bases.hpp>
7#include <polyfem/autogen/auto_q_bases.hpp>
8
9#include <polyfem/assembler/AssemblerUtils.hpp>
10
11#include <polyfem/utils/MaybeParallelFor.hpp>
12
13#include <cassert>
14#include <array>
16
17using namespace polyfem;
18using namespace polyfem::assembler;
19using namespace polyfem::basis;
20using namespace polyfem::mesh;
21using namespace polyfem::quadrature;
22using namespace polyfem::utils;
23
24/*
25
26Axes:
27 y
28 |
29 o──x
30
31Boundaries:
32X axis: left/right
33Y axis: bottom/top
34
35
36
37
38
39Corner nodes:
40 v2
41 │ \
42 │ \
43 │ \
44 x x
45 │ \
46 │ \
47 │ \
48v0──────x─────v1
49
50v0 = (0, 0)
51v1 = (1, 0)
52v2 = (0, 1)
53
54Edge nodes:
55 x
56 │ \
57 │ \
58 │ \
59e2 e1
60 │ \
61 │ \
62 │ \
63 x──────e0───x
64
65e0 = (0.5, 0)
66e1 = (0.5, 0.5)
67e2 = ( 0, 0.5)
68
69
70
71
72
73
74
75
76
77Corner nodes:
78v3──────x─────v2
79 │ ┆ │
80 │ ┆ │
81 │ ┆ │
82 x┄┄┄┄┄┄x┄┄┄┄┄┄x
83 │ ┆ │
84 │ ┆ │
85 │ ┆ │
86v0──────x─────v1
87
88v0 = (0, 0)
89v1 = (1, 0)
90v2 = (1, 1)
91v3 = (0, 1)
92
93Edge nodes:
94 x─────e2──────x
95 │ ┆ │
96 │ ┆ │
97 │ ┆ │
98e3┄┄┄┄┄┄x┄┄┄┄┄e1
99 │ ┆ │
100 │ ┆ │
101 │ ┆ │
102 x─────e0──────x
103
104e0 = (0.5, 0)
105e1 = ( 1, 0.5)
106e2 = (0.5, 1)
107e3 = ( 0, 0.5)
108
109*/
110
111namespace
112{
113
114 void sort(std::array<int, 2> &array)
115 {
116 if (array[0] > array[1])
117 {
118 std::swap(array[0], array[1]);
119 }
120
121 assert(array[0] <= array[1]);
122 }
123
124 template <class InputIterator, class T>
125 int find_index(InputIterator first, InputIterator last, const T &val)
126 {
127 return std::distance(first, std::find(first, last, val));
128 }
129
130 Navigation::Index find_edge(const Mesh2D &mesh, int f, int v1, int v2)
131 {
132 std::array<int, 2> v = {{v1, v2}};
133 std::array<int, 2> u;
134 // std::sort(v.begin(), v.end());
135 sort(v);
136 auto idx = mesh.get_index_from_face(f, 0);
137 for (int lv = 0; lv < mesh.n_face_vertices(idx.face); ++lv)
138 {
139 u[0] = idx.vertex;
140 u[1] = mesh.switch_vertex(idx).vertex;
141 // std::sort(u.begin(), u.end());
142 sort(u);
143 if (u == v)
144 {
145 if (idx.vertex != v1)
146 idx = mesh.switch_vertex(idx);
147 assert(idx.vertex == v1);
148 return idx;
149 }
150 idx = mesh.next_around_face(idx);
151 }
152 throw std::runtime_error("Edge not found");
153 }
154
155 std::array<int, 3> tri_vertices_local_to_global(const Mesh2D &mesh, int f)
156 {
157 assert(mesh.is_simplex(f));
158
159 // Vertex nodes
160 std::array<int, 3> l2g;
161 for (int lv = 0; lv < 3; ++lv)
162 {
163 l2g[lv] = mesh.face_vertex(f, lv);
164 }
165
166 return l2g;
167 }
168
169 std::array<int, 4> quad_vertices_local_to_global(const Mesh2D &mesh, int f)
170 {
171 assert(mesh.is_cube(f));
172
173 // Vertex nodes
174 std::array<int, 4> l2g;
175 for (int lv = 0; lv < 4; ++lv)
176 {
177 l2g[lv] = mesh.face_vertex(f, lv);
178 }
179
180 return l2g;
181 }
182
183 void tri_local_to_global(const bool is_geom_bases, const int p, const Mesh2D &mesh, int f, const Eigen::VectorXi &discr_order, const Eigen::VectorXi &edge_orders, std::vector<int> &res, MeshNodes &nodes, std::vector<std::vector<int>> &edge_virtual_nodes)
184 {
185 int edge_offset = mesh.n_vertices();
186 int face_offset = edge_offset + mesh.n_edges();
187
188 const int n_edge_nodes = p > 1 ? ((p - 1) * 3) : 0;
189 const int nn = p > 2 ? (p - 2) : 0;
190 const int n_face_nodes = nn * (nn + 1) / 2;
191
192 if (p == 0)
193 {
194 res.push_back(nodes.node_id_from_face(f));
195 return;
196 }
197
198 // std::vector<int> res;
199 res.reserve(3 + n_edge_nodes + n_face_nodes);
200
201 // Vertex nodes
202 auto v = tri_vertices_local_to_global(mesh, f);
203
204 // Edge nodes
205 Eigen::Matrix<Navigation::Index, 3, 1> e;
206 Eigen::Matrix<int, 3, 2> ev;
207 ev.row(0) << v[0], v[1];
208 ev.row(1) << v[1], v[2];
209 ev.row(2) << v[2], v[0];
210 for (int le = 0; le < e.rows(); ++le)
211 {
212 const auto index = find_edge(mesh, f, ev(le, 0), ev(le, 1));
213 e[le] = index;
214 }
215
216 for (size_t lv = 0; lv < v.size(); ++lv)
217 {
218 const auto index = e[lv];
219 const auto other_face = mesh.switch_face(index).face;
220
221 if (is_geom_bases)
222 res.push_back(nodes.node_id_from_primitive(v[lv]));
223 else
224 {
225 if (!mesh.is_conforming())
226 {
227 const auto &ncmesh = dynamic_cast<const NCMesh2D &>(mesh);
228 // hanging vertex
229 if (ncmesh.leader_edge_of_vertex(ncmesh.face_vertex(f, lv)) >= 0)
230 res.push_back(-lv - 1);
231 else
232 res.push_back(nodes.node_id_from_primitive(v[lv]));
233 }
234 else
235 {
236 if (discr_order.size() > 0 && other_face >= 0 && discr_order(f) > discr_order(other_face))
237 res.push_back(-lv - 1);
238 else
239 res.push_back(nodes.node_id_from_primitive(v[lv]));
240 }
241 }
242 }
243
244 for (int le = 0; le < e.rows(); ++le)
245 {
246 const auto index = e[le];
247 const auto other_face = mesh.switch_face(index).face;
248
249 if (is_geom_bases)
250 {
251 auto node_ids = nodes.node_ids_from_edge(index, p - 1);
252 res.insert(res.end(), node_ids.begin(), node_ids.end());
253 }
254 else
255 {
256 if (!mesh.is_conforming())
257 {
258 const auto &ncmesh = dynamic_cast<const NCMesh2D &>(mesh);
259 const int edge_id = ncmesh.face_edge(f, le);
260 // follower edge
261 if (ncmesh.leader_edge_of_edge(edge_id) >= 0)
262 {
263 for (int tmp = 0; tmp < p - 1; ++tmp)
264 res.push_back(-le - 1);
265 }
266 // leader or conforming edge with constrained order
267 else if (edge_orders[edge_id] < discr_order[f])
268 {
269 for (int tmp = 0; tmp < p - 1; ++tmp)
270 res.push_back(-le - 1);
271 // leader edge
272 if (ncmesh.n_follower_edges(edge_id) > 0)
273 edge_virtual_nodes[edge_id] = nodes.node_ids_from_edge(index, edge_orders[edge_id] - 1);
274 }
275 // leader or conforming edge with unconstrained order
276 else
277 {
278 auto node_ids = nodes.node_ids_from_edge(index, p - 1);
279 res.insert(res.end(), node_ids.begin(), node_ids.end());
280 }
281 }
282 else
283 {
284 bool skip_other = discr_order.size() > 0 && other_face >= 0 && discr_order(f) > discr_order(other_face);
285
286 if (discr_order.size() > 0 && other_face >= 0 && discr_order(f) > discr_order(other_face))
287 {
288 for (int tmp = 0; tmp < p - 1; ++tmp)
289 res.push_back(-le - 1);
290 }
291 else
292 {
293 auto node_ids = nodes.node_ids_from_edge(index, p - 1);
294 res.insert(res.end(), node_ids.begin(), node_ids.end());
295 }
296 }
297 }
298 }
299
300 if (n_face_nodes > 0)
301 {
302 const auto index = e[0];
303
304 auto node_ids = nodes.node_ids_from_face(index, p - 2);
305 res.insert(res.end(), node_ids.begin(), node_ids.end());
306 }
307
308 assert(res.size() == size_t(3 + n_edge_nodes + n_face_nodes));
309 // return res;
310 }
311
312 void quad_local_to_global(const bool serendipity, const int q, const Mesh2D &mesh, int f, const Eigen::VectorXi &discr_order, std::vector<int> &res, MeshNodes &nodes)
313 {
314 int edge_offset = mesh.n_vertices();
315 int face_offset = edge_offset + mesh.n_edges();
316
317 const int nn = q > 1 ? (q - 1) : 0;
318 const int n_edge_nodes = nn * 4;
319 const int n_face_nodes = serendipity ? 0 : nn * nn;
320
321 if (q == 0)
322 {
323 res.push_back(nodes.node_id_from_face(f));
324 return;
325 }
326
327 // std::vector<int> res;
328 res.reserve(4 + n_edge_nodes + n_face_nodes);
329
330 // Vertex nodes
331 auto v = quad_vertices_local_to_global(mesh, f);
332
333 // Edge nodes
334 Eigen::Matrix<Navigation::Index, 4, 1> e;
335 Eigen::Matrix<int, 4, 2> ev;
336 ev.row(0) << v[0], v[1];
337 ev.row(1) << v[1], v[2];
338 ev.row(2) << v[2], v[3];
339 ev.row(3) << v[3], v[0];
340 for (int le = 0; le < e.rows(); ++le)
341 {
342 const auto index = find_edge(mesh, f, ev(le, 0), ev(le, 1));
343 e[le] = index;
344 }
345
346 for (size_t lv = 0; lv < v.size(); ++lv)
347 {
348 const auto index = e[lv];
349 const auto other_face = mesh.switch_face(index).face;
350
351 if (discr_order.size() > 0 && other_face >= 0 && discr_order(f) > discr_order(other_face))
352 res.push_back(-lv - 1);
353 else
354 res.push_back(nodes.node_id_from_primitive(v[lv]));
355 }
356
357 for (int le = 0; le < e.rows(); ++le)
358 {
359 const auto index = e[le];
360 const auto other_face = mesh.switch_face(index).face;
361
362 bool skip_other = discr_order.size() > 0 && other_face >= 0 && discr_order(f) > discr_order(other_face);
363
364 if (discr_order.size() > 0 && other_face >= 0 && discr_order(f) > discr_order(other_face))
365 {
366 for (int tmp = 0; tmp < q - 1; ++tmp)
367 res.push_back(-le - 1);
368 }
369 else
370 {
371 auto node_ids = nodes.node_ids_from_edge(index, q - 1);
372 res.insert(res.end(), node_ids.begin(), node_ids.end());
373 }
374 }
375
376 if (n_face_nodes > 0)
377 {
378 const auto index = e[0];
379
380 auto node_ids = nodes.node_ids_from_face(index, q - 1);
381 res.insert(res.end(), node_ids.begin(), node_ids.end());
382 }
383
384 assert(res.size() == size_t(4 + n_edge_nodes + n_face_nodes));
385 // return res;
386 }
387
388 // -----------------------------------------------------------------------------
389
390 // -----------------------------------------------------------------------------
391
405 void compute_nodes(
406 const Mesh2D &mesh,
407 const Eigen::VectorXi &discr_orders,
408 const Eigen::VectorXi &edge_orders,
409 const bool serendipity,
410 const bool has_polys,
411 const bool is_geom_bases,
412 MeshNodes &nodes,
413 std::vector<std::vector<int>> &edge_virtual_nodes,
414 std::vector<std::vector<int>> &element_nodes_id,
415 std::vector<LocalBoundary> &local_boundary,
416 std::map<int, InterfaceData> &poly_edge_to_data)
417 {
418 // Step 1: Assign global node ids for each quads
419 local_boundary.clear();
420
421 element_nodes_id.resize(mesh.n_faces());
422
423 if (!mesh.is_conforming())
424 {
425 const auto &ncmesh = dynamic_cast<const NCMesh2D &>(mesh);
426 edge_virtual_nodes.resize(ncmesh.n_edges());
427 }
428
429 for (int f = 0; f < mesh.n_faces(); ++f)
430 {
431 const int discr_order = discr_orders(f);
432 if (mesh.is_cube(f))
433 {
434 quad_local_to_global(serendipity, discr_order, mesh, f, discr_orders, element_nodes_id[f], nodes);
435
436 LocalBoundary lb(f, BoundaryType::QUAD_LINE);
437
438 auto v = quad_vertices_local_to_global(mesh, f);
439 Eigen::Matrix<int, 4, 2> ev;
440 ev.row(0) << v[0], v[1];
441 ev.row(1) << v[1], v[2];
442 ev.row(2) << v[2], v[3];
443 ev.row(3) << v[3], v[0];
444
445 for (int i = 0; i < int(ev.rows()); ++i)
446 {
447 const auto index = find_edge(mesh, f, ev(i, 0), ev(i, 1));
448 const int edge = index.edge;
449
450 if (mesh.is_boundary_edge(edge) || mesh.get_boundary_id(edge) > 0)
451 {
452 lb.add_boundary_primitive(edge, i);
453 }
454 }
455
456 if (!lb.empty())
457 local_boundary.emplace_back(lb);
458 }
459 else if (mesh.is_simplex(f))
460 {
461 tri_local_to_global(is_geom_bases, discr_order, mesh, f, discr_orders, edge_orders, element_nodes_id[f], nodes, edge_virtual_nodes);
462
463 auto v = tri_vertices_local_to_global(mesh, f);
464
465 Eigen::Matrix<int, 3, 2> ev;
466 ev.row(0) << v[0], v[1];
467 ev.row(1) << v[1], v[2];
468 ev.row(2) << v[2], v[0];
469
470 LocalBoundary lb(f, BoundaryType::TRI_LINE);
471
472 for (int i = 0; i < int(ev.rows()); ++i)
473 {
474 const auto index = find_edge(mesh, f, ev(i, 0), ev(i, 1));
475 const int edge = index.edge;
476
477 if (mesh.is_boundary_edge(edge) || mesh.get_boundary_id(edge) > 0)
478 {
479 lb.add_boundary_primitive(edge, i);
480 }
481 }
482
483 if (!lb.empty())
484 local_boundary.emplace_back(lb);
485 }
486 }
487
488 if (!has_polys)
489 return;
490
491 // Step 2: Iterate over edges of polygons and compute interface weights
492 Eigen::VectorXi indices;
493 for (int f = 0; f < mesh.n_faces(); ++f)
494 {
495 // Skip non-polytopes
496 if (!mesh.is_polytope(f))
497 {
498 continue;
499 }
500
501 auto index = mesh.get_index_from_face(f, 0);
502 for (int lv = 0; lv < mesh.n_face_vertices(f); ++lv)
503 {
504 auto index2 = mesh.switch_face(index);
505 if (index2.face >= 0)
506 {
507 // Opposite face is a quad, we need to set interface data
508 int f2 = index2.face;
509 const int discr_order = discr_orders(f2);
510
511 if (mesh.is_cube(f2))
512 {
513 indices = LagrangeBasis2d::quad_edge_local_nodes(discr_order, mesh, index2);
514 }
515 else if (mesh.is_simplex(f2))
516 {
517 indices = LagrangeBasis2d::tri_edge_local_nodes(discr_order, mesh, index2);
518 }
519 else
520 {
521 // assert(false);
522 continue;
523 }
524 InterfaceData data;
525 assert(indices.size() == 2 + discr_order - 1);
526 data.local_indices.insert(data.local_indices.begin(), indices.data(), indices.data() + indices.size());
527 assert(indices.size() == data.local_indices.size());
528 poly_edge_to_data[index2.edge] = data;
529 }
530 index = mesh.next_around_face(index);
531 }
532 }
533 }
534
541 void local_to_global(const Eigen::MatrixXd &verts, const Eigen::MatrixXd &uv, Eigen::MatrixXd &pts)
542 {
543 const int dim = verts.cols();
544 const int N = uv.rows();
545 assert(dim == 2);
546 assert(uv.cols() == dim);
547 assert(verts.rows() == dim + 1);
548
549 pts.setZero(N, dim);
550 for (int i = 0; i < N; i++)
551 pts.row(i) = uv(i, 0) * verts.row(1) + uv(i, 1) * verts.row(2) + (1.0 - uv(i, 0) - uv(i, 1)) * verts.row(0);
552 }
553
560 void global_to_local(const Eigen::MatrixXd &verts, const Eigen::MatrixXd &pts, Eigen::MatrixXd &uv)
561 {
562 const int dim = verts.cols();
563 const int N = pts.rows();
564 assert(dim == 2);
565 assert(verts.rows() == dim + 1);
566 assert(pts.cols() == dim);
567
568 Eigen::Matrix2d J;
569 for (int i = 0; i < dim; i++)
570 J.col(i) = verts.row(i + 1) - verts.row(0);
571
572 double detJ = J(0, 0) * J(1, 1) - J(0, 1) * J(1, 0);
573 J /= detJ;
574
575 uv.setZero(N, dim);
576 maybe_parallel_for(N, [&](int start, int end, int thread_id) {
577 for (int i = start; i < end; i++)
578 {
579 auto point = pts.row(i) - verts.row(0);
580 uv(i, 0) = J(1, 1) * point(0) - J(0, 1) * point(1);
581 uv(i, 1) = J(0, 0) * point(1) - J(1, 0) * point(0);
582 }
583 });
584 }
585
592 void compute_edge_orders(const NCMesh2D &mesh, const Eigen::VectorXi &elem_orders, Eigen::VectorXi &edge_orders)
593 {
594 const int max_order = elem_orders.maxCoeff();
595 edge_orders.setConstant(mesh.n_edges(), max_order);
596
597 for (int i = 0; i < mesh.n_faces(); i++)
598 for (int j = 0; j < mesh.n_face_vertices(i); j++)
599 edge_orders[mesh.face_edge(i, j)] = std::min(edge_orders[mesh.face_edge(i, j)], elem_orders[i]);
600
601 for (int i = 0; i < mesh.n_edges(); i++)
602 if (mesh.leader_edge_of_edge(i) >= 0)
603 edge_orders[mesh.leader_edge_of_edge(i)] = std::min(edge_orders[mesh.leader_edge_of_edge(i)], edge_orders[i]);
604
605 for (int i = 0; i < mesh.n_edges(); i++)
606 if (mesh.leader_edge_of_edge(i) >= 0)
607 edge_orders[i] = std::min(edge_orders[mesh.leader_edge_of_edge(i)], edge_orders[i]);
608 }
609} // anonymous namespace
610
611Eigen::VectorXi LagrangeBasis2d::tri_edge_local_nodes(const int p, const Mesh2D &mesh, Navigation::Index index)
612{
613 int f = index.face;
614 assert(mesh.is_simplex(f));
615
616 // Local to global mapping of node indices
617 auto l2g = tri_vertices_local_to_global(mesh, f);
618
619 // Extract requested interface
620 Eigen::VectorXi result(2 + (p - 1));
621 result[0] = find_index(l2g.begin(), l2g.end(), index.vertex);
622 result[result.size() - 1] = find_index(l2g.begin(), l2g.end(), mesh.switch_vertex(index).vertex);
623
624 if ((result[0] == 0 && result[result.size() - 1] == 1) || (result[0] == 1 && result[result.size() - 1] == 2) || (result[0] == 2 && result[result.size() - 1] == 0))
625 {
626 for (int i = 0; i < p - 1; ++i)
627 {
628 result[i + 1] = 3 + result[0] * (p - 1) + i;
629 }
630 }
631 else
632 {
633 for (int i = 0; i < p - 1; ++i)
634 {
635 result[i + 1] = 3 + (result[0] + (result[0] == 0 ? 3 : 0)) * (p - 1) - i - 1;
636 }
637 }
638
639 return result;
640}
641
642Eigen::VectorXi LagrangeBasis2d::quad_edge_local_nodes(const int q, const Mesh2D &mesh, Navigation::Index index)
643{
644 int f = index.face;
645 assert(mesh.is_cube(f));
646
647 // Local to global mapping of node indices
648 auto l2g = quad_vertices_local_to_global(mesh, f);
649
650 // Extract requested interface
651 Eigen::VectorXi result(2 + (q - 1));
652 result[0] = find_index(l2g.begin(), l2g.end(), index.vertex);
653 result[result.size() - 1] = find_index(l2g.begin(), l2g.end(), mesh.switch_vertex(index).vertex);
654
655 if ((result[0] == 0 && result[result.size() - 1] == 1) || (result[0] == 1 && result[result.size() - 1] == 2) || (result[0] == 2 && result[result.size() - 1] == 3) || (result[0] == 3 && result[result.size() - 1] == 0))
656 {
657 for (int i = 0; i < q - 1; ++i)
658 {
659 result[i + 1] = 4 + result[0] * (q - 1) + i;
660 }
661 }
662 else
663 {
664 for (int i = 0; i < q - 1; ++i)
665 {
666 result[i + 1] = 4 + (result[0] + (result[0] == 0 ? 4 : 0)) * (q - 1) - i - 1;
667 }
668 }
669
670 return result;
671}
672
673// -----------------------------------------------------------------------------
674
675int LagrangeBasis2d::build_bases(
676 const Mesh2D &mesh,
677 const std::string &assembler,
678 const int quadrature_order,
679 const int mass_quadrature_order,
680 const int discr_order,
681 const bool serendipity,
682 const bool has_polys,
683 const bool is_geom_bases,
684 const bool use_corner_quadrature,
685 std::vector<ElementBases> &bases,
686 std::vector<LocalBoundary> &local_boundary,
687 std::map<int, InterfaceData> &poly_edge_to_data,
688 std::shared_ptr<MeshNodes> &mesh_nodes)
689{
690
691 Eigen::VectorXi discr_orders(mesh.n_faces());
692 discr_orders.setConstant(discr_order);
693
694 return build_bases(mesh, assembler, quadrature_order, mass_quadrature_order, discr_orders, serendipity, has_polys, is_geom_bases, use_corner_quadrature, bases, local_boundary, poly_edge_to_data, mesh_nodes);
695}
696
697int LagrangeBasis2d::build_bases(
698 const Mesh2D &mesh,
699 const std::string &assembler,
700 const int quadrature_order,
701 const int mass_quadrature_order,
702 const Eigen::VectorXi &discr_orders,
703 const bool serendipity,
704 const bool has_polys,
705 const bool is_geom_bases,
706 const bool use_corner_quadrature,
707 std::vector<ElementBases> &bases,
708 std::vector<LocalBoundary> &local_boundary,
709 std::map<int, InterfaceData> &poly_edge_to_data,
710 std::shared_ptr<MeshNodes> &mesh_nodes)
711{
712 assert(!mesh.is_volume());
713 assert(discr_orders.size() == mesh.n_faces());
714
715 const int max_p = discr_orders.maxCoeff();
716 const int nn = max_p > 1 ? (max_p - 1) : 0;
717 const int n_face_nodes = std::max(nn * nn, max_p == 1 ? 1 : 0);
718
719 Eigen::VectorXi edge_orders;
720 if (!mesh.is_conforming())
721 {
722 const auto &ncmesh = dynamic_cast<const NCMesh2D &>(mesh);
723 compute_edge_orders(ncmesh, discr_orders, edge_orders);
724 }
725
726 mesh_nodes = std::make_shared<MeshNodes>(mesh, has_polys, !is_geom_bases, (max_p > 1 ? (max_p - 1) : 0) * (is_geom_bases ? 2 : 1), max_p == 0 ? 1 : n_face_nodes);
727 MeshNodes &nodes = *mesh_nodes;
728 std::vector<std::vector<int>> element_nodes_id, edge_virtual_nodes;
729 compute_nodes(mesh, discr_orders, edge_orders, serendipity, has_polys, is_geom_bases, nodes, edge_virtual_nodes, element_nodes_id, local_boundary, poly_edge_to_data);
730 // boundary_nodes = nodes.boundary_nodes();
731
732 bases.resize(mesh.n_faces());
733 std::vector<int> interface_elements;
734 interface_elements.reserve(mesh.n_faces());
735
736 for (int e = 0; e < mesh.n_faces(); ++e)
737 {
738 ElementBases &b = bases[e];
739 const int discr_order = discr_orders(e);
740 const int n_el_bases = element_nodes_id[e].size();
741 b.bases.resize(n_el_bases);
742
743 bool skip_interface_element = false;
744
745 for (int j = 0; j < n_el_bases; ++j)
746 {
747 // mark interface between elements of different order
748 const int global_index = element_nodes_id[e][j];
749 if (global_index < 0)
750 {
751 skip_interface_element = true;
752 break;
753 }
754 }
755
756 if (skip_interface_element)
757 {
758 interface_elements.push_back(e);
759 }
760
761 if (mesh.is_cube(e))
762 {
763 const int real_order = quadrature_order > 0 ? quadrature_order : AssemblerUtils::quadrature_order(assembler, discr_order, AssemblerUtils::BasisType::CUBE_LAGRANGE, 2);
764 const int real_mass_order = mass_quadrature_order > 0 ? mass_quadrature_order : AssemblerUtils::quadrature_order("Mass", discr_order, AssemblerUtils::BasisType::CUBE_LAGRANGE, 2);
765 b.set_quadrature([real_order](Quadrature &quad) {
766 QuadQuadrature quad_quadrature;
767 quad_quadrature.get_quadrature(real_order, quad);
768 });
769 b.set_mass_quadrature([real_mass_order](Quadrature &quad) {
770 QuadQuadrature quad_quadrature;
771 quad_quadrature.get_quadrature(real_mass_order, quad);
772 });
773 // quad_quadrature.get_quadrature(real_order, b.quadrature);
774
775 b.set_local_node_from_primitive_func([discr_order, e](const int primitive_id, const Mesh &mesh) {
776 const auto &mesh2d = dynamic_cast<const Mesh2D &>(mesh);
777 auto index = mesh2d.get_index_from_face(e);
778
779 for (int le = 0; le < mesh2d.n_face_vertices(e); ++le)
780 {
781 if (index.edge == primitive_id)
782 break;
783 index = mesh2d.next_around_face(index);
784 }
785 assert(index.edge == primitive_id);
786 return quad_edge_local_nodes(discr_order, mesh2d, index);
787 });
788
789 for (int j = 0; j < n_el_bases; ++j)
790 {
791 const int global_index = element_nodes_id[e][j];
792
793 // if(!skip_interface_element)
794 b.bases[j].init(discr_order, global_index, j, nodes.node_position(global_index));
795
796 const int dtmp = serendipity ? -2 : discr_order;
797
798 b.bases[j].set_basis([dtmp, j](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { autogen::q_basis_value_2d(dtmp, j, uv, val); });
799 b.bases[j].set_grad([dtmp, j](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { autogen::q_grad_basis_value_2d(dtmp, j, uv, val); });
800 }
801 }
802 else if (mesh.is_simplex(e))
803 {
804 const int real_order = quadrature_order > 0 ? quadrature_order : AssemblerUtils::quadrature_order(assembler, discr_order, AssemblerUtils::BasisType::SIMPLEX_LAGRANGE, 2);
805 const int real_mass_order = mass_quadrature_order > 0 ? mass_quadrature_order : AssemblerUtils::quadrature_order("Mass", discr_order, AssemblerUtils::BasisType::SIMPLEX_LAGRANGE, 2);
806 b.set_quadrature([real_order, use_corner_quadrature](Quadrature &quad) {
807 TriQuadrature tri_quadrature(use_corner_quadrature);
808 tri_quadrature.get_quadrature(real_order, quad);
809 });
810 b.set_mass_quadrature([real_mass_order, use_corner_quadrature](Quadrature &quad) {
811 TriQuadrature tri_quadrature(use_corner_quadrature);
812 tri_quadrature.get_quadrature(real_mass_order, quad);
813 });
814
815 b.set_local_node_from_primitive_func([discr_order, e](const int primitive_id, const Mesh &mesh) {
816 const auto &mesh2d = dynamic_cast<const Mesh2D &>(mesh);
817 auto index = mesh2d.get_index_from_face(e);
818
819 for (int le = 0; le < mesh2d.n_face_vertices(e); ++le)
820 {
821 if (index.edge == primitive_id)
822 break;
823 index = mesh2d.next_around_face(index);
824 }
825 assert(index.edge == primitive_id);
826 return tri_edge_local_nodes(discr_order, mesh2d, index);
827 });
828
829 const bool rational = is_geom_bases && mesh.is_rational() && !mesh.cell_weights(e).empty();
830
831 for (int j = 0; j < n_el_bases; ++j)
832 {
833 const int global_index = element_nodes_id[e][j];
834
835 if (!skip_interface_element)
836 {
837 b.bases[j].init(discr_order, global_index, j, nodes.node_position(global_index));
838 }
839
840 if (rational)
841 {
842 const auto &w = mesh.cell_weights(e);
843 assert(discr_order == 2);
844 assert(w.size() == 6);
845
846 b.bases[j].set_basis([discr_order, j, w](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) {
847 autogen::p_basis_value_2d(discr_order, j, uv, val);
848 Eigen::MatrixXd denom = val;
849 denom.setZero();
850 Eigen::MatrixXd tmp;
851
852 for (int k = 0; k < 6; ++k)
853 {
854 autogen::p_basis_value_2d(discr_order, k, uv, tmp);
855 denom += w[k] * tmp;
856 }
857
858 val = (w[j] * val.array() / denom.array()).eval();
859 });
860
861 b.bases[j].set_grad([discr_order, j, w](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) {
862 Eigen::MatrixXd b;
863 autogen::p_basis_value_2d(discr_order, j, uv, b);
864 autogen::p_grad_basis_value_2d(discr_order, j, uv, val);
865 Eigen::MatrixXd denom = b;
866 denom.setZero();
867 Eigen::MatrixXd denom_prime = val;
868 denom_prime.setZero();
869 Eigen::MatrixXd tmp;
870
871 for (int k = 0; k < 6; ++k)
872 {
873 autogen::p_basis_value_2d(discr_order, k, uv, tmp);
874 denom += w[k] * tmp;
875
876 autogen::p_grad_basis_value_2d(discr_order, k, uv, tmp);
877 denom_prime += w[k] * tmp;
878 }
879
880 val.col(0) = ((w[j] * val.col(0).array() * denom.array() - w[j] * b.array() * denom_prime.col(0).array()) / (denom.array() * denom.array())).eval();
881 val.col(1) = ((w[j] * val.col(1).array() * denom.array() - w[j] * b.array() * denom_prime.col(1).array()) / (denom.array() * denom.array())).eval();
882 });
883 }
884 else
885 {
886 // pick out basis functions using autogenerated code
887 b.bases[j].set_basis([discr_order, j](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { autogen::p_basis_value_2d(discr_order, j, uv, val); });
888 b.bases[j].set_grad([discr_order, j](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { autogen::p_grad_basis_value_2d(discr_order, j, uv, val); });
889 }
890 }
891 }
892 else
893 {
894 // Polygon bases are built later on
895 }
896
897#ifndef NDEBUG
898 if (mesh.is_conforming())
899 {
900 Eigen::MatrixXd uv(4, 2);
901 uv << 0.1, 0.1, 0.3, 0.3, 0.9, 0.01, 0.01, 0.9;
902 Eigen::MatrixXd dx(4, 1);
903 dx.setConstant(1e-6);
904 Eigen::MatrixXd uvdx = uv;
905 uvdx.col(0) += dx;
906 Eigen::MatrixXd uvdy = uv;
907 uvdy.col(1) += dx;
908 Eigen::MatrixXd grad, val, vdx, vdy;
909
910 for (int j = 0; j < n_el_bases; ++j)
911 {
912 b.bases[j].eval_grad(uv, grad);
913
914 b.bases[j].eval_basis(uv, val);
915 b.bases[j].eval_basis(uvdx, vdx);
916 b.bases[j].eval_basis(uvdy, vdy);
917
918 assert((grad.col(0) - (vdx - val) / 1e-6).norm() < 1e-4);
919 assert((grad.col(1) - (vdy - val) / 1e-6).norm() < 1e-4);
920 }
921 }
922#endif
923 }
924
925 if (!is_geom_bases)
926 {
927 if (!mesh.is_conforming())
928 {
929 const auto &ncmesh = dynamic_cast<const NCMesh2D &>(mesh);
930
931 std::vector<std::vector<int>> elementOrder;
932 {
933 const int max_order = discr_orders.maxCoeff(), min_order = discr_orders.minCoeff();
934 int max_level = 0;
935 for (int e = 0; e < ncmesh.n_faces(); e++)
936 if (max_level < ncmesh.face_ref_level(e))
937 max_level = ncmesh.face_ref_level(e);
938
939 elementOrder.resize((max_level + 1) * (max_order - min_order + 1));
940 int N = 0;
941 int cur_level = 0;
942 while (cur_level <= max_level)
943 {
944 int order = min_order;
945 while (order <= max_order)
946 {
947 int cur_bucket = (max_order - min_order + 1) * cur_level + (order - min_order);
948 for (int i = 0; i < ncmesh.n_faces(); i++)
949 {
950 if (ncmesh.face_ref_level(i) != cur_level || discr_orders[i] != order)
951 continue;
952
953 N++;
954 elementOrder[cur_bucket].push_back(i);
955 }
956 order++;
957 }
958 cur_level++;
959 }
960 }
961
962 for (const auto &bucket : elementOrder)
963 {
964 if (bucket.size() == 0)
965 continue;
966 for (const int e : bucket)
967 {
968 ElementBases &b = bases[e];
969 const int discr_order = discr_orders(e);
970 const int n_el_bases = element_nodes_id[e].size();
971
972 for (int j = 0; j < n_el_bases; ++j)
973 {
974 const int global_index = element_nodes_id[e][j];
975
976 if (global_index >= 0)
977 b.bases[j].init(discr_order, global_index, j, nodes.node_position(global_index));
978 else
979 {
980 int large_edge = -1, local_edge = -1;
981 int opposite_element = -1;
982 bool is_on_leader_edge = false;
983
984 if (j < 3 + 3 * (discr_order - 1) && j >= 3)
985 {
986 local_edge = (j - 3) / (discr_order - 1);
987 large_edge = ncmesh.face_edge(e, local_edge);
988 if (ncmesh.n_follower_edges(large_edge) > 0)
989 is_on_leader_edge = true;
990 }
991
992 // this node is on a leader edge, but its order is constrained
993 if (is_on_leader_edge)
994 {
995 const int edge_order = edge_orders[large_edge];
996
997 // indices of fake nodes on this edge in the opposite element
998 Eigen::VectorXi indices;
999 indices.resize(discr_order + 1);
1000 for (int i = 0; i < discr_order - 1; i++)
1001 {
1002 indices[i + 1] = 3 + local_edge * (discr_order - 1) + i;
1003 }
1004 int index = indices[(j - 3) % (discr_order - 1) + 1];
1005
1006 // the position of node j in the opposite element
1007 Eigen::MatrixXd lnodes;
1008 autogen::p_nodes_2d(discr_order, lnodes);
1009 Eigen::Vector2d node_position = lnodes.row(index);
1010
1011 std::function<double(const int, const int, const double)> basis_1d = [](const int order, const int id, const double x) -> double {
1012 assert(id <= order && id >= 0);
1013 double y = 1;
1014 for (int o = 0; o <= order; o++)
1015 {
1016 if (o != id)
1017 y *= (x * order - o) / (id - o);
1018 }
1019 return y;
1020 };
1021
1022 // apply basis projection
1023 double x = NCMesh2D::element_weight_to_edge_weight(local_edge, node_position);
1024 for (int i = 0; i < edge_virtual_nodes[large_edge].size(); i++)
1025 {
1026 const int global_index = edge_virtual_nodes[large_edge][i];
1027 const double weight = basis_1d(edge_order, i + 1, x);
1028 if (std::abs(weight) < 1e-12)
1029 continue;
1030 b.bases[j].global().emplace_back(global_index, nodes.node_position(global_index), weight);
1031 }
1032
1033 const auto &global_1 = b.bases[local_edge].global();
1034 const auto &global_2 = b.bases[(local_edge + 1) % 3].global();
1035 double weight = basis_1d(edge_order, 0, x);
1036 if (std::abs(weight) > 1e-12)
1037 for (size_t ii = 0; ii < global_1.size(); ++ii)
1038 b.bases[j].global().emplace_back(global_1[ii].index, global_1[ii].node, weight * global_1[ii].val);
1039 weight = basis_1d(edge_order, edge_order, x);
1040 if (std::abs(weight) > 1e-12)
1041 for (size_t ii = 0; ii < global_2.size(); ++ii)
1042 b.bases[j].global().emplace_back(global_2[ii].index, global_2[ii].node, weight * global_2[ii].val);
1043 }
1044 else
1045 {
1046 Eigen::MatrixXd global_position;
1047 // this node is a hanging vertex, and it's not on a follower edge
1048 if (j < 3)
1049 {
1050 const int v_id = ncmesh.face_vertex(e, j);
1051 large_edge = ncmesh.leader_edge_of_vertex(v_id);
1052 opposite_element = ncmesh.face_neighbor(large_edge);
1053 global_position = ncmesh.point(v_id);
1054 }
1055 // this node is on a follower edge
1056 else
1057 {
1058 // find the local id of small edge
1059 int small_edge = ncmesh.face_edge(e, local_edge);
1060 large_edge = ncmesh.leader_edge_of_edge(small_edge);
1061
1062 // this edge is non-conforming
1063 if (ncmesh.n_face_neighbors(small_edge) == 1)
1064 opposite_element = ncmesh.face_neighbor(large_edge);
1065 // this edge is conforming, but the order on two sides is different
1066 else
1067 {
1068 Navigation::Index idx;
1069 idx.edge = small_edge;
1070 idx.vertex = ncmesh.edge_vertex(small_edge, 0);
1071 idx.face_corner = -1;
1072 idx.face = e;
1073 opposite_element = ncmesh.switch_face(idx).face;
1074 }
1075
1076 // the position of node j in the opposite element
1077 Eigen::MatrixXd lnodes;
1078 autogen::p_nodes_2d(discr_order, lnodes);
1079 global_position = lnodes.row(j);
1080
1081 Eigen::MatrixXd verts(ncmesh.n_face_vertices(e), 2);
1082 for (int i = 0; i < verts.rows(); i++)
1083 verts.row(i) = ncmesh.point(ncmesh.face_vertex(e, i));
1084
1085 Eigen::MatrixXd local_position = lnodes.row(j);
1086 local_to_global(verts, local_position, global_position);
1087 }
1088
1089 Eigen::MatrixXd verts(ncmesh.n_face_vertices(e), 2);
1090 for (int i = 0; i < verts.rows(); i++)
1091 verts.row(i) = ncmesh.point(ncmesh.face_vertex(opposite_element, i));
1092
1093 Eigen::MatrixXd node_position;
1094 global_to_local(verts, global_position, node_position);
1095
1096 // evaluate the basis of the opposite element at this node
1097 const auto &other_bases = bases[opposite_element];
1098 std::vector<AssemblyValues> w;
1099 other_bases.evaluate_bases(node_position, w);
1100
1101 // apply basis projection
1102 for (long i = 0; i < w.size(); ++i)
1103 {
1104 assert(w[i].val.size() == 1);
1105 if (std::abs(w[i].val(0)) < 1e-12)
1106 continue;
1107
1108 for (size_t ii = 0; ii < other_bases.bases[i].global().size(); ++ii)
1109 {
1110 const auto &other_global = other_bases.bases[i].global()[ii];
1111 b.bases[j].global().emplace_back(other_global.index, other_global.node, w[i].val(0) * other_global.val);
1112 }
1113 }
1114 }
1115
1116 auto &global_ = b.bases[j].global();
1117 if (global_.size() <= 1)
1118 continue;
1119
1120 std::map<int, Local2Global> list;
1121 for (size_t ii = 0; ii < global_.size(); ii++)
1122 {
1123 auto pair = list.insert({global_[ii].index, global_[ii]});
1124 if (!pair.second && pair.first != list.end())
1125 {
1126 assert((pair.first->second.node - global_[ii].node).norm() < 1e-12);
1127 pair.first->second.val += global_[ii].val;
1128 }
1129 }
1130
1131 global_.clear();
1132 for (auto it = list.begin(); it != list.end(); ++it)
1133 {
1134 if (std::abs(it->second.val) > 1e-12)
1135 global_.push_back(it->second);
1136 }
1137 }
1138 }
1139 }
1140 }
1141 }
1142 else
1143 {
1144 // loop over all potential element orders (since multiple loops may be needed to constrain interfaces between more than two elements)
1145 for (int pp = 2; pp <= autogen::MAX_P_BASES; ++pp)
1146 {
1147 // loop again over interface elements to address mismatched polynomial orders by constraining the higher order nodes
1148 for (int e : interface_elements)
1149 {
1150 ElementBases &b = bases[e];
1151 const int discr_order = discr_orders(e);
1152 const int n_el_bases = element_nodes_id[e].size();
1153
1154 assert(discr_order > 1);
1155 if (discr_order != pp)
1156 continue;
1157
1158 if (mesh.is_cube(e))
1159 {
1160 // TODO
1161 assert(false);
1162 }
1163 else if (mesh.is_simplex(e))
1164 {
1165 for (int j = 0; j < n_el_bases; ++j)
1166 {
1167 const int global_index = element_nodes_id[e][j];
1168
1169 if (global_index >= 0)
1170 b.bases[j].init(discr_order, global_index, j, nodes.node_position(global_index));
1171 else
1172 {
1173 const auto le = -(global_index + 1);
1174
1175 auto v = tri_vertices_local_to_global(mesh, e);
1176 Eigen::Matrix<int, 3, 2> ev;
1177 ev.row(0) << v[0], v[1];
1178 ev.row(1) << v[1], v[2];
1179 ev.row(2) << v[2], v[0];
1180
1181 const auto index = mesh.switch_face(find_edge(mesh, e, ev(le, 0), ev(le, 1)));
1182 const auto other_face = index.face;
1183
1184 Eigen::RowVector2d node_position;
1185 assert(discr_order > 1);
1186
1187 auto indices = tri_edge_local_nodes(discr_order, mesh, index);
1188 Eigen::MatrixXd lnodes;
1189 autogen::p_nodes_2d(discr_order, lnodes);
1190
1191 if (j < 3)
1192 node_position = lnodes.row(indices(0));
1193 else if (j < 3 + 3 * (discr_order - 1))
1194 node_position = lnodes.row(indices(((j - 3) % (discr_order - 1)) + 1));
1195 else
1196 assert(false);
1197
1198 const auto &other_bases = bases[other_face];
1199 // Eigen::MatrixXd w;
1200 std::vector<AssemblyValues> w;
1201 other_bases.evaluate_bases(node_position, w);
1202
1203 for (long i = 0; i < w.size(); ++i)
1204 {
1205 assert(w[i].val.size() == 1);
1206 if (std::abs(w[i].val(0)) < 1e-8)
1207 continue;
1208
1209 for (size_t ii = 0; ii < other_bases.bases[i].global().size(); ++ii)
1210 {
1211 const auto &other_global = other_bases.bases[i].global()[ii];
1212 // std::cout<<"e "<<e<<" " <<j << " gid "<<other_global.index<<std::endl;
1213 b.bases[j].global().emplace_back(other_global.index, other_global.node, w[i].val(0) * other_global.val);
1214 }
1215 }
1216 }
1217 }
1218 }
1219 else
1220 {
1221 // Polygon bases are built later on
1222 }
1223 }
1224 }
1225 }
1226 }
1227
1228 return nodes.n_nodes();
1229}
val
double val
Definition Assembler.cpp:86
y
int y
Definition SplineBasis3d.cpp:55
edge_id
int edge_id
Definition SplineBasis3d.cpp:56
x
int x
Definition SplineBasis3d.cpp:55
auto_p_bases.hpp
auto_q_bases.hpp
polyfem::assembler::AssemblerUtils::quadrature_order
static int quadrature_order(const std::string &assembler, const int basis_degree, const BasisType &b_type, const int dim)
utility for retrieving the needed quadrature order to precisely integrate the given form on the given...
Definition AssemblerUtils.cpp:182
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::set_quadrature
void set_quadrature(const QuadratureFunction &fun)
Definition ElementBases.hpp:66
polyfem::basis::ElementBases::set_local_node_from_primitive_func
void set_local_node_from_primitive_func(LocalNodeFromPrimitiveFunc fun)
sets mapping from local nodes to global nodes
Definition ElementBases.hpp:100
polyfem::basis::ElementBases::bases
std::vector< Basis > bases
one basis function per node in the element
Definition ElementBases.hpp:27
polyfem::basis::ElementBases::set_mass_quadrature
void set_mass_quadrature(const QuadratureFunction &fun)
Definition ElementBases.hpp:67
polyfem::basis::LagrangeBasis2d::build_bases
static int build_bases(const mesh::Mesh2D &mesh, const std::string &assembler, const int quadrature_order, const int mass_quadrature_order, const int discr_order, const bool serendipity, const bool has_polys, const bool is_geom_bases, const bool use_corner_quadrature, std::vector< ElementBases > &bases, std::vector< mesh::LocalBoundary > &local_boundary, std::map< int, InterfaceData > &poly_edge_to_data, std::shared_ptr< mesh::MeshNodes > &mesh_nodes)
Builds FE basis functions over the entire mesh (P1, P2 over triangles, Q1, Q2 over quads).
polyfem::basis::LagrangeBasis2d::tri_edge_local_nodes
static Eigen::VectorXi tri_edge_local_nodes(const int p, const mesh::Mesh2D &mesh, mesh::Navigation::Index index)
Definition LagrangeBasis2d.cpp:611
polyfem::basis::LagrangeBasis2d::quad_edge_local_nodes
static Eigen::VectorXi quad_edge_local_nodes(const int q, const mesh::Mesh2D &mesh, mesh::Navigation::Index index)
Definition LagrangeBasis2d.cpp:642
polyfem::mesh::LocalBoundary
Boundary primitive IDs for a single element.
Definition LocalBoundary.hpp:25
polyfem::mesh::Mesh2D::next_around_face
Navigation::Index next_around_face(Navigation::Index idx) const
Definition Mesh2D.hpp:61
polyfem::mesh::Mesh2D::is_volume
bool is_volume() const override
checks if mesh is volume
Definition Mesh2D.hpp:25
polyfem::mesh::Mesh2D::switch_vertex
virtual Navigation::Index switch_vertex(Navigation::Index idx) const =0
polyfem::mesh::Mesh2D::get_index_from_face
virtual Navigation::Index get_index_from_face(int f, int lv=0) const =0
polyfem::mesh::Mesh2D::switch_face
virtual Navigation::Index switch_face(Navigation::Index idx) 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_vertices
virtual int n_vertices() const =0
number of vertices
polyfem::mesh::Mesh::is_polytope
bool is_polytope(const int el_id) const
checks if element is polygon compatible
Definition Mesh.cpp:365
polyfem::mesh::Mesh::is_rational
bool is_rational() const
check if curved mesh has rational polynomials elements
Definition Mesh.hpp:287
polyfem::mesh::Mesh::is_conforming
virtual bool is_conforming() const =0
if the mesh is conforming
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::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:475
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::cell_weights
const std::vector< double > & cell_weights(const int cell_index) const
weights for rational polynomial meshes
Definition Mesh.hpp:550
polyfem::mesh::Mesh::is_boundary_edge
virtual bool is_boundary_edge(const int edge_global_id) const =0
is edge boundary
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::face_vertex
virtual int face_vertex(const int f_id, const int lv_id) const =0
id of the face vertex
polyfem::mesh::NCMesh2D::face_edge
int face_edge(const int f_id, const int le_id) const
Definition NCMesh2D.hpp:195
polyfem::mesh::NCMesh2D::element_weight_to_edge_weight
static double element_weight_to_edge_weight(const int l, const Eigen::Vector2d &pos)
Definition NCMesh2D.cpp:427
polyfem::mesh::NCMesh2D::n_faces
int n_faces() const override
number of faces
Definition NCMesh2D.hpp:168
polyfem::mesh::NCMesh2D::leader_edge_of_edge
int leader_edge_of_edge(const int e_id) const
Definition NCMesh2D.hpp:201
polyfem::mesh::NCMesh2D::face_ref_level
int face_ref_level(const int f_id) const
Definition NCMesh2D.hpp:189
polyfem::mesh::NCMesh2D::n_face_vertices
int n_face_vertices(const int f_id) const override
number of vertices of a face
Definition NCMesh2D.hpp:187
polyfem::mesh::NCMesh2D::n_edges
int n_edges() const override
number of edges
Definition NCMesh2D.hpp:169
polyfem::quadrature::QuadQuadrature::get_quadrature
void get_quadrature(const int order, Quadrature &quad)
Definition QuadQuadrature.cpp:16
JIXIE::sort
enable_if_t< t==0 > sort(Eigen::Matrix< T, 3, 3 > &U, Eigen::Matrix< T, 3, 1 > &sigma, Eigen::Matrix< T, 3, 3 > &V)
Helper function of 3X3 SVD for sorting singular values.
Definition ImplicitQRSVD.h:600
elasticity_rhs.f
f
Definition elasticity_rhs.py:173
generate_mooney_rivlin.grad
grad
Definition generate_mooney_rivlin.py:51
generate_mooney_rivlin.dim
int dim
Definition generate_mooney_rivlin.py:46
generate_rotation_mtx.e
list e
Definition generate_rotation_mtx.py:49
p_bases.tmp
list tmp
Definition p_bases.py:339
p_bases.nn
nn
Definition p_bases.py:315
p_bases.nodes
str nodes
Definition p_bases.py:372
p_bases.indices
list indices
Definition p_bases.py:232
p_bases.N
N
Definition p_bases.py:222
polyfem::assembler
Used for test only.
Definition AMIPSEnergy.cpp:6
polyfem::autogen::p_grad_basis_value_2d
void p_grad_basis_value_2d(const int p, const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val)
Definition auto_p_bases.cpp:882
polyfem::autogen::q_grad_basis_value_2d
void q_grad_basis_value_2d(const int q, const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val)
Definition auto_q_bases_2d_grad.cpp:182
polyfem::autogen::p_nodes_2d
void p_nodes_2d(const int p, Eigen::MatrixXd &val)
Definition auto_p_bases.cpp:835
polyfem::autogen::q_basis_value_2d
void q_basis_value_2d(const int q, const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val)
Definition auto_q_bases_2d_val.cpp:88
polyfem::autogen::p_basis_value_2d
void p_basis_value_2d(const int p, const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val)
Definition auto_p_bases.cpp:858
polyfem::utils::maybe_parallel_for
void maybe_parallel_for(int size, const std::function< void(int, int, int)> &partial_for)
polyfem::basis::InterfaceData::local_indices
std::vector< int > local_indices
Definition InterfaceData.hpp:31