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

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
libMesh::ElemCutter::~ElemCutter ( )

Destructor.

Definition at line 53 of file elem_cutter.C.

54 {}

Member Function Documentation

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

cutting algorithm in 1D.

Definition at line 210 of file elem_cutter.C.

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

212 {
213  libmesh_not_implemented();
214 }
void libMesh::ElemCutter::cut_2D ( const Elem elem,
const std::vector< Real > &  vertex_distance_func 
)
protected

cutting algorithm in 2D.

Definition at line 218 of file elem_cutter.C.

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

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

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

cutting algorithm in 3D.

Definition at line 307 of file elem_cutter.C.

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

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

309 {
310 #ifndef LIBMESH_HAVE_TETGEN
311 
312  // current implementation requires tetgen!
313  libMesh::err << "ERROR: current libMesh ElemCutter 3D implementation requires\n"
314  << " the \"tetgen\" library!\n"
315  << std::endl;
316  libmesh_not_implemented();
317 
318 #else // OK, LIBMESH_HAVE_TETGEN
319 
320  std::cout << "Inside cut cell element!\n";
321 
322  libmesh_assert (_inside_mesh_3D.get() != libmesh_nullptr);
323  libmesh_assert (_outside_mesh_3D.get() != libmesh_nullptr);
324 
325  _inside_mesh_3D->clear();
326  _outside_mesh_3D->clear();
327 
328  for (unsigned int v=0; v<elem.n_vertices(); v++)
329  {
330  if (vertex_distance_func[v] >= 0.)
331  _outside_mesh_3D->add_point (elem.point(v));
332 
333  if (vertex_distance_func[v] <= 0.)
334  _inside_mesh_3D->add_point (elem.point(v));
335  }
336 
337  for (std::vector<Point>::const_iterator it=_intersection_pts.begin();
338  it != _intersection_pts.end(); ++it)
339  {
340  _inside_mesh_3D->add_point(*it);
341  _outside_mesh_3D->add_point(*it);
342  }
343 
344 
345  // Triangulate!
346  _tetgen_inside->triangulate_pointset();
347  //_inside_mesh_3D->print_info();
348  _tetgen_outside->triangulate_pointset();
349  //_outside_mesh_3D->print_info();
350 
351 
352  // (below generates some horribly expensive meshes,
353  // but seems immune to the 0 volume problem).
354  // _tetgen_inside->pointset_convexhull();
355  // _inside_mesh_3D->find_neighbors();
356  // _inside_mesh_3D->print_info();
357  // _tetgen_inside->triangulate_conformingDelaunayMesh (1.e3, 100.);
358  // _inside_mesh_3D->print_info();
359 
360  // _tetgen_outside->pointset_convexhull();
361  // _outside_mesh_3D->find_neighbors();
362  // _outside_mesh_3D->print_info();
363  // _tetgen_outside->triangulate_conformingDelaunayMesh (1.e3, 100.);
364  // _outside_mesh_3D->print_info();
365 
366  std::ostringstream name;
367 
368  name << "cut_cell_"
369  << cut_cntr++
370  << ".dat";
371  _inside_mesh_3D->write ("in_" + name.str());
372  _outside_mesh_3D->write ("out_" + name.str());
373 
374  // finally, add the elements to our lists.
375  {
376  _inside_elem.clear(); _outside_elem.clear();
377 
378  MeshBase::const_element_iterator
379  it = _inside_mesh_3D->elements_begin(),
380  end = _inside_mesh_3D->elements_end();
381 
382  for (; it!=end; ++it)
383  if ((*it)->volume() > std::numeric_limits<Real>::epsilon())
384  _inside_elem.push_back (*it);
385 
386  it = _outside_mesh_3D->elements_begin();
387  end = _outside_mesh_3D->elements_end();
388 
389  for (; it!=end; ++it)
390  if ((*it)->volume() > std::numeric_limits<Real>::epsilon())
391  _outside_elem.push_back (*it);
392  }
393 
394 #endif
395 }
std::string name(const ElemQuality q)
Definition: elem_quality.C:39
std::unique_ptr< ReplicatedMesh > _inside_mesh_3D
Definition: elem_cutter.h:165
std::vector< Point > _intersection_pts
Definition: elem_cutter.h:171
const class libmesh_nullptr_t libmesh_nullptr
IterBase * end
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
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 155 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()(), and outside_elements().

157 {
158  _intersection_pts.clear();
159 
160  for (unsigned int e=0; e<elem.n_edges(); e++)
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 egde 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  std::cout << "adding cut point (d_star, x_star) = "
199  << d_star << " , " << x_star << std::endl;
200 
201  _intersection_pts.push_back (x_star);
202  }
203  }
204  }
205 }
std::vector< Point > _intersection_pts
Definition: elem_cutter.h:171
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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.

Referenced by libMesh::QComposite< QSubCell >::init().

117  { return _inside_elem; }
std::vector< Elem const * > _inside_elem
Definition: elem_cutter.h:154
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 libMesh::QComposite< QSubCell >::init(), and 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 (std::vector<Real>::const_iterator it=vertex_distance_func.begin();
98  it!=vertex_distance_func.end(); ++it)
99  {
100  vmin = std::min (vmin, *it);
101  vmax = std::max (vmax, *it);
102  }
103 
104  // if the distance function changes sign, we're cut.
105  return (vmin*vmax < 0.);
106 }
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)
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 (std::vector<Real>::const_iterator it=vertex_distance_func.begin();
64  it!=vertex_distance_func.end(); ++it)
65  if (*it > 0.) return false;
66 
67  // if the distance function is nonpositive, we are outside
68  return true;
69 }
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 (std::vector<Real>::const_iterator it=vertex_distance_func.begin();
79  it!=vertex_distance_func.end(); ++it)
80  if (*it < 0.) return false;
81 
82  // if the distance function is nonnegative, we are outside
83  return true;
84 }
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 110 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().

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 }
void find_intersection_points(const Elem &elem, const std::vector< Real > &vertex_distance_func)
Definition: elem_cutter.C:155
void cut_1D(const Elem &elem, const std::vector< Real > &vertex_distance_func)
Definition: elem_cutter.C:210
bool is_outside(const Elem &elem, const std::vector< Real > &vertex_distance_func) const
Definition: elem_cutter.C:73
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:218
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:307
std::vector< Elem const * > _outside_elem
Definition: elem_cutter.h:155
bool is_inside(const Elem &elem, const std::vector< Real > &vertex_distance_func) const
Definition: elem_cutter.C:58
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, cut_1D(), cut_2D(), cut_3D(), and find_intersection_points().

Referenced by libMesh::QComposite< QSubCell >::init().

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

Member Data Documentation

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

Definition at line 157 of file elem_cutter.h.

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()().

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

Definition at line 159 of file elem_cutter.h.

Referenced by cut_2D().

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

Definition at line 165 of file elem_cutter.h.

Referenced by cut_3D().

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().

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().

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

Definition at line 162 of file elem_cutter.h.

Referenced by cut_2D().

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

Definition at line 168 of file elem_cutter.h.

Referenced by cut_3D().

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

Definition at line 166 of file elem_cutter.h.

Referenced by cut_3D().

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

Definition at line 169 of file elem_cutter.h.

Referenced by cut_3D().

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

Definition at line 160 of file elem_cutter.h.

Referenced by cut_2D().

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: