libMesh
elem_cutter.C
Go to the documentation of this file.
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 
41  _inside_mesh_2D(std::make_unique<ReplicatedMesh>(_comm_self,2)),
42  _triangle_inside(std::make_unique<TriangleInterface>(*_inside_mesh_2D)),
43  _outside_mesh_2D(std::make_unique<ReplicatedMesh>(_comm_self,2)),
44  _triangle_outside(std::make_unique<TriangleInterface>(*_outside_mesh_2D)),
45  _inside_mesh_3D(std::make_unique<ReplicatedMesh>(_comm_self,3)),
46  _tetgen_inside(std::make_unique<TetGenMeshInterface>(*_inside_mesh_3D)),
47  _outside_mesh_3D(std::make_unique<ReplicatedMesh>(_comm_self,3)),
48  _tetgen_outside(std::make_unique<TetGenMeshInterface>(*_outside_mesh_3D))
49 {
50  cut_cntr = 0;
51 }
52 
53 
54 
55 ElemCutter::~ElemCutter() = default;
56 
57 
58 
59 bool ElemCutter::is_inside (const Elem & libmesh_dbg_var(elem),
60  const std::vector<Real> & vertex_distance_func) const
61 {
62  libmesh_assert_equal_to (elem.n_vertices(), vertex_distance_func.size());
63 
64  for (const auto & val : vertex_distance_func)
65  if (val > 0.)
66  return false;
67 
68  // if the distance function is nonpositive, we are outside
69  return true;
70 }
71 
72 
73 
74 bool ElemCutter::is_outside (const Elem & libmesh_dbg_var(elem),
75  const std::vector<Real> & vertex_distance_func) const
76 {
77  libmesh_assert_equal_to (elem.n_vertices(), vertex_distance_func.size());
78 
79  for (const auto & val : vertex_distance_func)
80  if (val < 0.)
81  return false;
82 
83  // if the distance function is nonnegative, we are outside
84  return true;
85 }
86 
87 
88 
89 bool ElemCutter::is_cut (const Elem & libmesh_dbg_var(elem),
90  const std::vector<Real> & vertex_distance_func) const
91 {
92  libmesh_assert_equal_to (elem.n_vertices(), vertex_distance_func.size());
93 
94  Real
95  vmin = vertex_distance_func.front(),
96  vmax = vmin;
97 
98  for (const auto & val : vertex_distance_func)
99  {
100  vmin = std::min (vmin, val);
101  vmax = std::max (vmax, val);
102  }
103 
104  // if the distance function changes sign, we're cut.
105  return (vmin*vmax < 0.);
106 }
107 
108 
109 
110 void ElemCutter::operator()(const Elem & elem,
111  const std::vector<Real> & vertex_distance_func)
112 
113 {
114  libmesh_assert_equal_to (vertex_distance_func.size(), elem.n_vertices());
115 
116  _inside_elem.clear();
117  _outside_elem.clear();
118 
119  // check for quick return?
120  {
121  // completely outside?
122  if (this->is_outside(elem, vertex_distance_func))
123  {
124  //std::cout << "element completely outside\n";
125  _outside_elem.push_back(& elem);
126  return;
127  }
128 
129  // completely inside?
130  else if (this->is_inside(elem, vertex_distance_func))
131  {
132  //std::cout << "element completely inside\n";
133  _inside_elem.push_back(&elem);
134  return;
135  }
136 
137  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  this->find_intersection_points (elem, vertex_distance_func);
142 
143  // and then dispatch the proper method
144  switch (elem.dim())
145  {
146  case 1: this->cut_1D(elem, vertex_distance_func); break;
147  case 2: this->cut_2D(elem, vertex_distance_func); break;
148  case 3: this->cut_3D(elem, vertex_distance_func); break;
149  default: libmesh_error_msg("Invalid element dimension: " << elem.dim());
150  }
151 }
152 
153 
154 
156  const std::vector<Real> & vertex_distance_func)
157 {
158  _intersection_pts.clear();
159 
160  for (auto e : elem.edge_index_range())
161  {
162  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  el0 = elem.get_node_index(edge->node_ptr(0)),
167  el1 = elem.get_node_index(edge->node_ptr(1));
168 
169  libmesh_assert (elem.is_vertex(el0));
170  libmesh_assert (elem.is_vertex(el1));
171  libmesh_assert_less (el0, vertex_distance_func.size());
172  libmesh_assert_less (el1, vertex_distance_func.size());
173 
174  const Real
175  d0 = vertex_distance_func[el0],
176  d1 = vertex_distance_func[el1];
177 
178  // if this edge has a 0 crossing
179  if (d0*d1 < 0.)
180  {
181  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  const Real d_star = d0 / (d0 - d1);
186 
187 
188  // Prevent adding nodes trivially close to existing
189  // nodes.
190  const Real endpoint_tol = 0.01;
191 
192  if ( (d_star > endpoint_tol) &&
193  (d_star < (1.-endpoint_tol)) )
194  {
195  const Point x_star = (edge->point(0)*(1-d_star) +
196  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  _intersection_pts.push_back (x_star);
204  }
205  }
206  }
207 }
208 
209 
210 
211 
212 void ElemCutter::cut_1D (const Elem & /*elem*/,
213  const std::vector<Real> &/*vertex_distance_func*/)
214 {
215  libmesh_not_implemented();
216 }
217 
218 
219 
220 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  libmesh_assert (_inside_mesh_2D.get() != nullptr);
236  libmesh_assert (_outside_mesh_2D.get() != nullptr);
237 
238  _inside_mesh_2D->clear();
239  _outside_mesh_2D->clear();
240 
241  for (auto v : make_range(elem.n_vertices()))
242  {
243  if (vertex_distance_func[v] >= 0.)
244  _outside_mesh_2D->add_point (elem.point(v));
245 
246  if (vertex_distance_func[v] <= 0.)
247  _inside_mesh_2D->add_point (elem.point(v));
248  }
249 
250  for (const auto & pt : _intersection_pts)
251  {
252  _inside_mesh_2D->add_point(pt);
253  _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  _triangle_inside->desired_area() = 100.;
263  _triangle_outside->desired_area() = 100.;
264 
265  // allow for small angles
266  _triangle_inside->minimum_angle() = 5.;
267  _triangle_outside->minimum_angle() = 5.;
268 
269  // Turn off Laplacian mesh smoothing after generation.
270  _triangle_inside->smooth_after_generating() = false;
271  _triangle_outside->smooth_after_generating() = false;
272 
273  // Triangulate!
274  _triangle_inside->triangulate();
275  _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  _inside_elem.clear();
287  _outside_elem.clear();
288 
289  for (const auto & el : _inside_mesh_2D->element_ptr_range())
290  _inside_elem.push_back (el);
291 
292  for (const auto & el : _outside_mesh_2D->element_ptr_range())
293  _outside_elem.push_back (el);
294 
295 #endif
296 }
297 
298 
299 
300 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  libmesh_assert (_inside_mesh_3D.get() != nullptr);
316  libmesh_assert (_outside_mesh_3D.get() != nullptr);
317 
318  _inside_mesh_3D->clear();
319  _outside_mesh_3D->clear();
320 
321  for (auto v : make_range(elem.n_vertices()))
322  {
323  if (vertex_distance_func[v] >= 0.)
324  _outside_mesh_3D->add_point (elem.point(v));
325 
326  if (vertex_distance_func[v] <= 0.)
327  _inside_mesh_3D->add_point (elem.point(v));
328  }
329 
330  for (const auto & pt : _intersection_pts)
331  {
332  _inside_mesh_3D->add_point(pt);
333  _outside_mesh_3D->add_point(pt);
334  }
335 
336 
337  // Triangulate!
338  _tetgen_inside->triangulate_pointset();
339  //_inside_mesh_3D->print_info();
340  _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  std::ostringstream name;
359 
360  name << "cut_cell_"
361  << cut_cntr++
362  << ".dat";
363  _inside_mesh_3D->write ("in_" + name.str());
364  _outside_mesh_3D->write ("out_" + name.str());
365 
366  // finally, add the elements to our lists.
367  _inside_elem.clear();
368  _outside_elem.clear();
369 
370  for (const auto & el : _inside_mesh_3D->element_ptr_range())
371  if (el->volume() > std::numeric_limits<Real>::epsilon())
372  _inside_elem.push_back (el);
373 
374  for (const auto & el : _outside_mesh_3D->element_ptr_range())
375  if (el->volume() > std::numeric_limits<Real>::epsilon())
376  _outside_elem.push_back (el);
377 
378 #endif
379 }
380 
381 
382 
383 } // namespace libMesh
384 
385 #endif // LIBMESH_HAVE_TRIANGLE && LIBMESH_HAVE_TETGEN
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:42
OStreamProxy err
std::unique_ptr< ReplicatedMesh > _inside_mesh_3D
Definition: elem_cutter.h:166
Class TetGenMeshInterface provides an interface for tetrahedralization of meshes using the TetGen lib...
The ReplicatedMesh class is derived from the MeshBase class, and is used to store identical copies of...
std::unique_ptr< ReplicatedMesh > _inside_mesh_2D
Definition: elem_cutter.h:160
void find_intersection_points(const Elem &elem, const std::vector< Real > &vertex_distance_func)
Finds the points where the cutting surface intersects the element edges.
Definition: elem_cutter.C:155
std::vector< Point > _intersection_pts
Definition: elem_cutter.h:172
std::unique_ptr< TriangleInterface > _triangle_outside
Definition: elem_cutter.h:164
unsigned int get_node_index(const Node *node_ptr) const
Definition: elem.h:2545
void cut_1D(const Elem &elem, const std::vector< Real > &vertex_distance_func)
cutting algorithm in 1D.
Definition: elem_cutter.C:212
ElemCutter()
Constructor.
Definition: elem_cutter.C:40
void operator()(const Elem &elem_in, const std::vector< Real > &vertex_distance_func)
This function implements cutting an element by a signed distance function.
Definition: elem_cutter.C:110
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
The libMesh namespace provides an interface to certain functionality in the library.
bool is_cut(const Elem &elem, const std::vector< Real > &vertex_distance_func) const
Definition: elem_cutter.C:89
void cut_2D(const Elem &elem, const std::vector< Real > &vertex_distance_func)
cutting algorithm in 2D.
Definition: elem_cutter.C:220
IntRange< unsigned short > edge_index_range() const
Definition: elem.h:2692
std::vector< Elem const * > _inside_elem
Definition: elem_cutter.h:155
std::unique_ptr< ReplicatedMesh > _outside_mesh_3D
Definition: elem_cutter.h:169
std::unique_ptr< TriangleInterface > _triangle_inside
Definition: elem_cutter.h:161
void cut_3D(const Elem &elem, const std::vector< Real > &vertex_distance_func)
cutting algorithm in 3D.
Definition: elem_cutter.C:300
libmesh_assert(ctx)
std::unique_ptr< TetGenMeshInterface > _tetgen_inside
Definition: elem_cutter.h:167
std::unique_ptr< TetGenMeshInterface > _tetgen_outside
Definition: elem_cutter.h:170
bool is_outside(const Elem &elem, const std::vector< Real > &vertex_distance_func) const
Definition: elem_cutter.C:74
bool is_inside(const Elem &elem, const std::vector< Real > &vertex_distance_func) const
Definition: elem_cutter.C:59
std::vector< Elem const * > _outside_elem
Definition: elem_cutter.h:156
~ElemCutter()
Destructor.
virtual unsigned int n_vertices() const =0
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual unsigned short dim() const =0
virtual bool is_vertex(const unsigned int i) const =0
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
std::unique_ptr< ReplicatedMesh > _outside_mesh_2D
Definition: elem_cutter.h:163
virtual std::unique_ptr< Elem > build_edge_ptr(const unsigned int i)=0
A C++ interface between LibMesh and the Triangle library written by J.R.
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
const Point & point(const unsigned int i) const
Definition: elem.h:2453