17 double merge_tolerance,
20 const std::string function_name = __FUNCTION__;
22 const auto& p = plane_point;
26 <<
"Ref. Point: " << p.PrintS()
27 <<
" Normal: " << n.PrintS();
34 size_t num_verts_snapped=0;
35 for (
auto& id_vertex : mesh.
vertices)
37 auto& vertex = id_vertex.second;
38 double d_from_plane = n.Dot(vertex - p);
40 if (std::fabs(d_from_plane) < merge_tolerance)
42 vertex -= n*d_from_plane;
47 Chi::log.
Log() <<
"Number of vertices snapped to plane: "
51 size_t num_bad_quality_cells = 0;
54 if (cell.Type() == CellType::POLYGON)
57 ++num_bad_quality_cells;
59 else if (cell.Type() == CellType::POLYHEDRON)
62 ++num_bad_quality_cells;
65 throw std::logic_error(function_name +
": Called for a mesh containing"
66 " an unsupported cell-type.");
68 if (num_bad_quality_cells > 0)
69 throw std::logic_error(function_name +
": Called for a mesh containing " +
70 std::to_string(num_bad_quality_cells) +
71 " bad quality cells.");
79 std::vector<chi_mesh::Cell*> cells_to_cut;
83 size_t num_neg_senses = 0;
84 size_t num_pos_senses = 0;
85 for (
auto vid : cell.vertex_ids_)
88 double new_sense = n.Dot(x-p);
90 if (new_sense < (0.0-float_compare)) ++num_neg_senses;
91 if (new_sense > (0.0+float_compare)) ++num_pos_senses;
93 if (num_neg_senses>0 && num_pos_senses>0)
95 cells_to_cut.emplace_back(&cell);
100 Chi::log.
Log() <<
"Number of cells to cut: " << cells_to_cut.size();
103 if (mesh.
local_cells[0].Type() == CellType::POLYGON)
106 std::set<uint64_t> cut_vertices;
108 for (
auto& cell_ptr : cells_to_cut)
109 for (uint64_t vid : cell_ptr->vertex_ids_)
111 const auto vertex = mesh.
vertices[vid];
112 double dv = std::fabs((vertex-p).
Dot(n));
113 if (dv<float_compare)
114 cut_vertices.insert(vid);
118 Chi::log.
Log() <<
"Number of cut vertices: " << cut_vertices.size();
121 size_t num_edges_cut=0;
122 std::set<Edge> edges_set;
123 for (
auto& cell_ptr : cells_to_cut)
125 const auto& cell = *cell_ptr;
126 const size_t num_edges = cell.vertex_ids_.size();
128 for (
size_t e=0; e<num_edges; ++e)
131 edges_set.insert(std::make_pair(std::min(edge.first,edge.second),
132 std::max(edge.first,edge.second)));
139 std::vector<ECI> cut_edges;
141 for (
auto& edge : edges_set)
143 const auto& v0 = mesh.
vertices[edge.first];
144 const auto& v1 = mesh.
vertices[edge.second];
150 double dv0 = std::fabs((v0-p).
Dot(n));
151 double dv1 = std::fabs((v1-p).
Dot(n));
152 if (dv0>float_compare and dv1>float_compare)
155 cut_edges.emplace_back(edge, new_vertex_address++);
162 Chi::log.
Log() <<
"Number of cut edges: " << num_edges_cut;
165 for (
auto& cell_ptr : cells_to_cut)
167 auto& cell = *cell_ptr;
169 CutPolygon(cut_edges,cut_vertices,p,n,mesh,cell);
174 if (mesh.
local_cells[0].Type() == CellType::POLYHEDRON)
177 std::set<uint64_t> cut_vertices;
179 for (
auto& cell_ptr : cells_to_cut)
180 for (uint64_t vid : cell_ptr->vertex_ids_)
182 const auto& vertex = mesh.
vertices[vid];
183 double dv = std::fabs((vertex-p).
Dot(n));
184 if (dv<float_compare)
185 cut_vertices.insert(vid);
190 size_t num_edges_cut=0;
191 std::set<Edge> edges_set;
192 for (
auto& cell_ptr : cells_to_cut)
194 const auto& cell = *cell_ptr;
196 for (
auto& face : cell.faces_)
198 const size_t num_edges = face.vertex_ids_.size();
200 for (
size_t e=0; e<num_edges; ++e)
203 edges_set.insert(std::make_pair(std::min(edge.first,edge.second),
204 std::max(edge.first,edge.second)));
212 std::vector<ECI> cut_edges;
214 for (
auto& edge : edges_set)
216 const auto& v0 = mesh.
vertices[edge.first];
217 const auto& v1 = mesh.
vertices[edge.second];
223 double dv0 = std::fabs((v0-p).
Dot(n));
224 double dv1 = std::fabs((v1-p).
Dot(n));
225 if (dv0>float_compare and dv1>float_compare)
228 cut_edges.emplace_back(edge, new_vertex_address++);
235 Chi::log.
Log() <<
"Number of cut edges: " << num_edges_cut;
238 for (
auto& cell_ptr : cells_to_cut)
241 plane_point,plane_normal,
247 Chi::log.
Log() <<
"Done cutting mesh with plane. Num cells = "
LogStream Log(LOG_LVL level=LOG_0)
LocalCellHandler local_cells
uint64_t GetGlobalVertexCount() const
void Insert(const uint64_t global_id, const chi_mesh::Vector3 &vec)
double Dot(const VecDbl &x, const VecDbl &y)
void CutMeshWithPlane(MeshContinuum &mesh, const Vector3 &plane_point, const Vector3 &plane_normal, double merge_tolerance=1.0e-3, double float_compare=1.0e-10)
void CutPolygon(const std::vector< ECI > &cut_edges, const std::set< uint64_t > &cut_vertices, const Vector3 &plane_point, const Vector3 &plane_normal, MeshContinuum &mesh, chi_mesh::Cell &cell)
Edge MakeEdgeFromPolygonEdgeIndex(const std::vector< uint64_t > &vertex_ids, size_t edge_index)
void Cut3DCell(const std::vector< ECI > &global_cut_edges, const std::set< uint64_t > &number, const Vector3 &plane_point, const Vector3 &plane_normal, double float_compare, MeshContinuum &mesh, chi_mesh::Cell &cell, bool verbose=false)
bool CheckPolyhedronQuality(const MeshContinuum &mesh, const chi_mesh::Cell &cell, bool check_convexity=false)
bool CheckPolygonQuality(const MeshContinuum &mesh, const chi_mesh::Cell &cell, bool check_convexity=false)
bool CheckPlaneLineIntersect(const chi_mesh::Normal &plane_normal, const chi_mesh::Vector3 &plane_point, const chi_mesh::Vector3 &line_point_0, const chi_mesh::Vector3 &line_point_1, chi_mesh::Vector3 &intersection_point, std::pair< double, double > *weights=nullptr)
Vector3 Normalized() const