libMesh
elem_cutter.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2013 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 
31 namespace
32 {
33 unsigned int cut_cntr;
34 }
35 
36 namespace libMesh
37 {
38 
39 // ------------------------------------------------------------
40 // ElemCutter implementation
42 {
45 
48 
49  cut_cntr = 0;
50 }
51 
52 
53 
55 {}
56 
57 
58 
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 (std::vector<Real>::const_iterator it=vertex_distance_func.begin();
65  it!=vertex_distance_func.end(); ++it)
66  if (*it > 0.) return false;
67 
68  // if the distance function is nonpositive, we are outside
69  return true;
70 }
71 
72 
73 
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 (std::vector<Real>::const_iterator it=vertex_distance_func.begin();
80  it!=vertex_distance_func.end(); ++it)
81  if (*it < 0.) return false;
82 
83  // if the distance function is nonnegative, we are outside
84  return true;
85 }
86 
87 
88 
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 (std::vector<Real>::const_iterator it=vertex_distance_func.begin();
99  it!=vertex_distance_func.end(); ++it)
100  {
101  vmin = std::min (vmin, *it);
102  vmax = std::max (vmax, *it);
103  }
104 
105  // if the distance function changes sign, we're cut.
106  return (vmin*vmax < 0.);
107 }
108 
109 
110 
111 void ElemCutter::operator()(const Elem & elem,
112  const std::vector<Real> & vertex_distance_func)
113 
114 {
115  libmesh_assert_equal_to (vertex_distance_func.size(), elem.n_vertices());
116 
117  _inside_elem.clear();
118  _outside_elem.clear();
119 
120  // check for quick return?
121  {
122  // completely outside?
123  if (this->is_outside(elem, vertex_distance_func))
124  {
125  //std::cout << "element completely outside\n";
126  _outside_elem.push_back(& elem);
127  return;
128  }
129 
130  // completely inside?
131  else if (this->is_inside(elem, vertex_distance_func))
132  {
133  //std::cout << "element completely inside\n";
134  _inside_elem.push_back(&elem);
135  return;
136  }
137 
138  libmesh_assert (this->is_cut (elem, vertex_distance_func));
139  }
140 
141  // we now know we are in a cut element, find the intersecting points.
142  this->find_intersection_points (elem, vertex_distance_func);
143 
144  // and then dispatch the proper method
145  switch (elem.dim())
146  {
147  case 1: this->cut_1D(elem, vertex_distance_func); break;
148  case 2: this->cut_2D(elem, vertex_distance_func); break;
149  case 3: this->cut_3D(elem, vertex_distance_func); break;
150  default: libmesh_error_msg("Invalid element dimension: " << elem.dim());
151  }
152 }
153 
154 
155 
157  const std::vector<Real> & vertex_distance_func)
158 {
159  _intersection_pts.clear();
160 
161  for (unsigned int e=0; e<elem.n_edges(); e++)
162  {
163  UniquePtr<const Elem> edge (elem.build_edge_ptr(e));
164 
165  // find the element nodes el0, el1 that map
166  unsigned int
167  el0 = elem.get_node_index(edge->node_ptr(0)),
168  el1 = elem.get_node_index(edge->node_ptr(1));
169 
170  libmesh_assert (elem.is_vertex(el0));
171  libmesh_assert (elem.is_vertex(el1));
172  libmesh_assert_less (el0, vertex_distance_func.size());
173  libmesh_assert_less (el1, vertex_distance_func.size());
174 
175  const Real
176  d0 = vertex_distance_func[el0],
177  d1 = vertex_distance_func[el1];
178 
179  // if this egde has a 0 crossing
180  if (d0*d1 < 0.)
181  {
182  libmesh_assert_not_equal_to (d0, d1);
183 
184  // then find d_star in [0,1], the
185  // distance from el0 to el1 where the 0 lives.
186  const Real d_star = d0 / (d0 - d1);
187 
188 
189  // Prevent adding nodes trivially close to existing
190  // nodes.
191  const Real endpoint_tol = 0.01;
192 
193  if ( (d_star > endpoint_tol) &&
194  (d_star < (1.-endpoint_tol)) )
195  {
196  const Point x_star = (edge->point(0)*(1-d_star) +
197  edge->point(1)*d_star);
198 
199  std::cout << "adding cut point (d_star, x_star) = "
200  << d_star << " , " << x_star << std::endl;
201 
202  _intersection_pts.push_back (x_star);
203  }
204  }
205  }
206 }
207 
208 
209 
210 
211 void ElemCutter::cut_1D (const Elem & /*elem*/,
212  const std::vector<Real> &/*vertex_distance_func*/)
213 {
214  libmesh_not_implemented();
215 }
216 
217 
218 
219 void ElemCutter::cut_2D (const Elem & elem,
220  const std::vector<Real> & vertex_distance_func)
221 {
222 #ifndef LIBMESH_HAVE_TRIANGLE
223 
224  // current implementation requires triangle!
225  libMesh::err << "ERROR: current libMesh ElemCutter 2D implementation requires\n"
226  << " the \"triangle\" library!\n"
227  << std::endl;
228  libmesh_not_implemented();
229 
230 #else // OK, LIBMESH_HAVE_TRIANGLE
231 
232  std::cout << "Inside cut face element!\n";
233 
236 
237  _inside_mesh_2D->clear();
238  _outside_mesh_2D->clear();
239 
240  for (unsigned int v=0; v<elem.n_vertices(); v++)
241  {
242  if (vertex_distance_func[v] >= 0.)
243  _outside_mesh_2D->add_point (elem.point(v));
244 
245  if (vertex_distance_func[v] <= 0.)
246  _inside_mesh_2D->add_point (elem.point(v));
247  }
248 
249  for (std::vector<Point>::const_iterator it=_intersection_pts.begin();
250  it != _intersection_pts.end(); ++it)
251  {
252  _inside_mesh_2D->add_point(*it);
253  _outside_mesh_2D->add_point(*it);
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  {
287  _inside_elem.clear(); _outside_elem.clear();
288 
290  it = _inside_mesh_2D->elements_begin(),
291  end = _inside_mesh_2D->elements_end();
292 
293  for (; it!=end; ++it)
294  _inside_elem.push_back (*it);
295 
296  it = _outside_mesh_2D->elements_begin();
297  end = _outside_mesh_2D->elements_end();
298 
299  for (; it!=end; ++it)
300  _outside_elem.push_back (*it);
301  }
302 
303 #endif
304 }
305 
306 
307 
308 void ElemCutter::cut_3D (const Elem & elem,
309  const std::vector<Real> & vertex_distance_func)
310 {
311 #ifndef LIBMESH_HAVE_TETGEN
312 
313  // current implementation requires tetgen!
314  libMesh::err << "ERROR: current libMesh ElemCutter 3D implementation requires\n"
315  << " the \"tetgen\" library!\n"
316  << std::endl;
317  libmesh_not_implemented();
318 
319 #else // OK, LIBMESH_HAVE_TETGEN
320 
321  std::cout << "Inside cut cell element!\n";
322 
325 
326  _inside_mesh_3D->clear();
327  _outside_mesh_3D->clear();
328 
329  for (unsigned int v=0; v<elem.n_vertices(); v++)
330  {
331  if (vertex_distance_func[v] >= 0.)
332  _outside_mesh_3D->add_point (elem.point(v));
333 
334  if (vertex_distance_func[v] <= 0.)
335  _inside_mesh_3D->add_point (elem.point(v));
336  }
337 
338  for (std::vector<Point>::const_iterator it=_intersection_pts.begin();
339  it != _intersection_pts.end(); ++it)
340  {
341  _inside_mesh_3D->add_point(*it);
342  _outside_mesh_3D->add_point(*it);
343  }
344 
345 
346  // Triangulate!
347  _tetgen_inside->triangulate_pointset();
348  //_inside_mesh_3D->print_info();
349  _tetgen_outside->triangulate_pointset();
350  //_outside_mesh_3D->print_info();
351 
352 
353  // (below generates some horribly expensive meshes,
354  // but seems immune to the 0 volume problem).
355  // _tetgen_inside->pointset_convexhull();
356  // _inside_mesh_3D->find_neighbors();
357  // _inside_mesh_3D->print_info();
358  // _tetgen_inside->triangulate_conformingDelaunayMesh (1.e3, 100.);
359  // _inside_mesh_3D->print_info();
360 
361  // _tetgen_outside->pointset_convexhull();
362  // _outside_mesh_3D->find_neighbors();
363  // _outside_mesh_3D->print_info();
364  // _tetgen_outside->triangulate_conformingDelaunayMesh (1.e3, 100.);
365  // _outside_mesh_3D->print_info();
366 
367  std::ostringstream name;
368 
369  name << "cut_cell_"
370  << cut_cntr++
371  << ".dat";
372  _inside_mesh_3D->write ("in_" + name.str());
373  _outside_mesh_3D->write ("out_" + name.str());
374 
375  // finally, add the elements to our lists.
376  {
377  _inside_elem.clear(); _outside_elem.clear();
378 
380  it = _inside_mesh_3D->elements_begin(),
381  end = _inside_mesh_3D->elements_end();
382 
383  for (; it!=end; ++it)
384  if ((*it)->volume() > std::numeric_limits<Real>::epsilon())
385  _inside_elem.push_back (*it);
386 
387  it = _outside_mesh_3D->elements_begin();
388  end = _outside_mesh_3D->elements_end();
389 
390  for (; it!=end; ++it)
391  if ((*it)->volume() > std::numeric_limits<Real>::epsilon())
392  _outside_elem.push_back (*it);
393  }
394 
395 #endif
396 }
397 
398 
399 
400 } // namespace libMesh
401 
402 #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:39
OStreamProxy err
Class TetGenMeshInterface provides an interface for tetrahedralization of meshes using the TetGen lib...
UniquePtr< ReplicatedMesh > _outside_mesh_3D
Definition: elem_cutter.h:159
UniquePtr< ReplicatedMesh > _inside_mesh_3D
Definition: elem_cutter.h:158
UniquePtr< TriangleInterface > _triangle_outside
Definition: elem_cutter.h:164
The ReplicatedMesh class is derived from the MeshBase class, and is used to store identical copies of...
UniquePtr< ReplicatedMesh > _inside_mesh_2D
Definition: elem_cutter.h:156
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:156
std::vector< Point > _intersection_pts
Definition: elem_cutter.h:168
void cut_1D(const Elem &elem, const std::vector< Real > &vertex_distance_func)
cutting algorithm in 1D.
Definition: elem_cutter.C:211
ElemCutter()
Constructor.
Definition: elem_cutter.C:41
The definition of the const_element_iterator struct.
Definition: mesh_base.h:1494
virtual unsigned int n_edges() const =0
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:111
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
bool is_outside(const Elem &elem, const std::vector< Real > &vertex_distance_func) const
Definition: elem_cutter.C:74
const class libmesh_nullptr_t libmesh_nullptr
IterBase * end
Also have a polymorphic pointer to the end object, this prevents iterating past the end...
bool is_cut(const Elem &elem, const std::vector< Real > &vertex_distance_func) const
Definition: elem_cutter.C:89
The libMesh namespace provides an interface to certain functionality in the library.
long double max(long double a, double b)
UniquePtr< ReplicatedMesh > _outside_mesh_2D
Definition: elem_cutter.h:157
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
void cut_2D(const Elem &elem, const std::vector< Real > &vertex_distance_func)
cutting algorithm in 2D.
Definition: elem_cutter.C:219
UniquePtr< TetGenMeshInterface > _tetgen_outside
Definition: elem_cutter.h:166
std::vector< Elem const * > _inside_elem
Definition: elem_cutter.h:153
void cut_3D(const Elem &elem, const std::vector< Real > &vertex_distance_func)
cutting algorithm in 3D.
Definition: elem_cutter.C:308
PetscErrorCode Vec Mat libmesh_dbg_var(j)
UniquePtr< TetGenMeshInterface > _tetgen_inside
Definition: elem_cutter.h:165
std::vector< Elem const * > _outside_elem
Definition: elem_cutter.h:154
~ElemCutter()
Destructor.
Definition: elem_cutter.C:54
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Point & point(const unsigned int i) const
Definition: elem.h:1809
virtual bool is_vertex(const unsigned int i) const =0
UniquePtr< TriangleInterface > _triangle_inside
Definition: elem_cutter.h:163
virtual unsigned int n_vertices() const =0
virtual unsigned int dim() const =0
virtual UniquePtr< Elem > build_edge_ptr(const unsigned int i)=0
bool is_inside(const Elem &elem, const std::vector< Real > &vertex_distance_func) const
Definition: elem_cutter.C:59
A C++ interface between LibMesh and the Triangle library written by J.R.
long double min(long double a, double b)
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
Parallel::Communicator _comm_self
Definition: elem_cutter.h:161
unsigned int get_node_index(const Node *node_ptr) const
Definition: elem.h:1929