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