Loading [MathJax]/extensions/tex2jax.js
PolyFEM
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
LagrangeBasis3d.cpp
Go to the documentation of this file.
1
3#include "LagrangeBasis3d.hpp"
4
5#include <polyfem/mesh/MeshNodes.hpp>
6#include <polyfem/quadrature/TetQuadrature.hpp>
7#include <polyfem/quadrature/HexQuadrature.hpp>
8
9#include <polyfem/assembler/AssemblerUtils.hpp>
10
11#include <polyfem/autogen/auto_p_bases.hpp>
12#include <polyfem/autogen/auto_q_bases.hpp>
13
14#include <polyfem/utils/MaybeParallelFor.hpp>
15
16#include <cassert>
17#include <array>
19
20using namespace polyfem;
21using namespace polyfem::assembler;
22using namespace polyfem::basis;
23using namespace polyfem::mesh;
24using namespace polyfem::quadrature;
25
26/*
27Axes:
28 z y
29 │╱
30 o──x
31Boundaries:
32X axis: left/right
33Y axis: front/back
34Z axis: bottom/top
35Corner nodes:
36v3──x──v2
37 │⋱ ╱ ╲
38 x ⋱╱ ╲
39 │ x x x
40 │ ╱ ⋱ ╲
41 │╱ ⋱╲
42v0─────x──────v1
43v0 = (0, 0, 0)
44v1 = (1, 0, 0)
45v2 = (0, 1, 0)
46v3 = (0, 0, 1)
47Edge nodes:
48e0 = (0.5, 0, 0)
49e1 = (0.5, 0.5, 0)
50e2 = ( 0, 0.5, 0)
51e3 = ( 0, 0, 0.5)
52e4 = (0.5, 0, 0.5)
53e5 = ( 0, 0.5, 0.5)
54Corner nodes:
55 v7──────x─────v6
56 ╱┆ ╱ ╱│
57 ╱ ┆ ╱ ╱ │
58 x┄┄┼┄┄┄x┄┄┄┄┄┄x │
59 ╱ x ╱ ╱┆ x
60 ╱ ┆ ╱ ╱ ┆ ╱│
61v4─────┼x─────v5 ┆╱ │
62 │ ┆┆ │ x │
63 │ v3┼┄┄┄┄┄x┼┄⌿┼┄v2
64 │ ╱ ┆ │╱ ┆ ╱
65 x┄┄┄⌿┄┄x┄┄┄┄┄┄x ┆╱
66 │ x ┆ │ x
67 │ ╱ ┆ │ ╱
68 │╱ ┆ │╱
69v0──────x─────v1
70v0 = (0, 0, 0)
71v1 = (1, 0, 0)
72v2 = (1, 1, 0)
73v3 = (0, 1, 0)
74v4 = (0, 0, 1)
75v5 = (1, 0, 1)
76v6 = (1, 1, 1)
77v7 = (0, 1, 1)
78Edge nodes:
79 x─────e10─────x
80 ╱┆ ╱ ╱│
81 ╱ ┆ ╱ ╱ │
82 e11┄┼┄┄┄x┄┄┄┄┄e9 │
83 ╱ e7 ╱ ╱┆ e6
84 ╱ ┆ ╱ ╱ ┆ ╱│
85 x─────e8──────x ┆╱ │
86 │ ┆┆ │ x │
87 │ x┼┄┄┄┄┄e2┄⌿┼┄┄x
88 │ ╱ ┆ │╱ ┆ ╱
89e4┄┄┄⌿┄┄x┄┄┄┄┄e5 ┆╱
90 │ e3 ┆ │ e1
91 │ ╱ ┆ │ ╱
92 │╱ ┆ │╱
93 x─────e0──────x
94e0 = (0.5, 0, 0)
95e1 = ( 1, 0.5, 0)
96e2 = (0.5, 1, 0)
97e3 = ( 0, 0.5, 0)
98e4 = ( 0, 0, 0.5)
99e5 = ( 1, 0, 0.5)
100e6 = ( 1, 1, 0.5)
101e7 = ( 0, 1, 0.5)
102e8 = (0.5, 0, 1)
103e9 = ( 1, 0.5, 1)
104e10 = (0.5, 1, 1)
105e11 = ( 0, 0.5, 1)
106Face nodes:
107 v7──────x─────v6
108 ╱┆ ╱ ╱│
109 ╱ ┆ ╱ ╱ │
110 x┄┄┼┄┄f5┄┄┄┄┄┄x │
111 ╱ x ╱ f3 ╱┆ x
112 ╱ ┆ ╱ ╱ ┆ ╱│
113v4─────┼x─────v5 ┆╱ │
114 │ f0 ┆┆ │ f1 │
115 │ v3┼┄┄┄┄┄x┼┄⌿┼┄v2
116 │ ╱ ┆ │╱ ┆ ╱
117 x┄┄┄⌿┄f2┄┄┄┄┄┄x ┆╱
118 │ x ┆ f4 │ x
119 │ ╱ ┆ │ ╱
120 │╱ ┆ │╱
121v0──────x─────v1
122f0 = ( 0, 0.5, 0.5)
123f1 = ( 1, 0.5, 0.5)
124f2 = (0.5, 0, 0.5)
125f3 = (0.5, 1, 0.5)
126f4 = (0.5, 0.5, 0)
127f5 = (0.5, 0.5, 1)
128*/
129
130namespace
131{
132
133 template <class InputIterator, class T>
134 int find_index(InputIterator first, InputIterator last, const T &val)
135 {
136 return std::distance(first, std::find(first, last, val));
137 }
138
139 Navigation3D::Index find_quad_face(const Mesh3D &mesh, int c, int v1, int v2, int v3, int v4)
140 {
141 std::array<int, 4> v = {{v1, v2, v3, v4}};
142 std::sort(v.begin(), v.end());
143 for (int lf = 0; lf < mesh.n_cell_faces(c); ++lf)
144 {
145 auto idx = mesh.get_index_from_element(c, lf, 0);
146 assert(mesh.n_face_vertices(idx.face) == 4);
147 std::array<int, 4> u;
148 for (int lv = 0; lv < mesh.n_face_vertices(idx.face); ++lv)
149 {
150 u[lv] = idx.vertex;
151 idx = mesh.next_around_face(idx);
152 }
153 std::sort(u.begin(), u.end());
154 if (u == v)
155 {
156 return idx;
157 }
158 }
159 assert(false);
160 return Navigation3D::Index();
161 }
162
163 std::array<int, 4> tet_vertices_local_to_global(const Mesh3D &mesh, int c)
164 {
165 // Vertex nodes
166 assert(mesh.is_simplex(c));
167 std::array<int, 4> l2g;
168 int lv = 0;
169 for (int vi : mesh.get_ordered_vertices_from_tet(c))
170 {
171 l2g[lv++] = vi;
172 }
173
174 return l2g;
175 }
176
177 std::array<int, 8> hex_vertices_local_to_global(const Mesh3D &mesh, int c)
178 {
179 assert(mesh.is_cube(c));
180
181 // Vertex nodes
182 std::array<int, 8> l2g;
183 int lv = 0;
184 for (int vi : mesh.get_ordered_vertices_from_hex(c))
185 {
186 l2g[lv++] = vi;
187 }
188
189 return l2g;
190 }
191
192 int lowest_order_elem_on_edge(const polyfem::mesh::NCMesh3D &mesh, const Eigen::VectorXi &discr_orders, const int eid)
193 {
194 auto elem_list = mesh.edge_neighs(eid);
195 int min = std::numeric_limits<int>::max();
196 int elem = -1;
197 for (const auto e : elem_list)
198 if (discr_orders[e] < min)
199 elem = e;
200 return elem;
201 }
202
203 void tet_local_to_global(const bool is_geom_bases, const int p, const Mesh3D &mesh, int c, const Eigen::VectorXi &discr_order, const Eigen::VectorXi &edge_orders, const Eigen::VectorXi &face_orders, std::vector<int> &res, polyfem::mesh::MeshNodes &nodes, std::vector<std::vector<int>> &edge_virtual_nodes, std::vector<std::vector<int>> &face_virtual_nodes)
204 {
205 const int n_edge_nodes = p > 1 ? ((p - 1) * 6) : 0;
206 const int nn = p > 2 ? (p - 2) : 0;
207 const int n_loc_f = (nn * (nn + 1) / 2);
208 const int n_face_nodes = n_loc_f * 4;
209 int n_cell_nodes = 0;
210 for (int pp = 4; pp <= p; ++pp)
211 n_cell_nodes += ((pp - 3) * ((pp - 3) + 1) / 2);
212
213 if (p == 0)
214 {
215 res.push_back(nodes.node_id_from_cell(c));
216 return;
217 }
218
219 // std::vector<int> res;
220 res.reserve(4 + n_edge_nodes + n_face_nodes + n_cell_nodes);
221
222 // Edge nodes
223 Eigen::Matrix<Navigation3D::Index, 4, 1> f;
224
225 auto v = tet_vertices_local_to_global(mesh, c);
226 Eigen::Matrix<int, 4, 3> fv;
227 fv.row(0) << v[0], v[1], v[2];
228 fv.row(1) << v[0], v[1], v[3];
229 fv.row(2) << v[1], v[2], v[3];
230 fv.row(3) << v[2], v[0], v[3];
231
232 for (long lf = 0; lf < fv.rows(); ++lf)
233 {
234 const auto index = mesh.get_index_from_element_face(c, fv(lf, 0), fv(lf, 1), fv(lf, 2));
235 f[lf] = index;
236 }
237
238 Eigen::Matrix<Navigation3D::Index, 6, 1> e;
239 Eigen::Matrix<int, 6, 2> ev;
240 ev.row(0) << v[0], v[1];
241 ev.row(1) << v[1], v[2];
242 ev.row(2) << v[2], v[0];
243
244 ev.row(3) << v[0], v[3];
245 ev.row(4) << v[1], v[3];
246 ev.row(5) << v[2], v[3];
247
248 for (int le = 0; le < e.rows(); ++le)
249 {
250 // const auto index = find_edge(mesh, c, ev(le, 0), ev(le, 1));
251 const auto index = mesh.get_index_from_element_edge(c, ev(le, 0), ev(le, 1));
252 e[le] = index;
253 }
254
255 // vertices
256 for (size_t lv = 0; lv < v.size(); ++lv)
257 {
258 if (!mesh.is_conforming() && !is_geom_bases)
259 {
260 const auto &ncmesh = dynamic_cast<const NCMesh3D &>(mesh);
261 // hanging vertex
262 if (ncmesh.leader_edge_of_vertex(v[lv]) >= 0 || ncmesh.leader_face_of_vertex(v[lv]) >= 0)
263 res.push_back(-lv - 1);
264 else
265 res.push_back(nodes.node_id_from_primitive(v[lv]));
266 }
267 else
268 res.push_back(nodes.node_id_from_primitive(v[lv]));
269 }
270
271 // Edges
272 for (int le = 0; le < e.rows(); ++le)
273 {
274 const auto index = e[le];
275 auto neighs = mesh.edge_neighs(index.edge);
276 int min_p = discr_order.size() > 0 ? discr_order(c) : 0;
277
278 if (is_geom_bases)
279 {
280 auto node_ids = nodes.node_ids_from_edge(index, p - 1);
281 res.insert(res.end(), node_ids.begin(), node_ids.end());
282 }
283 else
284 {
285 if (!mesh.is_conforming() && !is_geom_bases)
286 {
287 const auto &ncmesh = dynamic_cast<const NCMesh3D &>(mesh);
288 // slave edge
289 if (ncmesh.leader_edge_of_edge(index.edge) >= 0 || ncmesh.leader_face_of_edge(index.edge) >= 0)
290 {
291 for (int tmp = 0; tmp < p - 1; ++tmp)
292 res.push_back(-le - 1);
293 }
294 // master or conforming edge with constrained order
295 else if (edge_orders[index.edge] < discr_order(c))
296 {
297 for (int tmp = 0; tmp < p - 1; ++tmp)
298 res.push_back(-le - 1);
299
300 int min_order_elem = lowest_order_elem_on_edge(ncmesh, discr_order, index.edge);
301 // master edge, add extra nodes
302 if (min_order_elem == c)
303 edge_virtual_nodes[index.edge] = nodes.node_ids_from_edge(index, edge_orders[index.edge] - 1);
304 }
305 else
306 {
307 auto node_ids = nodes.node_ids_from_edge(index, p - 1);
308 res.insert(res.end(), node_ids.begin(), node_ids.end());
309 }
310 }
311 else
312 {
313 for (auto cid : neighs)
314 {
315 min_p = std::min(min_p, discr_order.size() > 0 ? discr_order(cid) : 0);
316 }
317
318 if (discr_order.size() > 0 && discr_order(c) > min_p)
319 {
320 for (int tmp = 0; tmp < p - 1; ++tmp)
321 res.push_back(-le - 10);
322 }
323 else
324 {
325 auto node_ids = nodes.node_ids_from_edge(index, p - 1);
326 res.insert(res.end(), node_ids.begin(), node_ids.end());
327 }
328 }
329 }
330 }
331
332 // faces
333 for (int lf = 0; lf < f.rows(); ++lf)
334 {
335 const auto index = f[lf];
336 const auto other_cell = mesh.switch_element(index).element;
337
338 const bool skip_other = discr_order.size() > 0 && other_cell >= 0 && discr_order(c) > discr_order(other_cell);
339
340 if (is_geom_bases)
341 {
342 auto node_ids = nodes.node_ids_from_face(index, p - 2);
343 res.insert(res.end(), node_ids.begin(), node_ids.end());
344 }
345 else
346 {
347 if (!mesh.is_conforming() && !is_geom_bases)
348 {
349 const auto &ncmesh = dynamic_cast<const NCMesh3D &>(mesh);
350 // slave face
351 if (ncmesh.leader_face_of_face(index.face) >= 0)
352 {
353 for (int tmp = 0; tmp < n_loc_f; ++tmp)
354 res.push_back(-lf - 1);
355 }
356 // master face or conforming face with constrained order
357 else if (face_orders[index.face] < discr_order[c])
358 {
359 for (int tmp = 0; tmp < n_loc_f; ++tmp)
360 res.push_back(-lf - 1);
361 // master face
362 if (ncmesh.n_follower_faces(index.face) > 0 && face_orders[index.face] > 2)
363 face_virtual_nodes[index.face] = nodes.node_ids_from_face(index, face_orders[index.face] - 2);
364 }
365 else
366 {
367 auto node_ids = nodes.node_ids_from_face(index, p - 2);
368 res.insert(res.end(), node_ids.begin(), node_ids.end());
369 }
370 }
371 else
372 {
373 if (skip_other)
374 {
375 for (int tmp = 0; tmp < n_loc_f; ++tmp)
376 res.push_back(-lf - 1);
377 }
378 else
379 {
380 auto node_ids = nodes.node_ids_from_face(index, p - 2);
381 res.insert(res.end(), node_ids.begin(), node_ids.end());
382 }
383 }
384 }
385 }
386
387 // cells
388 if (n_cell_nodes > 0)
389 {
390 const auto index = f[0];
391
392 auto node_ids = nodes.node_ids_from_cell(index, p - 3);
393 res.insert(res.end(), node_ids.begin(), node_ids.end());
394 }
395
396 assert(res.size() == size_t(4 + n_edge_nodes + n_face_nodes + n_cell_nodes));
397 }
398
399 void hex_local_to_global(const bool serendipity, const int q, const Mesh3D &mesh, int c, const Eigen::VectorXi &discr_order, std::vector<int> &res, MeshNodes &nodes)
400 {
401 assert(mesh.is_cube(c));
402
403 const int n_edge_nodes = ((q - 1) * 12);
404 const int nn = (q - 1);
405 const int n_loc_f = serendipity ? 0 : (nn * nn);
406 const int n_face_nodes = serendipity ? 0 : (n_loc_f * 6);
407 const int n_cell_nodes = serendipity ? 0 : (nn * nn * nn);
408
409 if (q == 0)
410 {
411 res.push_back(nodes.node_id_from_cell(c));
412 return;
413 }
414
415 // std::vector<int> res;
416 res.reserve(8 + n_edge_nodes + n_face_nodes + n_cell_nodes);
417
418 // Vertex nodes
419 auto v = hex_vertices_local_to_global(mesh, c);
420
421 // Edge nodes
422 Eigen::Matrix<Navigation3D::Index, 12, 1> e;
423 Eigen::Matrix<int, 12, 2> ev;
424 ev.row(0) << v[0], v[1];
425 ev.row(1) << v[1], v[2];
426 ev.row(2) << v[2], v[3];
427 ev.row(3) << v[3], v[0];
428 ev.row(4) << v[0], v[4];
429 ev.row(5) << v[1], v[5];
430 ev.row(6) << v[2], v[6];
431 ev.row(7) << v[3], v[7];
432 ev.row(8) << v[4], v[5];
433 ev.row(9) << v[5], v[6];
434 ev.row(10) << v[6], v[7];
435 ev.row(11) << v[7], v[4];
436 for (int le = 0; le < e.rows(); ++le)
437 {
438 // e[le] = find_edge(mesh, c, ev(le, 0), ev(le, 1)).edge;
439 e[le] = mesh.get_index_from_element_edge(c, ev(le, 0), ev(le, 1));
440 }
441
442 // Face nodes
443 Eigen::Matrix<Navigation3D::Index, 6, 1> f;
444 Eigen::Matrix<int, 6, 4> fv;
445 fv.row(0) << v[0], v[3], v[4], v[7];
446 fv.row(1) << v[1], v[2], v[5], v[6];
447 fv.row(2) << v[0], v[1], v[5], v[4];
448 fv.row(3) << v[3], v[2], v[6], v[7];
449 fv.row(4) << v[0], v[1], v[2], v[3];
450 fv.row(5) << v[4], v[5], v[6], v[7];
451 for (int lf = 0; lf < f.rows(); ++lf)
452 {
453 const auto index = find_quad_face(mesh, c, fv(lf, 0), fv(lf, 1), fv(lf, 2), fv(lf, 3));
454 f[lf] = index;
455 }
456
457 // vertices
458 for (size_t lv = 0; lv < v.size(); ++lv)
459 {
460 res.push_back(nodes.node_id_from_primitive(v[lv]));
461 }
462 assert(res.size() == size_t(8));
463
464 // Edges
465 for (int le = 0; le < e.rows(); ++le)
466 {
467 const auto index = e[le];
468 auto neighs = mesh.edge_neighs(index.edge);
469 int min_q = discr_order.size() > 0 ? discr_order(c) : 0;
470
471 for (auto cid : neighs)
472 {
473 min_q = std::min(min_q, discr_order.size() > 0 ? discr_order(cid) : 0);
474 }
475
476 if (discr_order.size() > 0 && discr_order(c) > min_q)
477 {
478 for (int tmp = 0; tmp < q - 1; ++tmp)
479 res.push_back(-le - 10);
480 }
481 else
482 {
483 auto node_ids = nodes.node_ids_from_edge(index, q - 1);
484 res.insert(res.end(), node_ids.begin(), node_ids.end());
485 }
486 }
487 assert(res.size() == size_t(8 + n_edge_nodes));
488
489 // faces
490 for (int lf = 0; lf < f.rows(); ++lf)
491 {
492 const auto index = f[lf];
493 const auto other_cell = mesh.switch_element(index).element;
494
495 const bool skip_other = discr_order.size() > 0 && other_cell >= 0 && discr_order(c) > discr_order(other_cell);
496
497 if (skip_other)
498 {
499 for (int tmp = 0; tmp < n_loc_f; ++tmp)
500 res.push_back(-lf - 1);
501 }
502 else
503 {
504 auto node_ids = nodes.node_ids_from_face(index, serendipity ? 0 : (q - 1));
505 assert(node_ids.size() == n_loc_f);
506 res.insert(res.end(), node_ids.begin(), node_ids.end());
507 }
508 }
509 assert(res.size() == size_t(8 + n_edge_nodes + n_face_nodes));
510
511 // cells
512 if (n_cell_nodes > 0)
513 {
514 const auto index = f[0];
515
516 auto node_ids = nodes.node_ids_from_cell(index, q - 1);
517 res.insert(res.end(), node_ids.begin(), node_ids.end());
518 }
519
520 assert(res.size() == size_t(8 + n_edge_nodes + n_face_nodes + n_cell_nodes));
521 }
522
523 // -----------------------------------------------------------------------------
524
538 void compute_nodes(
539 const Mesh3D &mesh,
540 const Eigen::VectorXi &discr_orders,
541 const Eigen::VectorXi &edge_orders,
542 const Eigen::VectorXi &face_orders,
543 const bool serendipity,
544 const bool has_polys,
545 const bool is_geom_bases,
546 MeshNodes &nodes,
547 std::vector<std::vector<int>> &edge_virtual_nodes,
548 std::vector<std::vector<int>> &face_virtual_nodes,
549 std::vector<std::vector<int>> &element_nodes_id,
550 std::vector<LocalBoundary> &local_boundary,
551 std::map<int, InterfaceData> &poly_face_to_data)
552 {
553 // Step 1: Assign global node ids for each quads
554 local_boundary.clear();
555 // local_boundary.resize(mesh.n_faces());
556 element_nodes_id.resize(mesh.n_faces());
557
558 if (!mesh.is_conforming())
559 {
560 const auto &ncmesh = dynamic_cast<const NCMesh3D &>(mesh);
561 edge_virtual_nodes.resize(ncmesh.n_edges());
562 face_virtual_nodes.resize(ncmesh.n_faces());
563 }
564
565 for (int c = 0; c < mesh.n_cells(); ++c)
566 {
567 const int discr_order = discr_orders(c);
568
569 if (mesh.is_cube(c))
570 {
571 hex_local_to_global(serendipity, discr_order, mesh, c, discr_orders, element_nodes_id[c], nodes);
572
573 auto v = hex_vertices_local_to_global(mesh, c);
574 Eigen::Matrix<int, 6, 4> fv;
575 fv.row(0) << v[0], v[3], v[4], v[7];
576 fv.row(1) << v[1], v[2], v[5], v[6];
577 fv.row(2) << v[0], v[1], v[5], v[4];
578 fv.row(3) << v[3], v[2], v[6], v[7];
579 fv.row(4) << v[0], v[1], v[2], v[3];
580 fv.row(5) << v[4], v[5], v[6], v[7];
581
582 LocalBoundary lb(c, BoundaryType::QUAD);
583 for (int i = 0; i < fv.rows(); ++i)
584 {
585 const int f = find_quad_face(mesh, c, fv(i, 0), fv(i, 1), fv(i, 2), fv(i, 3)).face;
586
587 if (mesh.is_boundary_face(f) || mesh.get_boundary_id(f) > 0)
588 {
589 lb.add_boundary_primitive(f, i);
590 }
591 }
592
593 if (!lb.empty())
594 local_boundary.emplace_back(lb);
595 }
596 else if (mesh.is_simplex(c))
597 {
598 // element_nodes_id[c] = polyfem::LagrangeBasis3d::tet_local_to_global(discr_order, mesh, c, discr_orders, nodes);
599 tet_local_to_global(is_geom_bases, discr_order, mesh, c, discr_orders, edge_orders, face_orders, element_nodes_id[c], nodes, edge_virtual_nodes, face_virtual_nodes);
600
601 auto v = tet_vertices_local_to_global(mesh, c);
602 Eigen::Matrix<int, 4, 3> fv;
603 fv.row(0) << v[0], v[1], v[2];
604 fv.row(1) << v[0], v[1], v[3];
605 fv.row(2) << v[1], v[2], v[3];
606 fv.row(3) << v[2], v[0], v[3];
607
608 LocalBoundary lb(c, BoundaryType::TRI);
609 for (long i = 0; i < fv.rows(); ++i)
610 {
611 const int f = mesh.get_index_from_element_face(c, fv(i, 0), fv(i, 1), fv(i, 2)).face;
612
613 if (mesh.is_boundary_face(f))
614 {
615 lb.add_boundary_primitive(f, i);
616 }
617 }
618
619 if (!lb.empty())
620 local_boundary.emplace_back(lb);
621 }
622 }
623
624 if (!has_polys)
625 return;
626
627 // Step 2: Iterate over edges of polygons and compute interface weights
628 Eigen::VectorXi indices;
629 for (int c = 0; c < mesh.n_cells(); ++c)
630 {
631 // Skip non-polytopes
632 if (!mesh.is_polytope(c))
633 {
634 continue;
635 }
636
637 for (int lf = 0; lf < mesh.n_cell_faces(c); ++lf)
638 {
639 auto index = mesh.get_index_from_element(c, lf, 0);
640 auto index2 = mesh.switch_element(index);
641 int c2 = index2.element;
642 assert(c2 >= 0);
643
644 const int discr_order = discr_orders(c2);
645 if (mesh.is_cube(c2))
646 {
647 indices = LagrangeBasis3d::hex_face_local_nodes(serendipity, discr_order, mesh, index2);
648 }
649 else if (mesh.is_simplex(c2))
650 {
651 indices = LagrangeBasis3d::tet_face_local_nodes(discr_order, mesh, index2);
652 }
653 else
654 continue;
655
656 InterfaceData data;
657 data.local_indices.insert(data.local_indices.begin(), indices.data(), indices.data() + indices.size());
658 assert(indices.size() == data.local_indices.size());
659 poly_face_to_data[index2.face] = data;
660 }
661 }
662 }
669 void local_to_global(const Eigen::MatrixXd &verts, const Eigen::MatrixXd &uv, Eigen::MatrixXd &pts)
670 {
671 const int dim = verts.cols();
672 const int N = uv.rows();
673 assert(dim == 3);
674 assert(uv.cols() == dim);
675 assert(verts.rows() == dim + 1);
676
677 pts.setZero(N, dim);
678 for (int i = 0; i < N; i++)
679 pts.row(i) = uv(i, 0) * verts.row(1) + uv(i, 1) * verts.row(2) + uv(i, 2) * verts.row(3) + (1.0 - uv(i, 0) - uv(i, 1) - uv(i, 2)) * verts.row(0);
680 }
681
682 void local_to_global_face(const Eigen::MatrixXd &verts, const Eigen::MatrixXd &uv, Eigen::MatrixXd &pts)
683 {
684 const int dim = verts.cols();
685 const int N = uv.rows();
686 assert(dim == 3);
687 assert(uv.cols() == 2);
688 assert(verts.rows() == 3);
689
690 pts.setZero(N, dim);
691 for (int i = 0; i < N; i++)
692 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);
693 }
694
701 void global_to_local(const Eigen::MatrixXd &verts, const Eigen::MatrixXd &pts, Eigen::MatrixXd &uv)
702 {
703 const int dim = verts.cols();
704 const int N = pts.rows();
705 assert(dim == 3);
706 assert(verts.rows() == dim + 1);
707 assert(pts.cols() == dim);
708
709 Eigen::Matrix3d J;
710 for (int i = 0; i < dim; i++)
711 J.col(i) = verts.row(i + 1) - verts.row(0);
712
713 Eigen::Matrix3d Jinv = J.inverse();
714
715 uv.setZero(N, dim);
716 polyfem::utils::maybe_parallel_for(N, [&](int start, int end, int thread_id) {
717 for (int i = start; i < end; i++)
718 {
719 auto point = pts.row(i) - verts.row(0);
720 uv.row(i) = Jinv * point.transpose();
721 }
722 });
723 }
724
725 void global_to_local_face(const Eigen::MatrixXd &verts, const Eigen::MatrixXd &pts, Eigen::MatrixXd &uv)
726 {
727 const int dim = verts.cols();
728 const int N = pts.rows();
729 assert(dim == 3);
730 assert(verts.rows() == 3);
731 assert(pts.cols() == dim);
732
733 Eigen::Matrix3d J;
734 for (int i = 0; i < 2; i++)
735 J.col(i) = verts.row(i + 1) - verts.row(0);
736
737 Eigen::Vector3d a = J.col(0);
738 Eigen::Vector3d b = J.col(1);
739 Eigen::Vector3d virtual_vert = a.cross(b);
740 J.col(2) = virtual_vert;
741
742 uv.setZero(N, 2);
743 polyfem::utils::maybe_parallel_for(N, [&](int start, int end, int thread_id) {
744 for (int i = start; i < end; i++)
745 {
746 auto point = pts.row(i) - verts.row(0);
747 Eigen::Vector3d x = J.colPivHouseholderQr().solve(point.transpose());
748 uv.row(i) = x.block(0, 0, 2, 1);
749 assert(std::abs(x(2)) < 1e-8);
750 }
751 });
752 }
753
754 void global_to_local_edge(const Eigen::MatrixXd &verts, const Eigen::MatrixXd &pts, Eigen::VectorXd &uv)
755 {
756 const int dim = verts.cols();
757 const int N = pts.rows();
758 assert(dim == 3);
759 assert(verts.rows() == 2);
760 assert(pts.cols() == dim);
761
762 auto edge = verts.row(1) - verts.row(0);
763 double squared_length = edge.squaredNorm();
764
765 uv.setZero(N);
766 polyfem::utils::maybe_parallel_for(N, [&](int start, int end, int thread_id) {
767 for (int i = start; i < end; i++)
768 {
769 auto vec = pts.row(i) - verts.row(0);
770 uv(i) = (vec.dot(edge)) / squared_length;
771 }
772 });
773 }
774
775 bool check_edge_face_orders(const polyfem::mesh::NCMesh3D &mesh, const Eigen::VectorXi &elem_orders, const Eigen::VectorXi &edge_orders, const Eigen::VectorXi &face_orders)
776 {
777 // same order for overlapping faces
778 for (int i = 0; i < mesh.n_faces(); i++)
779 if (mesh.leader_face_of_face(i) >= 0)
780 if (face_orders[mesh.leader_face_of_face(i)] != face_orders[i])
781 return false;
782
783 // face order no smaller than order of its edges
784 for (int i = 0; i < mesh.n_faces(); i++)
785 {
786 if (mesh.n_face_cells(i) == 0)
787 continue;
788 for (int j = 0; j < mesh.n_face_vertices(i); j++)
789 {
790 const int e_id = mesh.face_edge(i, j);
791 if (edge_orders[e_id] > face_orders[i])
792 return false;
793 }
794 }
795
796 // same order for overlapping edges
797 for (int i = 0; i < mesh.n_edges(); i++)
798 if (mesh.leader_edge_of_edge(i) >= 0)
799 if (edge_orders[mesh.leader_edge_of_edge(i)] != edge_orders[i])
800 return false;
801
802 // face order no larger than order of interior edges
803 for (int i = 0; i < mesh.n_edges(); i++)
804 {
805 if (mesh.n_edge_cells(i) == 0)
806 continue;
807 if (mesh.leader_face_of_edge(i) >= 0 && mesh.leader_edge_of_edge(i) < 0)
808 if (face_orders[mesh.leader_face_of_edge(i)] > edge_orders[i])
809 return false;
810 }
811 return true;
812 }
813
821 void compute_edge_face_orders(const polyfem::mesh::NCMesh3D &mesh, const Eigen::VectorXi &elem_orders, Eigen::VectorXi &edge_orders, Eigen::VectorXi &face_orders)
822 {
823 const int max_order = elem_orders.maxCoeff();
824 edge_orders.setConstant(mesh.n_edges(), max_order);
825 face_orders.setConstant(mesh.n_faces(), max_order);
826
827 for (int i = 0; i < mesh.n_cells(); i++)
828 for (int j = 0; j < mesh.n_cell_faces(i); j++)
829 face_orders[mesh.cell_face(i, j)] = std::min(face_orders[mesh.cell_face(i, j)], elem_orders[i]);
830
831 for (int i = 0; i < mesh.n_cells(); i++)
832 for (int j = 0; j < mesh.n_cell_edges(i); j++)
833 edge_orders[mesh.cell_edge(i, j)] = std::min(edge_orders[mesh.cell_edge(i, j)], elem_orders[i]);
834
835 while (!check_edge_face_orders(mesh, elem_orders, edge_orders, face_orders))
836 {
837 // same order for overlapping faces
838 for (int i = 0; i < mesh.n_faces(); i++)
839 if (mesh.leader_face_of_face(i) >= 0)
840 face_orders[mesh.leader_face_of_face(i)] = std::min(face_orders[mesh.leader_face_of_face(i)], face_orders[i]);
841
842 for (int i = 0; i < mesh.n_faces(); i++)
843 if (mesh.leader_face_of_face(i) >= 0)
844 face_orders[i] = std::min(face_orders[mesh.leader_face_of_face(i)], face_orders[i]);
845
846 // face order no smaller than order of its edges
847 for (int i = 0; i < mesh.n_faces(); i++)
848 {
849 if (mesh.n_face_cells(i) == 0)
850 continue;
851 for (int j = 0; j < mesh.n_face_vertices(i); j++)
852 {
853 const int e_id = mesh.face_edge(i, j);
854 edge_orders[e_id] = std::min(edge_orders[e_id], face_orders[i]);
855 }
856 }
857
858 // same order for overlapping edges
859 for (int i = 0; i < mesh.n_edges(); i++)
860 if (mesh.leader_edge_of_edge(i) >= 0)
861 edge_orders[mesh.leader_edge_of_edge(i)] = std::min(edge_orders[mesh.leader_edge_of_edge(i)], edge_orders[i]);
862
863 for (int i = 0; i < mesh.n_edges(); i++)
864 if (mesh.leader_edge_of_edge(i) >= 0)
865 edge_orders[i] = std::min(edge_orders[mesh.leader_edge_of_edge(i)], edge_orders[i]);
866
867 // face order no larger than order of interior edges
868 for (int i = 0; i < mesh.n_edges(); i++)
869 {
870 if (mesh.n_edge_cells(i) == 0)
871 continue;
872 if (mesh.leader_face_of_edge(i) >= 0 && mesh.leader_edge_of_edge(i) < 0)
873 face_orders[mesh.leader_face_of_edge(i)] = std::min(face_orders[mesh.leader_face_of_edge(i)], edge_orders[i]);
874 }
875 }
876 }
877} // anonymous namespace
878
879Eigen::VectorXi LagrangeBasis3d::tet_face_local_nodes(const int p, const Mesh3D &mesh, Navigation3D::Index index)
880{
881 const int nn = p > 2 ? (p - 2) : 0;
882 const int n_edge_nodes = (p - 1) * 6;
883 const int n_face_nodes = nn * (nn + 1) / 2;
884
885 const int c = index.element;
886 assert(mesh.is_simplex(c));
887
888 // Local to global mapping of node indices
889 const auto l2g = tet_vertices_local_to_global(mesh, c);
890
891 // Extract requested interface
892 Eigen::VectorXi result(3 + (p - 1) * 3 + n_face_nodes);
893 result[0] = find_index(l2g.begin(), l2g.end(), index.vertex);
894 result[1] = find_index(l2g.begin(), l2g.end(), mesh.next_around_face(index).vertex);
895 result[2] = find_index(l2g.begin(), l2g.end(), mesh.next_around_face(mesh.next_around_face(index)).vertex);
896
897 Eigen::Matrix<Navigation3D::Index, 6, 1> e;
898 Eigen::Matrix<int, 6, 2> ev;
899 ev.row(0) << l2g[0], l2g[1];
900 ev.row(1) << l2g[1], l2g[2];
901 ev.row(2) << l2g[2], l2g[0];
902
903 ev.row(3) << l2g[0], l2g[3];
904 ev.row(4) << l2g[1], l2g[3];
905 ev.row(5) << l2g[2], l2g[3];
906
907 Navigation3D::Index tmp = index;
908
909 for (int le = 0; le < e.rows(); ++le)
910 {
911 // const auto index = find_edge(mesh, c, ev(le, 0), ev(le, 1));
912 const auto l_index = mesh.get_index_from_element_edge(c, ev(le, 0), ev(le, 1));
913 e[le] = l_index;
914 }
915
916 int ii = 3;
917 for (int k = 0; k < 3; ++k)
918 {
919 bool reverse = false;
920 int le = 0;
921 for (; le < ev.rows(); ++le)
922 {
923 // const auto l_index = find_edge(mesh, c, ev(le, 0), ev(le, 1));
924 // const auto l_index = mesh.get_index_from_element_edge(c, ev(le, 0), ev(le, 1));
925 const auto l_index = e[le];
926 if (l_index.edge == tmp.edge)
927 {
928 if (l_index.vertex == tmp.vertex)
929 reverse = false;
930 else
931 {
932 reverse = true;
933 if (mesh.switch_vertex(tmp).vertex != l_index.vertex)
934 assert(false);
935 }
936
937 break;
938 }
939 }
940 assert(le < 6);
941
942 if (!reverse)
943 {
944
945 for (int i = 0; i < p - 1; ++i)
946 {
947 result[ii++] = 4 + le * (p - 1) + i;
948 }
949 }
950 else
951 {
952 for (int i = 0; i < p - 1; ++i)
953 {
954 result[ii++] = 4 + (le + 1) * (p - 1) - i - 1;
955 }
956 }
957
958 tmp = mesh.next_around_face(tmp);
959 }
960
961 // faces
962
963 Eigen::Matrix<int, 4, 3> fv;
964 fv.row(0) << l2g[0], l2g[1], l2g[2];
965 fv.row(1) << l2g[0], l2g[1], l2g[3];
966 fv.row(2) << l2g[1], l2g[2], l2g[3];
967 fv.row(3) << l2g[2], l2g[0], l2g[3];
968
969 long lf = 0;
970 for (; lf < fv.rows(); ++lf)
971 {
972 const auto l_index = mesh.get_index_from_element_face(c, fv(lf, 0), fv(lf, 1), fv(lf, 2));
973 if (l_index.face == index.face)
974 break;
975 }
976
977 assert(lf < fv.rows());
978
979 if (n_face_nodes == 0)
980 {
981 }
982 else if (n_face_nodes == 1)
983 result[ii++] = 4 + n_edge_nodes + lf;
984 else // if (n_face_nodes == 3)
985 {
986
987 const auto get_order = [&p, &nn, &n_face_nodes](const std::array<int, 3> &corners) {
988 int index;
989 int start;
990 int offset;
991
992 std::vector<int> order1(n_face_nodes); // A-> B
993 for (int k = 0; k < n_face_nodes; ++k)
994 order1[k] = k;
995
996 std::vector<int> order2(n_face_nodes); // B-> A
997 index = 0;
998 start = nn - 1;
999 for (int k = 0; k < nn; ++k)
1000 {
1001 for (int l = 0; l < nn - k; ++l)
1002 {
1003 order2[index] = start - l;
1004 index++;
1005 }
1006 start += (nn - 1) - k;
1007 }
1008
1009 std::vector<int> order3(n_face_nodes); // A->C
1010 index = 0;
1011 for (int k = 0; k < nn; ++k)
1012 {
1013 offset = k;
1014 for (int l = 0; l < nn - k; ++l)
1015 {
1016 order3[index] = offset;
1017 offset += nn - l;
1018
1019 index++;
1020 }
1021 }
1022
1023 std::vector<int> order4(n_face_nodes); // C-> A
1024 index = 0;
1025 start = n_face_nodes - 1;
1026 for (int k = 0; k < nn; ++k)
1027 {
1028 offset = 0;
1029 for (int l = 0; l < nn - k; ++l)
1030 {
1031 order4[index] = start - offset;
1032 offset += k + 2 + l;
1033 index++;
1034 }
1035
1036 start += -k - 1;
1037 }
1038
1039 std::vector<int> order5(n_face_nodes); // B-> C
1040 index = 0;
1041 start = nn - 1;
1042 for (int k = 0; k < nn; ++k)
1043 {
1044 offset = 0;
1045 for (int l = 0; l < nn - k; ++l)
1046 {
1047 order5[index] = start + offset;
1048 offset += nn - 1 - l;
1049 index++;
1050 }
1051
1052 start--;
1053 }
1054
1055 std::vector<int> order6(n_face_nodes); // C-> B
1056 index = 0;
1057 start = n_face_nodes;
1058 for (int k = 0; k < nn; ++k)
1059 {
1060 offset = 0;
1061 start = start - k - 1;
1062 for (int l = 0; l < nn - k; ++l)
1063 {
1064 order6[index] = start - offset;
1065 offset += l + 1 + k;
1066 index++;
1067 }
1068 }
1069
1070 if (corners[0] == order1[0] && corners[1] == order1[nn - 1])
1071 {
1072 assert(corners[2] == order1[n_face_nodes - 1]);
1073 return order1;
1074 }
1075
1076 if (corners[0] == order2[0] && corners[1] == order2[nn - 1])
1077 {
1078 assert(corners[2] == order2[n_face_nodes - 1]);
1079 return order2;
1080 }
1081
1082 if (corners[0] == order3[0] && corners[1] == order3[nn - 1])
1083 {
1084 assert(corners[2] == order3[n_face_nodes - 1]);
1085 return order3;
1086 }
1087
1088 if (corners[0] == order4[0] && corners[1] == order4[nn - 1])
1089 {
1090 assert(corners[2] == order4[n_face_nodes - 1]);
1091 return order4;
1092 }
1093
1094 if (corners[0] == order5[0] && corners[1] == order5[nn - 1])
1095 {
1096 assert(corners[2] == order5[n_face_nodes - 1]);
1097 return order5;
1098 }
1099
1100 if (corners[0] == order6[0] && corners[1] == order6[nn - 1])
1101 {
1102 assert(corners[2] == order6[n_face_nodes - 1]);
1103 return order6;
1104 }
1105
1106 assert(false);
1107 return order1;
1108 };
1109
1110 Eigen::MatrixXd nodes;
1111 autogen::p_nodes_3d(p, nodes);
1112 // auto pos = LagrangeBasis3d::linear_tet_face_local_nodes_coordinates(mesh, index);
1113 // Local to global mapping of node indices
1114
1115 // Extract requested interface
1116 std::array<int, 3> idx;
1117 for (int lv = 0; lv < 3; ++lv)
1118 {
1119 idx[lv] = find_index(l2g.begin(), l2g.end(), index.vertex);
1120 index = mesh.next_around_face(index);
1121 }
1122 Eigen::Matrix3d pos(3, 3);
1123 int cnt = 0;
1124 for (int i : idx)
1125 {
1126 pos.row(cnt++) = nodes.row(i);
1127 }
1128
1129 const Eigen::RowVector3d bary = pos.colwise().mean();
1130
1131 const int offset = 4 + n_edge_nodes;
1132 bool found = false;
1133 for (int lff = 0; lff < 4; ++lff)
1134 {
1135 Eigen::MatrixXd loc_nodes = nodes.block(offset + lff * n_face_nodes, 0, n_face_nodes, 3);
1136 Eigen::RowVector3d node_bary = loc_nodes.colwise().mean();
1137
1138 if ((node_bary - bary).norm() < 1e-10)
1139 {
1140 std::array<int, 3> corners;
1141 int sum = 0;
1142 for (int m = 0; m < 3; ++m)
1143 {
1144 auto t = pos.row(m);
1145 int min_n = -1;
1146 double min_dis = 10000;
1147
1148 for (int n = 0; n < n_face_nodes; ++n)
1149 {
1150 double dis = (loc_nodes.row(n) - t).squaredNorm();
1151 if (dis < min_dis)
1152 {
1153 min_dis = dis;
1154 min_n = n;
1155 }
1156 }
1157
1158 assert(min_n >= 0);
1159 assert(min_n < n_face_nodes);
1160 corners[m] = min_n;
1161 }
1162
1163 const auto indices = get_order(corners);
1164 for (int min_n : indices)
1165 {
1166 sum += min_n;
1167 result[ii++] = 4 + n_edge_nodes + min_n + lf * n_face_nodes;
1168 }
1169
1170 assert(sum == (n_face_nodes - 1) * n_face_nodes / 2);
1171
1172 found = true;
1173 assert(lff == lf);
1174
1175 break;
1176 }
1177 }
1178
1179 assert(found);
1180 }
1181 // else
1182 // {
1183 // assert(n_face_nodes == 0);
1184 // }
1185
1186 assert(ii == result.size());
1187 return result;
1188}
1189
1190Eigen::VectorXi LagrangeBasis3d::hex_face_local_nodes(const bool serendipity, const int q, const Mesh3D &mesh, Navigation3D::Index index)
1191{
1192 const int nn = q - 1;
1193 const int n_edge_nodes = nn * 12;
1194 const int n_face_nodes = serendipity ? 0 : nn * nn;
1195
1196 const int c = index.element;
1197 assert(mesh.is_cube(c));
1198
1199 // Local to global mapping of node indices
1200 const auto l2g = hex_vertices_local_to_global(mesh, c);
1201
1202 // Extract requested interface
1203 Eigen::VectorXi result(4 + nn * 4 + n_face_nodes);
1204 result[0] = find_index(l2g.begin(), l2g.end(), index.vertex);
1205 result[1] = find_index(l2g.begin(), l2g.end(), mesh.next_around_face(index).vertex);
1206 result[2] = find_index(l2g.begin(), l2g.end(), mesh.next_around_face(mesh.next_around_face(index)).vertex);
1207 result[3] = find_index(l2g.begin(), l2g.end(), mesh.next_around_face(mesh.next_around_face(mesh.next_around_face(index))).vertex);
1208
1209 Eigen::Matrix<Navigation3D::Index, 12, 1> e;
1210 Eigen::Matrix<int, 12, 2> ev;
1211 ev.row(0) << l2g[0], l2g[1];
1212 ev.row(1) << l2g[1], l2g[2];
1213 ev.row(2) << l2g[2], l2g[3];
1214 ev.row(3) << l2g[3], l2g[0];
1215 ev.row(4) << l2g[0], l2g[4];
1216 ev.row(5) << l2g[1], l2g[5];
1217 ev.row(6) << l2g[2], l2g[6];
1218 ev.row(7) << l2g[3], l2g[7];
1219 ev.row(8) << l2g[4], l2g[5];
1220 ev.row(9) << l2g[5], l2g[6];
1221 ev.row(10) << l2g[6], l2g[7];
1222 ev.row(11) << l2g[7], l2g[4];
1223
1224 Navigation3D::Index tmp = index;
1225
1226 for (int le = 0; le < e.rows(); ++le)
1227 {
1228 const auto l_index = mesh.get_index_from_element_edge(c, ev(le, 0), ev(le, 1));
1229 e[le] = l_index;
1230 }
1231
1232 int ii = 4;
1233 for (int k = 0; k < 4; ++k)
1234 {
1235 bool reverse = false;
1236 int le = 0;
1237 for (; le < ev.rows(); ++le)
1238 {
1239 // const auto l_index = find_edge(mesh, c, ev(le, 0), ev(le, 1));
1240 // const auto l_index = mesh.get_index_from_element_edge(c, ev(le, 0), ev(le, 1));
1241 const auto l_index = e[le];
1242 if (l_index.edge == tmp.edge)
1243 {
1244 if (l_index.vertex == tmp.vertex)
1245 reverse = false;
1246 else
1247 {
1248 reverse = true;
1249 assert(mesh.switch_vertex(tmp).vertex == l_index.vertex);
1250 }
1251
1252 break;
1253 }
1254 }
1255 assert(le < 12);
1256
1257 if (!reverse)
1258 {
1259
1260 for (int i = 0; i < q - 1; ++i)
1261 {
1262 result[ii++] = 8 + le * (q - 1) + i;
1263 }
1264 }
1265 else
1266 {
1267 for (int i = 0; i < q - 1; ++i)
1268 {
1269 result[ii++] = 8 + (le + 1) * (q - 1) - i - 1;
1270 }
1271 }
1272
1273 tmp = mesh.next_around_face(tmp);
1274 }
1275
1276 // faces
1277
1278 Eigen::Matrix<int, 6, 4> fv;
1279 fv.row(0) << l2g[0], l2g[3], l2g[4], l2g[7];
1280 fv.row(1) << l2g[1], l2g[2], l2g[5], l2g[6];
1281 fv.row(2) << l2g[0], l2g[1], l2g[5], l2g[4];
1282 fv.row(3) << l2g[3], l2g[2], l2g[6], l2g[7];
1283 fv.row(4) << l2g[0], l2g[1], l2g[2], l2g[3];
1284 fv.row(5) << l2g[4], l2g[5], l2g[6], l2g[7];
1285
1286 long lf = 0;
1287 for (; lf < fv.rows(); ++lf)
1288 {
1289 const auto l_index = find_quad_face(mesh, c, fv(lf, 0), fv(lf, 1), fv(lf, 2), fv(lf, 3));
1290 if (l_index.face == index.face)
1291 break;
1292 }
1293
1294 assert(lf < fv.rows());
1295
1296 if (n_face_nodes == 1)
1297 result[ii++] = 8 + n_edge_nodes + lf;
1298 else if (n_face_nodes != 0)
1299 {
1300 Eigen::MatrixXd nodes;
1301 autogen::q_nodes_3d(q, nodes);
1302 // auto pos = LagrangeBasis3d::linear_tet_face_local_nodes_coordinates(mesh, index);
1303 // Local to global mapping of node indices
1304
1305 // Extract requested interface
1306 std::array<int, 4> idx;
1307 for (int lv = 0; lv < 4; ++lv)
1308 {
1309 idx[lv] = find_index(l2g.begin(), l2g.end(), index.vertex);
1310 index = mesh.next_around_face(index);
1311 }
1312 Eigen::Matrix<double, 4, 3> pos(4, 3);
1313 int cnt = 0;
1314 for (int i : idx)
1315 {
1316 pos.row(cnt++) = nodes.row(i);
1317 }
1318
1319 const Eigen::RowVector3d bary = pos.colwise().mean();
1320
1321 const int offset = 8 + n_edge_nodes;
1322 bool found = false;
1323 for (int lff = 0; lff < 6; ++lff)
1324 {
1325 Eigen::Matrix<double, 4, 3> loc_nodes = nodes.block<4, 3>(offset + lff * n_face_nodes, 0);
1326 Eigen::RowVector3d node_bary = loc_nodes.colwise().mean();
1327
1328 if ((node_bary - bary).norm() < 1e-10)
1329 {
1330 int sum = 0;
1331 for (int m = 0; m < 4; ++m)
1332 {
1333 auto t = pos.row(m);
1334 int min_n = -1;
1335 double min_dis = 10000;
1336
1337 for (int n = 0; n < 4; ++n)
1338 {
1339 double dis = (loc_nodes.row(n) - t).squaredNorm();
1340 if (dis < min_dis)
1341 {
1342 min_dis = dis;
1343 min_n = n;
1344 }
1345 }
1346
1347 assert(min_n >= 0);
1348 assert(min_n < 4);
1349
1350 sum += min_n;
1351
1352 result[ii++] = 8 + n_edge_nodes + min_n + lf * n_face_nodes;
1353 }
1354
1355 assert(sum == 6); // 0 + 1 + 2 + 3
1356
1357 found = true;
1358 assert(lff == lf);
1359 }
1360
1361 if (found)
1362 break;
1363 }
1364
1365 assert(found);
1366 }
1367
1368 assert(ii == result.size());
1369 return result;
1370}
1371
1372int LagrangeBasis3d::build_bases(
1373 const Mesh3D &mesh,
1374 const std::string &assembler,
1375 const int quadrature_order,
1376 const int mass_quadrature_order,
1377 const int discr_order,
1378 const bool bernstein,
1379 const bool serendipity,
1380 const bool has_polys,
1381 const bool is_geom_bases,
1382 const bool use_corner_quadrature,
1383 std::vector<ElementBases> &bases,
1384 std::vector<LocalBoundary> &local_boundary,
1385 std::map<int, InterfaceData> &poly_face_to_data,
1386 std::shared_ptr<MeshNodes> &mesh_nodes)
1387{
1388 Eigen::VectorXi discr_orders(mesh.n_cells());
1389 discr_orders.setConstant(discr_order);
1390
1391 return build_bases(mesh, assembler, quadrature_order, mass_quadrature_order, discr_orders, bernstein, serendipity, has_polys, is_geom_bases, use_corner_quadrature, bases, local_boundary, poly_face_to_data, mesh_nodes);
1392}
1393
1394int LagrangeBasis3d::build_bases(
1395 const Mesh3D &mesh,
1396 const std::string &assembler,
1397 const int quadrature_order,
1398 const int mass_quadrature_order,
1399 const Eigen::VectorXi &discr_orders,
1400 const bool bernstein,
1401 const bool serendipity,
1402 const bool has_polys,
1403 const bool is_geom_bases,
1404 const bool use_corner_quadrature,
1405 std::vector<ElementBases> &bases,
1406 std::vector<LocalBoundary> &local_boundary,
1407 std::map<int, InterfaceData> &poly_face_to_data,
1408 std::shared_ptr<MeshNodes> &mesh_nodes)
1409{
1410 assert(mesh.is_volume());
1411 assert(discr_orders.size() == mesh.n_cells());
1412
1413 // Navigation3D::get_index_from_element_face_time = 0;
1414 // Navigation3D::switch_vertex_time = 0;
1415 // Navigation3D::switch_edge_time = 0;
1416 // Navigation3D::switch_face_time = 0;
1417 // Navigation3D::switch_element_time = 0;
1418
1419 const int max_p = discr_orders.maxCoeff();
1420 // assert(max_p < 5); //P5 not supported
1421
1422 const int nn = max_p > 1 ? (max_p - 1) : 0;
1423 const int n_face_nodes = nn * nn;
1424 const int n_cells_nodes = nn * nn * nn;
1425
1426 Eigen::VectorXi edge_orders, face_orders;
1427 if (!mesh.is_conforming())
1428 {
1429 const auto &ncmesh = dynamic_cast<const NCMesh3D &>(mesh);
1430 compute_edge_face_orders(ncmesh, discr_orders, edge_orders, face_orders);
1431 }
1432
1433 mesh_nodes = std::make_shared<MeshNodes>(mesh, has_polys, !is_geom_bases, nn, n_face_nodes * (is_geom_bases ? 2 : 1), max_p == 0 ? 1 : n_cells_nodes);
1434 MeshNodes &nodes = *mesh_nodes;
1435 std::vector<std::vector<int>> element_nodes_id, edge_virtual_nodes, face_virtual_nodes;
1436 compute_nodes(mesh, discr_orders, edge_orders, face_orders, serendipity, has_polys, is_geom_bases, nodes, edge_virtual_nodes, face_virtual_nodes, element_nodes_id, local_boundary, poly_face_to_data);
1437 // boundary_nodes = nodes.boundary_nodes();
1438
1439 bases.resize(mesh.n_cells());
1440 std::vector<int> interface_elements;
1441 interface_elements.reserve(mesh.n_faces());
1442
1443 for (int e = 0; e < mesh.n_cells(); ++e)
1444 {
1445 ElementBases &b = bases[e];
1446 const int discr_order = discr_orders(e);
1447 const int n_el_bases = (int)element_nodes_id[e].size();
1448 b.bases.resize(n_el_bases);
1449
1450 bool skip_interface_element = false;
1451
1452 for (int j = 0; j < n_el_bases; ++j)
1453 {
1454 const int global_index = element_nodes_id[e][j];
1455 if (global_index < 0)
1456 {
1457 skip_interface_element = true;
1458 break;
1459 }
1460 }
1461
1462 if (skip_interface_element)
1463 {
1464 interface_elements.push_back(e);
1465 }
1466
1467 if (mesh.is_cube(e))
1468 {
1469 const int real_order = quadrature_order > 0 ? quadrature_order : AssemblerUtils::quadrature_order(assembler, discr_order, AssemblerUtils::BasisType::CUBE_LAGRANGE, 3);
1470 const int real_mass_order = mass_quadrature_order > 0 ? mass_quadrature_order : AssemblerUtils::quadrature_order("Mass", discr_order, AssemblerUtils::BasisType::CUBE_LAGRANGE, 3);
1471 b.set_quadrature([real_order](Quadrature &quad) {
1472 HexQuadrature hex_quadrature;
1473 hex_quadrature.get_quadrature(real_order, quad);
1474 });
1475 b.set_mass_quadrature([real_mass_order](Quadrature &quad) {
1476 HexQuadrature hex_quadrature;
1477 hex_quadrature.get_quadrature(real_mass_order, quad);
1478 });
1479
1480 b.set_local_node_from_primitive_func([serendipity, discr_order, e](const int primitive_id, const Mesh &mesh) {
1481 const auto &mesh3d = dynamic_cast<const Mesh3D &>(mesh);
1482 Navigation3D::Index index;
1483
1484 for (int lf = 0; lf < 6; ++lf)
1485 {
1486 index = mesh3d.get_index_from_element(e, lf, 0);
1487 if (index.face == primitive_id)
1488 break;
1489 }
1490 assert(index.face == primitive_id);
1491 return hex_face_local_nodes(serendipity, discr_order, mesh3d, index);
1492 });
1493
1494 for (int j = 0; j < n_el_bases; ++j)
1495 {
1496 const int global_index = element_nodes_id[e][j];
1497
1498 b.bases[j].init(discr_order, global_index, j, nodes.node_position(global_index));
1499
1500 const int dtmp = serendipity ? -2 : discr_order;
1501
1502 b.bases[j].set_basis([dtmp, j](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { autogen::q_basis_value_3d(dtmp, j, uv, val); });
1503 b.bases[j].set_grad([dtmp, j](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { autogen::q_grad_basis_value_3d(dtmp, j, uv, val); });
1504 }
1505 }
1506 else if (mesh.is_simplex(e))
1507 {
1508 const int real_order = quadrature_order > 0 ? quadrature_order : AssemblerUtils::quadrature_order(assembler, discr_order, AssemblerUtils::BasisType::SIMPLEX_LAGRANGE, 3);
1509 const int real_mass_order = mass_quadrature_order > 0 ? mass_quadrature_order : AssemblerUtils::quadrature_order("Mass", discr_order, AssemblerUtils::BasisType::SIMPLEX_LAGRANGE, 3);
1510
1511 b.set_quadrature([real_order, use_corner_quadrature](Quadrature &quad) {
1512 TetQuadrature tet_quadrature(use_corner_quadrature);
1513 tet_quadrature.get_quadrature(real_order, quad);
1514 });
1515 b.set_mass_quadrature([real_mass_order, use_corner_quadrature](Quadrature &quad) {
1516 TetQuadrature tet_quadrature(use_corner_quadrature);
1517 tet_quadrature.get_quadrature(real_mass_order, quad);
1518 });
1519
1520 b.set_local_node_from_primitive_func([discr_order, e](const int primitive_id, const Mesh &mesh) {
1521 const auto &mesh3d = dynamic_cast<const Mesh3D &>(mesh);
1522 Navigation3D::Index index;
1523
1524 for (int lf = 0; lf < mesh3d.n_cell_faces(e); ++lf)
1525 {
1526 index = mesh3d.get_index_from_element(e, lf, 0);
1527 if (index.face == primitive_id)
1528 break;
1529 }
1530 assert(index.face == primitive_id);
1531 return tet_face_local_nodes(discr_order, mesh3d, index);
1532 });
1533
1534 const bool rational = is_geom_bases && mesh.is_rational() && !mesh.cell_weights(e).empty();
1535 assert(!rational);
1536
1537 for (int j = 0; j < n_el_bases; ++j)
1538 {
1539 const int global_index = element_nodes_id[e][j];
1540 if (!skip_interface_element)
1541 {
1542 b.bases[j].init(discr_order, global_index, j, nodes.node_position(global_index));
1543 }
1544
1545 b.bases[j].set_basis([bernstein, discr_order, j](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { autogen::p_basis_value_3d(bernstein, discr_order, j, uv, val); });
1546 b.bases[j].set_grad([bernstein, discr_order, j](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { autogen::p_grad_basis_value_3d(bernstein, discr_order, j, uv, val); });
1547 }
1548 }
1549 else
1550 {
1551 // Polyhedra bases are built later on
1552 // assert(false);
1553 }
1554 }
1555
1556 if (!is_geom_bases)
1557 {
1558 if (!mesh.is_conforming())
1559 {
1560 const auto &ncmesh = dynamic_cast<const NCMesh3D &>(mesh);
1561
1562 std::vector<std::vector<int>> elementOrder;
1563 {
1564 const int max_order = discr_orders.maxCoeff(), min_order = discr_orders.minCoeff();
1565 int max_level = 0;
1566 for (int e = 0; e < ncmesh.n_cells(); e++)
1567 if (max_level < ncmesh.cell_ref_level(e))
1568 max_level = ncmesh.cell_ref_level(e);
1569
1570 elementOrder.resize((max_level + 1) * (max_order - min_order + 1));
1571 int N = 0;
1572 int cur_level = 0;
1573 while (cur_level <= max_level)
1574 {
1575 int order = min_order;
1576 while (order <= max_order)
1577 {
1578 int cur_bucket = (max_order - min_order + 1) * cur_level + (order - min_order);
1579 for (int i = 0; i < ncmesh.n_cells(); i++)
1580 {
1581 if (ncmesh.cell_ref_level(i) != cur_level || discr_orders[i] != order)
1582 continue;
1583
1584 N++;
1585 elementOrder[cur_bucket].push_back(i);
1586 }
1587 order++;
1588 }
1589 cur_level++;
1590 }
1591 }
1592
1593 for (const auto &bucket : elementOrder)
1594 {
1595 if (bucket.size() == 0)
1596 continue;
1597 polyfem::utils::maybe_parallel_for((int)bucket.size(), [&](int start, int end, int thread_id) {
1598 for (int e_aux = start; e_aux < end; e_aux++)
1599 {
1600 const int e = bucket[e_aux];
1601 ElementBases &b = bases[e];
1602 const int discr_order = discr_orders(e);
1603 const int n_edge_nodes = discr_order - 1;
1604 const int n_face_nodes = (discr_order - 1) * (discr_order - 2) / 2;
1605 const int n_el_bases = element_nodes_id[e].size();
1606
1607 auto v = tet_vertices_local_to_global(mesh, e);
1608
1609 Eigen::Matrix<Navigation3D::Index, 4, 1> cell_faces;
1610 Eigen::Matrix<int, 4, 3> fv;
1611 fv.row(0) << v[0], v[1], v[2];
1612 fv.row(1) << v[0], v[1], v[3];
1613 fv.row(2) << v[1], v[2], v[3];
1614 fv.row(3) << v[2], v[0], v[3];
1615
1616 for (long lf = 0; lf < fv.rows(); ++lf)
1617 {
1618 const auto index = mesh.get_index_from_element_face(e, fv(lf, 0), fv(lf, 1), fv(lf, 2));
1619 cell_faces[lf] = index;
1620 }
1621
1622 Eigen::Matrix<Navigation3D::Index, 6, 1> cell_edges;
1623 Eigen::Matrix<int, 6, 2> ev;
1624 ev.row(0) << v[0], v[1];
1625 ev.row(1) << v[1], v[2];
1626 ev.row(2) << v[2], v[0];
1627
1628 ev.row(3) << v[0], v[3];
1629 ev.row(4) << v[1], v[3];
1630 ev.row(5) << v[2], v[3];
1631
1632 for (int le = 0; le < ev.rows(); ++le)
1633 {
1634 // const auto index = find_edge(mesh, c, ev(le, 0), ev(le, 1));
1635 const auto index = mesh.get_index_from_element_edge(e, ev(le, 0), ev(le, 1));
1636 cell_edges[le] = index;
1637 }
1638
1639 Eigen::MatrixXd verts(4, 3);
1640 for (int i = 0; i < ncmesh.n_cell_vertices(e); i++)
1641 verts.row(i) = ncmesh.point(v[i]);
1642
1643 for (int j = 0; j < n_el_bases; ++j)
1644 {
1645 const int global_index = element_nodes_id[e][j];
1646
1647 if (global_index >= 0)
1648 {
1649 b.bases[j].init(discr_order, global_index, j, nodes.node_position(global_index));
1650 }
1651 else
1652 {
1653 // vertex node - hanging vertex
1654 if (j < 4)
1655 {
1656 int large_elem = -1;
1657 if (ncmesh.leader_edge_of_vertex(v[j]) >= 0)
1658 {
1659 large_elem = lowest_order_elem_on_edge(ncmesh, discr_orders, ncmesh.leader_edge_of_vertex(v[j]));
1660 }
1661 else if (ncmesh.leader_face_of_vertex(v[j]) >= 0)
1662 {
1663 std::vector<int> ids;
1664 ncmesh.get_face_elements_neighs(ncmesh.leader_face_of_vertex(v[j]), ids);
1665 assert(ids.size() == 1);
1666 large_elem = ids[0];
1667 }
1668 else
1669 assert(false);
1670
1671 Eigen::MatrixXd large_elem_verts(4, 3);
1672 auto v_large = tet_vertices_local_to_global(mesh, large_elem);
1673 for (int i = 0; i < ncmesh.n_cell_vertices(large_elem); i++)
1674 large_elem_verts.row(i) = ncmesh.point(v_large[i]);
1675
1676 Eigen::MatrixXd node_position;
1677 global_to_local(large_elem_verts, verts.row(j), node_position);
1678
1679 // evaluate the basis of the large element at this node
1680 const auto &other_bases = bases[large_elem];
1681 std::vector<AssemblyValues> w;
1682 other_bases.evaluate_bases(node_position, w);
1683
1684 // apply basis projection
1685 for (long i = 0; i < w.size(); ++i)
1686 {
1687 assert(w[i].val.size() == 1);
1688 if (std::abs(w[i].val(0)) < 1e-12)
1689 continue;
1690
1691 assert(other_bases.bases[i].global().size() > 0);
1692 for (size_t ii = 0; ii < other_bases.bases[i].global().size(); ++ii)
1693 {
1694 const auto &other_global = other_bases.bases[i].global()[ii];
1695 assert(other_global.index >= 0);
1696 b.bases[j].global().emplace_back(other_global.index, other_global.node, w[i].val(0) * other_global.val);
1697 }
1698 }
1699 }
1700 // edge node - slave edge / edge on face / constrained order
1701 else if (j < 4 + 6 * n_edge_nodes)
1702 {
1703 const int local_edge_id = (j - 4) / n_edge_nodes;
1704 const int edge_id = cell_edges[local_edge_id].edge;
1705 bool need_extra_fake_nodes = false;
1706 int large_elem = -1;
1707
1708 // slave edge
1709 if (ncmesh.leader_edge_of_edge(edge_id) >= 0)
1710 {
1711 std::vector<int> ids;
1712 ncmesh.get_edge_elements_neighs(ncmesh.leader_edge_of_edge(edge_id), ids);
1713 large_elem = ids[0];
1714 }
1715 // edge on face
1716 else if (ncmesh.leader_face_of_edge(edge_id) >= 0)
1717 {
1718 std::vector<int> ids;
1719 ncmesh.get_face_elements_neighs(ncmesh.leader_face_of_edge(edge_id), ids);
1720 assert(ids.size() == 1);
1721 large_elem = ids[0];
1722 }
1723 // constrained order
1724 else if (discr_order > edge_orders[edge_id])
1725 {
1726 int min_order_elem = lowest_order_elem_on_edge(ncmesh, discr_orders, edge_id);
1727 // if haven't built min_order_elem? directly contribute to extra nodes
1728 if (discr_orders[min_order_elem] < discr_order)
1729 large_elem = min_order_elem;
1730
1731 // constrained order, master edge -- need extra fake nodes
1732 if (large_elem < 0)
1733 {
1734 // assert((edge.order < 2 || edge.global_ids.size() > 0) && edge.slaves.size() > 0);
1735 need_extra_fake_nodes = true;
1736 }
1737 }
1738 else
1739 assert(false);
1740
1741 assert(large_elem >= 0 || need_extra_fake_nodes);
1742 Eigen::MatrixXd lnodes;
1743 autogen::p_nodes_3d(discr_order, lnodes);
1744 Eigen::MatrixXd local_position = lnodes.row(j);
1745 if (need_extra_fake_nodes)
1746 {
1747 Eigen::MatrixXd global_position, edge_verts(2, 3);
1748 Eigen::VectorXd point_weight;
1749
1750 edge_verts.row(0) = ncmesh.point(ncmesh.edge_vertex(edge_id, 0));
1751 edge_verts.row(1) = ncmesh.point(ncmesh.edge_vertex(edge_id, 1));
1752
1753 local_to_global(verts, local_position, global_position);
1754 global_to_local_edge(edge_verts, global_position, point_weight);
1755
1756 std::function<double(const int, const int, const double)> basis_1d = [](const int order, const int id, const double x) -> double {
1757 assert(id <= order && id >= 0);
1758 double y = 1;
1759 for (int o = 0; o <= order; o++)
1760 {
1761 if (o != id)
1762 y *= (x * order - o) / (id - o);
1763 }
1764 return y;
1765 };
1766
1767 // contribution to edge nodes
1768 for (int i = 0; i < edge_virtual_nodes[edge_id].size(); i++)
1769 {
1770 const int global_index = edge_virtual_nodes[edge_id][i];
1771 // const double weight = basis_1d(edge_orders[edge_id], i+1, edge_weight);
1772 Eigen::VectorXd node_weight;
1773 global_to_local_edge(edge_verts, nodes.node_position(global_index), node_weight);
1774 const int basis_id = std::lround(node_weight(0) * edge_orders[edge_id]);
1775 const double weight = basis_1d(edge_orders[edge_id], basis_id, point_weight(0));
1776 if (std::abs(weight) < 1e-12)
1777 continue;
1778 b.bases[j].global().emplace_back(global_index, nodes.node_position(global_index), weight);
1779 }
1780
1781 // contribution to vertex nodes
1782 for (int i = 0; i < 2; i++)
1783 {
1784 const int lv = ev(local_edge_id, i);
1785 const auto &global_ = b.bases[lv].global();
1786 Eigen::VectorXd node_weight;
1787 global_to_local_edge(edge_verts, verts.row(lv), node_weight);
1788 const int basis_id = std::lround(node_weight(0) * edge_orders[edge_id]);
1789 const double weight = basis_1d(edge_orders[edge_id], basis_id, point_weight(0));
1790 if (std::abs(weight) > 1e-12)
1791 {
1792 assert(global_.size() > 0);
1793 for (size_t ii = 0; ii < global_.size(); ++ii)
1794 b.bases[j].global().emplace_back(global_[ii].index, global_[ii].node, weight * global_[ii].val);
1795 }
1796 }
1797 }
1798 else
1799 {
1800 Eigen::MatrixXd global_position, large_elem_verts(4, 3);
1801 auto v_large = tet_vertices_local_to_global(mesh, large_elem);
1802 for (int i = 0; i < ncmesh.n_cell_vertices(large_elem); i++)
1803 large_elem_verts.row(i) = ncmesh.point(v_large[i]);
1804 local_to_global(verts, local_position, global_position);
1805 global_to_local(large_elem_verts, global_position, local_position);
1806
1807 // evaluate the basis of the large element at this node
1808 const auto &other_bases = bases[large_elem];
1809 std::vector<AssemblyValues> w;
1810 other_bases.evaluate_bases(local_position, w);
1811
1812 // apply basis projection
1813 for (long i = 0; i < w.size(); ++i)
1814 {
1815 assert(w[i].val.size() == 1);
1816 if (std::abs(w[i].val(0)) < 1e-12)
1817 continue;
1818
1819 assert(other_bases.bases[i].global().size() > 0);
1820 for (size_t ii = 0; ii < other_bases.bases[i].global().size(); ++ii)
1821 {
1822 const auto &other_global = other_bases.bases[i].global()[ii];
1823 assert(other_global.index >= 0);
1824 b.bases[j].global().emplace_back(other_global.index, other_global.node, w[i].val(0) * other_global.val);
1825 }
1826 }
1827 }
1828 }
1829 // face node - slave face / constrained order
1830 else if (j < 4 + 6 * n_edge_nodes + 4 * n_face_nodes)
1831 {
1832 const int local_face_id = (j - (4 + 6 * n_edge_nodes)) / n_face_nodes;
1833 const int face_id = cell_faces(local_face_id).face;
1834 int large_elem = -1;
1835 bool need_extra_fake_nodes = false;
1836
1837 std::vector<int> ids;
1838 ncmesh.get_face_elements_neighs(ncmesh.leader_face_of_face(face_id), ids);
1839
1840 Eigen::MatrixXd face_verts(3, 3);
1841 for (int i = 0; i < ncmesh.n_face_vertices(face_id); i++)
1842 face_verts.row(i) = ncmesh.point(ncmesh.face_vertex(face_id, i));
1843
1844 // slave face
1845 if (ncmesh.leader_face_of_face(face_id) >= 0)
1846 {
1847 assert(ids.size() == 1);
1848 large_elem = ids[0];
1849 }
1850 // constrained order, conforming face
1851
1852 else if (face_orders[face_id] < discr_order && ids.size() == 2)
1853 {
1854 large_elem = ids[0] == e ? ids[1] : ids[0];
1855 }
1856 // constrained order, master face -- need extra fake nodes
1857 else if (face_orders[face_id] < discr_order && ncmesh.n_follower_faces(face_id) > 0)
1858 {
1859 // assert(ncface.global_ids.size() > 0 || ncface.order < 3);
1860 need_extra_fake_nodes = true;
1861 }
1862 else
1863 assert(false);
1864
1865 assert(large_elem >= 0 || need_extra_fake_nodes);
1866 Eigen::MatrixXd lnodes;
1867 autogen::p_nodes_3d(discr_order, lnodes);
1868 Eigen::MatrixXd local_position = lnodes.row(j);
1869 if (need_extra_fake_nodes)
1870 {
1871 Eigen::MatrixXd global_position;
1872 local_to_global(verts, local_position, global_position);
1873
1874 Eigen::MatrixXd tmp;
1875 global_to_local_face(face_verts, global_position, tmp);
1876 Eigen::VectorXd face_weight = tmp.transpose();
1877
1878 std::function<double(const int, const int, const double)> basis_aux = [](const int order, const int id, const double x) -> double {
1879 assert(id <= order && id >= 0);
1880 double y = 1;
1881 for (int o = 0; o < id; o++)
1882 y *= (x * order - o) / (id - o);
1883 return y;
1884 };
1885
1886 std::function<double(const int, const int, const int, const Eigen::Vector2d)> basis_2d = [&basis_aux](const int order, const int i, const int j, const Eigen::Vector2d uv) -> double {
1887 assert(i + j <= order && i >= 0 && j >= 0);
1888 double u = uv(0), v = uv(1);
1889 return basis_aux(order, i, u) * basis_aux(order, j, v) * basis_aux(order, order - i - j, 1 - u - v);
1890 };
1891
1892 // contribution to face nodes
1893 for (int global_ : face_virtual_nodes[face_id])
1894 {
1895 auto low_order_node = nodes.node_position(global_);
1896 Eigen::MatrixXd low_order_node_face_weight;
1897 global_to_local_face(face_verts, low_order_node, low_order_node_face_weight);
1898 int x = round(low_order_node_face_weight(0) * face_orders[face_id]), y = round(low_order_node_face_weight(1) * face_orders[face_id]);
1899 const double weight = basis_2d(face_orders[face_id], x, y, face_weight);
1900 if (std::abs(weight) < 1e-12)
1901 continue;
1902 b.bases[j].global().emplace_back(global_, nodes.node_position(global_), weight);
1903 }
1904
1905 // contribution to vertex nodes
1906 for (int i = 0; i < 3; i++)
1907 {
1908 const auto &global_ = b.bases[fv(local_face_id, i)].global();
1909 auto low_order_node = ncmesh.point(fv(local_face_id, i));
1910 Eigen::MatrixXd low_order_node_face_weight;
1911 global_to_local_face(face_verts, low_order_node, low_order_node_face_weight);
1912 int x = round(low_order_node_face_weight(0) * face_orders[face_id]), y = round(low_order_node_face_weight(1) * face_orders[face_id]);
1913 double weight = basis_2d(face_orders[face_id], x, y, face_weight);
1914 if (std::abs(weight) > 1e-12)
1915 {
1916 assert(global_.size() > 0);
1917 for (size_t ii = 0; ii < global_.size(); ++ii)
1918 b.bases[j].global().emplace_back(global_[ii].index, global_[ii].node, weight * global_[ii].val);
1919 }
1920 }
1921
1922 // contribution to edge nodes, two steps
1923 for (int x = 0, idx = 0; x <= face_orders[face_id]; x++)
1924 {
1925 for (int y = 0; x + y <= face_orders[face_id]; y++)
1926 {
1927 const int z = face_orders[face_id] - x - y;
1928 int flag = (int)(x == 0) + (int)(y == 0) + (int)(z == 0);
1929 if (flag != 1)
1930 continue;
1931
1932 // first step
1933 const double weight = basis_2d(face_orders[face_id], x, y, face_weight);
1934 if (std::abs(weight) < 1e-12)
1935 continue;
1936 Eigen::MatrixXd face_weight(1, 2);
1937 face_weight << (double)x / face_orders[face_id], (double)y / face_orders[face_id];
1938 Eigen::MatrixXd pos, local_pos;
1939 local_to_global_face(face_verts, face_weight, pos);
1940 global_to_local(verts, pos, local_pos);
1941 Local2Global step1(idx, local_pos, weight);
1942 idx++;
1943
1944 {
1945 // evaluate the basis of the large element at this node
1946 const auto &other_bases = bases[e];
1947 std::vector<AssemblyValues> w;
1948 other_bases.evaluate_bases(local_pos, w);
1949
1950 // apply basis projection
1951 for (long i = 0; i < w.size(); ++i)
1952 {
1953 assert(w[i].val.size() == 1);
1954 if (std::abs(w[i].val(0)) < 1e-12)
1955 continue;
1956
1957 assert(other_bases.bases[i].global().size() > 0);
1958 for (size_t ii = 0; ii < other_bases.bases[i].global().size(); ++ii)
1959 {
1960 const auto &other_global = other_bases.bases[i].global()[ii];
1961 assert(other_global.index >= 0);
1962 b.bases[j].global().emplace_back(other_global.index, other_global.node, step1.val * w[i].val(0) * other_global.val);
1963 }
1964 }
1965 }
1966 }
1967 }
1968 }
1969 else
1970 {
1971 Eigen::MatrixXd global_position, large_elem_verts(4, 3);
1972 auto v_large = tet_vertices_local_to_global(mesh, large_elem);
1973 for (int i = 0; i < ncmesh.n_cell_vertices(large_elem); i++)
1974 large_elem_verts.row(i) = ncmesh.point(v_large[i]);
1975 local_to_global(verts, local_position, global_position);
1976 global_to_local(large_elem_verts, global_position, local_position);
1977
1978 // evaluate the basis of the large element at this node
1979 const auto &other_bases = bases[large_elem];
1980 std::vector<AssemblyValues> w;
1981 other_bases.evaluate_bases(local_position, w);
1982
1983 // apply basis projection
1984 for (long i = 0; i < w.size(); ++i)
1985 {
1986 assert(w[i].val.size() == 1);
1987 if (std::abs(w[i].val(0)) < 1e-12)
1988 continue;
1989
1990 assert(other_bases.bases[i].global().size() > 0);
1991 for (size_t ii = 0; ii < other_bases.bases[i].global().size(); ++ii)
1992 {
1993 const auto &other_global = other_bases.bases[i].global()[ii];
1994 assert(other_global.index >= 0);
1995 b.bases[j].global().emplace_back(other_global.index, other_global.node, w[i].val(0) * other_global.val);
1996 }
1997 }
1998 }
1999 }
2000 else
2001 assert(false);
2002
2003 auto &global_ = b.bases[j].global();
2004 if (global_.size() <= 1)
2005 continue;
2006
2007 std::map<int, Local2Global> list;
2008 for (size_t ii = 0; ii < global_.size(); ii++)
2009 {
2010 auto pair = list.insert({global_[ii].index, global_[ii]});
2011 if (!pair.second && pair.first != list.end())
2012 {
2013 assert((pair.first->second.node - global_[ii].node).norm() < 1e-12);
2014 pair.first->second.val += global_[ii].val;
2015 }
2016 }
2017
2018 global_.clear();
2019 for (auto it = list.begin(); it != list.end(); ++it)
2020 {
2021 if (std::abs(it->second.val) > 1e-12)
2022 {
2023 global_.push_back(it->second);
2024 }
2025 }
2026 }
2027 }
2028 }
2029 });
2030 }
2031 }
2032 else
2033 {
2034 for (int pp = 2; pp <= autogen::MAX_P_BASES; ++pp)
2035 {
2036 for (int e : interface_elements)
2037 {
2038 ElementBases &b = bases[e];
2039 const int discr_order = discr_orders(e);
2040 const int n_el_bases = element_nodes_id[e].size();
2041 assert(discr_order > 1);
2042 if (discr_order != pp)
2043 continue;
2044
2045 if (mesh.is_cube(e))
2046 {
2047 // TODO
2048 assert(false);
2049 }
2050 else if (mesh.is_simplex(e))
2051 {
2052 for (int j = 0; j < n_el_bases; ++j)
2053 {
2054 const int global_index = element_nodes_id[e][j];
2055
2056 if (global_index >= 0)
2057 b.bases[j].init(discr_order, global_index, j, nodes.node_position(global_index));
2058 else
2059 {
2060 const int lnn = max_p > 2 ? (discr_order - 2) : 0;
2061 const int ln_edge_nodes = discr_order - 1;
2062 const int ln_face_nodes = lnn * (lnn + 1) / 2;
2063
2064 const auto v = tet_vertices_local_to_global(mesh, e);
2065 Navigation3D::Index index;
2066 if (global_index <= -30)
2067 {
2068 assert(false);
2069 // const auto lv = -(global_index + 30);
2070 // assert(lv>=0 && lv < 4);
2071 // assert(j < 4);
2072
2073 // if(lv == 3)
2074 // {
2075 // index = mesh.switch_element(find_edge(mesh, e, v[lv], v[0]));
2076 // if(index.element < 0)
2077 // index = mesh.switch_element(find_edge(mesh, e, v[lv], v[1]));
2078 // if(index.element < 0)
2079 // index = mesh.switch_element(find_edge(mesh, e, v[lv], v[2]));
2080 // }
2081 // else
2082 // {
2083 // index = mesh.switch_element(find_edge(mesh, e, v[lv], v[(lv+1)%3]));
2084 // if(index.element < 0)
2085 // index = mesh.switch_element(find_edge(mesh, e, v[lv], v[(lv+2)%3]));
2086 // if(index.element < 0)
2087 // index = mesh.switch_element(find_edge(mesh, e, v[lv], v[3]));
2088 // }
2089 }
2090 else if (global_index <= -10)
2091 {
2092 const auto le = -(global_index + 10);
2093 assert(le >= 0 && le < 6);
2094 assert(j >= 4 && j < 4 + 6 * ln_edge_nodes);
2095
2096 Eigen::Matrix<int, 6, 2> ev;
2097 ev.row(0) << v[0], v[1];
2098 ev.row(1) << v[1], v[2];
2099 ev.row(2) << v[2], v[0];
2100
2101 ev.row(3) << v[0], v[3];
2102 ev.row(4) << v[1], v[3];
2103 ev.row(5) << v[2], v[3];
2104
2105 // const auto edge_index = find_edge(mesh, e, ev(le, 0), ev(le, 1));
2106 const auto edge_index = mesh.get_index_from_element_edge(e, ev(le, 0), ev(le, 1));
2107 auto neighs = mesh.edge_neighs(edge_index.edge);
2108 int min_p = discr_order;
2109 int min_cell = edge_index.element;
2110
2111 for (auto cid : neighs)
2112 {
2113 if (discr_orders[cid] < min_p)
2114 {
2115 min_p = discr_orders[cid];
2116 min_cell = cid;
2117 }
2118 }
2119
2120 bool found = false;
2121 for (int lf = 0; lf < 4; ++lf)
2122 {
2123 for (int lv = 0; lv < 4; ++lv)
2124 {
2125 index = mesh.get_index_from_element(min_cell, lf, lv);
2126
2127 if (index.vertex == edge_index.vertex)
2128 {
2129 if (index.edge != edge_index.edge)
2130 {
2131 auto tmp = index;
2132 index = mesh.switch_edge(tmp);
2133
2134 if (index.edge != edge_index.edge)
2135 {
2136 index = mesh.switch_edge(mesh.switch_face(tmp));
2137 }
2138 }
2139 found = true;
2140 break;
2141 }
2142 }
2143
2144 if (found)
2145 break;
2146 }
2147
2148 assert(found);
2149 assert(index.vertex == edge_index.vertex && index.edge == edge_index.edge);
2150 assert(index.element != edge_index.element);
2151 }
2152 else
2153 {
2154 const auto lf = -(global_index + 1);
2155 assert(lf >= 0 && lf < 4);
2156 assert(j >= 4 + 6 * ln_edge_nodes && j < 4 + 6 * ln_edge_nodes + 4 * ln_face_nodes);
2157
2158 Eigen::Matrix<int, 4, 3> fv;
2159 fv.row(0) << v[0], v[1], v[2];
2160 fv.row(1) << v[0], v[1], v[3];
2161 fv.row(2) << v[1], v[2], v[3];
2162 fv.row(3) << v[2], v[0], v[3];
2163
2164 index = mesh.switch_element(mesh.get_index_from_element_face(e, fv(lf, 0), fv(lf, 1), fv(lf, 2)));
2165 }
2166
2167 const auto other_cell = index.element;
2168 assert(other_cell >= 0);
2169 assert(discr_order > discr_orders(other_cell));
2170
2171 auto indices = tet_face_local_nodes(discr_order, mesh, index);
2172 Eigen::MatrixXd lnodes;
2173 autogen::p_nodes_3d(discr_order, lnodes);
2174 Eigen::RowVector3d node_position; // = lnodes.row(indices(ii));
2175
2176 if (j < 4)
2177 node_position = lnodes.row(indices(0));
2178 else if (j < 4 + 6 * ln_edge_nodes)
2179 node_position = lnodes.row(indices(((j - 4) % ln_edge_nodes) + 3));
2180 else if (j < 4 + 6 * ln_edge_nodes + 4 * ln_face_nodes)
2181 {
2182 // node_position = lnodes.row(indices(((j - 4 - 6*ln_edge_nodes) % ln_face_nodes) + 3 + 3*ln_edge_nodes));
2183 auto me_indices = tet_face_local_nodes(discr_order, mesh, mesh.switch_element(index));
2184 int ii;
2185 for (ii = 0; ii < me_indices.size(); ++ii)
2186 {
2187 if (me_indices(ii) == j)
2188 break;
2189 }
2190
2191 assert(ii >= 3 + 3 * ln_edge_nodes);
2192 assert(ii < me_indices.size());
2193
2194 node_position = lnodes.row(indices(ii));
2195 }
2196 else
2197 assert(false);
2198
2199 const auto &other_bases = bases[other_cell];
2200 // Eigen::MatrixXd w;
2201 std::vector<AssemblyValues> w;
2202 other_bases.evaluate_bases(node_position, w);
2203
2204 assert(b.bases[j].global().size() == 0);
2205
2206 for (long i = 0; i < w.size(); ++i)
2207 {
2208 assert(w[i].val.size() == 1);
2209 if (std::abs(w[i].val(0)) < 1e-8)
2210 continue;
2211
2212 // assert(other_bases.bases[i].global().size() == 1);
2213 for (size_t ii = 0; ii < other_bases.bases[i].global().size(); ++ii)
2214 {
2215 const auto &other_global = other_bases.bases[i].global()[ii];
2216 // logger().trace("e {} j {} gid {}", e, j, other_global.index);
2217 b.bases[j].global().emplace_back(other_global.index, other_global.node, w[i].val(0) * other_global.val);
2218 }
2219 }
2220 }
2221 }
2222 }
2223 else
2224 {
2225 // Polygon bases are built later on
2226 }
2227 }
2228 }
2229 }
2230 }
2231
2232 return nodes.n_nodes();
2233}
vec
Eigen::MatrixXd vec
Definition Assembler.cpp:72
val
double val
Definition Assembler.cpp:86
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:189
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::LagrangeBasis3d::build_bases
static int build_bases(const mesh::Mesh3D &mesh, const std::string &assembler, const int quadrature_order, const int mass_quadrature_order, const int discr_order, const bool bernstein, 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_face_to_data, std::shared_ptr< mesh::MeshNodes > &mesh_nodes)
Builds FE basis functions over the entire mesh (P1, P2 over tets, Q1, Q2 over hes).
polyfem::basis::LagrangeBasis3d::hex_face_local_nodes
static Eigen::VectorXi hex_face_local_nodes(const bool serendipity, const int q, const mesh::Mesh3D &mesh, mesh::Navigation3D::Index index)
Definition LagrangeBasis3d.cpp:1190
polyfem::basis::LagrangeBasis3d::tet_face_local_nodes
static Eigen::VectorXi tet_face_local_nodes(const int p, const mesh::Mesh3D &mesh, mesh::Navigation3D::Index index)
Definition LagrangeBasis3d.cpp:879
polyfem::mesh::LocalBoundary
Boundary primitive IDs for a single element.
Definition LocalBoundary.hpp:25
polyfem::mesh::Mesh3D::get_index_from_element
virtual Navigation3D::Index get_index_from_element(int hi, int lf, int lv) const =0
polyfem::mesh::Mesh3D::switch_element
virtual Navigation3D::Index switch_element(Navigation3D::Index idx) const =0
polyfem::mesh::Mesh3D::get_index_from_element_edge
virtual Navigation3D::Index get_index_from_element_edge(int hi, int v0, int v1) const =0
polyfem::mesh::Mesh3D::n_cell_faces
virtual int n_cell_faces(const int c_id) const =0
polyfem::mesh::Mesh3D::edge_neighs
virtual std::vector< uint32_t > edge_neighs(const int e_gid) const =0
polyfem::mesh::Mesh3D::is_volume
bool is_volume() const override
checks if mesh is volume
Definition Mesh3D.hpp:28
polyfem::mesh::Mesh3D::next_around_face
virtual Navigation3D::Index next_around_face(Navigation3D::Index idx) const =0
polyfem::mesh::Mesh3D::get_index_from_element_face
virtual Navigation3D::Index get_index_from_element_face(int hi, int v0, int v1, int v2) const =0
polyfem::mesh::Mesh3D::switch_vertex
virtual Navigation3D::Index switch_vertex(Navigation3D::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::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::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: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::n_cells
virtual int n_cells() const =0
number of cells
polyfem::mesh::Mesh::n_faces
virtual int n_faces() const =0
number of faces
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::NCMesh3D::cell_ref_level
int cell_ref_level(const int c_id) const
Definition NCMesh3D.hpp:227
polyfem::mesh::NCMesh3D::n_faces
int n_faces() const override
number of faces
Definition NCMesh3D.hpp:189
polyfem::mesh::NCMesh3D::face_edge
int face_edge(const int f_id, const int le_id) const
Definition NCMesh3D.cpp:18
polyfem::mesh::NCMesh3D::n_cell_faces
int n_cell_faces(const int c_id) const override
Definition NCMesh3D.hpp:219
polyfem::mesh::NCMesh3D::n_edge_cells
int n_edge_cells(const int e_id) const
Definition NCMesh3D.hpp:216
polyfem::mesh::NCMesh3D::n_cells
int n_cells() const override
number of cells
Definition NCMesh3D.hpp:188
polyfem::mesh::NCMesh3D::n_edges
int n_edges() const override
number of edges
Definition NCMesh3D.hpp:197
polyfem::mesh::NCMesh3D::n_face_vertices
int n_face_vertices(const int f_id) const override
number of vertices of a face
Definition NCMesh3D.hpp:214
polyfem::mesh::NCMesh3D::edge_neighs
std::vector< uint32_t > edge_neighs(const int e_gid) const override
Definition NCMesh3D.cpp:446
polyfem::mesh::NCMesh3D::leader_face_of_edge
int leader_face_of_edge(const int e_id) const
Definition NCMesh3D.hpp:298
polyfem::mesh::NCMesh3D::cell_face
int cell_face(const int c_id, const int lf_id) const override
Definition NCMesh3D.hpp:221
polyfem::mesh::NCMesh3D::leader_edge_of_edge
int leader_edge_of_edge(const int e_id) const
Definition NCMesh3D.hpp:286
polyfem::mesh::NCMesh3D::n_cell_edges
int n_cell_edges(const int c_id) const override
Definition NCMesh3D.hpp:218
polyfem::mesh::NCMesh3D::cell_edge
int cell_edge(const int c_id, const int le_id) const override
Definition NCMesh3D.hpp:222
polyfem::mesh::NCMesh3D::leader_face_of_face
int leader_face_of_face(const int f_id) const
Definition NCMesh3D.hpp:304
polyfem::mesh::NCMesh3D::n_face_cells
int n_face_cells(const int f_id) const
Definition NCMesh3D.hpp:215
polyfem::quadrature::HexQuadrature::get_quadrature
void get_quadrature(const int order, Quadrature &quad)
Definition HexQuadrature.cpp:17
elasticity_rhs.f
f
Definition elasticity_rhs.py:173
generate_mooney_rivlin.c2
c2
Definition generate_mooney_rivlin.py:5
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:365
p_bases.nn
nn
Definition p_bases.py:341
p_bases.nodes
str nodes
Definition p_bases.py:398
p_bases.indices
list indices
Definition p_bases.py:258
p_bases.N
N
Definition p_bases.py:248
polyfem::assembler
Used for test only.
Definition AMIPSEnergy.cpp:6
polyfem::autogen::q_grad_basis_value_3d
void q_grad_basis_value_3d(const int q, const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val)
Definition auto_q_bases_3d_grad.cpp:12
polyfem::autogen::p_grad_basis_value_3d
void p_grad_basis_value_3d(const bool bernstein, const int p, const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val)
Definition auto_p_bases.cpp:852
polyfem::autogen::p_nodes_3d
void p_nodes_3d(const int p, Eigen::MatrixXd &val)
Definition auto_p_bases.cpp:830
polyfem::autogen::p_basis_value_3d
void p_basis_value_3d(const bool bernstein, const int p, const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val)
Definition auto_p_bases.cpp:839
polyfem::autogen::q_nodes_3d
void q_nodes_3d(const int q, Eigen::MatrixXd &val)
Definition auto_q_bases_3d_nodes.cpp:149
polyfem::autogen::q_basis_value_3d
void q_basis_value_3d(const int q, const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val)
Definition auto_q_bases_3d_val.cpp:175
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