libMesh::ElemCutter Class Reference

Subdivides a single element using a mesh generator. More...

#include <elem_cutter.h>

Public Member Functions

 ElemCutter ()
 
 ~ElemCutter ()
 
bool is_inside (const Elem &elem, const std::vector< Real > &vertex_distance_func) const
 
bool is_outside (const Elem &elem, const std::vector< Real > &vertex_distance_func) const
 
bool is_cut (const Elem &elem, const std::vector< Real > &vertex_distance_func) const
 
void operator() (const Elem &elem_in, const std::vector< Real > &vertex_distance_func)
 
const std::vector< Elem const * > & inside_elements () const
 
const std::vector< Elem const * > & outside_elements () const
 

Protected Member Functions

void find_intersection_points (const Elem &elem, const std::vector< Real > &vertex_distance_func)
 
void cut_1D (const Elem &elem, const std::vector< Real > &vertex_distance_func)
 
void cut_2D (const Elem &elem, const std::vector< Real > &vertex_distance_func)
 
void cut_3D (const Elem &elem, const std::vector< Real > &vertex_distance_func)
 

Protected Attributes

std::vector< Elem const * > _inside_elem
 
std::vector< Elem const * > _outside_elem
 
Parallel::Communicator _comm_self
 
std::unique_ptr< ReplicatedMesh_inside_mesh_2D
 
std::unique_ptr< TriangleInterface_triangle_inside
 
std::unique_ptr< ReplicatedMesh_outside_mesh_2D
 
std::unique_ptr< TriangleInterface_triangle_outside
 
std::unique_ptr< ReplicatedMesh_inside_mesh_3D
 
std::unique_ptr< TetGenMeshInterface_tetgen_inside
 
std::unique_ptr< ReplicatedMesh_outside_mesh_3D
 
std::unique_ptr< TetGenMeshInterface_tetgen_outside
 
std::vector< Point_intersection_pts
 

Detailed Description

Subdivides a single element using a mesh generator.

This class implements cutting a single element into a collection of subelements. This class depends on libmesh's Triangle and Tetgen interfaces, the former of which is only defined if libmesh is configured with –disable-strict-lgpl.

Author
Benjamin S. Kirk
Date
2013

Definition at line 60 of file elem_cutter.h.

Constructor & Destructor Documentation

◆ ElemCutter()

libMesh::ElemCutter::ElemCutter ( )

Constructor. Initializes pointer data without requiring a full ReplicatedMesh in this header file.

Definition at line 38 of file elem_cutter.C.

38  :
39  _inside_mesh_2D(libmesh_make_unique<ReplicatedMesh>(_comm_self,2)),
40  _triangle_inside(libmesh_make_unique<TriangleInterface>(*_inside_mesh_2D)),
41  _outside_mesh_2D(libmesh_make_unique<ReplicatedMesh>(_comm_self,2)),
42  _triangle_outside(libmesh_make_unique<TriangleInterface>(*_outside_mesh_2D)),
43  _inside_mesh_3D(libmesh_make_unique<ReplicatedMesh>(_comm_self,3)),
44  _tetgen_inside(libmesh_make_unique<TetGenMeshInterface>(*_inside_mesh_3D)),
45  _outside_mesh_3D(libmesh_make_unique<ReplicatedMesh>(_comm_self,3)),
46  _tetgen_outside(libmesh_make_unique<TetGenMeshInterface>(*_outside_mesh_3D))
47 {
48  cut_cntr = 0;
49 }
std::unique_ptr< ReplicatedMesh > _inside_mesh_3D
Definition: elem_cutter.h:165
std::unique_ptr< ReplicatedMesh > _inside_mesh_2D
Definition: elem_cutter.h:159
std::unique_ptr< TriangleInterface > _triangle_outside
Definition: elem_cutter.h:163
std::unique_ptr< ReplicatedMesh > _outside_mesh_3D
Definition: elem_cutter.h:168
std::unique_ptr< TriangleInterface > _triangle_inside
Definition: elem_cutter.h:160
std::unique_ptr< TetGenMeshInterface > _tetgen_inside
Definition: elem_cutter.h:166
std::unique_ptr< TetGenMeshInterface > _tetgen_outside
Definition: elem_cutter.h:169
std::unique_ptr< ReplicatedMesh > _outside_mesh_2D
Definition: elem_cutter.h:162
Parallel::Communicator _comm_self
Definition: elem_cutter.h:157

