11#include <polysolve/linear/Solver.hpp>
18#include <Eigen/Sparse>
31 using namespace Eigen;
32 using namespace assembler;
40 typedef Matrix<std::vector<int>, 3, 3> SpaceMatrix;
41 typedef Matrix<RowVectorNd, 3, 3> NodeMatrix;
43 void print_local_space(
const SpaceMatrix &space)
47 for (
int j = 2; j >= 0; --j)
49 for (
int i = 0; i < 3; ++i)
51 if (space(i, j).size() > 0)
53 for (std::size_t l = 0; l < space(i, j).size(); ++l)
54 ss << space(i, j)[l] <<
",";
64 logger().trace(
"Local space:\n{}", ss.str());
67 int node_id_from_edge_index(
const Mesh2D &mesh, MeshNodes &mesh_nodes,
const Navigation::Index &index)
69 const int face_id = mesh.switch_face(index).face;
70 if (face_id >= 0 && mesh.is_cube(face_id))
71 return mesh_nodes.node_id_from_face(face_id);
73 return mesh_nodes.node_id_from_edge(index.edge);
76 void explore_direction(
const Navigation::Index &index,
const Mesh2D &mesh, MeshNodes &mesh_nodes,
const int x,
const int y,
const bool is_x,
const bool invert, LocalBoundary &lb, SpaceMatrix &space, NodeMatrix &node, std::map<int, InterfaceData> &poly_edge_to_data)
78 int node_id = node_id_from_edge_index(mesh, mesh_nodes, index);
81 assert(std::find(space(
x,
y).begin(), space(
x,
y).end(), node_id) == space(
x,
y).end());
82 space(
x,
y).push_back(node_id);
83 node(
x,
y) = mesh_nodes.node_position(node_id);
84 assert(node(
x,
y).size() == 2);
86 const int x1 = is_x ?
x : (invert ? 2 : 0);
87 const int y1 = !is_x ?
y : (invert ? 2 : 0);
89 const int x2 = is_x ?
x : (invert ? 0 : 2);
90 const int y2 = !is_x ?
y : (invert ? 0 : 2);
92 const bool is_boundary = mesh_nodes.is_boundary(node_id);
93 const bool is_interface = mesh_nodes.is_interface(node_id);
101 else if (is_interface)
103 InterfaceData &data = poly_edge_to_data[index.edge];
104 data.local_indices.push_back(
y * 3 +
x);
108 assert(!is_boundary && !is_interface);
110 Navigation::Index start_index = mesh.switch_face(index);
111 assert(start_index.vertex == index.vertex);
112 assert(start_index.face >= 0);
114 Navigation::Index edge1 = mesh.switch_edge(start_index);
115 node_id = node_id_from_edge_index(mesh, mesh_nodes, edge1);
119 if (std::find(space(x1, y1).begin(), space(x1, y1).end(), node_id) == space(x1, y1).end())
121 space(x1, y1).push_back(node_id);
122 node(x1, y1) = mesh_nodes.node_position(node_id);
125 Navigation::Index edge2 = mesh.switch_edge(mesh.switch_vertex(start_index));
126 node_id = node_id_from_edge_index(mesh, mesh_nodes, edge2);
129 if (std::find(space(x2, y2).begin(), space(x2, y2).end(), node_id) == space(x2, y2).end())
131 space(x2, y2).push_back(node_id);
132 node(x2, y2) = mesh_nodes.node_position(node_id);
138 void add_id_for_poly(
const Navigation::Index &index,
const int x1,
const int y1,
const int x2,
const int y2,
const SpaceMatrix &space, std::map<int, InterfaceData> &poly_edge_to_data)
140 auto it = poly_edge_to_data.find(index.edge);
141 if (it != poly_edge_to_data.end())
143 InterfaceData &data = it->second;
145 assert(space(x1, y1).size() == 1);
146 data.local_indices.push_back(y1 * 3 + x1);
148 assert(space(x2, y2).size() == 1);
149 data.local_indices.push_back(y2 * 3 + x2);
153 void build_local_space(
const Mesh2D &mesh, MeshNodes &mesh_nodes,
const int el_index, SpaceMatrix &space, NodeMatrix &node, std::vector<LocalBoundary> &local_boundary, std::map<int, InterfaceData> &poly_edge_to_data)
155 assert(!mesh.is_volume());
157 Navigation::Index index;
160 const int el_node_id = mesh_nodes.node_id_from_face(el_index);
161 space(1, 1).push_back(el_node_id);
162 node(1, 1) = mesh_nodes.node_position(el_node_id);
165 LocalBoundary lb(el_index, BoundaryType::QUAD_LINE);
168 index = mesh.get_index_from_face(el_index);
169 explore_direction(index, mesh, mesh_nodes, 1, 0,
false,
false, lb, space, node, poly_edge_to_data);
172 index = mesh.next_around_face(index);
173 explore_direction(index, mesh, mesh_nodes, 2, 1,
true,
false, lb, space, node, poly_edge_to_data);
176 index = mesh.next_around_face(index);
177 explore_direction(index, mesh, mesh_nodes, 1, 2,
false,
true, lb, space, node, poly_edge_to_data);
180 index = mesh.next_around_face(index);
181 explore_direction(index, mesh, mesh_nodes, 0, 1,
true,
true, lb, space, node, poly_edge_to_data);
184 if (mesh_nodes.is_boundary_or_interface(space(1, 2).front()) && mesh_nodes.is_boundary_or_interface(space(2, 1).front()))
186 assert(space(2, 2).empty());
188 Navigation::Index start_index = mesh.get_index_from_face(el_index);
189 start_index = mesh.next_around_face(start_index);
190 start_index = mesh.next_around_face(start_index);
192 const int node_id = mesh_nodes.node_id_from_vertex(start_index.vertex);
194 space(2, 2).push_back(node_id);
195 node(2, 2) = mesh_nodes.node_position(node_id);
201 if (mesh_nodes.is_boundary_or_interface(space(1, 0).front()) && mesh_nodes.is_boundary_or_interface(space(2, 1).front()))
203 assert(space(2, 0).empty());
205 Navigation::Index start_index = mesh.get_index_from_face(el_index);
206 start_index = mesh.next_around_face(start_index);
208 const int node_id = mesh_nodes.node_id_from_vertex(start_index.vertex);
210 space(2, 0).push_back(node_id);
211 node(2, 0) = mesh_nodes.node_position(node_id);
217 if (mesh_nodes.is_boundary_or_interface(space(1, 2).front()) && mesh_nodes.is_boundary_or_interface(space(0, 1).front()))
219 assert(space(0, 2).empty());
221 Navigation::Index start_index = mesh.get_index_from_face(el_index);
222 start_index = mesh.next_around_face(start_index);
223 start_index = mesh.next_around_face(start_index);
224 start_index = mesh.next_around_face(start_index);
227 const int node_id = mesh_nodes.node_id_from_vertex(start_index.vertex);
228 space(0, 2).push_back(node_id);
229 node(0, 2) = mesh_nodes.node_position(node_id);
235 if (mesh_nodes.is_boundary_or_interface(space(1, 0).front()) && mesh_nodes.is_boundary_or_interface(space(0, 1).front()))
237 Navigation::Index start_index = mesh.get_index_from_face(el_index);
240 const int node_id = mesh_nodes.node_id_from_vertex(start_index.vertex);
241 space(0, 0).push_back(node_id);
242 node(0, 0) = mesh_nodes.node_position(node_id);
251 index = mesh.get_index_from_face(el_index);
252 add_id_for_poly(index, 0, 0, 2, 0, space, poly_edge_to_data);
254 index = mesh.next_around_face(index);
255 add_id_for_poly(index, 2, 0, 2, 2, space, poly_edge_to_data);
257 index = mesh.next_around_face(index);
258 add_id_for_poly(index, 2, 2, 0, 2, space, poly_edge_to_data);
260 index = mesh.next_around_face(index);
261 add_id_for_poly(index, 0, 2, 0, 0, space, poly_edge_to_data);
264 local_boundary.emplace_back(lb);
267 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)
270 if (mesh_nodes.is_boundary_or_interface(space(0, 1).front()) && mesh_nodes.is_boundary_or_interface(space(2, 1).front()))
272 h_knots[0] = {{0, 0, 0, 1}};
273 h_knots[1] = {{0, 0, 1, 1}};
274 h_knots[2] = {{0, 1, 1, 1}};
277 else if (mesh_nodes.is_boundary_or_interface(space(0, 1).front()))
279 h_knots[0] = {{0, 0, 0, 1}};
280 h_knots[1] = {{0, 0, 1, 2}};
281 h_knots[2] = {{0, 1, 2, 3}};
284 else if (mesh_nodes.is_boundary_or_interface(space(2, 1).front()))
286 h_knots[0] = {{-2, -1, 0, 1}};
287 h_knots[1] = {{-1, 0, 1, 1}};
288 h_knots[2] = {{0, 1, 1, 1}};
292 h_knots[0] = {{-2, -1, 0, 1}};
293 h_knots[1] = {{-1, 0, 1, 2}};
294 h_knots[2] = {{0, 1, 2, 3}};
298 if (mesh_nodes.is_boundary_or_interface(space(1, 0).front()) && mesh_nodes.is_boundary_or_interface(space(1, 2).front()))
300 v_knots[0] = {{0, 0, 0, 1}};
301 v_knots[1] = {{0, 0, 1, 1}};
302 v_knots[2] = {{0, 1, 1, 1}};
305 else if (mesh_nodes.is_boundary_or_interface(space(1, 0).front()))
307 v_knots[0] = {{0, 0, 0, 1}};
308 v_knots[1] = {{0, 0, 1, 2}};
309 v_knots[2] = {{0, 1, 2, 3}};
312 else if (mesh_nodes.is_boundary_or_interface(space(1, 2).front()))
314 v_knots[0] = {{-2, -1, 0, 1}};
315 v_knots[1] = {{-1, 0, 1, 1}};
316 v_knots[2] = {{0, 1, 1, 1}};
320 v_knots[0] = {{-2, -1, 0, 1}};
321 v_knots[1] = {{-1, 0, 1, 2}};
322 v_knots[2] = {{0, 1, 2, 3}};
326 void basis_for_regular_quad(
const SpaceMatrix &space,
const NodeMatrix &loc_nodes,
const std::array<std::array<double, 4>, 3> &h_knots,
const std::array<std::array<double, 4>, 3> &v_knots, ElementBases &b)
328 for (
int y = 0;
y < 3; ++
y)
330 for (
int x = 0;
x < 3; ++
x)
332 if (space(
x,
y).size() == 1)
334 const int global_index = space(
x,
y).front();
335 const Eigen::MatrixXd &node = loc_nodes(
x,
y);
336 assert(node.size() == 2);
338 const int local_index =
y * 3 +
x;
339 b.bases[local_index].init(2, global_index, local_index, node);
341 const QuadraticBSpline2d spline(h_knots[
x], v_knots[
y]);
342 b.bases[local_index].set_basis([spline](
const Eigen::MatrixXd &uv, Eigen::MatrixXd &
val) { spline.interpolate(uv,
val); });
343 b.bases[local_index].set_grad([spline](
const Eigen::MatrixXd &uv, Eigen::MatrixXd &
val) { spline.derivative(uv,
val); });
349 void basis_for_irregulard_quad(
const int el_id,
const Mesh2D &mesh, MeshNodes &mesh_nodes,
const SpaceMatrix &space,
const NodeMatrix &loc_nodes,
const std::array<std::array<double, 4>, 3> &h_knots,
const std::array<std::array<double, 4>, 3> &v_knots, ElementBases &b)
351 for (
int y = 0;
y < 3; ++
y)
353 for (
int x = 0;
x < 3; ++
x)
355 if (space(
x,
y).size() > 1)
363 std::vector<int> other_indices;
364 const auto ¢er = b.bases[1 * 3 + 1].global().front();
366 const auto &el1 = b.bases[mpy * 3 + mpx].global().front();
367 const auto &el2 = b.bases[mmy * 3 + mmx].global().front();
369 Navigation::Index start_index = mesh.get_index_from_face(el_id);
371 for (
int i = 0; i < 4; ++i)
373 other_indices.clear();
375 Navigation::Index index = start_index;
378 const int f_index = mesh_nodes.node_id_from_face(index.face);
379 if (f_index != el1.index && f_index != el2.index && f_index != center.index)
380 other_indices.push_back(f_index);
383 index = mesh.next_around_vertex(index);
384 }
while (index.face != start_index.face);
391 start_index = mesh.next_around_face(start_index);
395 const int local_index =
y * 3 +
x;
396 auto &
base = b.bases[local_index];
398 const int k = int(other_indices.size()) + 3;
400 base.global().resize(k);
402 base.global()[0].index = center.index;
403 base.global()[0].val = (4. - k) / k;
404 base.global()[0].node = center.node;
406 base.global()[1].index = el1.index;
407 base.global()[1].val = (4. - k) / k;
408 base.global()[1].node = el1.node;
410 base.global()[2].index = el2.index;
411 base.global()[2].val = (4. - k) / k;
412 base.global()[2].node = el2.node;
414 for (std::size_t n = 0; n < other_indices.size(); ++n)
416 base.global()[3 + n].index = other_indices[n];
417 base.global()[3 + n].val = 4. / k;
418 base.global()[3 + n].node = mesh_nodes.node_position(other_indices[n]);
421 const QuadraticBSpline2d spline(h_knots[
x], v_knots[
y]);
422 b.bases[local_index].set_basis([spline](
const Eigen::MatrixXd &uv, Eigen::MatrixXd &
val) { spline.interpolate(uv,
val); });
423 b.bases[local_index].set_grad([spline](
const Eigen::MatrixXd &uv, Eigen::MatrixXd &
val) { spline.derivative(uv,
val); });
429 void create_q2_nodes(
const Mesh2D &mesh,
const int el_index, std::set<int> &vertex_id, std::set<int> &
edge_id, ElementBases &b, std::vector<LocalBoundary> &local_boundary,
int &n_bases)
433 LocalBoundary lb(el_index, BoundaryType::QUAD_LINE);
435 Navigation::Index index = mesh.get_index_from_face(el_index);
436 for (
int j = 0; j < 4; ++j)
438 int current_vertex_node_id = -1;
439 int current_edge_node_id = -1;
440 Eigen::Matrix<double, 1, 2> current_edge_node;
441 Eigen::MatrixXd current_vertex_node;
446 int vertex_basis_id = e2l[0];
447 int edge_basis_id = e2l[1];
449 const int opposite_face = mesh.switch_face(index).face;
452 bool is_vertex_q2 =
true;
453 bool is_vertex_boundary =
false;
455 Navigation::Index vindex = index;
461 is_vertex_boundary =
true;
464 if (mesh.is_spline_compatible(vindex.face))
466 is_vertex_q2 =
false;
469 vindex = mesh.next_around_vertex(vindex);
470 }
while (vindex.edge != index.edge);
474 vindex = mesh.switch_face(index);
479 is_vertex_boundary =
true;
483 if (mesh.is_spline_compatible(vindex.face))
485 is_vertex_q2 =
false;
488 vindex = mesh.next_around_vertex(vindex);
489 }
while (vindex.edge != index.edge);
492 const bool is_edge_q2 = opposite_face < 0 || !mesh.is_spline_compatible(opposite_face);
496 const bool is_new_edge =
edge_id.insert(index.edge).second;
500 current_edge_node_id = n_bases++;
501 current_edge_node = mesh.edge_barycenter(index.edge);
503 if (opposite_face < 0)
506 lb.add_boundary_primitive(index.edge, edge_basis_id - 4);
514 const bool is_new_vertex = vertex_id.insert(index.vertex).second;
518 current_vertex_node_id = n_bases++;
519 current_vertex_node = mesh.point(index.vertex);
527 if (current_vertex_node_id >= 0)
528 b.bases[vertex_basis_id].init(2, current_vertex_node_id, vertex_basis_id, current_vertex_node);
530 if (current_edge_node_id >= 0)
531 b.bases[edge_basis_id].init(2, current_edge_node_id, edge_basis_id, current_edge_node);
540 index = mesh.next_around_face(index);
544 const int face_basis_id = 8;
545 b.bases[face_basis_id].init(2, n_bases++, face_basis_id, mesh.face_barycenter(el_index));
550 local_boundary.emplace_back(lb);
553 void insert_into_global(
const Local2Global &data, std::vector<Local2Global> &
vec)
556 if (fabs(data.val) < 1e-10)
561 for (std::size_t i = 0; i <
vec.size(); ++i)
563 if (
vec[i].index == data.index)
566 assert(fabs(
vec[i].
val - data.val) < 1e-10);
567 assert((
vec[i].node - data.node).norm() < 1e-10);
577 void assign_q2_weights(
const Mesh2D &mesh,
const int el_index, std::vector<ElementBases> &bases)
580 std::vector<AssemblyValues> eval_p;
581 Navigation::Index index = mesh.get_index_from_face(el_index);
582 ElementBases &b = bases[el_index];
584 for (
int j = 0; j < 4; ++j)
586 const int opposite_face = mesh.switch_face(index).face;
588 if (opposite_face < 0 || !mesh.is_cube(opposite_face))
590 index = mesh.next_around_face(index);
598 Eigen::Matrix<double, 3, 2> param_p;
601 Eigen::MatrixXd quad_loc_nodes;
604 for (
int k = 0; k < 3; ++k)
605 param_p.row(k) = quad_loc_nodes.row(opposite_indices[k]);
612 const auto &other_bases = bases[opposite_face];
613 other_bases.evaluate_bases(param_p, eval_p);
615 for (std::size_t i = 0; i < other_bases.bases.size(); ++i)
617 const auto &other_b = other_bases.bases[i];
619 if (other_b.global().empty())
622 assert(eval_p[i].
val.size() == 3);
625 if (eval_p[i].
val.cwiseAbs().maxCoeff() <= 1e-10)
628 for (std::size_t k = 0; k < other_b.global().size(); ++k)
634 auto glob0 = other_b.global()[k];
635 glob0.val *= eval_p[i].val(0);
636 auto glob1 = other_b.global()[k];
637 glob1.val *= eval_p[i].val(1);
638 auto glob2 = other_b.global()[k];
639 glob2.val *= eval_p[i].val(2);
641 insert_into_global(glob0, b.bases[i0].global());
642 insert_into_global(glob1, b.bases[i1].global());
643 insert_into_global(glob2, b.bases[i2].global());
647 index = mesh.next_around_face(index);
651 void setup_data_for_polygons(
const Mesh2D &mesh,
const int el_index,
const ElementBases &b, std::map<int, InterfaceData> &poly_edge_to_data)
653 Navigation::Index index = mesh.get_index_from_face(el_index);
654 for (
int j = 0; j < 4; ++j)
656 const int opposite_face = mesh.switch_face(index).face;
657 const bool is_neigh_poly = opposite_face >= 0 && mesh.is_polytope(opposite_face);
663 const int vertex_basis_id = e2l[0];
664 const int edge_basis_id = e2l[1];
665 const int vertex_basis_id2 = e2l[2];
667 InterfaceData &data = poly_edge_to_data[index.edge];
669 data.local_indices.push_back(edge_basis_id);
670 data.local_indices.push_back(vertex_basis_id);
671 data.local_indices.push_back(vertex_basis_id2);
674 index = mesh.next_around_face(index);
680 const std::string &assembler,
681 const int quadrature_order,
const int mass_quadrature_order, std::vector<ElementBases> &bases, std::vector<LocalBoundary> &local_boundary, std::map<int, InterfaceData> &poly_edge_to_data)
686 MeshNodes mesh_nodes(mesh,
true,
true, 1, 1);
691 local_boundary.clear();
695 for (
int e = 0; e < n_els; ++e)
701 NodeMatrix loc_nodes;
704 build_local_space(mesh, mesh_nodes, e, space, loc_nodes, local_boundary, poly_edge_to_data);
710 const int real_mass_order = mass_quadrature_order > 0 ? mass_quadrature_order :
AssemblerUtils::quadrature_order(
"Mass", 2, AssemblerUtils::BasisType::SPLINE, 2);
723 Eigen::VectorXi res(3);
724 const auto &mesh2d =
dynamic_cast<const Mesh2D &
>(mesh);
727 for (le = 0; le < mesh2d.n_face_vertices(e); ++le)
729 if (index.edge == primitive_id)
731 index = mesh2d.next_around_face(index);
733 assert(index.edge == primitive_id);
738 res << (3 * 0 + 0), (3 * 1 + 0), (3 * 2 + 0);
741 res << (3 * 0 + 0), (3 * 0 + 1), (3 * 0 + 2);
744 res << (3 * 0 + 2), (3 * 1 + 2), (3 * 2 + 2);
747 res << (3 * 2 + 0), (3 * 2 + 1), (3 * 2 + 2);
756 std::array<std::array<double, 4>, 3> h_knots;
757 std::array<std::array<double, 4>, 3> v_knots;
759 setup_knots_vectors(mesh_nodes, space, h_knots, v_knots);
763 basis_for_regular_quad(space, loc_nodes, h_knots, v_knots, b);
764 basis_for_irregulard_quad(e, mesh, mesh_nodes, space, loc_nodes, h_knots, v_knots, b);
768 std::set<int> vertex_id;
770 int n_bases = mesh_nodes.
n_nodes();
772 for (
int e = 0; e < n_els; ++e)
779 const int real_order = quadrature_order > 0 ? quadrature_order :
AssemblerUtils::quadrature_order(assembler, 2, AssemblerUtils::BasisType::CUBE_LAGRANGE, 2);
780 const int real_mass_order = mass_quadrature_order > 0 ? mass_quadrature_order :
AssemblerUtils::quadrature_order(
"Mass", 2, AssemblerUtils::BasisType::CUBE_LAGRANGE, 2);
792 const auto &mesh2d =
dynamic_cast<const Mesh2D &
>(mesh);
795 for (
int le = 0; le < mesh2d.n_face_vertices(e); ++le)
797 if (index.edge == primitive_id)
799 index = mesh2d.next_around_face(index);
801 assert(index.edge == primitive_id);
805 Eigen::VectorXi res(indices.size());
807 for (
size_t i = 0; i < indices.size(); ++i)
813 create_q2_nodes(mesh, e, vertex_id,
edge_id, b, local_boundary, n_bases);
816 bool missing_bases =
false;
819 missing_bases =
false;
820 for (
int e = 0; e < n_els; ++e)
829 assign_q2_weights(mesh, e, bases);
831 missing_bases = missing_bases || b.is_complete();
833 }
while (missing_bases);
835 for (
int e = 0; e < n_els; ++e)
840 setup_data_for_polygons(mesh, e, b, poly_edge_to_data);
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 quad_edge_local_nodes(const int q, const mesh::Mesh2D &mesh, mesh::Navigation::Index index)
static void fit_nodes(const mesh::Mesh2D &mesh, const int n_bases, std::vector< ElementBases > &gbases)
static int build_bases(const mesh::Mesh2D &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_edge_to_data)
bool is_volume() const override
checks if mesh is volume
virtual Navigation::Index get_index_from_face(int f, int lv=0) const =0
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_nodes_2d(const int q, Eigen::MatrixXd &val)
void q_grad_basis_value_2d(const int q, const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val)
void q_basis_value_2d(const int q, const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val)
spdlog::logger & logger()
Retrieves the current logger.