9#include <polysolve/linear/Solver.hpp>
17#include <Eigen/Sparse>
30 using namespace Eigen;
32 using namespace assembler;
44 inline const int &operator()(
const int i,
const int j,
const int k)
const
49 inline int &operator()(
const int i,
const int j,
const int k)
58 bool is_regular(
const int xx,
const int yy,
const int zz)
const
63 return !((x == 1 && y == yy && z == zz) || (x == xx && y == 1 && z == zz) || (x == xx && y == yy && z == 1));
68 std::array<Matrix<int, 3, 3>, 3>
space_;
71 bool is_edge_singular(
const Navigation3D::Index &index,
const Mesh3D &mesh)
74 mesh.get_edge_elements_neighs(index.edge, ids);
76 if (ids.size() == 4 || mesh.is_boundary_edge(index.edge))
81 if (mesh.is_polytope(idx))
88 void print_local_space(
const SpaceMatrix &space)
91 for (
int k = 2; k >= 0; --k)
93 for (
int j = 2; j >= 0; --j)
95 for (
int i = 0; i < 3; ++i)
99 ss << space(i, j, k) <<
"\t";
111 logger().trace(
"Local space:\n{}", ss.str());
114 int node_id_from_face_index(
const Mesh3D &mesh, MeshNodes &mesh_nodes,
const Navigation3D::Index &index)
116 int el_id = mesh.switch_element(index).element;
117 if (el_id >= 0 && mesh.is_cube(el_id))
119 return mesh_nodes.node_id_from_cell(el_id);
122 return mesh_nodes.node_id_from_face(index.face);
125 int node_id_from_edge_index(
const Mesh3D &mesh, MeshNodes &mesh_nodes,
const Navigation3D::Index &index)
127 Navigation3D::Index new_index = mesh.switch_element(index);
128 int el_id = new_index.element;
130 if (el_id < 0 || mesh.is_polytope(el_id))
132 new_index = mesh.switch_element(mesh.switch_face(index));
133 el_id = new_index.element;
135 if (el_id < 0 || mesh.is_polytope(el_id))
137 return mesh_nodes.node_id_from_edge(index.edge);
140 return node_id_from_face_index(mesh, mesh_nodes, mesh.switch_face(new_index));
143 return node_id_from_face_index(mesh, mesh_nodes, mesh.switch_face(new_index));
146 int node_id_from_vertex_index_explore(
const Mesh3D &mesh,
const MeshNodes &mesh_nodes,
const Navigation3D::Index &index,
int &node_id)
148 Navigation3D::Index new_index = mesh.switch_element(index);
150 int id = new_index.element;
152 if (
id < 0 || mesh.is_polytope(
id))
156 node_id = mesh_nodes.primitive_from_vertex(index.vertex);
160 new_index = mesh.switch_element(mesh.switch_face(new_index));
161 id = new_index.element;
163 if (
id < 0 || mesh.is_polytope(
id))
167 node_id = mesh_nodes.primitive_from_edge(mesh.switch_edge(new_index).edge);
171 new_index = mesh.switch_element(mesh.switch_face(mesh.switch_edge(new_index)));
172 id = new_index.element;
174 if (
id < 0 || mesh.is_polytope(
id))
178 node_id = mesh_nodes.primitive_from_face(new_index.face);
183 node_id = mesh_nodes.primitive_from_cell(
id);
187 int node_id_from_vertex_index(
const Mesh3D &mesh, MeshNodes &mesh_nodes,
const Navigation3D::Index &index)
189 std::array<int, 6>
path;
190 std::array<int, 6> primitive_ids;
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]);
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]);
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]);
201 const int min_path = *std::min_element(
path.begin(),
path.end());
203 int primitive_id = 0;
204 for (
int i = 0; i < 6; ++i)
206 if (path[i] == min_path)
208 primitive_id = primitive_ids[i];
213 return mesh_nodes.node_id_from_primitive(primitive_id);
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)
218 std::array<std::function<Navigation3D::Index(Navigation3D::Index)>, 12> to_edge;
219 mesh.to_edge_functions(to_edge);
221 Navigation3D::Index index;
222 for (
int i = 0; i < 12; ++i)
224 index = to_edge[i](mesh.get_index_from_element(element_id));
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);
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));
252 index = mesh.next_around_edge(index);
253 }
while (index.element != element_id);
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));
266 index = mesh.next_around_edge(index);
267 }
while (index.element != element_id);
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)
277 const int f1 = index.face;
278 const int f2 = mesh.switch_face(index).face;
280 const int id1 = mesh_nodes.primitive_from_face(f1);
281 const int id2 = mesh_nodes.primitive_from_face(f2);
283 if (mesh_nodes.is_primitive_interface(id1))
285 InterfaceData &data = poly_face_to_data[f1];
286 data.local_indices.push_back(global_index);
289 if (mesh_nodes.is_primitive_interface(id2))
291 InterfaceData &data = poly_face_to_data[f2];
292 data.local_indices.push_back(global_index);
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)
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;
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);
306 if (mesh_nodes.is_primitive_interface(id1))
308 InterfaceData &data = poly_face_to_data[f1];
309 data.local_indices.push_back(global_index);
312 if (mesh_nodes.is_primitive_interface(id2))
314 InterfaceData &data = poly_face_to_data[f2];
315 data.local_indices.push_back(global_index);
318 if (mesh_nodes.is_primitive_interface(id3))
320 InterfaceData &data = poly_face_to_data[f3];
321 data.local_indices.push_back(global_index);
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)
327 int node_id = node_id_from_edge_index(mesh, mesh_nodes, index);
328 space(
x,
y,
z) = node_id;
331 if (is_edge_singular(index, mesh))
333 std::vector<int> ids;
334 mesh.get_edge_elements_neighs(index.edge, ids);
337 assert(!space.is_k_regular);
338 space.is_k_regular =
true;
342 space.edge_id = index.edge;
348 add_edge_id_for_poly(index, mesh, mesh_nodes, 9 *
z + 3 *
y +
x, poly_face_to_data);
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)
353 int node_id = node_id_from_vertex_index(mesh, mesh_nodes, index);
354 space(
x,
y,
z) = node_id;
360 add_vertex_id_for_poly(index, mesh, mesh_nodes, 9 *
z + 3 *
y +
x, poly_face_to_data);
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)
365 int node_id = node_id_from_face_index(mesh, mesh_nodes, index);
366 space(
x,
y,
z) = node_id;
369 if (mesh_nodes.is_boundary(node_id))
375 else if (mesh_nodes.is_interface(node_id))
377 InterfaceData &data = poly_face_to_data[index.face];
378 data.local_indices.push_back(9 *
z + 3 *
y +
x);
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)
386 assert(mesh.is_volume());
388 Navigation3D::Index start_index = mesh.get_index_from_element(el_index);
389 Navigation3D::Index index;
391 std::array<std::function<Navigation3D::Index(Navigation3D::Index)>, 6> to_face;
392 mesh.to_face_functions(to_face);
394 std::array<std::function<Navigation3D::Index(Navigation3D::Index)>, 12> to_edge;
395 mesh.to_edge_functions(to_edge);
397 std::array<std::function<Navigation3D::Index(Navigation3D::Index)>, 8> to_vertex;
398 mesh.to_vertex_functions(to_vertex);
400 const int node_id = mesh_nodes.node_id_from_cell(el_index);
401 space(1, 1, 1) = node_id;
404 LocalBoundary lb(el_index, BoundaryType::QUAD);
407 index = to_face[1](start_index);
408 explore_face(index, mesh, mesh_nodes, 1, 1, 0, space, lb, poly_face_to_data);
410 index = to_face[0](start_index);
411 explore_face(index, mesh, mesh_nodes, 1, 1, 2, space, lb, poly_face_to_data);
413 index = to_face[3](start_index);
414 explore_face(index, mesh, mesh_nodes, 0, 1, 1, space, lb, poly_face_to_data);
416 index = to_face[2](start_index);
417 explore_face(index, mesh, mesh_nodes, 2, 1, 1, space, lb, poly_face_to_data);
419 index = to_face[5](start_index);
420 explore_face(index, mesh, mesh_nodes, 1, 0, 1, space, lb, poly_face_to_data);
422 index = to_face[4](start_index);
423 explore_face(index, mesh, mesh_nodes, 1, 2, 1, space, lb, poly_face_to_data);
426 index = to_edge[0](start_index);
427 explore_edge(index, mesh, mesh_nodes, 1, 0, 0, space, lb, poly_face_to_data);
429 index = to_edge[1](start_index);
430 explore_edge(index, mesh, mesh_nodes, 2, 1, 0, space, lb, poly_face_to_data);
432 index = to_edge[2](start_index);
433 explore_edge(index, mesh, mesh_nodes, 1, 2, 0, space, lb, poly_face_to_data);
435 index = to_edge[3](start_index);
436 explore_edge(index, mesh, mesh_nodes, 0, 1, 0, space, lb, poly_face_to_data);
438 index = to_edge[4](start_index);
439 explore_edge(index, mesh, mesh_nodes, 0, 0, 1, space, lb, poly_face_to_data);
441 index = to_edge[5](start_index);
442 explore_edge(index, mesh, mesh_nodes, 2, 0, 1, space, lb, poly_face_to_data);
444 index = to_edge[6](start_index);
445 explore_edge(index, mesh, mesh_nodes, 2, 2, 1, space, lb, poly_face_to_data);
447 index = to_edge[7](start_index);
448 explore_edge(index, mesh, mesh_nodes, 0, 2, 1, space, lb, poly_face_to_data);
450 index = to_edge[8](start_index);
451 explore_edge(index, mesh, mesh_nodes, 1, 0, 2, space, lb, poly_face_to_data);
453 index = to_edge[9](start_index);
454 explore_edge(index, mesh, mesh_nodes, 2, 1, 2, space, lb, poly_face_to_data);
456 index = to_edge[10](start_index);
457 explore_edge(index, mesh, mesh_nodes, 1, 2, 2, space, lb, poly_face_to_data);
459 index = to_edge[11](start_index);
460 explore_edge(index, mesh, mesh_nodes, 0, 1, 2, space, lb, poly_face_to_data);
463 index = to_vertex[0](start_index);
464 explore_vertex(index, mesh, mesh_nodes, 0, 0, 0, space, lb, poly_face_to_data);
466 index = to_vertex[1](start_index);
467 explore_vertex(index, mesh, mesh_nodes, 2, 0, 0, space, lb, poly_face_to_data);
469 index = to_vertex[2](start_index);
470 explore_vertex(index, mesh, mesh_nodes, 2, 2, 0, space, lb, poly_face_to_data);
472 index = to_vertex[3](start_index);
473 explore_vertex(index, mesh, mesh_nodes, 0, 2, 0, space, lb, poly_face_to_data);
475 index = to_vertex[4](start_index);
476 explore_vertex(index, mesh, mesh_nodes, 0, 0, 2, space, lb, poly_face_to_data);
478 index = to_vertex[5](start_index);
479 explore_vertex(index, mesh, mesh_nodes, 2, 0, 2, space, lb, poly_face_to_data);
481 index = to_vertex[6](start_index);
482 explore_vertex(index, mesh, mesh_nodes, 2, 2, 2, space, lb, poly_face_to_data);
484 index = to_vertex[7](start_index);
485 explore_vertex(index, mesh, mesh_nodes, 0, 2, 2, space, lb, poly_face_to_data);
488 local_boundary.emplace_back(lb);
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)
494 if (mesh_nodes.is_boundary_or_interface(space(0, 1, 1)) && mesh_nodes.is_boundary_or_interface(space(2, 1, 1)))
496 h_knots[0] = {{0, 0, 0, 1}};
497 h_knots[1] = {{0, 0, 1, 1}};
498 h_knots[2] = {{0, 1, 1, 1}};
501 else if (mesh_nodes.is_boundary_or_interface(space(0, 1, 1)))
503 h_knots[0] = {{0, 0, 0, 1}};
504 h_knots[1] = {{0, 0, 1, 2}};
505 h_knots[2] = {{0, 1, 2, 3}};
508 else if (mesh_nodes.is_boundary_or_interface(space(2, 1, 1)))
510 h_knots[0] = {{-2, -1, 0, 1}};
511 h_knots[1] = {{-1, 0, 1, 1}};
512 h_knots[2] = {{0, 1, 1, 1}};
516 h_knots[0] = {{-2, -1, 0, 1}};
517 h_knots[1] = {{-1, 0, 1, 2}};
518 h_knots[2] = {{0, 1, 2, 3}};
522 if (mesh_nodes.is_boundary_or_interface(space(1, 0, 1)) && mesh_nodes.is_boundary_or_interface(space(1, 2, 1)))
524 v_knots[0] = {{0, 0, 0, 1}};
525 v_knots[1] = {{0, 0, 1, 1}};
526 v_knots[2] = {{0, 1, 1, 1}};
529 else if (mesh_nodes.is_boundary_or_interface(space(1, 0, 1)))
531 v_knots[0] = {{0, 0, 0, 1}};
532 v_knots[1] = {{0, 0, 1, 2}};
533 v_knots[2] = {{0, 1, 2, 3}};
536 else if (mesh_nodes.is_boundary_or_interface(space(1, 2, 1)))
538 v_knots[0] = {{-2, -1, 0, 1}};
539 v_knots[1] = {{-1, 0, 1, 1}};
540 v_knots[2] = {{0, 1, 1, 1}};
544 v_knots[0] = {{-2, -1, 0, 1}};
545 v_knots[1] = {{-1, 0, 1, 2}};
546 v_knots[2] = {{0, 1, 2, 3}};
550 if (mesh_nodes.is_boundary_or_interface(space(1, 1, 0)) && mesh_nodes.is_boundary_or_interface(space(1, 1, 2)))
552 w_knots[0] = {{0, 0, 0, 1}};
553 w_knots[1] = {{0, 0, 1, 1}};
554 w_knots[2] = {{0, 1, 1, 1}};
557 else if (mesh_nodes.is_boundary_or_interface(space(1, 1, 0)))
559 w_knots[0] = {{0, 0, 0, 1}};
560 w_knots[1] = {{0, 0, 1, 2}};
561 w_knots[2] = {{0, 1, 2, 3}};
564 else if (mesh_nodes.is_boundary_or_interface(space(1, 1, 2)))
566 w_knots[0] = {{-2, -1, 0, 1}};
567 w_knots[1] = {{-1, 0, 1, 1}};
568 w_knots[2] = {{0, 1, 1, 1}};
572 w_knots[0] = {{-2, -1, 0, 1}};
573 w_knots[1] = {{-1, 0, 1, 2}};
574 w_knots[2] = {{0, 1, 2, 3}};
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)
580 for (
int z = 0;
z < 3; ++
z)
582 for (
int y = 0;
y < 3; ++
y)
584 for (
int x = 0;
x < 3; ++
x)
586 if (space.is_regular(
x,
y,
z))
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);
593 b.bases[local_index].init(2, global_index, local_index, node);
595 const QuadraticBSpline3d spline(h_knots[
x], v_knots[
y], w_knots[
z]);
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); });
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)
607 for (
int z = 0;
z < 3; ++
z)
609 for (
int y = 0;
y < 3; ++
y)
611 for (
int x = 0;
x < 3; ++
x)
613 if (!space.is_regular(
x,
y,
z))
615 const int local_index = 9 *
z + 3 *
y +
x;
629 const int edge_id = space.edge_id;
632 if (space.x ==
x && space.y ==
y && space.z == 1)
645 else if (space.x ==
x && space.y == 1 && space.z ==
z)
658 else if (space.x == 1 && space.y ==
y && space.z ==
z)
674 const auto ¢er = b.bases[zz * 9 + yy * 3 + xx].global().front();
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();
679 std::vector<int> ids;
680 get_edge_elements_neighs(mesh, mesh_nodes, el_index,
edge_id, dir, ids);
682 if (ids.front() != center.index)
686 get_edge_elements_neighs(mesh, mesh_nodes, el_index,
edge_id, dir == 2 ? 0 : 2, ids);
689 assert(ids.front() == center.index);
691 std::vector<int> other_indices;
692 std::vector<Eigen::MatrixXd> other_nodes;
693 for (
size_t i = 0; i < ids.size(); ++i)
695 const int node_id = ids[i];
696 if (node_id != center.index && node_id != el1.index && node_id != el2.index)
698 other_indices.push_back(node_id);
702 auto &
base = b.bases[local_index];
704 const int k = int(other_indices.size()) + 3;
709 base.global().resize(k);
711 base.global()[0].index = center.index;
712 base.global()[0].val = (4. - k) / k;
713 base.global()[0].node = center.node;
715 base.global()[1].index = el1.index;
716 base.global()[1].val = (4. - k) / k;
717 base.global()[1].node = el1.node;
719 base.global()[2].index = el2.index;
720 base.global()[2].val = (4. - k) / k;
721 base.global()[2].node = el2.node;
727 for (std::size_t n = 0; n < other_indices.size(); ++n)
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]);
734 const QuadraticBSpline3d spline(h_knots[
x], v_knots[
y], w_knots[
z]);
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); });
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)
748 std::array<std::function<Navigation3D::Index(Navigation3D::Index)>, 6> to_face;
749 mesh.to_face_functions(to_face);
751 std::array<std::function<Navigation3D::Index(Navigation3D::Index)>, 12> to_edge;
752 mesh.to_edge_functions(to_edge);
754 std::array<std::function<Navigation3D::Index(Navigation3D::Index)>, 8> to_vertex;
755 mesh.to_vertex_functions(to_vertex);
757 LocalBoundary lb(el_index, BoundaryType::QUAD);
759 const Navigation3D::Index start_index = mesh.get_index_from_element(el_index);
760 for (
int j = 0; j < 8; ++j)
762 const Navigation3D::Index index = to_vertex[j](start_index);
766 int current_vertex_node_id = -1;
767 Eigen::MatrixXd current_vertex_node;
770 bool is_vertex_q2 =
true;
772 std::vector<int> vertex_neighs;
773 mesh.get_vertex_elements_neighs(index.vertex, vertex_neighs);
775 for (
size_t i = 0; i < vertex_neighs.size(); ++i)
777 if (mesh.is_spline_compatible(vertex_neighs[i]))
779 is_vertex_q2 =
false;
783 const bool is_vertex_boundary = mesh.is_boundary_vertex(index.vertex);
787 const bool is_new_vertex = vertex_id.insert(index.vertex).second;
791 current_vertex_node_id = n_bases++;
792 current_vertex_node = mesh.point(index.vertex);
800 if (current_vertex_node_id >= 0)
801 b.bases[loc_index].init(2, current_vertex_node_id, loc_index, current_vertex_node);
807 for (
int j = 0; j < 12; ++j)
809 Navigation3D::Index index = to_edge[j](start_index);
811 int current_edge_node_id = -1;
812 Eigen::Matrix<double, 1, 3> current_edge_node;
816 bool is_edge_q2 =
true;
818 std::vector<int> edge_neighs;
819 mesh.get_edge_elements_neighs(index.edge, edge_neighs);
821 for (
size_t i = 0; i < edge_neighs.size(); ++i)
823 if (mesh.is_spline_compatible(edge_neighs[i]))
829 const bool is_edge_boundary = mesh.is_boundary_edge(index.edge);
833 const bool is_new_edge =
edge_id.insert(index.edge).second;
837 current_edge_node_id = n_bases++;
838 current_edge_node = mesh.edge_barycenter(index.edge);
846 if (current_edge_node_id >= 0)
847 b.bases[loc_index].init(2, current_edge_node_id, loc_index, current_edge_node);
853 for (
int j = 0; j < 6; ++j)
855 Navigation3D::Index index = to_face[j](start_index);
857 int current_face_node_id = -1;
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);
867 const bool is_new_face = face_id.insert(index.face).second;
871 current_face_node_id = n_bases++;
872 current_face_node = mesh.face_barycenter(index.face);
874 const int b_index = loc_index - 20;
876 if (opposite_element < 0)
879 lb.add_boundary_primitive(index.face, b_index);
885 if (current_face_node_id >= 0)
886 b.bases[loc_index].init(2, current_face_node_id, loc_index, current_face_node);
893 b.bases[26].init(2, n_bases++, 26, mesh.cell_barycenter(el_index));
898 local_boundary.emplace_back(lb);
901 void insert_into_global(
const int el_index,
const Local2Global &data, std::vector<Local2Global> &
vec,
const int size)
904 if (fabs(data.val) < 1e-10)
909 for (
int i = 0; i < size; ++i)
911 if (
vec[i].index == data.index)
914 assert((
vec[i].node - data.node).norm() < 1e-10);
924 void assign_q2_weights(
const Mesh3D &mesh,
const int el_index, std::vector<ElementBases> &bases)
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];
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)
935 const Navigation3D::Index index = to_face[
f](start_index);
936 const int opposite_element = mesh.switch_element(index).element;
938 if (opposite_element < 0 || !mesh.is_cube(opposite_element))
942 Eigen::Matrix<double, 9, 3> param_p;
944 Eigen::MatrixXd hex_loc_nodes;
947 for (
int k = 0; k < 9; ++k)
948 param_p.row(k) = hex_loc_nodes.row(opposite_indices[k]);
950 const auto &other_bases = bases[opposite_element];
955 std::array<int, 9> sizes;
957 for (
int l = 0; l < 9; ++l)
958 sizes[l] = b.bases[indices[l]].global().size();
960 other_bases.evaluate_bases(param_p, eval_p);
961 for (std::size_t i = 0; i < other_bases.bases.size(); ++i)
963 const auto &other_b = other_bases.bases[i];
965 if (other_b.global().empty())
969 assert(eval_p[i].
val.size() == 9);
972 if (eval_p[i].
val.cwiseAbs().maxCoeff() <= 1e-10)
975 for (std::size_t k = 0; k < other_b.global().size(); ++k)
977 for (
int l = 0; l < 9; ++l)
979 Local2Global glob = other_b.global()[k];
980 glob.val *= eval_p[i].val(l);
982 insert_into_global(el_index, glob, b.bases[indices[l]].global(), sizes[l]);
989 void setup_data_for_polygons(
const Mesh3D &mesh,
const int el_index,
const ElementBases &b, std::map<int, InterfaceData> &poly_face_to_data)
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)
996 const Navigation3D::Index index = to_face[
f](start_index);
998 const int opposite_element = mesh.switch_element(index).element;
999 const bool is_neigh_poly = opposite_element >= 0 && mesh.is_polytope(opposite_element);
1006 InterfaceData &data = poly_face_to_data[index.face];
1008 for (
int kk = 0; kk < e2l.size(); ++kk)
1010 const auto idx = e2l(kk);
1011 data.local_indices.push_back(idx);
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)
1025 MeshNodes mesh_nodes(mesh,
true,
true, 1, 1, 1);
1028 bases.resize(n_els);
1029 local_boundary.clear();
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;
1039 for (
int e = 0; e < n_els; ++e)
1046 build_local_space(mesh, mesh_nodes, e, space, local_boundary, poly_face_to_data);
1050 const int real_mass_order = mass_quadrature_order > 0 ? mass_quadrature_order :
AssemblerUtils::quadrature_order(
"Mass", 2, AssemblerUtils::BasisType::SPLINE, 3);
1064 const auto &mesh3d =
dynamic_cast<const Mesh3D &
>(mesh);
1067 mesh3d.to_face_functions(to_face);
1069 auto start_index = mesh3d.get_index_from_element(e);
1070 auto index = start_index;
1073 for (lf = 0; lf < mesh3d.n_cell_faces(e); ++lf)
1075 index = to_face[lf](start_index);
1076 if (index.face == primitive_id)
1079 assert(index.face == primitive_id);
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}},
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}},
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}},
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}},
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}},
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}},
1092 Eigen::VectorXi res(9);
1094 for (
int i = 0; i < 9; ++i)
1095 res(i) = face_to_index[lf][i];
1100 setup_knots_vectors(mesh_nodes, space, h_knots, v_knots, w_knots);
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);
1107 int n_bases = mesh_nodes.
n_nodes();
1109 std::set<int> face_id;
1111 std::set<int> vertex_id;
1113 for (
int e = 0; e < n_els; ++e)
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);
1134 const auto &mesh3d =
dynamic_cast<const Mesh3D &
>(mesh);
1137 for (
int lf = 0; lf < 6; ++lf)
1139 index = mesh3d.get_index_from_element(e, lf, 0);
1140 if (index.
face == primitive_id)
1143 assert(index.
face == primitive_id);
1147 Eigen::VectorXi res(indices.size());
1149 for (
size_t i = 0; i < indices.size(); ++i)
1150 res(i) = indices[i];
1155 create_q2_nodes(mesh, e, vertex_id,
edge_id, face_id, b, local_boundary, n_bases);
1158 bool missing_bases =
false;
1161 missing_bases =
false;
1162 for (
int e = 0; e < n_els; ++e)
1168 if (b.is_complete())
1171 assign_q2_weights(mesh, e, bases);
1173 missing_bases = missing_bases || b.is_complete();
1175 }
while (missing_bases);
1177 for (
int e = 0; e < n_els; ++e)
1183 setup_data_for_polygons(mesh, e, b, poly_face_to_data);
1204 for (
auto &k : poly_face_to_data)
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));
std::array< Matrix< int, 3, 3 >, 3 > space_
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...
Stores the basis functions for a given element in a mesh (facet in 2d, cell in 3d).
void set_quadrature(const QuadratureFunction &fun)
void set_local_node_from_primitive_func(LocalNodeFromPrimitiveFunc fun)
sets mapping from local nodes to global nodes
std::vector< Basis > bases
one basis function per node in the element
void set_mass_quadrature(const QuadratureFunction &fun)
static Eigen::VectorXi hex_face_local_nodes(const bool serendipity, const int q, const mesh::Mesh3D &mesh, mesh::Navigation3D::Index index)
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)
static void fit_nodes(const mesh::Mesh3D &mesh, const int n_bases, std::vector< ElementBases > &gbases)
bool is_volume() const override
checks if mesh is volume
Abstract mesh class to capture 2d/3d conforming and non-conforming meshes.
int n_elements() const
utitlity to return the number of elements, cells or faces in 3d and 2d
bool is_polytope(const int el_id) const
checks if element is polygon compatible
bool is_spline_compatible(const int el_id) const
checks if element is spline compatible
void get_quadrature(const int order, Quadrature &quad)
void q_grad_basis_value_3d(const int q, const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val)
void q_nodes_3d(const int q, Eigen::MatrixXd &val)
void q_basis_value_3d(const int q, const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val)
spdlog::logger & logger()
Retrieves the current logger.