◆ ~ElemCutter()

libMesh::ElemCutter::~ElemCutter ( )

Destructor.

Definition at line 53 of file elem_cutter.C.

54 {}

Member Function Documentation

◆ cut_1D()

void libMesh::ElemCutter::cut_1D ( const Elem elem,
const std::vector< Real > &  vertex_distance_func 
)
protected

cutting algorithm in 1D.

Definition at line 209 of file elem_cutter.C.

Referenced by operator()().

211 {
212  libmesh_not_implemented();
213 }

◆ cut_2D()

void libMesh::ElemCutter::cut_2D ( const Elem elem,
const std::vector< Real > &  vertex_distance_func 
)
protected

cutting algorithm in 2D.

Definition at line 217 of file elem_cutter.C.

References _inside_elem, _inside_mesh_2D, _intersection_pts, _outside_elem, _outside_mesh_2D, _triangle_inside, _triangle_outside, libMesh::err, libMesh::Elem::n_vertices(), and libMesh::Elem::point().

Referenced by operator()().

219 {
220 #ifndef LIBMESH_HAVE_TRIANGLE
221 
222  // current implementation requires triangle!
223  libMesh::err << "ERROR: current libMesh ElemCutter 2D implementation requires\n"
224  << " the \"triangle\" library!\n"
225  << std::endl;
226  libmesh_not_implemented();
227 
228 #else // OK, LIBMESH_HAVE_TRIANGLE
229 
230  std::cout << "Inside cut face element!\n";
231 
232  libmesh_assert (_inside_mesh_2D.get() != nullptr);
233  libmesh_assert (_outside_mesh_2D.get() != nullptr);
234 
235  _inside_mesh_2D->clear();
236  _outside_mesh_2D->clear();
237 
238  for (unsigned int v=0; v<elem.n_vertices(); v++)
239  {
240  if (vertex_distance_func[v] >= 0.)
241  _outside_mesh_2D->add_point (elem.point(v));
242 
243  if (vertex_distance_func[v] <= 0.)
244  _inside_mesh_2D->add_point (elem.point(v));
245  }
246 
247  for (const auto & pt : _intersection_pts)
248  {
249  _inside_mesh_2D->add_point(pt);
250  _outside_mesh_2D->add_point(pt);
251  }
252 
253 
254  // Customize the variables for the triangulation
255  // we will be cutting reference cell, and want as few triangles
256  // as possible, so jack this up larger than the area we will be
257  // triangulating so we are governed only by accurately defining
258  // the boundaries.
259  _triangle_inside->desired_area() = 100.;
260  _triangle_outside->desired_area() = 100.;
261 
262  // allow for small angles
263  _triangle_inside->minimum_angle() = 5.;
264  _triangle_outside->minimum_angle() = 5.;
265 
266  // Turn off Laplacian mesh smoothing after generation.
267  _triangle_inside->smooth_after_generating() = false;
268  _triangle_outside->smooth_after_generating() = false;
269 
270  // Triangulate!
271  _triangle_inside->triangulate();
272  _triangle_outside->triangulate();
273 
274  // std::ostringstream name;
275 
276  // name << "cut_face_"
277  // << cut_cntr++
278  // << ".dat";
279  // _inside_mesh_2D->write ("in_" + name.str());
280  // _outside_mesh_2D->write ("out_" + name.str());
281 
282  // finally, add the elements to our lists.
283  _inside_elem.clear();
284  _outside_elem.clear();
285 
286  for (const auto & elem : _inside_mesh_2D->element_ptr_range())
287  _inside_elem.push_back (elem);
288 
289  for (const auto & elem : _outside_mesh_2D->element_ptr_range())
290  _outside_elem.push_back (elem);
291 
292 #endif
293 }
std::unique_ptr< ReplicatedMesh > _inside_mesh_2D
Definition: elem_cutter.h:159
std::vector< Point > _intersection_pts
Definition: elem_cutter.h:171
std::unique_ptr< TriangleInterface > _triangle_outside
Definition: elem_cutter.h:163
std::vector< Elem const * > _inside_elem
Definition: elem_cutter.h:154
std::unique_ptr< TriangleInterface > _triangle_inside
Definition: elem_cutter.h:160
OStreamProxy err(std::cerr)
std::vector< Elem const * > _outside_elem
Definition: elem_cutter.h:155
std::unique_ptr< ReplicatedMesh > _outside_mesh_2D
Definition: elem_cutter.h:162

