10#include <geogram/basic/file_system.h>
11#include <geogram/mesh/mesh_io.h>
12#include <geogram/mesh/mesh_geometry.h>
13#include <geogram/mesh/mesh_repair.h>
21 using namespace utils;
28 if (n_refinement <= 0)
35 bool all_simplicial =
true;
41 for (
int i = 0; i < n_refinement; ++i)
50 mesh_.clear(
false,
false);
67 c2e_ = std::make_unique<GEO::Attribute<GEO::index_t>>(
mesh_.facet_corners.attributes(),
"edge_id");
69 boundary_edges_ = std::make_unique<GEO::Attribute<bool>>(
mesh_.edges.attributes(),
"boundary_edge");
83 for (
int e = 0; e < (int)
mesh_.edges.nb(); ++e)
85 for (
int lv = 0; lv < 2; ++lv)
152 mesh_.clear(
false,
false);
157 c2e_ = std::make_unique<GEO::Attribute<GEO::index_t>>(
mesh_.facet_corners.attributes(),
"edge_id");
159 boundary_edges_ = std::make_unique<GEO::Attribute<bool>>(
mesh_.edges.attributes(),
"boundary_edge");
167 if (!mesh_save(
mesh_, path))
183 mesh_.clear(
false,
false);
189 c2e_ = std::make_unique<GEO::Attribute<GEO::index_t>>(
mesh_.facet_corners.attributes(),
"edge_id");
191 boundary_edges_ = std::make_unique<GEO::Attribute<bool>>(
mesh_.edges.attributes(),
"boundary_edge");
208 assert(nodes.size() ==
n_faces());
210 for (
int f = 0; f <
n_faces(); ++f)
214 const auto &nodes_ids = nodes[f];
216 if (nodes_ids.size() == 3)
222 else if (nodes_ids.size() == 6)
226 for (
int le = 0; le < 3; ++le)
231 if (n.nodes.size() <= 0)
237 if ((n.v1 == nodes_ids[0] && n.v2 == nodes_ids[1]) || (n.v2 == nodes_ids[0] && n.v1 == nodes_ids[1]))
239 else if ((n.v1 == nodes_ids[1] && n.v2 == nodes_ids[2]) || (n.v2 == nodes_ids[1] && n.v1 == nodes_ids[2]))
244 n.nodes.resize(1, 2);
245 n.nodes <<
V(nodes_ids[node_index], 0),
V(nodes_ids[node_index], 1);
251 else if (nodes_ids.size() == 10)
255 for (
int le = 0; le < 3; ++le)
260 if (n.nodes.size() <= 0)
267 if (n.v1 == nodes_ids[0] && n.v2 == nodes_ids[1])
272 else if (n.v2 == nodes_ids[0] && n.v1 == nodes_ids[1])
277 else if (n.v1 == nodes_ids[1] && n.v2 == nodes_ids[2])
282 else if (n.v2 == nodes_ids[1] && n.v1 == nodes_ids[2])
287 else if (n.v1 == nodes_ids[2] && n.v2 == nodes_ids[0])
294 assert(n.v2 == nodes_ids[2] && n.v1 == nodes_ids[0]);
299 n.nodes.resize(2, 2);
300 n.nodes.row(0) <<
V(nodes_ids[node_index1], 0),
V(nodes_ids[node_index1], 1);
301 n.nodes.row(1) <<
V(nodes_ids[node_index2], 0),
V(nodes_ids[node_index2], 1);
308 n.v1 =
mesh_.facets.vertex(f, 0);
309 n.v2 =
mesh_.facets.vertex(f, 1);
310 n.v3 =
mesh_.facets.vertex(f, 2);
311 n.nodes.resize(1, 2);
312 n.nodes <<
V(nodes_ids[9], 0),
V(nodes_ids[9], 1);
316 else if (nodes_ids.size() == 15)
340 const double t = i / (n_new_nodes + 1.0);
342 return (1 - t) * v1 + t * v2;
347 return n.nodes.row(i - 1);
350 assert(n.v2 == index.
vertex);
351 return n.nodes.row(n.nodes.rows() - i);
365 const double b2 = i / (n_new_nodes + 2.0);
366 const double b3 = j / (n_new_nodes + 2.0);
367 const double b1 = 1 - b3 - b2;
371 return b1 * v1 + b2 * v2 + b3 * v3;
377 return n.nodes.row(0);
389 const double b1 = i / (n_new_nodes + 1.0);
390 const double b2 = j / (n_new_nodes + 1.0);
392 return v1 * (1 - b1) * (1 - b2) + v2 * b1 * (1 - b2) + v3 * b1 * b2 + v4 * (1 - b1) * b2;
401 GEO::vec3 min_corner, max_corner;
402 GEO::get_bbox(
mesh_, &min_corner[0], &max_corner[0]);
406 min(0) = min_corner.x;
407 min(1) = min_corner.y;
409 max(0) = max_corner.x;
410 max(1) = max_corner.y;
416 GEO::vec3 min_corner, max_corner;
417 GEO::get_bbox(
mesh_, &min_corner[0], &max_corner[0]);
418 GEO::vec3 extent = max_corner - min_corner;
419 double scaling = std::max(extent[0], std::max(extent[1], extent[2]));
421 const GEO::vec3 origin = min_corner;
422 for (GEO::index_t v = 0; v <
mesh_.vertices.nb(); ++v)
424 mesh_.vertices.point(v) = (
mesh_.vertices.point(v) - origin) / scaling;
426 Eigen::RowVector2d shift;
427 shift << origin[0], origin[1];
430 if (n.nodes.size() > 0)
431 n.nodes = (n.nodes.rowwise() - shift) / scaling;
435 if (n.nodes.size() > 0)
436 n.nodes = (n.nodes.rowwise() - shift) / scaling;
439 logger().debug(
"-- bbox before normalization:");
440 logger().debug(
" min : {} {}", min_corner[0], min_corner[1]);
441 logger().debug(
" max : {} {}", max_corner[0], max_corner[1]);
442 logger().debug(
" extent: {} {}", max_corner[0] - min_corner[0], max_corner[1] - min_corner[1]);
443 GEO::get_bbox(
mesh_, &min_corner[0], &max_corner[0]);
444 logger().debug(
"-- bbox after normalization:");
445 logger().debug(
" min : {} {}", min_corner[0], min_corner[1]);
446 logger().debug(
" max : {} {}", max_corner[0], max_corner[1]);
447 logger().debug(
" extent: {} {}", max_corner[0] - min_corner[0], max_corner[1] - min_corner[1]);
449 Eigen::MatrixXd p0, p1, p;
452 logger().debug(
"-- edge length after normalization:");
453 logger().debug(
" min: {}", p.rowwise().norm().minCoeff());
454 logger().debug(
" max: {}", p.rowwise().norm().maxCoeff());
455 logger().debug(
" avg: {}", p.rowwise().norm().mean());
460 const int v0 =
mesh_.edges.vertex(gid, 0);
461 const int v1 =
mesh_.edges.vertex(gid, 1);
468 mesh_.vertices.point(global_index).x = p(0);
469 mesh_.vertices.point(global_index).y = p(1);
474 const double *ptr =
mesh_.vertices.point_ptr(global_index);
500 std::vector<Eigen::MatrixXi> local_tris(
mesh_.facets.nb());
501 std::vector<Eigen::MatrixXd> local_pts(
mesh_.facets.nb());
508 for (GEO::index_t f = 0; f <
mesh_.facets.nb(); ++f)
518 const int vertex =
mesh_.facets.vertex(f, i);
519 const double *pt =
mesh_.vertices.point_ptr(vertex);
520 face_pts(i, 0) = pt[0];
521 face_pts(i, 1) = pt[1];
529 local_tris[f].row(i - 1) << 0, i, i + 1;
532 local_pts[f] = face_pts;
534 total_tris += local_tris[f].rows();
535 total_pts += local_pts[f].rows();
537 ranges.push_back(total_tris);
539 assert(local_pts[f].rows() == face_pts.rows());
542 tris.resize(total_tris, 3);
543 pts.resize(total_pts, 2);
547 for (std::size_t i = 0; i < local_tris.size(); ++i)
549 tris.block(tri_index, 0, local_tris[i].rows(), local_tris[i].cols()) = local_tris[i].array() + pts_index;
550 tri_index += local_tris[i].rows();
552 pts.block(pts_index, 0, local_pts[i].rows(), local_pts[i].cols()) = local_pts[i];
553 pts_index += local_pts[i].rows();
570 const int v0 =
mesh_.edges.vertex(index, 0);
571 const int v1 =
mesh_.edges.vertex(index, 1);
597 for (
int e = 0; e <
n_edges(); ++e)
604 if (fabs(p(0) - min_corner[0]) < eps)
606 else if (fabs(p(1) - min_corner[1]) < eps)
608 else if (fabs(p(0) - max_corner[0]) < eps)
610 else if (fabs(p(1) - max_corner[1]) < eps)
624 for (
int e = 0; e <
n_edges(); ++e)
639 for (
int e = 0; e <
n_edges(); ++e)
651 for (
int e = 0; e <
n_edges(); ++e)
663 for (
int e = 0; e <
n_edges(); ++e)
667 std::sort(vs.begin(), vs.end());
676 for (
int e = 0; e <
n_edges(); ++e)
681 std::sort(vs.begin(), vs.end());
688 assert(
typeid(mesh) ==
typeid(
CMesh2D));
697 for (
int i = n_v; i < (int)
mesh_.vertices.nb(); ++i)
699 GEO::vec3 &p =
mesh_.vertices.point(i);
703 std::vector<GEO::index_t> indices;
704 for (
int i = 0; i < mesh2d.
n_faces(); ++i)
707 for (
int j = 0; j < mesh2d.
mesh_.facets.nb_vertices(i); ++j)
708 indices.push_back(mesh2d.
mesh_.facets.vertex(i, j) + n_v);
710 mesh_.facets.create_polygon(indices.size(), &indices[0]);
720 c2e_ = std::make_unique<GEO::Attribute<GEO::index_t>>(
mesh_.facet_corners.attributes(),
"edge_id");
722 boundary_edges_ = std::make_unique<GEO::Attribute<bool>>(
mesh_.edges.attributes(),
"boundary_edge");
727 std::unique_ptr<CMesh2D> copy_mesh = std::make_unique<CMesh2D>();
728 copy_mesh->load(this->
mesh_);
735 copy_mesh->orders_ = this->
orders_;
virtual RowVectorNd edge_barycenter(const int index) const override
edge barycenter
Navigation::Index switch_edge(Navigation::Index idx) const override
bool load(const std::string &path) override
loads a mesh from the path
void normalize() override
normalize the mesh
void compute_elements_tag() override
compute element types, see ElementType
RowVectorNd edge_node(const Navigation::Index &index, const int n_new_nodes, const int i) const override
int n_faces() const override
number of faces
int n_face_vertices(const int f_id) const override
number of vertices of a face
int n_vertices() const override
number of vertices
void compute_boundary_ids(const double eps) override
computes the selection based on the bbx of the mesh.
int n_edges() const override
number of edges
bool is_boundary_element(const int element_global_id) const override
is cell boundary
double edge_length(const int gid) const override
edge length
std::unique_ptr< GEO::Attribute< bool > > boundary_edges_
std::unique_ptr< Mesh > copy() const override
Create a copy of the mesh.
Navigation::Index switch_vertex(Navigation::Index idx) const override
void set_point(const int global_index, const RowVectorNd &p) override
Set the point.
bool is_boundary_edge(const int edge_global_id) const override
is edge boundary
std::unique_ptr< GEO::Attribute< GEO::index_t > > c2e_
void triangulate_faces(Eigen::MatrixXi &tris, Eigen::MatrixXd &pts, std::vector< int > &ranges) const override
generate a triangular representation of every face
std::unique_ptr< GEO::Attribute< bool > > boundary_vertices_
void compute_body_ids(const std::function< int(const size_t, const RowVectorNd &)> &marker) override
computes boundary selections based on a function
void append(const Mesh &mesh) override
appends a new mesh to the end of this
int edge_vertex(const int e_id, const int lv_id) const override
id of the edge vertex
Navigation::Index get_index_from_face(int f, int lv=0) const override
virtual void bounding_box(RowVectorNd &min, RowVectorNd &max) const override
computes the bbox of the mesh
void attach_higher_order_nodes(const Eigen::MatrixXd &V, const std::vector< std::vector< int > > &nodes) override
attach high order nodes
bool build_from_matrices(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F) override
build a mesh from matrices
void refine(const int n_refinement, const double t) override
refine the mesh
bool save(const std::string &path) const override
virtual RowVectorNd point(const int global_index) const override
point coordinates
RowVectorNd face_node(const Navigation::Index &index, const int n_new_nodes, const int i, const int j) const override
virtual void update_elements_tag() override
Update elements types.
RowVectorNd face_barycenter(const int index) const override
face barycenter
Navigation::Index next_around_face(Navigation::Index idx) const
void get_edges(Eigen::MatrixXd &p0, Eigen::MatrixXd &p1) const override
Get all the edges.
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
Eigen::MatrixXi orders_
list of geometry orders, one per cell
std::vector< ElementType > elements_tag_
list of element types
bool is_rational_
stores if the mesh is rational
bool is_cube(const int el_id) const
checks if element is cube compatible
Eigen::MatrixXi in_ordered_faces_
Order of the input faces, TODO: change to std::vector of Eigen::Vector.
bool is_simplex(const int el_id) const
checks if element is simples compatible
std::vector< int > boundary_ids_
list of surface labels
std::vector< int > node_ids_
list of node labels
std::vector< CellNodes > cell_nodes_
high-order nodes associates to cells
std::vector< std::vector< double > > cell_weights_
weights associates to cells for rational polynomail meshes
std::vector< int > body_ids_
list of volume labels
std::vector< FaceNodes > face_nodes_
high-order nodes associates to faces
std::vector< EdgeNodes > edge_nodes_
high-order nodes associates to edges
Eigen::MatrixXi in_ordered_edges_
Order of the input edges.
virtual void append(const Mesh &mesh)
appends a new mesh to the end of this
Eigen::VectorXi in_ordered_vertices_
Order of the input vertices.
void prepare_mesh(GEO::Mesh &M)
SplitFunction catmul_clark_split_func()
SplitFunction polar_split_func(double t)
Helper function.
void refine_triangle_mesh(const GEO::Mesh &M_in, GEO::Mesh &M_out)
Refine a triangle mesh.
void orient_normals_2d(GEO::Mesh &M)
Orient facets of a 2D mesh so that each connected component has positive volume.
void to_geogram_mesh(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, GEO::Mesh &M)
Converts a triangle mesh to a Geogram mesh.
void generate_edges(GEO::Mesh &M)
assing edges to M
void refine_polygonal_mesh(const GEO::Mesh &M_in, GEO::Mesh &M_out, Polygons::SplitFunction split_func)
Refine a polygonal mesh.
void compute_element_tags(const GEO::Mesh &M, std::vector< ElementType > &element_tags)
Compute the type of each facet in a surface mesh.
spdlog::logger & logger()
Retrieves the current logger.
Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor, 1, 3 > RowVectorNd