7#include <igl/barycentric_coordinates.h>
9#include <geogram/mesh/mesh_io.h>
14 using namespace utils;
27 const Eigen::Vector3d e0 = (v2 - v1).transpose();
28 const Eigen::Vector3d e1 = (v3 - v1).transpose();
30 return e0.cross(e1).norm() / 2;
36 p1.resize(p0.rows(), p0.cols());
38 for (std::size_t e = 0; e <
n_edges(); ++e)
43 p0.row(e) =
point(v0);
44 p1.row(e) =
point(v1);
48 void Mesh3D::get_edges(Eigen::MatrixXd &p0, Eigen::MatrixXd &p1,
const std::vector<bool> &valid_elements)
const
51 for (
size_t i = 0; i < valid_elements.size(); ++i)
53 if (valid_elements[i])
64 for (
size_t i = 0; i < valid_elements.size(); ++i)
66 if (!valid_elements[i])
87 const double t = i / (n_new_nodes + 1.0);
89 return std::make_pair((1 - t) * v1 + t * v2, -1);
94 return std::make_pair(n.nodes.row(i - 1), n.nodes_ids[i - 1]);
96 return std::make_pair(n.nodes.row(n.nodes.rows() - i), n.nodes_ids[n.nodes_ids.size() - i]);
101 const int tmp = n_new_nodes == 2 ? 3 : n_new_nodes;
117 const double b2 = i / (n_new_nodes + 2.0);
118 const double b3 = j / (n_new_nodes + 2.0);
119 const double b1 = 1 - b3 - b2;
123 return std::make_pair(b1 * v1 + b2 * v2 + b3 * v3, -1);
126 const int ii = i - 1;
127 const int jj = j - 1;
131 int lindex = jj * n_new_nodes + ii;
135 static const std::array<int, 3> remapping = {{0, 2, 1}};
140 lindex = remapping[lindex];
148 else if (n.v2 == index.
vertex)
153 lindex = remapping[lindex];
161 lindex = (lindex + 1) % 3;
163 else if (n.v3 == index.
vertex)
168 lindex = remapping[lindex];
176 lindex = (lindex + 2) % 3;
184 return std::make_pair(n.nodes.row(lindex), n.nodes_ids[lindex]);
186 else if (
is_cube(index.
element) || is_prism_quad || is_pyramid_quad)
196 const double b1 = i / (n_new_nodes + 1.0);
197 const double b2 = j / (n_new_nodes + 1.0);
199 return std::make_pair(v1 * (1 - b1) * (1 - b2) + v2 * b1 * (1 - b2) + v3 * b1 * b2 + v4 * (1 - b1) * b2, -1);
210 assert(n_new_nodes == 1);
212 assert(n.nodes.rows() == 1);
213 return std::make_pair(n.nodes, n.nodes_ids[0]);
216 if (n_new_nodes == 1)
221 if (n_new_nodes == 1)
230 const double w1 = double(i) / (n_new_nodes + 3);
231 const double w2 = double(j) / (n_new_nodes + 3);
232 const double w3 = double(k) / (n_new_nodes + 3);
233 const double w4 = 1 - w1 - w2 - w3;
240 return std::make_pair(w4 * v1 + w1 * v2 + w2 * v3 + w3 * v4, -1);
259 const double b1 = i / (n_new_nodes + 1.0);
260 const double b2 = j / (n_new_nodes + 1.0);
262 const double b3 = k / (n_new_nodes + 1.0);
264 RowVectorNd blin1 = v1 * (1 - b1) * (1 - b2) + v2 * b1 * (1 - b2) + v3 * b1 * b2 + v4 * (1 - b1) * b2;
265 RowVectorNd blin2 = v5 * (1 - b1) * (1 - b2) + v6 * b1 * (1 - b2) + v7 * b1 * b2 + v8 * (1 - b1) * b2;
267 return std::make_pair((1 - b3) * blin1 + b3 * blin2, -1);
272 const double n = n_new_nodes;
273 const double zv = 1.0 - k / n;
278 xv = 0.5 * (1.0 - zv);
279 yv = 0.5 * (1.0 - zv);
283 const double hk = 1.0 / (k + 1);
284 xv = i * hk * (1.0 - zv);
285 yv = j * hk * (1.0 - zv);
287 return std::make_pair(
RowVectorNd{{xv, yv, zv}}, -1);
351 assert(
is_cube(element_index));
353 std::array<int, 8> v;
357 for (
int i = 0; i < 8; ++i)
358 v[i] = to_vertex[i](idx).vertex;
376 std::array<int, 4> v;
378 for (
int lv = 0; lv < 3; ++lv)
400 std::array<int, 6> v;
404 for (
int i = 0; i < 3; ++i)
407 if (idx.vertex != start.
vertex)
410 for (
int i = 0; i < 3; ++i)
415 assert(start.
vertex == v[0]);
418 for (
int i = 0; i < 3; ++i)
423 assert(start.
vertex == v[3]);
430 std::array<int, 5> v;
441 for (
int k = 0; k < 3; ++k)
456 for (
int lv = 0; lv < 4; ++lv)
474 auto &box = boxes[i];
475 box[0].setConstant(std::numeric_limits<double>::max());
476 box[1].setConstant(std::numeric_limits<double>::min());
481 for (
int d = 0; d < 3; ++d)
483 box[0][d] = std::min(box[0][d],
point(v_id)[d]);
484 box[1][d] = std::max(box[1][d],
point(v_id)[d]);
496 const auto A =
point(indices[0]);
497 const auto B =
point(indices[1]);
498 const auto C =
point(indices[2]);
499 const auto D =
point(indices[3]);
501 igl::barycentric_coordinates(p, A, B, C, D, coord);
510 const auto A =
point(indices[0]);
511 const auto B =
point(indices[1]);
512 const auto C =
point(indices[2]);
513 const auto D =
point(indices[3]);
515 Eigen::MatrixXd coords(4, 4);
516 coords << A, 1, B, 1, C, 1, D, 1;
517 coords.transposeInPlace();
519 jacobian = coords * reference_map;
521 assert(jacobian.determinant() > 0);
void to_vertex_functions(std::array< std::function< Navigation3D::Index(Navigation3D::Index)>, 8 > &to_vertex) const
virtual Navigation3D::Index get_index_from_element(int hi, int lf, int lv) const =0
double tri_area(const int gid) const override
area of a tri face of a tet mesh
void compute_cell_jacobian(const int el_id, const Eigen::MatrixXd &reference_map, Eigen::MatrixXd &jacobian) const
virtual int n_cell_edges(const int c_id) const =0
virtual int cell_edge(const int c_id, const int le_id) const =0
std::array< int, 5 > get_ordered_vertices_from_pyramid(const int element_index) const
void elements_boxes(std::vector< std::array< Eigen::Vector3d, 2 > > &boxes) const override
constructs a box around every element (3d cell, 2d face)
std::array< int, 8 > get_ordered_vertices_from_hex(const int element_index) const
std::pair< RowVectorNd, int > cell_node(const Navigation3D::Index &index, const int n_new_nodes, const int i, const int j, const int k) const
void get_edges(Eigen::MatrixXd &p0, Eigen::MatrixXd &p1) const override
Get all the edges.
virtual Navigation3D::Index switch_edge(Navigation3D::Index idx) const =0
void barycentric_coords(const RowVectorNd &p, const int el_id, Eigen::MatrixXd &coord) const override
constructs barycentric coodiantes for a point p.
std::pair< RowVectorNd, int > face_node(const Navigation3D::Index &index, const int n_new_nodes, const int i, const int j) const
std::pair< RowVectorNd, int > edge_node(const Navigation3D::Index &index, const int n_new_nodes, const int i) const
virtual std::array< int, 4 > get_ordered_vertices_from_tet(const int element_index) const
void to_edge_functions(std::array< std::function< Navigation3D::Index(Navigation3D::Index)>, 12 > &to_edge) const
std::array< int, 6 > get_ordered_vertices_from_prism(const int element_index) const
virtual Navigation3D::Index next_around_face(Navigation3D::Index idx) const =0
virtual Navigation3D::Index switch_face(Navigation3D::Index idx) const =0
void to_face_functions(std::array< std::function< Navigation3D::Index(Navigation3D::Index)>, 6 > &to_face) const
virtual Navigation3D::Index switch_vertex(Navigation3D::Index idx) const =0
int n_elements() const
utitlity to return the number of elements, cells or faces in 3d and 2d
virtual int n_vertices() const =0
number of vertices
Eigen::MatrixXi orders_
list of geometry orders, one per cell
bool is_cube(const int el_id) const
checks if element is cube compatible
virtual RowVectorNd point(const int global_index) const =0
point coordinates
bool is_simplex(const int el_id) const
checks if element is simplex
bool is_prism(const int el_id) const
checks if element is a prism
std::vector< CellNodes > cell_nodes_
high-order nodes associates to cells
virtual int edge_vertex(const int e_id, const int lv_id) const =0
id of the edge vertex
std::vector< FaceNodes > face_nodes_
high-order nodes associates to faces
std::vector< EdgeNodes > edge_nodes_
high-order nodes associates to edges
bool is_pyramid(const int el_id) const
checks if element is a pyramid
virtual RowVectorNd cell_barycenter(const int c) const =0
cell barycenter
virtual int n_edges() const =0
number of edges
virtual int cell_vertex(const int f_id, const int lv_id) const =0
id of the vertex of a cell
virtual int n_face_vertices(const int f_id) const =0
number of vertices of a face
virtual int n_cell_vertices(const int c_id) const =0
number of vertices of a cell
virtual int face_vertex(const int f_id, const int lv_id) const =0
id of the face vertex
Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor, 1, 3 > RowVectorNd