◆ cut_3D()

void libMesh::ElemCutter::cut_3D ( const Elem elem,
const std::vector< Real > &  vertex_distance_func 
)
protected

cutting algorithm in 3D.

Definition at line 297 of file elem_cutter.C.

References _inside_elem, _inside_mesh_3D, _intersection_pts, _outside_elem, _outside_mesh_3D, _tetgen_inside, _tetgen_outside, libMesh::err, libMesh::Elem::n_vertices(), libMesh::Quality::name(), libMesh::Elem::point(), and libMesh::Elem::volume().

Referenced by operator()().

299 {
300 #ifndef LIBMESH_HAVE_TETGEN
301 
302  // current implementation requires tetgen!
303  libMesh::err << "ERROR: current libMesh ElemCutter 3D implementation requires\n"
304  << " the \"tetgen\" library!\n"
305  << std::endl;
306  libmesh_not_implemented();
307 
308 #else // OK, LIBMESH_HAVE_TETGEN
309 
310  std::cout << "Inside cut cell element!\n";
311 
312  libmesh_assert (_inside_mesh_3D.get() != nullptr);
313  libmesh_assert (_outside_mesh_3D.get() != nullptr);
314 
315  _inside_mesh_3D->clear();
316  _outside_mesh_3D->clear();
317 
318  for (unsigned int v=0; v<elem.n_vertices(); v++)
319  {
320  if (vertex_distance_func[v] >= 0.)
321  _outside_mesh_3D->add_point (elem.point(v));
322 
323  if (vertex_distance_func[v] <= 0.)
324  _inside_mesh_3D->add_point (elem.point(v));
325  }
326 
327  for (const auto & pt : _intersection_pts)
328  {
329  _inside_mesh_3D->add_point(pt);
330  _outside_mesh_3D->add_point(pt);
331  }
332 
333 
334  // Triangulate!
335  _tetgen_inside->triangulate_pointset();
336  //_inside_mesh_3D->print_info();
337  _tetgen_outside->triangulate_pointset();
338  //_outside_mesh_3D->print_info();
339 
340 
341  // (below generates some horribly expensive meshes,
342  // but seems immune to the 0 volume problem).
343  // _tetgen_inside->pointset_convexhull();
344  // _inside_mesh_3D->find_neighbors();
345  // _inside_mesh_3D->print_info();
346  // _tetgen_inside->triangulate_conformingDelaunayMesh (1.e3, 100.);
347  // _inside_mesh_3D->print_info();
348 
349  // _tetgen_outside->pointset_convexhull();
350  // _outside_mesh_3D->find_neighbors();
351  // _outside_mesh_3D->print_info();
352  // _tetgen_outside->triangulate_conformingDelaunayMesh (1.e3, 100.);
353  // _outside_mesh_3D->print_info();
354 
355  std::ostringstream name;
356 
357  name << "cut_cell_"
358  << cut_cntr++
359  << ".dat";
360  _inside_mesh_3D->write ("in_" + name.str());
361  _outside_mesh_3D->write ("out_" + name.str());
362 
363  // finally, add the elements to our lists.
364  _inside_elem.clear();
365  _outside_elem.clear();
366 
367  for (const auto & elem : _inside_mesh_3D->element_ptr_range())
368  if (elem->volume() > std::numeric_limits<Real>::epsilon())
369  _inside_elem.push_back (elem);
370 
371  for (const auto & elem : _outside_mesh_3D->element_ptr_range())
372  if (elem->volume() > std::numeric_limits<Real>::epsilon())
373  _outside_elem.push_back (elem);
374 
375 #endif
376 }
std::string name(const ElemQuality q)
Definition: elem_quality.C:42
std::unique_ptr< ReplicatedMesh > _inside_mesh_3D
Definition: elem_cutter.h:165
std::vector< Point > _intersection_pts
Definition: elem_cutter.h:171
std::vector< Elem const * > _inside_elem
Definition: elem_cutter.h:154
std::unique_ptr< ReplicatedMesh > _outside_mesh_3D
Definition: elem_cutter.h:168
OStreamProxy err(std::cerr)
std::unique_ptr< TetGenMeshInterface > _tetgen_inside
Definition: elem_cutter.h:166
std::unique_ptr< TetGenMeshInterface > _tetgen_outside
Definition: elem_cutter.h:169
std::vector< Elem const * > _outside_elem
Definition: elem_cutter.h:155

