PolyFEM
Loading...
Searching...
No Matches
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 serendipity,
1379 const bool has_polys,
1380 const bool is_geom_bases,
1381 std::vector<ElementBases> &bases,
1382 std::vector<LocalBoundary> &local_boundary,
1383 std::map<int, InterfaceData> &poly_face_to_data,
1384 std::shared_ptr<MeshNodes> &mesh_nodes)
1385{
1386 Eigen::VectorXi discr_orders(mesh.n_cells());
1387 discr_orders.setConstant(discr_order);
1388
1389 return build_bases(mesh, assembler, quadrature_order, mass_quadrature_order, discr_orders, serendipity, has_polys, is_geom_bases, bases, local_boundary, poly_face_to_data, mesh_nodes);
1390}
1391
1392int LagrangeBasis3d::build_bases(
1393 const Mesh3D &mesh,
1394 const std::string &assembler,
1395 const int quadrature_order,
1396 const int mass_quadrature_order,
1397 const Eigen::VectorXi &discr_orders,
1398 const bool serendipity,
1399 const bool has_polys,
1400 const bool is_geom_bases,
1401 std::vector<ElementBases> &bases,
1402 std::vector<LocalBoundary> &local_boundary,
1403 std::map<int, InterfaceData> &poly_face_to_data,
1404 std::shared_ptr<MeshNodes> &mesh_nodes)
1405{
1406 assert(mesh.is_volume());
1407 assert(discr_orders.size() == mesh.n_cells());
1408
1409 // Navigation3D::get_index_from_element_face_time = 0;
1410 // Navigation3D::switch_vertex_time = 0;
1411 // Navigation3D::switch_edge_time = 0;
1412 // Navigation3D::switch_face_time = 0;
1413 // Navigation3D::switch_element_time = 0;
1414
1415 const int max_p = discr_orders.maxCoeff();
1416 // assert(max_p < 5); //P5 not supported
1417
1418 const int nn = max_p > 1 ? (max_p - 1) : 0;
1419 const int n_face_nodes = nn * nn;
1420 const int n_cells_nodes = nn * nn * nn;
1421
1422 Eigen::VectorXi edge_orders, face_orders;
1423 if (!mesh.is_conforming())
1424 {
1425 const auto &ncmesh = dynamic_cast<const NCMesh3D &>(mesh);
1426 compute_edge_face_orders(ncmesh, discr_orders, edge_orders, face_orders);
1427 }
1428
1429 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);
1430 MeshNodes &nodes = *mesh_nodes;
1431 std::vector<std::vector<int>> element_nodes_id, edge_virtual_nodes, face_virtual_nodes;
1432 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);
1433 // boundary_nodes = nodes.boundary_nodes();
1434
1435 // std::cout<<"get_index_from_element_face_time " << Navigation3D::get_index_from_element_face_time <<std::endl;
1436 // std::cout<<"switch_vertex_time " << Navigation3D::switch_vertex_time <<std::endl;
1437 // std::cout<<"switch_edge_time " << Navigation3D::switch_edge_time <<std::endl;
1438 // std::cout<<"switch_face_time " << Navigation3D::switch_face_time <<std::endl;
1439 // std::cout<<"switch_element_time " << Navigation3D::switch_element_time <<std::endl;
1440
1441 bases.resize(mesh.n_cells());
1442 std::vector<int> interface_elements;
1443 interface_elements.reserve(mesh.n_faces());
1444
1445 for (int e = 0; e < mesh.n_cells(); ++e)
1446 {
1447 ElementBases &b = bases[e];
1448 const int discr_order = discr_orders(e);
1449 const int n_el_bases = (int)element_nodes_id[e].size();
1450 b.bases.resize(n_el_bases);
1451
1452 bool skip_interface_element = false;
1453
1454 for (int j = 0; j < n_el_bases; ++j)
1455 {
1456 const int global_index = element_nodes_id[e][j];
1457 if (global_index < 0)
1458 {
1459 skip_interface_element = true;
1460 break;
1461 }
1462 }
1463
1464 if (skip_interface_element)
1465 {
1466 interface_elements.push_back(e);
1467 }
1468
1469 if (mesh.is_cube(e))
1470 {
1471 const int real_order = quadrature_order > 0 ? quadrature_order : AssemblerUtils::quadrature_order(assembler, discr_order, AssemblerUtils::BasisType::CUBE_LAGRANGE, 3);
1472 const int real_mass_order = mass_quadrature_order > 0 ? mass_quadrature_order : AssemblerUtils::quadrature_order("Mass", discr_order, AssemblerUtils::BasisType::CUBE_LAGRANGE, 3);
1473 b.set_quadrature([real_order](Quadrature &quad) {
1474 HexQuadrature hex_quadrature;
1475 hex_quadrature.get_quadrature(real_order, quad);
1476 });
1477 b.set_mass_quadrature([real_mass_order](Quadrature &quad) {
1478 HexQuadrature hex_quadrature;
1479 hex_quadrature.get_quadrature(real_mass_order, quad);
1480 });
1481
1482 b.set_local_node_from_primitive_func([serendipity, discr_order, e](const int primitive_id, const Mesh &mesh) {
1483 const auto &mesh3d = dynamic_cast<const Mesh3D &>(mesh);
1484 Navigation3D::Index index;
1485
1486 for (int lf = 0; lf < 6; ++lf)
1487 {
1488 index = mesh3d.get_index_from_element(e, lf, 0);
1489 if (index.face == primitive_id)
1490 break;
1491 }
1492 assert(index.face == primitive_id);
1493 return hex_face_local_nodes(serendipity, discr_order, mesh3d, index);
1494 });
1495
1496 for (int j = 0; j < n_el_bases; ++j)
1497 {
1498 const int global_index = element_nodes_id[e][j];
1499
1500 b.bases[j].init(discr_order, global_index, j, nodes.node_position(global_index));
1501
1502 const int dtmp = serendipity ? -2 : discr_order;
1503
1504 b.bases[j].set_basis([dtmp, j](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { autogen::q_basis_value_3d(dtmp, j, uv, val); });
1505 b.bases[j].set_grad([dtmp, j](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { autogen::q_grad_basis_value_3d(dtmp, j, uv, val); });
1506 }
1507 }
1508 else if (mesh.is_simplex(e))
1509 {
1510 const int real_order = quadrature_order > 0 ? quadrature_order : AssemblerUtils::quadrature_order(assembler, discr_order, AssemblerUtils::BasisType::SIMPLEX_LAGRANGE, 3);
1511 const int real_mass_order = mass_quadrature_order > 0 ? mass_quadrature_order : AssemblerUtils::quadrature_order("Mass", discr_order, AssemblerUtils::BasisType::SIMPLEX_LAGRANGE, 3);
1512
1513 b.set_quadrature([real_order](Quadrature &quad) {
1514 TetQuadrature tet_quadrature;
1515 tet_quadrature.get_quadrature(real_order, quad);
1516 });
1517 b.set_mass_quadrature([real_mass_order](Quadrature &quad) {
1518 TetQuadrature tet_quadrature;
1519 tet_quadrature.get_quadrature(real_mass_order, quad);
1520 });
1521
1522 b.set_local_node_from_primitive_func([discr_order, e](const int primitive_id, const Mesh &mesh) {
1523 const auto &mesh3d = dynamic_cast<const Mesh3D &>(mesh);
1524 Navigation3D::Index index;
1525
1526 for (int lf = 0; lf < mesh3d.n_cell_faces(e); ++lf)
1527 {
1528 index = mesh3d.get_index_from_element(e, lf, 0);
1529 if (index.face == primitive_id)
1530 break;
1531 }
1532 assert(index.face == primitive_id);
1533 return tet_face_local_nodes(discr_order, mesh3d, index);
1534 });
1535
1536 const bool rational = is_geom_bases && mesh.is_rational() && !mesh.cell_weights(e).empty();
1537 assert(!rational);
1538
1539 for (int j = 0; j < n_el_bases; ++j)
1540 {
1541 const int global_index = element_nodes_id[e][j];
1542 if (!skip_interface_element)
1543 {
1544 b.bases[j].init(discr_order, global_index, j, nodes.node_position(global_index));
1545 }
1546
1547 b.bases[j].set_basis([discr_order, j](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { autogen::p_basis_value_3d(discr_order, j, uv, val); });
1548 b.bases[j].set_grad([discr_order, j](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { autogen::p_grad_basis_value_3d(discr_order, j, uv, val); });
1549 }
1550 }
1551 else
1552 {
1553 // Polyhedra bases are built later on
1554 // assert(false);
1555 }
1556 }
1557
1558 if (!is_geom_bases)
1559 {
1560 if (!mesh.is_conforming())
1561 {
1562 const auto &ncmesh = dynamic_cast<const NCMesh3D &>(mesh);
1563
1564 std::vector<std::vector<int>> elementOrder;
1565 {
1566 const int max_order = discr_orders.maxCoeff(), min_order = discr_orders.minCoeff();
1567 int max_level = 0;
1568 for (int e = 0; e < ncmesh.n_cells(); e++)
1569 if (max_level < ncmesh.cell_ref_level(e))
1570 max_level = ncmesh.cell_ref_level(e);
1571
1572 elementOrder.resize((max_level + 1) * (max_order - min_order + 1));
1573 int N = 0;
1574 int cur_level = 0;
1575 while (cur_level <= max_level)
1576 {
1577 int order = min_order;
1578 while (order <= max_order)
1579 {
1580 int cur_bucket = (max_order - min_order + 1) * cur_level + (order - min_order);
1581 for (int i = 0; i < ncmesh.n_cells(); i++)
1582 {
1583 if (ncmesh.cell_ref_level(i) != cur_level || discr_orders[i] != order)
1584 continue;
1585
1586 N++;
1587 elementOrder[cur_bucket].push_back(i);
1588 }
1589 order++;
1590 }
1591 cur_level++;
1592 }
1593 }
1594
1595 for (const auto &bucket : elementOrder)
1596 {
1597 if (bucket.size() == 0)
1598 continue;
1599 polyfem::utils::maybe_parallel_for((int)bucket.size(), [&](int start, int end, int thread_id) {
1600 for (int e_aux = start; e_aux < end; e_aux++)
1601 {
1602 const int e = bucket[e_aux];
1603 ElementBases &b = bases[e];
1604 const int discr_order = discr_orders(e);
1605 const int n_edge_nodes = discr_order - 1;
1606 const int n_face_nodes = (discr_order - 1) * (discr_order - 2) / 2;
1607 const int n_el_bases = element_nodes_id[e].size();
1608
1609 auto v = tet_vertices_local_to_global(mesh, e);
1610
1611 Eigen::Matrix<Navigation3D::Index, 4, 1> cell_faces;
1612 Eigen::Matrix<int, 4, 3> fv;
1613 fv.row(0) << v[0], v[1], v[2];
1614 fv.row(1) << v[0], v[1], v[3];
1615 fv.row(2) << v[1], v[2], v[3];
1616 fv.row(3) << v[2], v[0], v[3];
1617
1618 for (long lf = 0; lf < fv.rows(); ++lf)
1619 {
1620 const auto index = mesh.get_index_from_element_face(e, fv(lf, 0), fv(lf, 1), fv(lf, 2));
1621 cell_faces[lf] = index;
1622 }
1623
1624 Eigen::Matrix<Navigation3D::Index, 6, 1> cell_edges;
1625 Eigen::Matrix<int, 6, 2> ev;
1626 ev.row(0) << v[0], v[1];
1627 ev.row(1) << v[1], v[2];
1628 ev.row(2) << v[2], v[0];
1629
1630 ev.row(3) << v[0], v[3];
1631 ev.row(4) << v[1], v[3];
1632 ev.row(5) << v[2], v[3];
1633
1634 for (int le = 0; le < ev.rows(); ++le)
1635 {
1636 // const auto index = find_edge(mesh, c, ev(le, 0), ev(le, 1));
1637 const auto index = mesh.get_index_from_element_edge(e, ev(le, 0), ev(le, 1));
1638 cell_edges[le] = index;
1639 }
1640
1641 Eigen::MatrixXd verts(4, 3);
1642 for (int i = 0; i < ncmesh.n_cell_vertices(e); i++)
1643 verts.row(i) = ncmesh.point(v[i]);
1644
1645 for (int j = 0; j < n_el_bases; ++j)
1646 {
1647 const int global_index = element_nodes_id[e][j];
1648
1649 if (global_index >= 0)
1650 {
1651 b.bases[j].init(discr_order, global_index, j, nodes.node_position(global_index));
1652 }
1653 else
1654 {
1655 // vertex node - hanging vertex
1656 if (j < 4)
1657 {
1658 int large_elem = -1;
1659 if (ncmesh.leader_edge_of_vertex(v[j]) >= 0)
1660 {
1661 large_elem = lowest_order_elem_on_edge(ncmesh, discr_orders, ncmesh.leader_edge_of_vertex(v[j]));
1662 }
1663 else if (ncmesh.leader_face_of_vertex(v[j]) >= 0)
1664 {
1665 std::vector<int> ids;
1666 ncmesh.get_face_elements_neighs(ncmesh.leader_face_of_vertex(v[j]), ids);
1667 assert(ids.size() == 1);
1668 large_elem = ids[0];
1669 }
1670 else
1671 assert(false);
1672
1673 Eigen::MatrixXd large_elem_verts(4, 3);
1674 auto v_large = tet_vertices_local_to_global(mesh, large_elem);
1675 for (int i = 0; i < ncmesh.n_cell_vertices(large_elem); i++)
1676 large_elem_verts.row(i) = ncmesh.point(v_large[i]);
1677
1678 Eigen::MatrixXd node_position;
1679 global_to_local(large_elem_verts, verts.row(j), node_position);
1680
1681 // evaluate the basis of the large element at this node
1682 const auto &other_bases = bases[large_elem];
1683 std::vector<AssemblyValues> w;
1684 other_bases.evaluate_bases(node_position, w);
1685
1686 // apply basis projection
1687 for (long i = 0; i < w.size(); ++i)
1688 {
1689 assert(w[i].val.size() == 1);
1690 if (std::abs(w[i].val(0)) < 1e-12)
1691 continue;
1692
1693 assert(other_bases.bases[i].global().size() > 0);
1694 for (size_t ii = 0; ii < other_bases.bases[i].global().size(); ++ii)
1695 {
1696 const auto &other_global = other_bases.bases[i].global()[ii];
1697 assert(other_global.index >= 0);
1698 b.bases[j].global().emplace_back(other_global.index, other_global.node, w[i].val(0) * other_global.val);
1699 }
1700 }
1701 }
1702 // edge node - slave edge / edge on face / constrained order
1703 else if (j < 4 + 6 * n_edge_nodes)
1704 {
1705 const int local_edge_id = (j - 4) / n_edge_nodes;
1706 const int edge_id = cell_edges[local_edge_id].edge;
1707 bool need_extra_fake_nodes = false;
1708 int large_elem = -1;
1709
1710 // slave edge
1711 if (ncmesh.leader_edge_of_edge(edge_id) >= 0)
1712 {
1713 std::vector<int> ids;
1714 ncmesh.get_edge_elements_neighs(ncmesh.leader_edge_of_edge(edge_id), ids);
1715 large_elem = ids[0];
1716 }
1717 // edge on face
1718 else if (ncmesh.leader_face_of_edge(edge_id) >= 0)
1719 {
1720 std::vector<int> ids;
1721 ncmesh.get_face_elements_neighs(ncmesh.leader_face_of_edge(edge_id), ids);
1722 assert(ids.size() == 1);
1723 large_elem = ids[0];
1724 }
1725 // constrained order
1726 else if (discr_order > edge_orders[edge_id])
1727 {
1728 int min_order_elem = lowest_order_elem_on_edge(ncmesh, discr_orders, edge_id);
1729 // if haven't built min_order_elem? directly contribute to extra nodes
1730 if (discr_orders[min_order_elem] < discr_order)
1731 large_elem = min_order_elem;
1732
1733 // constrained order, master edge -- need extra fake nodes
1734 if (large_elem < 0)
1735 {
1736 // assert((edge.order < 2 || edge.global_ids.size() > 0) && edge.slaves.size() > 0);
1737 need_extra_fake_nodes = true;
1738 }
1739 }
1740 else
1741 assert(false);
1742
1743 assert(large_elem >= 0 || need_extra_fake_nodes);
1744 Eigen::MatrixXd lnodes;
1745 autogen::p_nodes_3d(discr_order, lnodes);
1746 Eigen::MatrixXd local_position = lnodes.row(j);
1747 if (need_extra_fake_nodes)
1748 {
1749 Eigen::MatrixXd global_position, edge_verts(2, 3);
1750 Eigen::VectorXd point_weight;
1751
1752 edge_verts.row(0) = ncmesh.point(ncmesh.edge_vertex(edge_id, 0));
1753 edge_verts.row(1) = ncmesh.point(ncmesh.edge_vertex(edge_id, 1));
1754
1755 local_to_global(verts, local_position, global_position);
1756 global_to_local_edge(edge_verts, global_position, point_weight);
1757
1758 std::function<double(const int, const int, const double)> basis_1d = [](const int order, const int id, const double x) -> double {
1759 assert(id <= order && id >= 0);
1760 double y = 1;
1761 for (int o = 0; o <= order; o++)
1762 {
1763 if (o != id)
1764 y *= (x * order - o) / (id - o);
1765 }
1766 return y;
1767 };
1768
1769 // contribution to edge nodes
1770 for (int i = 0; i < edge_virtual_nodes[edge_id].size(); i++)
1771 {
1772 const int global_index = edge_virtual_nodes[edge_id][i];
1773 // const double weight = basis_1d(edge_orders[edge_id], i+1, edge_weight);
1774 Eigen::VectorXd node_weight;
1775 global_to_local_edge(edge_verts, nodes.node_position(global_index), node_weight);
1776 const int basis_id = std::lround(node_weight(0) * edge_orders[edge_id]);
1777 const double weight = basis_1d(edge_orders[edge_id], basis_id, point_weight(0));
1778 if (std::abs(weight) < 1e-12)
1779 continue;
1780 b.bases[j].global().emplace_back(global_index, nodes.node_position(global_index), weight);
1781 }
1782
1783 // contribution to vertex nodes
1784 for (int i = 0; i < 2; i++)
1785 {
1786 const int lv = ev(local_edge_id, i);
1787 const auto &global_ = b.bases[lv].global();
1788 Eigen::VectorXd node_weight;
1789 global_to_local_edge(edge_verts, verts.row(lv), node_weight);
1790 const int basis_id = std::lround(node_weight(0) * edge_orders[edge_id]);
1791 const double weight = basis_1d(edge_orders[edge_id], basis_id, point_weight(0));
1792 if (std::abs(weight) > 1e-12)
1793 {
1794 assert(global_.size() > 0);
1795 for (size_t ii = 0; ii < global_.size(); ++ii)
1796 b.bases[j].global().emplace_back(global_[ii].index, global_[ii].node, weight * global_[ii].val);
1797 }
1798 }
1799 }
1800 else
1801 {
1802 Eigen::MatrixXd global_position, large_elem_verts(4, 3);
1803 auto v_large = tet_vertices_local_to_global(mesh, large_elem);
1804 for (int i = 0; i < ncmesh.n_cell_vertices(large_elem); i++)
1805 large_elem_verts.row(i) = ncmesh.point(v_large[i]);
1806 local_to_global(verts, local_position, global_position);
1807 global_to_local(large_elem_verts, global_position, local_position);
1808
1809 // evaluate the basis of the large element at this node
1810 const auto &other_bases = bases[large_elem];
1811 std::vector<AssemblyValues> w;
1812 other_bases.evaluate_bases(local_position, w);
1813
1814 // apply basis projection
1815 for (long i = 0; i < w.size(); ++i)
1816 {
1817 assert(w[i].val.size() == 1);
1818 if (std::abs(w[i].val(0)) < 1e-12)
1819 continue;
1820
1821 assert(other_bases.bases[i].global().size() > 0);
1822 for (size_t ii = 0; ii < other_bases.bases[i].global().size(); ++ii)
1823 {
1824 const auto &other_global = other_bases.bases[i].global()[ii];
1825 assert(other_global.index >= 0);
1826 b.bases[j].global().emplace_back(other_global.index, other_global.node, w[i].val(0) * other_global.val);
1827 }
1828 }
1829 }
1830 }
1831 // face node - slave face / constrained order
1832 else if (j < 4 + 6 * n_edge_nodes + 4 * n_face_nodes)
1833 {
1834 const int local_face_id = (j - (4 + 6 * n_edge_nodes)) / n_face_nodes;
1835 const int face_id = cell_faces(local_face_id).face;
1836 int large_elem = -1;
1837 bool need_extra_fake_nodes = false;
1838
1839 std::vector<int> ids;
1840 ncmesh.get_face_elements_neighs(ncmesh.leader_face_of_face(face_id), ids);
1841
1842 Eigen::MatrixXd face_verts(3, 3);
1843 for (int i = 0; i < ncmesh.n_face_vertices(face_id); i++)
1844 face_verts.row(i) = ncmesh.point(ncmesh.face_vertex(face_id, i));
1845
1846 // slave face
1847 if (ncmesh.leader_face_of_face(face_id) >= 0)
1848 {
1849 assert(ids.size() == 1);
1850 large_elem = ids[0];
1851 }
1852 // constrained order, conforming face
1853
1854 else if (face_orders[face_id] < discr_order && ids.size() == 2)
1855 {
1856 large_elem = ids[0] == e ? ids[1] : ids[0];
1857 }
1858 // constrained order, master face -- need extra fake nodes
1859 else if (face_orders[face_id] < discr_order && ncmesh.n_follower_faces(face_id) > 0)
1860 {
1861 // assert(ncface.global_ids.size() > 0 || ncface.order < 3);
1862 need_extra_fake_nodes = true;
1863 }
1864 else
1865 assert(false);
1866
1867 assert(large_elem >= 0 || need_extra_fake_nodes);
1868 Eigen::MatrixXd lnodes;
1869 autogen::p_nodes_3d(discr_order, lnodes);
1870 Eigen::MatrixXd local_position = lnodes.row(j);
1871 if (need_extra_fake_nodes)
1872 {
1873 Eigen::MatrixXd global_position;
1874 local_to_global(verts, local_position, global_position);
1875
1876 Eigen::MatrixXd tmp;
1877 global_to_local_face(face_verts, global_position, tmp);
1878 Eigen::VectorXd face_weight = tmp.transpose();
1879
1880 std::function<double(const int, const int, const double)> basis_aux = [](const int order, const int id, const double x) -> double {
1881 assert(id <= order && id >= 0);
1882 double y = 1;
1883 for (int o = 0; o < id; o++)
1884 y *= (x * order - o) / (id - o);
1885 return y;
1886 };
1887
1888 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 {
1889 assert(i + j <= order && i >= 0 && j >= 0);
1890 double u = uv(0), v = uv(1);
1891 return basis_aux(order, i, u) * basis_aux(order, j, v) * basis_aux(order, order - i - j, 1 - u - v);
1892 };
1893
1894 // contribution to face nodes
1895 for (int global_ : face_virtual_nodes[face_id])
1896 {
1897 auto low_order_node = nodes.node_position(global_);
1898 Eigen::MatrixXd low_order_node_face_weight;
1899 global_to_local_face(face_verts, low_order_node, low_order_node_face_weight);
1900 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]);
1901 const double weight = basis_2d(face_orders[face_id], x, y, face_weight);
1902 if (std::abs(weight) < 1e-12)
1903 continue;
1904 b.bases[j].global().emplace_back(global_, nodes.node_position(global_), weight);
1905 }
1906
1907 // contribution to vertex nodes
1908 for (int i = 0; i < 3; i++)
1909 {
1910 const auto &global_ = b.bases[fv(local_face_id, i)].global();
1911 auto low_order_node = ncmesh.point(fv(local_face_id, i));
1912 Eigen::MatrixXd low_order_node_face_weight;
1913 global_to_local_face(face_verts, low_order_node, low_order_node_face_weight);
1914 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]);
1915 double weight = basis_2d(face_orders[face_id], x, y, face_weight);
1916 if (std::abs(weight) > 1e-12)
1917 {
1918 assert(global_.size() > 0);
1919 for (size_t ii = 0; ii < global_.size(); ++ii)
1920 b.bases[j].global().emplace_back(global_[ii].index, global_[ii].node, weight * global_[ii].val);
1921 }
1922 }
1923
1924 // contribution to edge nodes, two steps
1925 for (int x = 0, idx = 0; x <= face_orders[face_id]; x++)
1926 {
1927 for (int y = 0; x + y <= face_orders[face_id]; y++)
1928 {
1929 const int z = face_orders[face_id] - x - y;
1930 int flag = (int)(x == 0) + (int)(y == 0) + (int)(z == 0);
1931 if (flag != 1)
1932 continue;
1933
1934 // first step
1935 const double weight = basis_2d(face_orders[face_id], x, y, face_weight);
1936 if (std::abs(weight) < 1e-12)
1937 continue;
1938 Eigen::MatrixXd face_weight(1, 2);
1939 face_weight << (double)x / face_orders[face_id], (double)y / face_orders[face_id];
1940 Eigen::MatrixXd pos, local_pos;
1941 local_to_global_face(face_verts, face_weight, pos);
1942 global_to_local(verts, pos, local_pos);
1943 Local2Global step1(idx, local_pos, weight);
1944 idx++;
1945
1946 {
1947 // evaluate the basis of the large element at this node
1948 const auto &other_bases = bases[e];
1949 std::vector<AssemblyValues> w;
1950 other_bases.evaluate_bases(local_pos, w);
1951
1952 // apply basis projection
1953 for (long i = 0; i < w.size(); ++i)
1954 {
1955 assert(w[i].val.size() == 1);
1956 if (std::abs(w[i].val(0)) < 1e-12)
1957 continue;
1958
1959 assert(other_bases.bases[i].global().size() > 0);
1960 for (size_t ii = 0; ii < other_bases.bases[i].global().size(); ++ii)
1961 {
1962 const auto &other_global = other_bases.bases[i].global()[ii];
1963 assert(other_global.index >= 0);
1964 b.bases[j].global().emplace_back(other_global.index, other_global.node, step1.val * w[i].val(0) * other_global.val);
1965 }
1966 }
1967 }
1968 }
1969 }
1970 }
1971 else
1972 {
1973 Eigen::MatrixXd global_position, large_elem_verts(4, 3);
1974 auto v_large = tet_vertices_local_to_global(mesh, large_elem);
1975 for (int i = 0; i < ncmesh.n_cell_vertices(large_elem); i++)
1976 large_elem_verts.row(i) = ncmesh.point(v_large[i]);
1977 local_to_global(verts, local_position, global_position);
1978 global_to_local(large_elem_verts, global_position, local_position);
1979
1980 // evaluate the basis of the large element at this node
1981 const auto &other_bases = bases[large_elem];
1982 std::vector<AssemblyValues> w;
1983 other_bases.evaluate_bases(local_position, w);
1984
1985 // apply basis projection
1986 for (long i = 0; i < w.size(); ++i)
1987 {
1988 assert(w[i].val.size() == 1);
1989 if (std::abs(w[i].val(0)) < 1e-12)
1990 continue;
1991
1992 assert(other_bases.bases[i].global().size() > 0);
1993 for (size_t ii = 0; ii < other_bases.bases[i].global().size(); ++ii)
1994 {
1995 const auto &other_global = other_bases.bases[i].global()[ii];
1996 assert(other_global.index >= 0);
1997 b.bases[j].global().emplace_back(other_global.index, other_global.node, w[i].val(0) * other_global.val);
1998 }
1999 }
2000 }
2001 }
2002 else
2003 assert(false);
2004
2005 auto &global_ = b.bases[j].global();
2006 if (global_.size() <= 1)
2007 continue;
2008
2009 std::map<int, Local2Global> list;
2010 for (size_t ii = 0; ii < global_.size(); ii++)
2011 {
2012 auto pair = list.insert({global_[ii].index, global_[ii]});
2013 if (!pair.second && pair.first != list.end())
2014 {
2015 assert((pair.first->second.node - global_[ii].node).norm() < 1e-12);
2016 pair.first->second.val += global_[ii].val;
2017 }
2018 }
2019
2020 global_.clear();
2021 for (auto it = list.begin(); it != list.end(); ++it)
2022 {
2023 if (std::abs(it->second.val) > 1e-12)
2024 {
2025 global_.push_back(it->second);
2026 }
2027 }
2028 }
2029 }
2030 }
2031 });
2032 }
2033 }
2034 else
2035 {
2036 for (int pp = 2; pp <= autogen::MAX_P_BASES; ++pp)
2037 {
2038 for (int e : interface_elements)
2039 {
2040 ElementBases &b = bases[e];
2041 const int discr_order = discr_orders(e);
2042 const int n_el_bases = element_nodes_id[e].size();
2043 assert(discr_order > 1);
2044 if (discr_order != pp)
2045 continue;
2046
2047 if (mesh.is_cube(e))
2048 {
2049 // TODO
2050 assert(false);
2051 }
2052 else if (mesh.is_simplex(e))
2053 {
2054 for (int j = 0; j < n_el_bases; ++j)
2055 {
2056 const int global_index = element_nodes_id[e][j];
2057
2058 if (global_index >= 0)
2059 b.bases[j].init(discr_order, global_index, j, nodes.node_position(global_index));
2060 else
2061 {
2062 const int lnn = max_p > 2 ? (discr_order - 2) : 0;
2063 const int ln_edge_nodes = discr_order - 1;
2064 const int ln_face_nodes = lnn * (lnn + 1) / 2;
2065
2066 const auto v = tet_vertices_local_to_global(mesh, e);
2067 Navigation3D::Index index;
2068 if (global_index <= -30)
2069 {
2070 assert(false);
2071 // const auto lv = -(global_index + 30);
2072 // assert(lv>=0 && lv < 4);
2073 // assert(j < 4);
2074
2075 // if(lv == 3)
2076 // {
2077 // index = mesh.switch_element(find_edge(mesh, e, v[lv], v[0]));
2078 // if(index.element < 0)
2079 // index = mesh.switch_element(find_edge(mesh, e, v[lv], v[1]));
2080 // if(index.element < 0)
2081 // index = mesh.switch_element(find_edge(mesh, e, v[lv], v[2]));
2082 // }
2083 // else
2084 // {
2085 // index = mesh.switch_element(find_edge(mesh, e, v[lv], v[(lv+1)%3]));
2086 // if(index.element < 0)
2087 // index = mesh.switch_element(find_edge(mesh, e, v[lv], v[(lv+2)%3]));
2088 // if(index.element < 0)
2089 // index = mesh.switch_element(find_edge(mesh, e, v[lv], v[3]));
2090 // }
2091 }
2092 else if (global_index <= -10)
2093 {
2094 const auto le = -(global_index + 10);
2095 assert(le >= 0 && le < 6);
2096 assert(j >= 4 && j < 4 + 6 * ln_edge_nodes);
2097
2098 Eigen::Matrix<int, 6, 2> ev;
2099 ev.row(0) << v[0], v[1];
2100 ev.row(1) << v[1], v[2];
2101 ev.row(2) << v[2], v[0];
2102
2103 ev.row(3) << v[0], v[3];
2104 ev.row(4) << v[1], v[3];
2105 ev.row(5) << v[2], v[3];
2106
2107 // const auto edge_index = find_edge(mesh, e, ev(le, 0), ev(le, 1));
2108 const auto edge_index = mesh.get_index_from_element_edge(e, ev(le, 0), ev(le, 1));
2109 auto neighs = mesh.edge_neighs(edge_index.edge);
2110 int min_p = discr_order;
2111 int min_cell = index.element;
2112
2113 for (auto cid : neighs)
2114 {
2115 if (discr_orders[cid] < min_p)
2116 {
2117 min_p = discr_orders[cid];
2118 min_cell = cid;
2119 }
2120 }
2121
2122 bool found = false;
2123 for (int lf = 0; lf < 4; ++lf)
2124 {
2125 for (int lv = 0; lv < 4; ++lv)
2126 {
2127 index = mesh.get_index_from_element(min_cell, lf, lv);
2128
2129 if (index.vertex == edge_index.vertex)
2130 {
2131 if (index.edge != edge_index.edge)
2132 {
2133 auto tmp = index;
2134 index = mesh.switch_edge(tmp);
2135
2136 if (index.edge != edge_index.edge)
2137 {
2138 index = mesh.switch_edge(mesh.switch_face(tmp));
2139 }
2140 }
2141 found = true;
2142 break;
2143 }
2144 }
2145
2146 if (found)
2147 break;
2148 }
2149
2150 assert(found);
2151 assert(index.vertex == edge_index.vertex && index.edge == edge_index.edge);
2152 assert(index.element != edge_index.element);
2153 }
2154 else
2155 {
2156 const auto lf = -(global_index + 1);
2157 assert(lf >= 0 && lf < 4);
2158 assert(j >= 4 + 6 * ln_edge_nodes && j < 4 + 6 * ln_edge_nodes + 4 * ln_face_nodes);
2159
2160 Eigen::Matrix<int, 4, 3> fv;
2161 fv.row(0) << v[0], v[1], v[2];
2162 fv.row(1) << v[0], v[1], v[3];
2163 fv.row(2) << v[1], v[2], v[3];
2164 fv.row(3) << v[2], v[0], v[3];
2165
2166 index = mesh.switch_element(mesh.get_index_from_element_face(e, fv(lf, 0), fv(lf, 1), fv(lf, 2)));
2167 }
2168
2169 const auto other_cell = index.element;
2170 assert(other_cell >= 0);
2171 assert(discr_order > discr_orders(other_cell));
2172
2173 auto indices = tet_face_local_nodes(discr_order, mesh, index);
2174 Eigen::MatrixXd lnodes;
2175 autogen::p_nodes_3d(discr_order, lnodes);
2176 Eigen::RowVector3d node_position; // = lnodes.row(indices(ii));
2177
2178 if (j < 4)
2179 node_position = lnodes.row(indices(0));
2180 else if (j < 4 + 6 * ln_edge_nodes)
2181 node_position = lnodes.row(indices(((j - 4) % ln_edge_nodes) + 3));
2182 else if (j < 4 + 6 * ln_edge_nodes + 4 * ln_face_nodes)
2183 {
2184 // node_position = lnodes.row(indices(((j - 4 - 6*ln_edge_nodes) % ln_face_nodes) + 3 + 3*ln_edge_nodes));
2185 auto me_indices = tet_face_local_nodes(discr_order, mesh, mesh.switch_element(index));
2186 int ii;
2187 for (ii = 0; ii < me_indices.size(); ++ii)
2188 {
2189 if (me_indices(ii) == j)
2190 break;
2191 }
2192
2193 assert(ii >= 3 + 3 * ln_edge_nodes);
2194 assert(ii < me_indices.size());
2195
2196 node_position = lnodes.row(indices(ii));
2197 }
2198 else
2199 assert(false);
2200
2201 // std::cout<<indices.transpose()<<std::endl;
2202 // auto asd = quadr_tri_edge_local_nodes(mesh, index);
2203 // std::cout<<asd[0]<<" "<<asd[1]<<" "<<asd[2]<<std::endl;
2204
2205 // std::cout<<"\n"<<lnodes<<"\nnewp\n"<<node_position<<"\n"<<std::endl;
2206 // const auto param_p = quadr_tri_edge_local_nodes_coordinates(mesh, index);
2207
2208 // if( j < 3)
2209 // node_position = param_p.row(0);
2210 // else if( j < 3 + 3*(discr_order-1)){
2211 // node_position = param_p.row( (j-3) % (discr_order-1) + 1);
2212 // }
2213 // else
2214 // assert(false);
2215 // std::cout<<node_position<<"\n\n----\n"<<std::endl;
2216
2217 const auto &other_bases = bases[other_cell];
2218 // Eigen::MatrixXd w;
2219 std::vector<AssemblyValues> w;
2220 other_bases.evaluate_bases(node_position, w);
2221
2222 assert(b.bases[j].global().size() == 0);
2223
2224 for (long i = 0; i < w.size(); ++i)
2225 {
2226 assert(w[i].val.size() == 1);
2227 if (std::abs(w[i].val(0)) < 1e-8)
2228 continue;
2229
2230 // assert(other_bases.bases[i].global().size() == 1);
2231 for (size_t ii = 0; ii < other_bases.bases[i].global().size(); ++ii)
2232 {
2233 const auto &other_global = other_bases.bases[i].global()[ii];
2234 // std::cout<<"e "<<e<<" " <<j << " gid "<<other_global.index<<std::endl;
2235 b.bases[j].global().emplace_back(other_global.index, other_global.node, w[i].val(0) * other_global.val);
2236 }
2237 }
2238 }
2239 }
2240 }
2241 else
2242 {
2243 // Polygon bases are built later on
2244 }
2245 }
2246 }
2247 }
2248 }
2249
2250 return nodes.n_nodes();
2251}
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:177
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::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::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 serendipity, const bool has_polys, const bool is_geom_bases, 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::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:361
polyfem::mesh::Mesh::is_rational
bool is_rational() const
check if curved mesh has rational polynomials elements
Definition Mesh.hpp:271
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:348
polyfem::mesh::Mesh::is_boundary_face
virtual bool is_boundary_face(const int face_global_id) const =0
is face boundary
polyfem::mesh::Mesh::get_boundary_id
virtual int get_boundary_id(const int primitive) const
Get the boundary selection of an element (face in 3d, edge in 2d)
Definition Mesh.hpp:478
polyfem::mesh::Mesh::is_simplex
bool is_simplex(const int el_id) const
checks if element is simples compatible
Definition Mesh.cpp:418
polyfem::mesh::Mesh::cell_weights
const std::vector< double > & cell_weights(const int cell_index) const
weights for rational polynomial meshes
Definition Mesh.hpp:552
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:543
polyfem::mesh::NCMesh3D::leader_face_of_edge
int leader_face_of_edge(const int e_id) const
Definition NCMesh3D.hpp:303
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:291
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:309
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
polyfem::quadrature::TetQuadrature::get_quadrature
void get_quadrature(const int order, Quadrature &quad)
Definition TetQuadrature.cpp:31
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
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::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_nodes_3d
void p_nodes_3d(const int p, Eigen::MatrixXd &val)
Definition auto_p_bases.cpp:2805
polyfem::autogen::p_grad_basis_value_3d
void p_grad_basis_value_3d(const int p, const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val)
Definition auto_p_bases.cpp:2852
polyfem::autogen::p_basis_value_3d
void p_basis_value_3d(const int p, const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val)
Definition auto_p_bases.cpp:2828
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