PolyFEM
Loading...
Searching...
No Matches
SplineBasis3d.cpp
Go to the documentation of this file.
1#include "SplineBasis3d.hpp"
2
3#include "LagrangeBasis3d.hpp"
4#include "function/QuadraticBSpline3d.hpp"
5#include <polyfem/quadrature/HexQuadrature.hpp>
6
7#include <polyfem/assembler/AssemblerUtils.hpp>
8
9#include <polysolve/linear/Solver.hpp>
10#include <polyfem/mesh/MeshNodes.hpp>
11
12#include <polyfem/utils/Types.hpp>
13
14#include <polyfem/Common.hpp>
15#include <polyfem/autogen/auto_q_bases.hpp>
16
17#include <Eigen/Sparse>
18
19#include <cassert>
20#include <iostream>
21#include <vector>
22#include <array>
23#include <map>
24#include <numeric>
25
26// TODO carefull with simplices
27
28namespace polyfem
29{
30 using namespace Eigen;
31 using namespace polysolve;
32 using namespace assembler;
33 using namespace mesh;
34 using namespace quadrature;
35
36 namespace basis
37 {
38
39 namespace
40 {
41 class SpaceMatrix
42 {
43 public:
44 inline const int &operator()(const int i, const int j, const int k) const
45 {
46 return space_[k](i, j);
47 }
48
49 inline int &operator()(const int i, const int j, const int k)
50 {
51 return space_[k](i, j);
52 }
53
54 bool is_k_regular = false;
55 int x, y, z;
56 int edge_id;
57
58 bool is_regular(const int xx, const int yy, const int zz) const
59 {
60 if (!is_k_regular)
61 return true;
62
63 return !((x == 1 && y == yy && z == zz) || (x == xx && y == 1 && z == zz) || (x == xx && y == yy && z == 1));
64 // space(1, y, z).size() <= 1 && space(x, 1, z).size() <= 1 && space(x, y, 1).size() <= 1
65 }
66
67 private:
68 std::array<Matrix<int, 3, 3>, 3> space_;
69 };
70
71 bool is_edge_singular(const Navigation3D::Index &index, const Mesh3D &mesh)
72 {
73 std::vector<int> ids;
74 mesh.get_edge_elements_neighs(index.edge, ids);
75
76 if (ids.size() == 4 || mesh.is_boundary_edge(index.edge))
77 return false;
78
79 for (auto idx : ids)
80 {
81 if (mesh.is_polytope(idx))
82 return false;
83 }
84
85 return true;
86 }
87
88 void print_local_space(const SpaceMatrix &space)
89 {
90 for (int k = 2; k >= 0; --k)
91 {
92 for (int j = 2; j >= 0; --j)
93 {
94 for (int i = 0; i < 3; ++i)
95 {
96 // if(space(i, j, k).size() > 0){
97 // for(std::size_t l = 0; l < space(i, j, k).size(); ++l)
98 std::cout << space(i, j, k) << "\t";
99 // }
100 // else
101 // std::cout<<"x\t";
102 }
103 std::cout << std::endl;
104 }
105
106 std::cout << "\n"
107 << std::endl;
108 }
109 }
110
111 int node_id_from_face_index(const Mesh3D &mesh, MeshNodes &mesh_nodes, const Navigation3D::Index &index)
112 {
113 int el_id = mesh.switch_element(index).element;
114 if (el_id >= 0 && mesh.is_cube(el_id))
115 {
116 return mesh_nodes.node_id_from_cell(el_id);
117 }
118
119 return mesh_nodes.node_id_from_face(index.face);
120 }
121
122 int node_id_from_edge_index(const Mesh3D &mesh, MeshNodes &mesh_nodes, const Navigation3D::Index &index)
123 {
124 Navigation3D::Index new_index = mesh.switch_element(index);
125 int el_id = new_index.element;
126
127 if (el_id < 0 || mesh.is_polytope(el_id))
128 {
129 new_index = mesh.switch_element(mesh.switch_face(index));
130 el_id = new_index.element;
131
132 if (el_id < 0 || mesh.is_polytope(el_id))
133 {
134 return mesh_nodes.node_id_from_edge(index.edge);
135 }
136
137 return node_id_from_face_index(mesh, mesh_nodes, mesh.switch_face(new_index));
138 }
139
140 return node_id_from_face_index(mesh, mesh_nodes, mesh.switch_face(new_index));
141 }
142
143 int node_id_from_vertex_index_explore(const Mesh3D &mesh, const MeshNodes &mesh_nodes, const Navigation3D::Index &index, int &node_id)
144 {
145 Navigation3D::Index new_index = mesh.switch_element(index);
146
147 int id = new_index.element;
148
149 if (id < 0 || mesh.is_polytope(id))
150 {
151 // id = vertex_node_id(index.vertex);
152 // node = node_from_vertex(index.vertex);
153 node_id = mesh_nodes.primitive_from_vertex(index.vertex);
154 return 3;
155 }
156
157 new_index = mesh.switch_element(mesh.switch_face(new_index));
158 id = new_index.element;
159
160 if (id < 0 || mesh.is_polytope(id))
161 {
162 // id = edge_node_id(switch_edge(new_index).edge);
163 // node = node_from_edge(switch_edge(new_index).edge);
164 node_id = mesh_nodes.primitive_from_edge(mesh.switch_edge(new_index).edge);
165 return 2;
166 }
167
168 new_index = mesh.switch_element(mesh.switch_face(mesh.switch_edge(new_index)));
169 id = new_index.element;
170
171 if (id < 0 || mesh.is_polytope(id))
172 {
173 // id = face_node_id(new_index.face);
174 // node = node_from_face(new_index.face);
175 node_id = mesh_nodes.primitive_from_face(new_index.face);
176 return 1;
177 }
178
179 // node = node_from_element(id);
180 node_id = mesh_nodes.primitive_from_cell(id);
181 return 0;
182 }
183
184 int node_id_from_vertex_index(const Mesh3D &mesh, MeshNodes &mesh_nodes, const Navigation3D::Index &index)
185 {
186 std::array<int, 6> path;
187 std::array<int, 6> primitive_ids;
188
189 path[0] = node_id_from_vertex_index_explore(mesh, mesh_nodes, index, primitive_ids[0]);
190 path[1] = node_id_from_vertex_index_explore(mesh, mesh_nodes, mesh.switch_face(index), primitive_ids[1]);
191
192 path[2] = node_id_from_vertex_index_explore(mesh, mesh_nodes, mesh.switch_edge(index), primitive_ids[2]);
193 path[3] = node_id_from_vertex_index_explore(mesh, mesh_nodes, mesh.switch_face(mesh.switch_edge(index)), primitive_ids[3]);
194
195 path[4] = node_id_from_vertex_index_explore(mesh, mesh_nodes, mesh.switch_edge(mesh.switch_face(index)), primitive_ids[4]);
196 path[5] = node_id_from_vertex_index_explore(mesh, mesh_nodes, mesh.switch_face(mesh.switch_edge(mesh.switch_face(index))), primitive_ids[5]);
197
198 const int min_path = *std::min_element(path.begin(), path.end());
199
200 int primitive_id = 0;
201 for (int i = 0; i < 6; ++i)
202 {
203 if (path[i] == min_path)
204 {
205 primitive_id = primitive_ids[i];
206 break;
207 }
208 }
209
210 return mesh_nodes.node_id_from_primitive(primitive_id);
211 }
212
213 void get_edge_elements_neighs(const Mesh3D &mesh, MeshNodes &mesh_nodes, const int element_id, const int edge_id, int dir, std::vector<int> &ids)
214 {
215 std::array<std::function<Navigation3D::Index(Navigation3D::Index)>, 12> to_edge;
216 mesh.to_edge_functions(to_edge);
217
218 Navigation3D::Index index;
219 for (int i = 0; i < 12; ++i)
220 {
221 index = to_edge[i](mesh.get_index_from_element(element_id));
222
223 if (index.edge == edge_id)
224 break;
225 }
226
227 assert(index.edge == edge_id);
228
229 if (dir == 1)
230 {
231 int id;
232 do
233 {
234 ids.push_back(mesh_nodes.node_id_from_cell(index.element));
235 index = mesh.next_around_edge(index);
236 } while (index.element != element_id);
237
238 return;
239 }
240
241 if (dir == 0)
242 {
243 int id;
244 do
245 {
246 const Navigation3D::Index f_index = mesh.switch_face(mesh.switch_edge(index));
247 ids.push_back(node_id_from_face_index(mesh, mesh_nodes, f_index));
248
249 index = mesh.next_around_edge(index);
250 } while (index.element != element_id);
251
252 return;
253 }
254
255 if (dir == 2)
256 {
257 int id;
258 do
259 {
260 const Navigation3D::Index f_index = mesh.switch_face(mesh.switch_edge(mesh.switch_vertex(index)));
261 ids.push_back(node_id_from_face_index(mesh, mesh_nodes, f_index));
262
263 index = mesh.next_around_edge(index);
264 } while (index.element != element_id);
265
266 return;
267 }
268
269 assert(false);
270 }
271
272 void add_edge_id_for_poly(const Navigation3D::Index &index, const Mesh3D &mesh, MeshNodes &mesh_nodes, const int global_index, std::map<int, InterfaceData> &poly_face_to_data)
273 {
274 const int f1 = index.face;
275 const int f2 = mesh.switch_face(index).face;
276
277 const int id1 = mesh_nodes.primitive_from_face(f1);
278 const int id2 = mesh_nodes.primitive_from_face(f2);
279
280 if (mesh_nodes.is_primitive_interface(id1))
281 {
282 InterfaceData &data = poly_face_to_data[f1];
283 data.local_indices.push_back(global_index);
284 }
285
286 if (mesh_nodes.is_primitive_interface(id2))
287 {
288 InterfaceData &data = poly_face_to_data[f2];
289 data.local_indices.push_back(global_index);
290 }
291 }
292
293 void add_vertex_id_for_poly(const Navigation3D::Index &index, const Mesh3D &mesh, MeshNodes &mesh_nodes, const int global_index, std::map<int, InterfaceData> &poly_face_to_data)
294 {
295 const int f1 = index.face;
296 const int f2 = mesh.switch_face(index).face;
297 const int f3 = mesh.switch_face(mesh.switch_edge(index)).face;
298
299 const int id1 = mesh_nodes.primitive_from_face(f1);
300 const int id2 = mesh_nodes.primitive_from_face(f2);
301 const int id3 = mesh_nodes.primitive_from_face(f3);
302
303 if (mesh_nodes.is_primitive_interface(id1))
304 {
305 InterfaceData &data = poly_face_to_data[f1];
306 data.local_indices.push_back(global_index);
307 }
308
309 if (mesh_nodes.is_primitive_interface(id2))
310 {
311 InterfaceData &data = poly_face_to_data[f2];
312 data.local_indices.push_back(global_index);
313 }
314
315 if (mesh_nodes.is_primitive_interface(id3))
316 {
317 InterfaceData &data = poly_face_to_data[f3];
318 data.local_indices.push_back(global_index);
319 }
320 }
321
322 void explore_edge(const Navigation3D::Index &index, const Mesh3D &mesh, MeshNodes &mesh_nodes, const int x, const int y, const int z, SpaceMatrix &space, LocalBoundary &local_boundary, std::map<int, InterfaceData> &poly_face_to_data)
323 {
324 int node_id = node_id_from_edge_index(mesh, mesh_nodes, index);
325 space(x, y, z) = node_id;
326 // node(x, y, z) = mesh_nodes.node_position(node_id);
327
328 if (is_edge_singular(index, mesh))
329 {
330 std::vector<int> ids;
331 mesh.get_edge_elements_neighs(index.edge, ids);
332 // irregular edge
333
334 assert(!space.is_k_regular);
335 space.is_k_regular = true;
336 space.x = x;
337 space.y = y;
338 space.z = z;
339 space.edge_id = index.edge;
340 }
341
342 // if(mesh_nodes.is_boundary(node_id))
343 // bounday_nodes.push_back(node_id);
344
345 add_edge_id_for_poly(index, mesh, mesh_nodes, 9 * z + 3 * y + x, poly_face_to_data);
346 }
347
348 void explore_vertex(const Navigation3D::Index &index, const Mesh3D &mesh, MeshNodes &mesh_nodes, const int x, const int y, const int z, SpaceMatrix &space, LocalBoundary &local_boundary, std::map<int, InterfaceData> &poly_face_to_data)
349 {
350 int node_id = node_id_from_vertex_index(mesh, mesh_nodes, index);
351 space(x, y, z) = node_id;
352 // node(x, y, z) = mesh_nodes.node_position(node_id);
353
354 // if(mesh_nodes.is_boundary(node_id))
355 // bounday_nodes.push_back(node_id);
356
357 add_vertex_id_for_poly(index, mesh, mesh_nodes, 9 * z + 3 * y + x, poly_face_to_data);
358 }
359
360 void explore_face(const Navigation3D::Index &index, const Mesh3D &mesh, MeshNodes &mesh_nodes, const int x, const int y, const int z, SpaceMatrix &space, LocalBoundary &local_boundary, std::map<int, InterfaceData> &poly_face_to_data)
361 {
362 int node_id = node_id_from_face_index(mesh, mesh_nodes, index);
363 space(x, y, z) = node_id;
364 // node(x, y, z) = mesh_nodes.node_position(node_id);
365
366 if (mesh_nodes.is_boundary(node_id))
367 {
368 // local_boundary.add_boundary_primitive(index.face, LagrangeBasis3d::quadr_hex_face_local_nodes(mesh, index)[8]-20);
369 local_boundary.add_boundary_primitive(index.face, LagrangeBasis3d::hex_face_local_nodes(false, 2, mesh, index)[8] - 20);
370 // bounday_nodes.push_back(node_id);
371 }
372 else if (mesh_nodes.is_interface(node_id))
373 {
374 InterfaceData &data = poly_face_to_data[index.face];
375 data.local_indices.push_back(9 * z + 3 * y + x);
376 // igl::viewer::Viewer &viewer = UIState::ui_state().viewer;
377 // viewer.data.add_points(mesh_nodes.node_position(node_id), Eigen::MatrixXd::Constant(1, 3, 0));
378 }
379 }
380
381 void build_local_space(const Mesh3D &mesh, MeshNodes &mesh_nodes, const int el_index, SpaceMatrix &space, std::vector<LocalBoundary> &local_boundary, std::map<int, InterfaceData> &poly_face_to_data)
382 {
383 assert(mesh.is_volume());
384
385 Navigation3D::Index start_index = mesh.get_index_from_element(el_index);
386 Navigation3D::Index index;
387
388 std::array<std::function<Navigation3D::Index(Navigation3D::Index)>, 6> to_face;
389 mesh.to_face_functions(to_face);
390
391 std::array<std::function<Navigation3D::Index(Navigation3D::Index)>, 12> to_edge;
392 mesh.to_edge_functions(to_edge);
393
394 std::array<std::function<Navigation3D::Index(Navigation3D::Index)>, 8> to_vertex;
395 mesh.to_vertex_functions(to_vertex);
396
397 const int node_id = mesh_nodes.node_id_from_cell(el_index);
398 space(1, 1, 1) = node_id;
399 // node(1, 1, 1) = mesh_nodes.node_position(node_id);
400
401 LocalBoundary lb(el_index, BoundaryType::QUAD);
402
404 index = to_face[1](start_index);
405 explore_face(index, mesh, mesh_nodes, 1, 1, 0, space, lb, poly_face_to_data);
406
407 index = to_face[0](start_index);
408 explore_face(index, mesh, mesh_nodes, 1, 1, 2, space, lb, poly_face_to_data);
409
410 index = to_face[3](start_index);
411 explore_face(index, mesh, mesh_nodes, 0, 1, 1, space, lb, poly_face_to_data);
412
413 index = to_face[2](start_index);
414 explore_face(index, mesh, mesh_nodes, 2, 1, 1, space, lb, poly_face_to_data);
415
416 index = to_face[5](start_index);
417 explore_face(index, mesh, mesh_nodes, 1, 0, 1, space, lb, poly_face_to_data);
418
419 index = to_face[4](start_index);
420 explore_face(index, mesh, mesh_nodes, 1, 2, 1, space, lb, poly_face_to_data);
421
423 index = to_edge[0](start_index);
424 explore_edge(index, mesh, mesh_nodes, 1, 0, 0, space, lb, poly_face_to_data);
425
426 index = to_edge[1](start_index);
427 explore_edge(index, mesh, mesh_nodes, 2, 1, 0, space, lb, poly_face_to_data);
428
429 index = to_edge[2](start_index);
430 explore_edge(index, mesh, mesh_nodes, 1, 2, 0, space, lb, poly_face_to_data);
431
432 index = to_edge[3](start_index);
433 explore_edge(index, mesh, mesh_nodes, 0, 1, 0, space, lb, poly_face_to_data);
434
435 index = to_edge[4](start_index);
436 explore_edge(index, mesh, mesh_nodes, 0, 0, 1, space, lb, poly_face_to_data);
437
438 index = to_edge[5](start_index);
439 explore_edge(index, mesh, mesh_nodes, 2, 0, 1, space, lb, poly_face_to_data);
440
441 index = to_edge[6](start_index);
442 explore_edge(index, mesh, mesh_nodes, 2, 2, 1, space, lb, poly_face_to_data);
443
444 index = to_edge[7](start_index);
445 explore_edge(index, mesh, mesh_nodes, 0, 2, 1, space, lb, poly_face_to_data);
446
447 index = to_edge[8](start_index);
448 explore_edge(index, mesh, mesh_nodes, 1, 0, 2, space, lb, poly_face_to_data);
449
450 index = to_edge[9](start_index);
451 explore_edge(index, mesh, mesh_nodes, 2, 1, 2, space, lb, poly_face_to_data);
452
453 index = to_edge[10](start_index);
454 explore_edge(index, mesh, mesh_nodes, 1, 2, 2, space, lb, poly_face_to_data);
455
456 index = to_edge[11](start_index);
457 explore_edge(index, mesh, mesh_nodes, 0, 1, 2, space, lb, poly_face_to_data);
458
460 index = to_vertex[0](start_index);
461 explore_vertex(index, mesh, mesh_nodes, 0, 0, 0, space, lb, poly_face_to_data);
462
463 index = to_vertex[1](start_index);
464 explore_vertex(index, mesh, mesh_nodes, 2, 0, 0, space, lb, poly_face_to_data);
465
466 index = to_vertex[2](start_index);
467 explore_vertex(index, mesh, mesh_nodes, 2, 2, 0, space, lb, poly_face_to_data);
468
469 index = to_vertex[3](start_index);
470 explore_vertex(index, mesh, mesh_nodes, 0, 2, 0, space, lb, poly_face_to_data);
471
472 index = to_vertex[4](start_index);
473 explore_vertex(index, mesh, mesh_nodes, 0, 0, 2, space, lb, poly_face_to_data);
474
475 index = to_vertex[5](start_index);
476 explore_vertex(index, mesh, mesh_nodes, 2, 0, 2, space, lb, poly_face_to_data);
477
478 index = to_vertex[6](start_index);
479 explore_vertex(index, mesh, mesh_nodes, 2, 2, 2, space, lb, poly_face_to_data);
480
481 index = to_vertex[7](start_index);
482 explore_vertex(index, mesh, mesh_nodes, 0, 2, 2, space, lb, poly_face_to_data);
483
484 if (!lb.empty())
485 local_boundary.emplace_back(lb);
486 }
487
488 void setup_knots_vectors(MeshNodes &mesh_nodes, const SpaceMatrix &space, std::array<std::array<double, 4>, 3> &h_knots, std::array<std::array<double, 4>, 3> &v_knots, std::array<std::array<double, 4>, 3> &w_knots)
489 {
490 // left and right neigh are absent
491 if (mesh_nodes.is_boundary_or_interface(space(0, 1, 1)) && mesh_nodes.is_boundary_or_interface(space(2, 1, 1)))
492 {
493 h_knots[0] = {{0, 0, 0, 1}};
494 h_knots[1] = {{0, 0, 1, 1}};
495 h_knots[2] = {{0, 1, 1, 1}};
496 }
497 // left neigh is absent
498 else if (mesh_nodes.is_boundary_or_interface(space(0, 1, 1)))
499 {
500 h_knots[0] = {{0, 0, 0, 1}};
501 h_knots[1] = {{0, 0, 1, 2}};
502 h_knots[2] = {{0, 1, 2, 3}};
503 }
504 // right neigh is absent
505 else if (mesh_nodes.is_boundary_or_interface(space(2, 1, 1)))
506 {
507 h_knots[0] = {{-2, -1, 0, 1}};
508 h_knots[1] = {{-1, 0, 1, 1}};
509 h_knots[2] = {{0, 1, 1, 1}};
510 }
511 else
512 {
513 h_knots[0] = {{-2, -1, 0, 1}};
514 h_knots[1] = {{-1, 0, 1, 2}};
515 h_knots[2] = {{0, 1, 2, 3}};
516 }
517
518 // top and bottom neigh are absent
519 if (mesh_nodes.is_boundary_or_interface(space(1, 0, 1)) && mesh_nodes.is_boundary_or_interface(space(1, 2, 1)))
520 {
521 v_knots[0] = {{0, 0, 0, 1}};
522 v_knots[1] = {{0, 0, 1, 1}};
523 v_knots[2] = {{0, 1, 1, 1}};
524 }
525 // bottom neigh is absent
526 else if (mesh_nodes.is_boundary_or_interface(space(1, 0, 1)))
527 {
528 v_knots[0] = {{0, 0, 0, 1}};
529 v_knots[1] = {{0, 0, 1, 2}};
530 v_knots[2] = {{0, 1, 2, 3}};
531 }
532 // top neigh is absent
533 else if (mesh_nodes.is_boundary_or_interface(space(1, 2, 1)))
534 {
535 v_knots[0] = {{-2, -1, 0, 1}};
536 v_knots[1] = {{-1, 0, 1, 1}};
537 v_knots[2] = {{0, 1, 1, 1}};
538 }
539 else
540 {
541 v_knots[0] = {{-2, -1, 0, 1}};
542 v_knots[1] = {{-1, 0, 1, 2}};
543 v_knots[2] = {{0, 1, 2, 3}};
544 }
545
546 // front and back neigh are absent
547 if (mesh_nodes.is_boundary_or_interface(space(1, 1, 0)) && mesh_nodes.is_boundary_or_interface(space(1, 1, 2)))
548 {
549 w_knots[0] = {{0, 0, 0, 1}};
550 w_knots[1] = {{0, 0, 1, 1}};
551 w_knots[2] = {{0, 1, 1, 1}};
552 }
553 // back neigh is absent
554 else if (mesh_nodes.is_boundary_or_interface(space(1, 1, 0)))
555 {
556 w_knots[0] = {{0, 0, 0, 1}};
557 w_knots[1] = {{0, 0, 1, 2}};
558 w_knots[2] = {{0, 1, 2, 3}};
559 }
560 // front neigh is absent
561 else if (mesh_nodes.is_boundary_or_interface(space(1, 1, 2)))
562 {
563 w_knots[0] = {{-2, -1, 0, 1}};
564 w_knots[1] = {{-1, 0, 1, 1}};
565 w_knots[2] = {{0, 1, 1, 1}};
566 }
567 else
568 {
569 w_knots[0] = {{-2, -1, 0, 1}};
570 w_knots[1] = {{-1, 0, 1, 2}};
571 w_knots[2] = {{0, 1, 2, 3}};
572 }
573 }
574
575 void basis_for_regular_hex(MeshNodes &mesh_nodes, const SpaceMatrix &space, const std::array<std::array<double, 4>, 3> &h_knots, const std::array<std::array<double, 4>, 3> &v_knots, const std::array<std::array<double, 4>, 3> &w_knots, ElementBases &b)
576 {
577 for (int z = 0; z < 3; ++z)
578 {
579 for (int y = 0; y < 3; ++y)
580 {
581 for (int x = 0; x < 3; ++x)
582 {
583 if (space.is_regular(x, y, z)) // space(1, y, z).size() <= 1 && space(x, 1, z).size() <= 1 && space(x, y, 1).size() <= 1)
584 {
585 const int local_index = 9 * z + 3 * y + x;
586 const int global_index = space(x, y, z);
587 const auto node = mesh_nodes.node_position(global_index);
588 // loc_nodes(x, y, z);
589
590 b.bases[local_index].init(2, global_index, local_index, node);
591
592 const QuadraticBSpline3d spline(h_knots[x], v_knots[y], w_knots[z]);
593
594 b.bases[local_index].set_basis([spline](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { spline.interpolate(uv, val); });
595 b.bases[local_index].set_grad([spline](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { spline.derivative(uv, val); });
596 }
597 }
598 }
599 }
600 }
601
602 void basis_for_irregulard_hex(const int el_index, const Mesh3D &mesh, MeshNodes &mesh_nodes, const SpaceMatrix &space, const std::array<std::array<double, 4>, 3> &h_knots, const std::array<std::array<double, 4>, 3> &v_knots, const std::array<std::array<double, 4>, 3> &w_knots, ElementBases &b, std::map<int, InterfaceData> &poly_face_to_data)
603 {
604 for (int z = 0; z < 3; ++z)
605 {
606 for (int y = 0; y < 3; ++y)
607 {
608 for (int x = 0; x < 3; ++x)
609 {
610 if (!space.is_regular(x, y, z)) // space(1, y, z).size() > 1 || space(x, 1, z).size() > 1 || space(x, y, 1).size() > 1)
611 {
612 const int local_index = 9 * z + 3 * y + x;
613
614 int mpx = -1;
615 int mpy = -1;
616 int mpz = -1;
617
618 int mmx = -1;
619 int mmy = -1;
620 int mmz = -1;
621
622 int xx = 1;
623 int yy = 1;
624 int zz = 1;
625
626 const int edge_id = space.edge_id;
627 int dir = -1;
628
629 if (space.x == x && space.y == y && space.z == 1)
630 {
631 mpx = 1;
632 mpy = y;
633 mpz = z;
634
635 mmx = x;
636 mmy = 1;
637 mmz = z;
638
639 zz = z;
640 dir = z;
641 }
642 else if (space.x == x && space.y == 1 && space.z == z)
643 {
644 mpx = 1;
645 mpy = y;
646 mpz = z;
647
648 mmx = x;
649 mmy = y;
650 mmz = 1;
651
652 yy = y;
653 dir = y;
654 }
655 else if (space.x == 1 && space.y == y && space.z == z)
656 {
657 mpx = x;
658 mpy = y;
659 mpz = 1;
660
661 mmx = x;
662 mmy = 1;
663 mmz = z;
664
665 xx = x;
666 dir = x;
667 }
668 else
669 assert(false);
670
671 const auto &center = b.bases[zz * 9 + yy * 3 + xx].global().front();
672
673 const auto &el1 = b.bases[mpz * 9 + mpy * 3 + mpx].global().front();
674 const auto &el2 = b.bases[mmz * 9 + mmy * 3 + mmx].global().front();
675
676 // std::cout<<"el_index "<<el_index<<std::endl;
677 // std::cout<<"center "<<center.index<<std::endl;
678 // std::cout<<"el1 "<<el1.index<<std::endl;
679 // std::cout<<"el2 "<<el2.index<<std::endl;
680
681 std::vector<int> ids;
682 get_edge_elements_neighs(mesh, mesh_nodes, el_index, edge_id, dir, ids);
683
684 if (ids.front() != center.index)
685 {
686 assert(dir != 1);
687 ids.clear();
688 get_edge_elements_neighs(mesh, mesh_nodes, el_index, edge_id, dir == 2 ? 0 : 2, ids);
689 }
690
691 assert(ids.front() == center.index);
692
693 std::vector<int> other_indices;
694 std::vector<Eigen::MatrixXd> other_nodes;
695 for (size_t i = 0; i < ids.size(); ++i)
696 {
697 const int node_id = ids[i];
698 if (node_id != center.index && node_id != el1.index && node_id != el2.index)
699 {
700 other_indices.push_back(node_id);
701 }
702 }
703
704 // std::cout<<ids.size()<< " " << other_indices.size()<<std::endl;
705
706 auto &base = b.bases[local_index];
707
708 const int k = int(other_indices.size()) + 3;
709
710 // const bool is_interface = mesh_nodes.is_interface(center.index);
711 // const int face_id = is_interface ? mesh_nodes.face_from_node_id(center.index) : -1;
712
713 base.global().resize(k);
714
715 base.global()[0].index = center.index;
716 base.global()[0].val = (4. - k) / k;
717 base.global()[0].node = center.node;
718
719 base.global()[1].index = el1.index;
720 base.global()[1].val = (4. - k) / k;
721 base.global()[1].node = el1.node;
722
723 base.global()[2].index = el2.index;
724 base.global()[2].val = (4. - k) / k;
725 base.global()[2].node = el2.node;
726
727 // if(is_interface){
728 // poly_face_to_data[face_id].local_indices.push_back(local_index);
729 // }
730
731 for (std::size_t n = 0; n < other_indices.size(); ++n)
732 {
733 base.global()[3 + n].index = other_indices[n];
734 base.global()[3 + n].val = 4. / k;
735 base.global()[3 + n].node = mesh_nodes.node_position(other_indices[n]);
736 }
737
738 const QuadraticBSpline3d spline(h_knots[x], v_knots[y], w_knots[z]);
739
740 b.bases[local_index].set_basis([spline](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { spline.interpolate(uv, val); });
741 b.bases[local_index].set_grad([spline](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { spline.derivative(uv, val); });
742 }
743 }
744 }
745 }
746 }
747
748 void create_q2_nodes(const Mesh3D &mesh, const int el_index, std::set<int> &vertex_id, std::set<int> &edge_id, std::set<int> &face_id, ElementBases &b, std::vector<LocalBoundary> &local_boundary, int &n_bases)
749 {
750 b.bases.resize(27);
751
752 std::array<std::function<Navigation3D::Index(Navigation3D::Index)>, 6> to_face;
753 mesh.to_face_functions(to_face);
754
755 std::array<std::function<Navigation3D::Index(Navigation3D::Index)>, 12> to_edge;
756 mesh.to_edge_functions(to_edge);
757
758 std::array<std::function<Navigation3D::Index(Navigation3D::Index)>, 8> to_vertex;
759 mesh.to_vertex_functions(to_vertex);
760
761 LocalBoundary lb(el_index, BoundaryType::QUAD);
762
763 const Navigation3D::Index start_index = mesh.get_index_from_element(el_index);
764 for (int j = 0; j < 8; ++j)
765 {
766 const Navigation3D::Index index = to_vertex[j](start_index);
767 // const int loc_index = LagrangeBasis3d::quadr_hex_face_local_nodes(mesh, index)[0];
768 const int loc_index = LagrangeBasis3d::hex_face_local_nodes(false, 2, mesh, index)[0];
769
770 int current_vertex_node_id = -1;
771 Eigen::MatrixXd current_vertex_node;
772
773 // if the edge/vertex is boundary the it is a Q2 edge
774 bool is_vertex_q2 = true;
775
776 std::vector<int> vertex_neighs;
777 mesh.get_vertex_elements_neighs(index.vertex, vertex_neighs);
778
779 for (size_t i = 0; i < vertex_neighs.size(); ++i)
780 {
781 if (mesh.is_spline_compatible(vertex_neighs[i]))
782 {
783 is_vertex_q2 = false;
784 break;
785 }
786 }
787 const bool is_vertex_boundary = mesh.is_boundary_vertex(index.vertex);
788
789 if (is_vertex_q2)
790 {
791 const bool is_new_vertex = vertex_id.insert(index.vertex).second;
792
793 if (is_new_vertex)
794 {
795 current_vertex_node_id = n_bases++;
796 current_vertex_node = mesh.point(index.vertex);
797
798 // if(is_vertex_boundary)//mesh.is_vertex_boundary(index.vertex))
799 // bounday_nodes.push_back(current_vertex_node_id);
800 }
801 }
802
803 // init new Q2 nodes
804 if (current_vertex_node_id >= 0)
805 b.bases[loc_index].init(2, current_vertex_node_id, loc_index, current_vertex_node);
806
807 b.bases[loc_index].set_basis([loc_index](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { autogen::q_basis_value_3d(2, loc_index, uv, val); });
808 b.bases[loc_index].set_grad([loc_index](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { autogen::q_grad_basis_value_3d(2, loc_index, uv, val); });
809 }
810
811 for (int j = 0; j < 12; ++j)
812 {
813 Navigation3D::Index index = to_edge[j](start_index);
814
815 int current_edge_node_id = -1;
816 Eigen::Matrix<double, 1, 3> current_edge_node;
817 // const int loc_index = LagrangeBasis3d::quadr_hex_face_local_nodes(mesh, index)[1];
818 const int loc_index = LagrangeBasis3d::hex_face_local_nodes(false, 2, mesh, index)[4];
819
820 bool is_edge_q2 = true;
821
822 std::vector<int> edge_neighs;
823 mesh.get_edge_elements_neighs(index.edge, edge_neighs);
824
825 for (size_t i = 0; i < edge_neighs.size(); ++i)
826 {
827 if (mesh.is_spline_compatible(edge_neighs[i]))
828 {
829 is_edge_q2 = false;
830 break;
831 }
832 }
833 const bool is_edge_boundary = mesh.is_boundary_edge(index.edge);
834
835 if (is_edge_q2)
836 {
837 const bool is_new_edge = edge_id.insert(index.edge).second;
838
839 if (is_new_edge)
840 {
841 current_edge_node_id = n_bases++;
842 current_edge_node = mesh.edge_barycenter(index.edge);
843
844 // if(is_edge_boundary)
845 // bounday_nodes.push_back(current_edge_node_id);
846 }
847 }
848
849 // init new Q2 nodes
850 if (current_edge_node_id >= 0)
851 b.bases[loc_index].init(2, current_edge_node_id, loc_index, current_edge_node);
852
853 b.bases[loc_index].set_basis([loc_index](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { autogen::q_basis_value_3d(2, loc_index, uv, val); });
854 b.bases[loc_index].set_grad([loc_index](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { autogen::q_grad_basis_value_3d(2, loc_index, uv, val); });
855 }
856
857 for (int j = 0; j < 6; ++j)
858 {
859 Navigation3D::Index index = to_face[j](start_index);
860
861 int current_face_node_id = -1;
862
863 Eigen::Matrix<double, 1, 3> current_face_node;
864 const int opposite_element = mesh.switch_element(index).element;
865 const bool is_face_q2 = opposite_element < 0 || !mesh.is_spline_compatible(opposite_element);
866 // const int loc_index = LagrangeBasis3d::quadr_hex_face_local_nodes(mesh, index)[8];
867 const int loc_index = LagrangeBasis3d::hex_face_local_nodes(false, 2, mesh, index)[8];
868
869 if (is_face_q2)
870 {
871 const bool is_new_face = face_id.insert(index.face).second;
872
873 if (is_new_face)
874 {
875 current_face_node_id = n_bases++;
876 current_face_node = mesh.face_barycenter(index.face);
877
878 const int b_index = loc_index - 20;
879
880 if (opposite_element < 0) // && mesh.n_element_faces(opposite_element) == 6 && mesh.n_element_vertices(opposite_element) == 8)
881 {
882 // bounday_nodes.push_back(current_face_node_id);
883 lb.add_boundary_primitive(index.face, b_index);
884 }
885 }
886 }
887
888 // init new Q2 nodes
889 if (current_face_node_id >= 0)
890 b.bases[loc_index].init(2, current_face_node_id, loc_index, current_face_node);
891
892 b.bases[loc_index].set_basis([loc_index](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { autogen::q_basis_value_3d(2, loc_index, uv, val); });
893 b.bases[loc_index].set_grad([loc_index](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { autogen::q_grad_basis_value_3d(2, loc_index, uv, val); });
894 }
895
896 // //central node always present
897 b.bases[26].init(2, n_bases++, 26, mesh.cell_barycenter(el_index));
898 b.bases[26].set_basis([](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { autogen::q_basis_value_3d(2, 26, uv, val); });
899 b.bases[26].set_grad([](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { autogen::q_grad_basis_value_3d(2, 26, uv, val); });
900
901 if (!lb.empty())
902 local_boundary.emplace_back(lb);
903 }
904
905 void insert_into_global(const int el_index, const Local2Global &data, std::vector<Local2Global> &vec, const int size)
906 {
907 // ignore small weights
908 if (fabs(data.val) < 1e-10)
909 return;
910
911 bool found = false;
912
913 for (int i = 0; i < size; ++i)
914 {
915 if (vec[i].index == data.index)
916 {
917 // if(fabs(vec[i].val - data.val) > 1e-10){
918 // std::cout<<el_index <<" "<<vec[i].val <<" "<< data.val<<" "<<fabs(vec[i].val - data.val)<<std::endl;
919 // // vec[i].val += data.val;
920 // }
921 // assert(fabs(vec[i].val - data.val) < 1e-10);
922 assert((vec[i].node - data.node).norm() < 1e-10);
923 found = true;
924 break;
925 }
926 }
927
928 if (!found)
929 vec.push_back(data);
930 }
931
932 void assign_q2_weights(const Mesh3D &mesh, const int el_index, std::vector<ElementBases> &bases)
933 {
934 // Eigen::MatrixXd eval_p;
935 std::vector<AssemblyValues> eval_p;
936 const Navigation3D::Index start_index = mesh.get_index_from_element(el_index);
937 ElementBases &b = bases[el_index];
938
939 std::array<std::function<Navigation3D::Index(Navigation3D::Index)>, 6> to_face;
940 mesh.to_face_functions(to_face);
941 for (int f = 0; f < 6; ++f)
942 {
943 const Navigation3D::Index index = to_face[f](start_index);
944 const int opposite_element = mesh.switch_element(index).element;
945
946 if (opposite_element < 0 || !mesh.is_cube(opposite_element))
947 continue;
948
949 // const auto &param_p = LagrangeBasis3d::quadr_hex_face_local_nodes_coordinates(mesh, mesh.switch_element(index));
950 Eigen::Matrix<double, 9, 3> param_p;
951 {
952 Eigen::MatrixXd hex_loc_nodes;
953 polyfem::autogen::q_nodes_3d(2, hex_loc_nodes);
954 const auto opposite_indices = LagrangeBasis3d::hex_face_local_nodes(false, 2, mesh, mesh.switch_element(index));
955 for (int k = 0; k < 9; ++k)
956 param_p.row(k) = hex_loc_nodes.row(opposite_indices[k]);
957 }
958 const auto &other_bases = bases[opposite_element];
959
960 // const auto &indices = LagrangeBasis3d::quadr_hex_face_local_nodes(mesh, index);
961 const auto &indices = LagrangeBasis3d::hex_face_local_nodes(false, 2, mesh, index);
962
963 std::array<int, 9> sizes;
964
965 for (int l = 0; l < 9; ++l)
966 sizes[l] = b.bases[indices[l]].global().size();
967
968 other_bases.evaluate_bases(param_p, eval_p);
969 for (std::size_t i = 0; i < other_bases.bases.size(); ++i)
970 {
971 const auto &other_b = other_bases.bases[i];
972
973 if (other_b.global().empty())
974 continue;
975
976 // other_b.basis(param_p, eval_p);
977 assert(eval_p[i].val.size() == 9);
978
979 // basis i of element opposite element is zero on this elements
980 if (eval_p[i].val.cwiseAbs().maxCoeff() <= 1e-10)
981 continue;
982
983 for (std::size_t k = 0; k < other_b.global().size(); ++k)
984 {
985 for (int l = 0; l < 9; ++l)
986 {
987 Local2Global glob = other_b.global()[k];
988 glob.val *= eval_p[i].val(l);
989
990 insert_into_global(el_index, glob, b.bases[indices[l]].global(), sizes[l]);
991 }
992 }
993 }
994 }
995 }
996
997 void setup_data_for_polygons(const Mesh3D &mesh, const int el_index, const ElementBases &b, std::map<int, InterfaceData> &poly_face_to_data)
998 {
999 const Navigation3D::Index start_index = mesh.get_index_from_element(el_index);
1000 std::array<std::function<Navigation3D::Index(Navigation3D::Index)>, 6> to_face;
1001 mesh.to_face_functions(to_face);
1002 for (int f = 0; f < 6; ++f)
1003 {
1004 const Navigation3D::Index index = to_face[f](start_index);
1005
1006 const int opposite_element = mesh.switch_element(index).element;
1007 const bool is_neigh_poly = opposite_element >= 0 && mesh.is_polytope(opposite_element);
1008
1009 if (is_neigh_poly)
1010 {
1011 // auto e2l = LagrangeBasis3d::quadr_hex_face_local_nodes(mesh, index);
1012 auto e2l = LagrangeBasis3d::hex_face_local_nodes(false, 2, mesh, index);
1013
1014 InterfaceData &data = poly_face_to_data[index.face];
1015
1016 for (int kk = 0; kk < e2l.size(); ++kk)
1017 {
1018 const auto idx = e2l(kk);
1019 data.local_indices.push_back(idx);
1020 }
1021 }
1022 }
1023 }
1024 } // namespace
1025
1026 int SplineBasis3d::build_bases(const Mesh3D &mesh,
1027 const std::string &assembler,
1028 const int quadrature_order, const int mass_quadrature_order, std::vector<ElementBases> &bases, std::vector<LocalBoundary> &local_boundary, std::map<int, InterfaceData> &poly_face_to_data)
1029 {
1030 using std::max;
1031 assert(mesh.is_volume());
1032
1033 MeshNodes mesh_nodes(mesh, true, true, 1, 1, 1);
1034
1035 const int n_els = mesh.n_elements();
1036 bases.resize(n_els);
1037 local_boundary.clear();
1038
1039 // bounday_nodes.clear();
1040
1041 // HexQuadrature hex_quadrature;
1042
1043 std::array<std::array<double, 4>, 3> h_knots;
1044 std::array<std::array<double, 4>, 3> v_knots;
1045 std::array<std::array<double, 4>, 3> w_knots;
1046
1047 for (int e = 0; e < n_els; ++e)
1048 {
1049 if (!mesh.is_spline_compatible(e))
1050 continue;
1051
1052 SpaceMatrix space;
1053
1054 build_local_space(mesh, mesh_nodes, e, space, local_boundary, poly_face_to_data);
1055
1056 ElementBases &b = bases[e];
1057 const int real_order = quadrature_order > 0 ? quadrature_order : AssemblerUtils::quadrature_order(assembler, 2, AssemblerUtils::BasisType::SPLINE, 3);
1058 const int real_mass_order = mass_quadrature_order > 0 ? mass_quadrature_order : AssemblerUtils::quadrature_order("Mass", 2, AssemblerUtils::BasisType::SPLINE, 3);
1059
1060 b.set_quadrature([real_order](Quadrature &quad) {
1061 HexQuadrature hex_quadrature;
1062 hex_quadrature.get_quadrature(real_order, quad);
1063 });
1064 b.set_mass_quadrature([real_mass_order](Quadrature &quad) {
1065 HexQuadrature hex_quadrature;
1066 hex_quadrature.get_quadrature(real_mass_order, quad);
1067 });
1068 // hex_quadrature.get_quadrature(quadrature_order, b.quadrature);
1069 b.bases.resize(27);
1070
1071 b.set_local_node_from_primitive_func([e](const int primitive_id, const Mesh &mesh) {
1072 const auto &mesh3d = dynamic_cast<const Mesh3D &>(mesh);
1073
1074 std::array<std::function<Navigation3D::Index(Navigation3D::Index)>, 6> to_face;
1075 mesh3d.to_face_functions(to_face);
1076
1077 auto start_index = mesh3d.get_index_from_element(e);
1078 auto index = start_index;
1079
1080 int lf;
1081 for (lf = 0; lf < mesh3d.n_cell_faces(e); ++lf)
1082 {
1083 index = to_face[lf](start_index);
1084 if (index.face == primitive_id)
1085 break;
1086 }
1087 assert(index.face == primitive_id);
1088
1089 static constexpr std::array<std::array<int, 9>, 6> face_to_index = {{
1090 {{2 * 9 + 0 * 3 + 0, 2 * 9 + 1 * 3 + 0, 2 * 9 + 2 * 3 + 0, 2 * 9 + 0 * 3 + 1, 2 * 9 + 1 * 3 + 1, 2 * 9 + 2 * 3 + 1, 2 * 9 + 0 * 3 + 2, 2 * 9 + 1 * 3 + 2, 2 * 9 + 2 * 3 + 2}}, // 0
1091 {{0 * 9 + 0 * 3 + 0, 0 * 9 + 1 * 3 + 0, 0 * 9 + 2 * 3 + 0, 0 * 9 + 0 * 3 + 1, 0 * 9 + 1 * 3 + 1, 0 * 9 + 2 * 3 + 1, 0 * 9 + 0 * 3 + 2, 0 * 9 + 1 * 3 + 2, 0 * 9 + 2 * 3 + 2}}, // 1
1092
1093 {{0 * 9 + 0 * 3 + 2, 0 * 9 + 1 * 3 + 2, 0 * 9 + 2 * 3 + 2, 1 * 9 + 0 * 3 + 2, 1 * 9 + 1 * 3 + 2, 1 * 9 + 2 * 3 + 2, 2 * 9 + 0 * 3 + 2, 2 * 9 + 1 * 3 + 2, 2 * 9 + 2 * 3 + 2}}, // 2
1094 {{0 * 9 + 0 * 3 + 0, 0 * 9 + 1 * 3 + 0, 0 * 9 + 2 * 3 + 0, 1 * 9 + 0 * 3 + 0, 1 * 9 + 1 * 3 + 0, 1 * 9 + 2 * 3 + 0, 2 * 9 + 0 * 3 + 0, 2 * 9 + 1 * 3 + 0, 2 * 9 + 2 * 3 + 0}}, // 3
1095
1096 {{0 * 9 + 2 * 3 + 0, 0 * 9 + 2 * 3 + 1, 0 * 9 + 2 * 3 + 2, 1 * 9 + 2 * 3 + 0, 1 * 9 + 2 * 3 + 1, 1 * 9 + 2 * 3 + 2, 2 * 9 + 2 * 3 + 0, 2 * 9 + 2 * 3 + 1, 2 * 9 + 2 * 3 + 2}}, // 4
1097 {{0 * 9 + 0 * 3 + 0, 0 * 9 + 0 * 3 + 1, 0 * 9 + 0 * 3 + 2, 1 * 9 + 0 * 3 + 0, 1 * 9 + 0 * 3 + 1, 1 * 9 + 0 * 3 + 2, 2 * 9 + 0 * 3 + 0, 2 * 9 + 0 * 3 + 1, 2 * 9 + 0 * 3 + 2}}, // 5
1098 }};
1099
1100 Eigen::VectorXi res(9);
1101
1102 for (int i = 0; i < 9; ++i)
1103 res(i) = face_to_index[lf][i];
1104
1105 return res;
1106 });
1107
1108 setup_knots_vectors(mesh_nodes, space, h_knots, v_knots, w_knots);
1109 // print_local_space(space);
1110
1111 basis_for_regular_hex(mesh_nodes, space, h_knots, v_knots, w_knots, b);
1112 basis_for_irregulard_hex(e, mesh, mesh_nodes, space, h_knots, v_knots, w_knots, b, poly_face_to_data);
1113 }
1114
1115 int n_bases = mesh_nodes.n_nodes();
1116
1117 std::set<int> face_id;
1118 std::set<int> edge_id;
1119 std::set<int> vertex_id;
1120
1121 for (int e = 0; e < n_els; ++e)
1122 {
1123 if (mesh.is_polytope(e) || mesh.is_spline_compatible(e))
1124 continue;
1125
1126 ElementBases &b = bases[e];
1127
1128 const int real_order = quadrature_order > 0 ? quadrature_order : AssemblerUtils::quadrature_order(assembler, 2, AssemblerUtils::BasisType::CUBE_LAGRANGE, 3);
1129 const int real_mass_order = mass_quadrature_order > 0 ? mass_quadrature_order : AssemblerUtils::quadrature_order("Mass", 2, AssemblerUtils::BasisType::CUBE_LAGRANGE, 3);
1130
1131 // hex_quadrature.get_quadrature(quadrature_order, b.quadrature);
1132 b.set_quadrature([real_order](Quadrature &quad) {
1133 HexQuadrature hex_quadrature;
1134 hex_quadrature.get_quadrature(real_order, quad);
1135 });
1136 b.set_mass_quadrature([real_mass_order](Quadrature &quad) {
1137 HexQuadrature hex_quadrature;
1138 hex_quadrature.get_quadrature(real_mass_order, quad);
1139 });
1140
1141 b.set_local_node_from_primitive_func([e](const int primitive_id, const Mesh &mesh) {
1142 const auto &mesh3d = dynamic_cast<const Mesh3D &>(mesh);
1143 Navigation3D::Index index;
1144
1145 for (int lf = 0; lf < 6; ++lf)
1146 {
1147 index = mesh3d.get_index_from_element(e, lf, 0);
1148 if (index.face == primitive_id)
1149 break;
1150 }
1151 assert(index.face == primitive_id);
1152
1153 // const auto indices = LagrangeBasis3d::quadr_hex_face_local_nodes(mesh3d, index);
1154 const auto indices = LagrangeBasis3d::hex_face_local_nodes(false, 2, mesh3d, index);
1155 Eigen::VectorXi res(indices.size());
1156
1157 for (size_t i = 0; i < indices.size(); ++i)
1158 res(i) = indices[i];
1159
1160 return res;
1161 });
1162
1163 create_q2_nodes(mesh, e, vertex_id, edge_id, face_id, b, local_boundary, n_bases);
1164 }
1165
1166 bool missing_bases = false;
1167 do
1168 {
1169 missing_bases = false;
1170 for (int e = 0; e < n_els; ++e)
1171 {
1172 if (mesh.is_polytope(e) || mesh.is_spline_compatible(e))
1173 continue;
1174
1175 auto &b = bases[e];
1176 if (b.is_complete())
1177 continue;
1178
1179 assign_q2_weights(mesh, e, bases);
1180
1181 missing_bases = missing_bases || b.is_complete();
1182 }
1183 } while (missing_bases);
1184
1185 for (int e = 0; e < n_els; ++e)
1186 {
1187 if (mesh.is_polytope(e) || mesh.is_spline_compatible(e))
1188 continue;
1189
1190 const ElementBases &b = bases[e];
1191 setup_data_for_polygons(mesh, e, b, poly_face_to_data);
1192 }
1193
1194 // for(int e = 0; e < n_els; ++e)
1195 // {
1196 // if(!mesh.is_polytope(e))
1197 // continue;
1198
1199 // for (int lf = 0; lf < mesh.n_cell_faces(e); ++lf)
1200 // {
1201 // auto index = mesh.get_index_from_element(e, lf, 0);
1202 // auto index2 = mesh.switch_element(index);
1203 // if (index2.element >= 0) {
1204 // auto &array = poly_face_to_data[index.face].local_indices;
1205 // auto &b = bases[index2.element];
1206 // array.resize(b.bases.size());
1207 // std::iota(array.begin(), array.end(), 0);
1208 // }
1209 // }
1210 // }
1211
1212 for (auto &k : poly_face_to_data)
1213 {
1214 auto &array = k.second.local_indices;
1215 std::sort(array.begin(), array.end());
1216 auto it = std::unique(array.begin(), array.end());
1217 array.resize(std::distance(array.begin(), it));
1218 }
1219
1220 return n_bases;
1221 }
1222
1223 void SplineBasis3d::fit_nodes(const Mesh3D &mesh, const int n_bases, std::vector<ElementBases> &gbases)
1224 {
1225 assert(false);
1226 // const int dim = 3;
1227 // const int n_constraints = 27;
1228 // const int n_elements = mesh.n_elements();
1229
1230 // std::vector< Eigen::Triplet<double> > entries, entries_t;
1231
1232 // MeshNodes nodes(mesh, true, true, 1, 1, 1);
1233 // // Eigen::MatrixXd tmp;
1234 // std::vector<AssemblyValues> tmp_val;
1235
1236 // Eigen::MatrixXd node_rhs(n_constraints*n_elements, dim);
1237 // Eigen::MatrixXd samples(n_constraints, dim);
1238
1239 // for(int i = 0; i < n_constraints; ++i)
1240 // samples.row(i) = LagrangeBasis3d::quadr_hex_local_node_coordinates(i);
1241
1242 // for(int i = 0; i < n_elements; ++i)
1243 // {
1244 // auto &base = gbases[i];
1245
1246 // if(!mesh.is_cube(i))
1247 // continue;
1248
1249 // auto global_ids = LagrangeBasis3d::quadr_hex_local_to_global(mesh, i);
1250 // assert(global_ids.size() == n_constraints);
1251
1252 // for(int j = 0; j < n_constraints; ++j)
1253 // {
1254 // auto n_id = nodes.node_id_from_primitive(global_ids[j]);
1255 // auto n = nodes.node_position(n_id);
1256 // for(int d = 0; d < dim; ++d)
1257 // node_rhs(n_constraints*i + j, d) = n(d);
1258 // }
1259
1260 // base.evaluate_bases(samples, tmp_val);
1261 // const auto &lbs = base.bases;
1262
1263 // const int n_local_bases = int(lbs.size());
1264 // for(int j = 0; j < n_local_bases; ++j)
1265 // {
1266 // const Basis &b = lbs[j];
1267 // const auto &tmp = tmp_val[j].val;
1268
1269 // for(std::size_t ii = 0; ii < b.global().size(); ++ii)
1270 // {
1271 // for (long k = 0; k < tmp.size(); ++k)
1272 // {
1273 // entries.emplace_back(n_constraints*i + k, b.global()[ii].index, tmp(k)*b.global()[ii].val);
1274 // entries_t.emplace_back(b.global()[ii].index, n_constraints*i + k, tmp(k)*b.global()[ii].val);
1275 // }
1276 // }
1277 // }
1278 // }
1279
1280 // Eigen::MatrixXd new_nodes(n_bases, dim);
1281 // {
1282 // StiffnessMatrix mat(n_constraints*n_elements, n_bases);
1283 // StiffnessMatrix mat_t(n_bases, n_constraints*n_elements);
1284
1285 // mat.setFromTriplets(entries.begin(), entries.end());
1286 // mat_t.setFromTriplets(entries_t.begin(), entries_t.end());
1287
1288 // StiffnessMatrix A = mat_t * mat;
1289 // Eigen::MatrixXd b = mat_t * node_rhs;
1290
1291 // json params = {
1292 // {"mtype", -2}, // matrix type for Pardiso (2 = SPD)
1293 // // {"max_iter", 0}, // for iterative solvers
1294 // // {"tolerance", 1e-9}, // for iterative solvers
1295 // };
1296 // auto solver = LinearSolver::create("", "");
1297 // solver->setParameters(params);
1298 // solver->analyzePattern(A);
1299 // solver->factorize(A);
1300
1301 // for(int d = 0; d < dim; ++d)
1302 // solver->solve(b.col(d), new_nodes.col(d));
1303 // }
1304
1305 // for(int i = 0; i < n_elements; ++i)
1306 // {
1307 // auto &base = gbases[i];
1308
1309 // if(!mesh.is_cube(i))
1310 // continue;
1311
1312 // auto &lbs = base.bases;
1313 // const int n_local_bases = int(lbs.size());
1314 // for(int j = 0; j < n_local_bases; ++j)
1315 // {
1316 // Basis &b = lbs[j];
1317
1318 // for(std::size_t ii = 0; ii < b.global().size(); ++ii)
1319 // {
1320 // // if(nodes.is_primitive_boundary(b.global()[ii].index))
1321 // // continue;
1322
1323 // for(int d = 0; d < dim; ++d)
1324 // {
1325 // b.global()[ii].node(d) = new_nodes(b.global()[ii].index, d);
1326 // }
1327 // }
1328 // }
1329 // }
1330 }
1331 } // namespace basis
1332} // namespace polyfem
vec
Eigen::MatrixXd vec
Definition Assembler.cpp:72
val
double val
Definition Assembler.cpp:86
quadrature
Quadrature quadrature
Definition MassMatrixAssembler.cpp:30
y
int y
Definition SplineBasis3d.cpp:55
z
int z
Definition SplineBasis3d.cpp:55
edge_id
int edge_id
Definition SplineBasis3d.cpp:56
is_k_regular
bool is_k_regular
Definition SplineBasis3d.cpp:54
x
int x
Definition SplineBasis3d.cpp:55
space_
std::array< Matrix< int, 3, 3 >, 3 > space_
Definition SplineBasis3d.cpp:68
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::SplineBasis3d::build_bases
static int build_bases(const mesh::Mesh3D &mesh, const std::string &assembler, const int quadrature_order, const int mass_quadrature_order, std::vector< ElementBases > &bases, std::vector< mesh::LocalBoundary > &local_boundary, std::map< int, InterfaceData > &poly_face_to_data)
Definition SplineBasis3d.cpp:1026
polyfem::basis::SplineBasis3d::fit_nodes
static void fit_nodes(const mesh::Mesh3D &mesh, const int n_bases, std::vector< ElementBases > &gbases)
Definition SplineBasis3d.cpp:1223
polyfem::mesh::Mesh3D::is_volume
bool is_volume() const override
checks if mesh is volume
Definition Mesh3D.hpp:28
polyfem::mesh::Mesh
Abstract mesh class to capture 2d/3d conforming and non-conforming meshes.
Definition Mesh.hpp:39
polyfem::mesh::Mesh::n_elements
int n_elements() const
utitlity to return the number of elements, cells or faces in 3d and 2d
Definition Mesh.hpp:161
polyfem::mesh::Mesh::is_polytope
bool is_polytope(const int el_id) const
checks if element is polygon compatible
Definition Mesh.cpp:365
polyfem::mesh::Mesh::is_spline_compatible
bool is_spline_compatible(const int el_id) const
checks if element is spline compatible
Definition Mesh.cpp:332
polyfem::quadrature::HexQuadrature::get_quadrature
void get_quadrature(const int order, Quadrature &quad)
Definition HexQuadrature.cpp:17
eigs.path
path
Definition eigs.py:128
elasticity_rhs.f
f
Definition elasticity_rhs.py:173
p_bases.base
str base
Definition p_bases.py:402
p_bases.indices
list indices
Definition p_bases.py:232
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::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