◆ find_intersection_points()

void libMesh::ElemCutter::find_intersection_points ( const Elem elem,
const std::vector< Real > &  vertex_distance_func 
)
protected

Finds the points where the cutting surface intersects the element edges.

Definition at line 154 of file elem_cutter.C.

References _intersection_pts, libMesh::Elem::build_edge_ptr(), libMesh::Elem::get_node_index(), libMesh::Elem::is_vertex(), libMesh::Elem::n_edges(), and libMesh::Real.

Referenced by operator()().

156 {
157  _intersection_pts.clear();
158 
159  for (unsigned int e=0; e<elem.n_edges(); e++)
160  {
161  std::unique_ptr<const Elem> edge (elem.build_edge_ptr(e));
162 
163  // find the element nodes el0, el1 that map
164  unsigned int
165  el0 = elem.get_node_index(edge->node_ptr(0)),
166  el1 = elem.get_node_index(edge->node_ptr(1));
167 
168  libmesh_assert (elem.is_vertex(el0));
169  libmesh_assert (elem.is_vertex(el1));
170  libmesh_assert_less (el0, vertex_distance_func.size());
171  libmesh_assert_less (el1, vertex_distance_func.size());
172 
173  const Real
174  d0 = vertex_distance_func[el0],
175  d1 = vertex_distance_func[el1];
176 
177  // if this egde has a 0 crossing
178  if (d0*d1 < 0.)
179  {
180  libmesh_assert_not_equal_to (d0, d1);
181 
182  // then find d_star in [0,1], the
183  // distance from el0 to el1 where the 0 lives.
184  const Real d_star = d0 / (d0 - d1);
185 
186 
187  // Prevent adding nodes trivially close to existing
188  // nodes.
189  const Real endpoint_tol = 0.01;
190 
191  if ( (d_star > endpoint_tol) &&
192  (d_star < (1.-endpoint_tol)) )
193  {
194  const Point x_star = (edge->point(0)*(1-d_star) +
195  edge->point(1)*d_star);
196 
197  std::cout << "adding cut point (d_star, x_star) = "
198  << d_star << " , " << x_star << std::endl;
199 
200  _intersection_pts.push_back (x_star);
201  }
202  }
203  }
204 }
std::vector< Point > _intersection_pts
Definition: elem_cutter.h:171
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

◆ inside_elements()

const std::vector<Elem const *>& libMesh::ElemCutter::inside_elements ( ) const
inline
Returns
A list of general element pieces considered inside the cutting surface. These are subelements whose geometric union defines the spatial domain of the inside portion of the cut element.

Definition at line 116 of file elem_cutter.h.

References _inside_elem.

117  { return _inside_elem; }
std::vector< Elem const * > _inside_elem
Definition: elem_cutter.h:154

◆ is_cut()

bool libMesh::ElemCutter::is_cut ( const Elem elem,
const std::vector< Real > &  vertex_distance_func 
) const
Returns
true if the element is cut by the interface defined implicitly by the vertex values of the signed vertex_distance_func.

