Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 : #include "libmesh/libmesh_config.h"
20 : #if defined(LIBMESH_HAVE_TRIANGLE) && defined(LIBMESH_HAVE_TETGEN)
21 :
22 : // Local includes
23 : #include "libmesh/elem_cutter.h"
24 : #include "libmesh/elem.h"
25 : #include "libmesh/replicated_mesh.h"
26 : #include "libmesh/mesh_triangle_interface.h"
27 : #include "libmesh/mesh_tetgen_interface.h"
28 :
29 : // C++ Includes
30 : #include <memory>
31 :
32 : namespace
33 : {
34 : unsigned int cut_cntr;
35 : }
36 :
37 : namespace libMesh
38 : {
39 :
40 8 : ElemCutter::ElemCutter() :
41 4 : _inside_mesh_2D(std::make_unique<ReplicatedMesh>(_comm_self,2)),
42 4 : _triangle_inside(std::make_unique<TriangleInterface>(*_inside_mesh_2D)),
43 0 : _outside_mesh_2D(std::make_unique<ReplicatedMesh>(_comm_self,2)),
44 4 : _triangle_outside(std::make_unique<TriangleInterface>(*_outside_mesh_2D)),
45 0 : _inside_mesh_3D(std::make_unique<ReplicatedMesh>(_comm_self,3)),
46 4 : _tetgen_inside(std::make_unique<TetGenMeshInterface>(*_inside_mesh_3D)),
47 0 : _outside_mesh_3D(std::make_unique<ReplicatedMesh>(_comm_self,3)),
48 20 : _tetgen_outside(std::make_unique<TetGenMeshInterface>(*_outside_mesh_3D))
49 : {
50 8 : cut_cntr = 0;
51 8 : }
52 :
53 :
54 :
55 : ElemCutter::~ElemCutter() = default;
56 :
57 :
58 :
59 770 : bool ElemCutter::is_inside (const Elem & libmesh_dbg_var(elem),
60 : const std::vector<Real> & vertex_distance_func) const
61 : {
62 385 : libmesh_assert_equal_to (elem.n_vertices(), vertex_distance_func.size());
63 :
64 1236 : for (const auto & val : vertex_distance_func)
65 1236 : if (val > 0.)
66 385 : return false;
67 :
68 : // if the distance function is nonpositive, we are outside
69 0 : return true;
70 : }
71 :
72 :
73 :
74 770 : bool ElemCutter::is_outside (const Elem & libmesh_dbg_var(elem),
75 : const std::vector<Real> & vertex_distance_func) const
76 : {
77 385 : libmesh_assert_equal_to (elem.n_vertices(), vertex_distance_func.size());
78 :
79 1670 : for (const auto & val : vertex_distance_func)
80 1670 : if (val < 0.)
81 385 : return false;
82 :
83 : // if the distance function is nonnegative, we are outside
84 0 : return true;
85 : }
86 :
87 :
88 :
89 38273 : bool ElemCutter::is_cut (const Elem & libmesh_dbg_var(elem),
90 : const std::vector<Real> & vertex_distance_func) const
91 : {
92 19329 : libmesh_assert_equal_to (elem.n_vertices(), vertex_distance_func.size());
93 :
94 : Real
95 38273 : vmin = vertex_distance_func.front(),
96 38273 : vmax = vmin;
97 :
98 220173 : for (const auto & val : vertex_distance_func)
99 : {
100 181900 : vmin = std::min (vmin, val);
101 181900 : vmax = std::max (vmax, val);
102 : }
103 :
104 : // if the distance function changes sign, we're cut.
105 38273 : return (vmin*vmax < 0.);
106 : }
107 :
108 :
109 :
110 770 : void ElemCutter::operator()(const Elem & elem,
111 : const std::vector<Real> & vertex_distance_func)
112 :
113 : {
114 385 : libmesh_assert_equal_to (vertex_distance_func.size(), elem.n_vertices());
115 :
116 770 : _inside_elem.clear();
117 770 : _outside_elem.clear();
118 :
119 : // check for quick return?
120 : {
121 : // completely outside?
122 770 : if (this->is_outside(elem, vertex_distance_func))
123 : {
124 : //std::cout << "element completely outside\n";
125 0 : _outside_elem.push_back(& elem);
126 0 : return;
127 : }
128 :
129 : // completely inside?
130 770 : else if (this->is_inside(elem, vertex_distance_func))
131 : {
132 : //std::cout << "element completely inside\n";
133 0 : _inside_elem.push_back(&elem);
134 0 : return;
135 : }
136 :
137 385 : libmesh_assert (this->is_cut (elem, vertex_distance_func));
138 : }
139 :
140 : // we now know we are in a cut element, find the intersecting points.
141 770 : this->find_intersection_points (elem, vertex_distance_func);
142 :
143 : // and then dispatch the proper method
144 770 : switch (elem.dim())
145 : {
146 0 : case 1: this->cut_1D(elem, vertex_distance_func); break;
147 230 : case 2: this->cut_2D(elem, vertex_distance_func); break;
148 540 : case 3: this->cut_3D(elem, vertex_distance_func); break;
149 0 : default: libmesh_error_msg("Invalid element dimension: " << elem.dim());
150 : }
151 : }
152 :
153 :
154 :
155 770 : void ElemCutter::find_intersection_points(const Elem & elem,
156 : const std::vector<Real> & vertex_distance_func)
157 : {
158 770 : _intersection_pts.clear();
159 :
160 5422 : for (auto e : elem.edge_index_range())
161 : {
162 6978 : std::unique_ptr<const Elem> edge (elem.build_edge_ptr(e));
163 :
164 : // find the element nodes el0, el1 that map
165 : unsigned int
166 4652 : el0 = elem.get_node_index(edge->node_ptr(0)),
167 4652 : el1 = elem.get_node_index(edge->node_ptr(1));
168 :
169 2326 : libmesh_assert (elem.is_vertex(el0));
170 2326 : libmesh_assert (elem.is_vertex(el1));
171 2326 : libmesh_assert_less (el0, vertex_distance_func.size());
172 2326 : libmesh_assert_less (el1, vertex_distance_func.size());
173 :
174 : const Real
175 4652 : d0 = vertex_distance_func[el0],
176 4652 : d1 = vertex_distance_func[el1];
177 :
178 : // if this edge has a 0 crossing
179 4652 : if (d0*d1 < 0.)
180 : {
181 1135 : libmesh_assert_not_equal_to (d0, d1);
182 :
183 : // then find d_star in [0,1], the
184 : // distance from el0 to el1 where the 0 lives.
185 2270 : const Real d_star = d0 / (d0 - d1);
186 :
187 :
188 : // Prevent adding nodes trivially close to existing
189 : // nodes.
190 1135 : const Real endpoint_tol = 0.01;
191 :
192 2270 : if ( (d_star > endpoint_tol) &&
193 1135 : (d_star < (1.-endpoint_tol)) )
194 : {
195 2270 : const Point x_star = (edge->point(0)*(1-d_star) +
196 3405 : edge->point(1)*d_star);
197 :
198 : /*
199 : std::cout << "adding cut point (d_star, x_star) = "
200 : << d_star << " , " << x_star << std::endl;
201 : */
202 :
203 2270 : _intersection_pts.push_back (x_star);
204 : }
205 : }
206 : }
207 770 : }
208 :
209 :
210 :
211 :
212 0 : void ElemCutter::cut_1D (const Elem & /*elem*/,
213 : const std::vector<Real> &/*vertex_distance_func*/)
214 : {
215 0 : libmesh_not_implemented();
216 : }
217 :
218 :
219 :
220 230 : void ElemCutter::cut_2D (const Elem & elem,
221 : const std::vector<Real> & vertex_distance_func)
222 : {
223 : #ifndef LIBMESH_HAVE_TRIANGLE
224 :
225 : // current implementation requires triangle!
226 : libMesh::err << "ERROR: current libMesh ElemCutter 2D implementation requires\n"
227 : << " the \"triangle\" library!\n"
228 : << std::endl;
229 : libmesh_not_implemented();
230 :
231 : #else // OK, LIBMESH_HAVE_TRIANGLE
232 :
233 : // std::cout << "Inside cut face element!\n";
234 :
235 115 : libmesh_assert (_inside_mesh_2D.get() != nullptr);
236 115 : libmesh_assert (_outside_mesh_2D.get() != nullptr);
237 :
238 230 : _inside_mesh_2D->clear();
239 230 : _outside_mesh_2D->clear();
240 :
241 982 : for (auto v : make_range(elem.n_vertices()))
242 : {
243 1128 : if (vertex_distance_func[v] >= 0.)
244 579 : _outside_mesh_2D->add_point (elem.point(v));
245 :
246 1128 : if (vertex_distance_func[v] <= 0.)
247 555 : _inside_mesh_2D->add_point (elem.point(v));
248 : }
249 :
250 686 : for (const auto & pt : _intersection_pts)
251 : {
252 456 : _inside_mesh_2D->add_point(pt);
253 456 : _outside_mesh_2D->add_point(pt);
254 : }
255 :
256 :
257 : // Customize the variables for the triangulation
258 : // we will be cutting reference cell, and want as few triangles
259 : // as possible, so jack this up larger than the area we will be
260 : // triangulating so we are governed only by accurately defining
261 : // the boundaries.
262 230 : _triangle_inside->desired_area() = 100.;
263 230 : _triangle_outside->desired_area() = 100.;
264 :
265 : // allow for small angles
266 230 : _triangle_inside->minimum_angle() = 5.;
267 230 : _triangle_outside->minimum_angle() = 5.;
268 :
269 : // Turn off Laplacian mesh smoothing after generation.
270 230 : _triangle_inside->smooth_after_generating() = false;
271 230 : _triangle_outside->smooth_after_generating() = false;
272 :
273 : // Triangulate!
274 230 : _triangle_inside->triangulate();
275 230 : _triangle_outside->triangulate();
276 :
277 : // std::ostringstream name;
278 :
279 : // name << "cut_face_"
280 : // << cut_cntr++
281 : // << ".dat";
282 : // _inside_mesh_2D->write ("in_" + name.str());
283 : // _outside_mesh_2D->write ("out_" + name.str());
284 :
285 : // finally, add the elements to our lists.
286 230 : _inside_elem.clear();
287 230 : _outside_elem.clear();
288 :
289 721 : for (const auto & el : _inside_mesh_2D->element_ptr_range())
290 376 : _inside_elem.push_back (el);
291 :
292 749 : for (const auto & el : _outside_mesh_2D->element_ptr_range())
293 404 : _outside_elem.push_back (el);
294 :
295 : #endif
296 230 : }
297 :
298 :
299 :
300 540 : void ElemCutter::cut_3D (const Elem & elem,
301 : const std::vector<Real> & vertex_distance_func)
302 : {
303 : #ifndef LIBMESH_HAVE_TETGEN
304 :
305 : // current implementation requires tetgen!
306 : libMesh::err << "ERROR: current libMesh ElemCutter 3D implementation requires\n"
307 : << " the \"tetgen\" library!\n"
308 : << std::endl;
309 : libmesh_not_implemented();
310 :
311 : #else // OK, LIBMESH_HAVE_TETGEN
312 :
313 : // std::cout << "Inside cut cell element!\n";
314 :
315 270 : libmesh_assert (_inside_mesh_3D.get() != nullptr);
316 270 : libmesh_assert (_outside_mesh_3D.get() != nullptr);
317 :
318 540 : _inside_mesh_3D->clear();
319 540 : _outside_mesh_3D->clear();
320 :
321 3140 : for (auto v : make_range(elem.n_vertices()))
322 : {
323 3900 : if (vertex_distance_func[v] >= 0.)
324 2139 : _outside_mesh_3D->add_point (elem.point(v));
325 :
326 3900 : if (vertex_distance_func[v] <= 0.)
327 1791 : _inside_mesh_3D->add_point (elem.point(v));
328 : }
329 :
330 2354 : for (const auto & pt : _intersection_pts)
331 : {
332 1814 : _inside_mesh_3D->add_point(pt);
333 1814 : _outside_mesh_3D->add_point(pt);
334 : }
335 :
336 :
337 : // Triangulate!
338 540 : _tetgen_inside->triangulate_pointset();
339 : //_inside_mesh_3D->print_info();
340 540 : _tetgen_outside->triangulate_pointset();
341 : //_outside_mesh_3D->print_info();
342 :
343 :
344 : // (below generates some horribly expensive meshes,
345 : // but seems immune to the 0 volume problem).
346 : // _tetgen_inside->pointset_convexhull();
347 : // _inside_mesh_3D->find_neighbors();
348 : // _inside_mesh_3D->print_info();
349 : // _tetgen_inside->triangulate_conformingDelaunayMesh (1.e3, 100.);
350 : // _inside_mesh_3D->print_info();
351 :
352 : // _tetgen_outside->pointset_convexhull();
353 : // _outside_mesh_3D->find_neighbors();
354 : // _outside_mesh_3D->print_info();
355 : // _tetgen_outside->triangulate_conformingDelaunayMesh (1.e3, 100.);
356 : // _outside_mesh_3D->print_info();
357 :
358 1080 : std::ostringstream name;
359 :
360 270 : name << "cut_cell_"
361 540 : << cut_cntr++
362 540 : << ".dat";
363 810 : _inside_mesh_3D->write ("in_" + name.str());
364 810 : _outside_mesh_3D->write ("out_" + name.str());
365 :
366 : // finally, add the elements to our lists.
367 540 : _inside_elem.clear();
368 540 : _outside_elem.clear();
369 :
370 2366 : for (const auto & el : _inside_mesh_3D->element_ptr_range())
371 1556 : if (el->volume() > std::numeric_limits<Real>::epsilon())
372 1480 : _inside_elem.push_back (el);
373 :
374 2728 : for (const auto & el : _outside_mesh_3D->element_ptr_range())
375 1918 : if (el->volume() > std::numeric_limits<Real>::epsilon())
376 1770 : _outside_elem.push_back (el);
377 :
378 : #endif
379 540 : }
380 :
381 :
382 :
383 : } // namespace libMesh
384 :
385 : #endif // LIBMESH_HAVE_TRIANGLE && LIBMESH_HAVE_TETGEN
|
