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