Definition at line 88 of file elem_cutter.C.

References std::max(), std::min(), and libMesh::Real.

Referenced by operator()().

90 {
91  libmesh_assert_equal_to (elem.n_vertices(), vertex_distance_func.size());
92 
93  Real
94  vmin = vertex_distance_func.front(),
95  vmax = vmin;
96 
97  for (const auto & val : vertex_distance_func)
98  {
99  vmin = std::min (vmin, val);
100  vmax = std::max (vmax, val);
101  }
102 
103  // if the distance function changes sign, we're cut.
104  return (vmin*vmax < 0.);
105 }
long double max(long double a, double b)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
long double min(long double a, double b)

◆ is_inside()

bool libMesh::ElemCutter::is_inside ( const Elem elem,
const std::vector< Real > &  vertex_distance_func 
) const
Returns
true if the element is completely inside the interface defined implicitly by the vertex values of the signed vertex_distance_func.

Definition at line 58 of file elem_cutter.C.

Referenced by operator()().

60 {
61  libmesh_assert_equal_to (elem.n_vertices(), vertex_distance_func.size());
62 
63  for (const auto & val : vertex_distance_func)
64  if (val > 0.)
65  return false;
66 
67  // if the distance function is nonpositive, we are outside
68  return true;
69 }

◆ is_outside()

bool libMesh::ElemCutter::is_outside ( const Elem elem,
const std::vector< Real > &  vertex_distance_func 
) const
Returns
true if the element is completely outside the interface defined implicitly by the vertex values of the signed vertex_distance_func.

Definition at line 73 of file elem_cutter.C.

Referenced by operator()().

75 {
76  libmesh_assert_equal_to (elem.n_vertices(), vertex_distance_func.size());
77 
78  for (const auto & val : vertex_distance_func)
79  if (val < 0.)
80  return false;
81 
82  // if the distance function is nonnegative, we are outside
83  return true;
84 }

◆ operator()()

void libMesh::ElemCutter::operator() ( const Elem elem_in,
const std::vector< Real > &  vertex_distance_func 
)

This function implements cutting an element by a signed distance function. The input array vertex_distance_func contains the vertex values of a signed distance function, from which the cutting interface is inferred from the 0 level set. If all vertex values are positive, the element is outside the cutting surface and is not cut. Likewise if all vertex values are negative, the element is inside the cutting surface and is not cut.

Definition at line 109 of file elem_cutter.C.

References _inside_elem, _outside_elem, cut_1D(), cut_2D(), cut_3D(), libMesh::Elem::dim(), find_intersection_points(), is_cut(), is_inside(), is_outside(), and libMesh::Elem::n_vertices().

112 {
113  libmesh_assert_equal_to (vertex_distance_func.size(), elem.n_vertices());
114 
115  _inside_elem.clear();
116  _outside_elem.clear();
117 
118  // check for quick return?
119  {
120  // completely outside?
121  if (this->is_outside(elem, vertex_distance_func))
122  {
123  //std::cout << "element completely outside\n";
124  _outside_elem.push_back(& elem);
125  return;
126  }
127 
128  // completely inside?
129  else if (this->is_inside(elem, vertex_distance_func))
130  {
131  //std::cout << "element completely inside\n";
132  _inside_elem.push_back(&elem);
133  return;
134  }
135 
136  libmesh_assert (this->is_cut (elem, vertex_distance_func));
137  }
138 
139  // we now know we are in a cut element, find the intersecting points.
140  this->find_intersection_points (elem, vertex_distance_func);
141 
142  // and then dispatch the proper method
143  switch (elem.dim())
144  {
145  case 1: this->cut_1D(elem, vertex_distance_func); break;
146  case 2: this->cut_2D(elem, vertex_distance_func); break;
147  case 3: this->cut_3D(elem, vertex_distance_func); break;
148  default: libmesh_error_msg("Invalid element dimension: " << elem.dim());
149  }
150 }
void find_intersection_points(const Elem &elem, const std::vector< Real > &vertex_distance_func)
Definition: elem_cutter.C:154
void cut_1D(const Elem &elem, const std::vector< Real > &vertex_distance_func)
Definition: elem_cutter.C:209
bool is_cut(const Elem &elem, const std::vector< Real > &vertex_distance_func) const
Definition: elem_cutter.C:88
void cut_2D(const Elem &elem, const std::vector< Real > &vertex_distance_func)
Definition: elem_cutter.C:217
std::vector< Elem const * > _inside_elem
Definition: elem_cutter.h:154
void cut_3D(const Elem &elem, const std::vector< Real > &vertex_distance_func)
Definition: elem_cutter.C:297
bool is_outside(const Elem &elem, const std::vector< Real > &vertex_distance_func) const
Definition: elem_cutter.C:73
bool is_inside(const Elem &elem, const std::vector< Real > &vertex_distance_func) const
Definition: elem_cutter.C:58
std::vector< Elem const * > _outside_elem
Definition: elem_cutter.h:155

