Chi-Tech
cutmesh_plane.cc
Go to the documentation of this file.
1#include "meshcutting.h"
2
5
6#include "chi_runtime.h"
7#include "chi_log.h"
8
9#include <algorithm>
10
11//###################################################################
12/**Cuts a mesh with a plane.*/
15 const Vector3 &plane_point,
16 const Vector3 &plane_normal,
17 double merge_tolerance/*=1.0e-3*/,
18 double float_compare/*=1.0e-10*/)
19{
20 const std::string function_name = __FUNCTION__;
21
22 const auto& p = plane_point;
23 const auto& n = plane_normal.Normalized();
24
25 Chi::log.Log() << "Cutting mesh with plane. "
26 << "Ref. Point: " << p.PrintS()
27 << " Normal: " << n.PrintS();
28
29 //============================================= Snap vertices to plane within
30 // merge tolerance to avoid
31 // creating small cells or cutting
32 // parallel faces
33 // Order N, num_vertices
34 size_t num_verts_snapped=0;
35 for (auto& id_vertex : mesh.vertices)
36 {
37 auto& vertex = id_vertex.second;
38 double d_from_plane = n.Dot(vertex - p);
39
40 if (std::fabs(d_from_plane) < merge_tolerance)
41 {
42 vertex -= n*d_from_plane;
43 ++num_verts_snapped;
44 }
45 }
46
47 Chi::log.Log() << "Number of vertices snapped to plane: "
48 << num_verts_snapped;
49
50 //============================================= Perform quality checks
51 size_t num_bad_quality_cells = 0;
52 for (const auto& cell : mesh.local_cells)
53 {
54 if (cell.Type() == CellType::POLYGON)
55 {
56 if (not CheckPolygonQuality(mesh,cell))
57 ++num_bad_quality_cells;
58 }
59 else if (cell.Type() == CellType::POLYHEDRON)
60 {
61 if (not CheckPolyhedronQuality(mesh,cell,/*check_convexity*/true))
62 ++num_bad_quality_cells;
63 }
64 else
65 throw std::logic_error(function_name + ": Called for a mesh containing"
66 " an unsupported cell-type.");
67 }
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.");
72
73 //============================================= Determine cells to cut
74 // Order N, num_cells
75 // A cell is a candidate for cutting if its
76 // vertices lay on both sides of the plane.
77 // So this algorithm just checks the sense wrt
78 // the plane. Works for both 2D and 3D
79 std::vector<chi_mesh::Cell*> cells_to_cut;
80
81 for (auto& cell : mesh.local_cells)
82 {
83 size_t num_neg_senses = 0;
84 size_t num_pos_senses = 0;
85 for (auto vid : cell.vertex_ids_)
86 {
87 const auto& x = mesh.vertices[vid];
88 double new_sense = n.Dot(x-p);
89
90 if (new_sense < (0.0-float_compare)) ++num_neg_senses;
91 if (new_sense > (0.0+float_compare)) ++num_pos_senses;
92
93 if (num_neg_senses>0 && num_pos_senses>0)
94 {
95 cells_to_cut.emplace_back(&cell);
96 break;
97 }
98 }//for vid
99 }//for cell
100 Chi::log.Log() << "Number of cells to cut: " << cells_to_cut.size();
101
102 //============================================= Two-D algorithm
103 if (mesh.local_cells[0].Type() == CellType::POLYGON)
104 {
105 //====================================== Determine cut vertices
106 std::set<uint64_t> cut_vertices;
107 {
108 for (auto& cell_ptr : cells_to_cut)
109 for (uint64_t vid : cell_ptr->vertex_ids_)
110 {
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);
115 }//for vid
116 }//populate cut_vertices
117
118 Chi::log.Log() << "Number of cut vertices: " << cut_vertices.size();
119
120 //====================================== Build unique edges
121 size_t num_edges_cut=0;
122 std::set<Edge> edges_set;
123 for (auto& cell_ptr : cells_to_cut)
124 {
125 const auto& cell = *cell_ptr;
126 const size_t num_edges = cell.vertex_ids_.size();
127
128 for (size_t e=0; e<num_edges; ++e)
129 {
130 auto edge = MakeEdgeFromPolygonEdgeIndex(cell.vertex_ids_, e);
131 edges_set.insert(std::make_pair(std::min(edge.first,edge.second),
132 std::max(edge.first,edge.second)));
133 }
134 }//for cell - built edges_set
135
136 uint64_t new_vertex_address = mesh.GetGlobalVertexCount();
137
138 //====================================== Determine cut edges
139 std::vector<ECI> cut_edges;
140 {
141 for (auto& edge : edges_set)
142 {
143 const auto& v0 = mesh.vertices[edge.first];
144 const auto& v1 = mesh.vertices[edge.second];
145
146 chi_mesh::Vector3 cut_point;
147
148 if (CheckPlaneLineIntersect(n,p,v0,v1,cut_point))
149 {
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)
153 {
154 mesh.vertices.Insert(new_vertex_address,cut_point);
155 cut_edges.emplace_back(edge, new_vertex_address++);
156 ++num_edges_cut;
157 }
158 }
159 }//for edge - determine cut
160 }//populate edges cut
161
162 Chi::log.Log() << "Number of cut edges: " << num_edges_cut;
163
164 //====================================== Process cells that are cut
165 for (auto& cell_ptr : cells_to_cut)
166 {
167 auto& cell = *cell_ptr;
168
169 CutPolygon(cut_edges,cut_vertices,p,n,mesh,cell);
170 }//for cell_ptr
171 }//two-D
172
173 //============================================= Three-D algorithm
174 if (mesh.local_cells[0].Type() == CellType::POLYHEDRON)
175 {
176 //====================================== Determine cut vertices
177 std::set<uint64_t> cut_vertices;
178 {
179 for (auto& cell_ptr : cells_to_cut)
180 for (uint64_t vid : cell_ptr->vertex_ids_)
181 {
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);
186 }//for vid
187 }//populate cut_vertices
188
189 //====================================== Build unique edges
190 size_t num_edges_cut=0;
191 std::set<Edge> edges_set;
192 for (auto& cell_ptr : cells_to_cut)
193 {
194 const auto& cell = *cell_ptr;
195
196 for (auto& face : cell.faces_)
197 {
198 const size_t num_edges = face.vertex_ids_.size();
199
200 for (size_t e=0; e<num_edges; ++e)
201 {
202 auto edge = MakeEdgeFromPolygonEdgeIndex(face.vertex_ids_, e);
203 edges_set.insert(std::make_pair(std::min(edge.first,edge.second),
204 std::max(edge.first,edge.second)));
205 }
206 }//for face
207 }//for cell - built edges_set
208
209 uint64_t new_vertex_address = mesh.GetGlobalVertexCount();
210
211 //====================================== Determine cut edges
212 std::vector<ECI> cut_edges;
213 {
214 for (auto& edge : edges_set)
215 {
216 const auto& v0 = mesh.vertices[edge.first];
217 const auto& v1 = mesh.vertices[edge.second];
218
219 chi_mesh::Vector3 cut_point;
220
221 if (CheckPlaneLineIntersect(n,p,v0,v1,cut_point))
222 {
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)
226 {
227 mesh.vertices.Insert(new_vertex_address,cut_point);
228 cut_edges.emplace_back(edge, new_vertex_address++);
229 ++num_edges_cut;
230 }
231 }
232 }//for edge - determine cut
233 }//populate edges cut
234
235 Chi::log.Log() << "Number of cut edges: " << num_edges_cut;
236
237 //====================================== Process cells that are cut
238 for (auto& cell_ptr : cells_to_cut)
239 {
240 Cut3DCell(cut_edges, cut_vertices,
241 plane_point,plane_normal,
242 float_compare,
243 mesh,*cell_ptr);
244 }//for cell_ptr
245 }
246
247 Chi::log.Log() << "Done cutting mesh with plane. Num cells = "
248 << mesh.local_cells.size();
249}
static chi::ChiLog & log
Definition: chi_runtime.h:81
LogStream Log(LOG_LVL level=LOG_0)
Definition: chi_log.cc:35
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