◆ outside_elements()

const std::vector<Elem const *>& libMesh::ElemCutter::outside_elements ( ) const
inline
Returns
A list of general element pieces considered outside the cutting surface. These are subelements whose geometric union defines the spatial domain of the outside portion of the cut element.

Definition at line 124 of file elem_cutter.h.

References _outside_elem.

125  { return _outside_elem; }
std::vector< Elem const * > _outside_elem
Definition: elem_cutter.h:155

Member Data Documentation

◆ _comm_self

Parallel::Communicator libMesh::ElemCutter::_comm_self
protected

Definition at line 157 of file elem_cutter.h.

◆ _inside_elem

std::vector<Elem const *> libMesh::ElemCutter::_inside_elem
protected

Definition at line 154 of file elem_cutter.h.

Referenced by cut_2D(), cut_3D(), inside_elements(), and operator()().

◆ _inside_mesh_2D

std::unique_ptr<ReplicatedMesh> libMesh::ElemCutter::_inside_mesh_2D
protected

Definition at line 159 of file elem_cutter.h.

Referenced by cut_2D().

◆ _inside_mesh_3D

std::unique_ptr<ReplicatedMesh> libMesh::ElemCutter::_inside_mesh_3D
protected

Definition at line 165 of file elem_cutter.h.

Referenced by cut_3D().

◆ _intersection_pts

std::vector<Point> libMesh::ElemCutter::_intersection_pts
protected

Definition at line 171 of file elem_cutter.h.

Referenced by cut_2D(), cut_3D(), and find_intersection_points().

◆ _outside_elem

std::vector<Elem const *> libMesh::ElemCutter::_outside_elem
protected

Definition at line 155 of file elem_cutter.h.

Referenced by cut_2D(), cut_3D(), operator()(), and outside_elements().

◆ _outside_mesh_2D

std::unique_ptr<ReplicatedMesh> libMesh::ElemCutter::_outside_mesh_2D
protected

Definition at line 162 of file elem_cutter.h.

Referenced by cut_2D().

◆ _outside_mesh_3D

std::unique_ptr<ReplicatedMesh> libMesh::ElemCutter::_outside_mesh_3D
protected

Definition at line 168 of file elem_cutter.h.

Referenced by cut_3D().

◆ _tetgen_inside

std::unique_ptr<TetGenMeshInterface> libMesh::ElemCutter::_tetgen_inside
protected

Definition at line 166 of file elem_cutter.h.

Referenced by cut_3D().

◆ _tetgen_outside

std::unique_ptr<TetGenMeshInterface> libMesh::ElemCutter::_tetgen_outside
protected

Definition at line 169 of file elem_cutter.h.

Referenced by cut_3D().

◆ _triangle_inside

std::unique_ptr<TriangleInterface> libMesh::ElemCutter::_triangle_inside
protected

Definition at line 160 of file elem_cutter.h.

Referenced by cut_2D().

◆ _triangle_outside

std::unique_ptr<TriangleInterface> libMesh::ElemCutter::_triangle_outside
protected

Definition at line 163 of file elem_cutter.h.

Referenced by cut_2D().


The documentation for this class was generated from the following files: