#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
UniquePtr< SerialMesh >	_inside_mesh_2D
UniquePtr< SerialMesh >	_outside_mesh_2D
UniquePtr< SerialMesh >	_inside_mesh_3D
UniquePtr< SerialMesh >	_outside_mesh_3D
Parallel::Communicator	_comm_self
UniquePtr< TriangleInterface >	_triangle_inside
UniquePtr< TriangleInterface >	_triangle_outside
UniquePtr< TetGenMeshInterface >	_tetgen_inside
UniquePtr< TetGenMeshInterface >	_tetgen_outside
std::vector< Point >	_intersection_pts

Detailed Description

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, 2013

Definition at line 58 of file elem_cutter.h.

Constructor & Destructor Documentation

libMesh::ElemCutter::ElemCutter ( )

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

Definition at line 41 of file elem_cutter.C.

References _comm_self, _inside_mesh_2D, _inside_mesh_3D, _outside_mesh_2D, _outside_mesh_3D, _tetgen_inside, _tetgen_outside, _triangle_inside, and _triangle_outside.

{
  _inside_mesh_2D.reset  (new SerialMesh(_comm_self,2));  _triangle_inside.reset  (new TriangleInterface (*_inside_mesh_2D));
  _outside_mesh_2D.reset (new SerialMesh(_comm_self,2));  _triangle_outside.reset (new TriangleInterface (*_outside_mesh_2D));

  _inside_mesh_3D.reset  (new SerialMesh(_comm_self,3));  _tetgen_inside.reset  (new TetGenMeshInterface (*_inside_mesh_3D));
  _outside_mesh_3D.reset (new SerialMesh(_comm_self,3));  _tetgen_outside.reset (new TetGenMeshInterface (*_outside_mesh_3D));

  cut_cntr = 0;
}

libMesh::ElemCutter::~ElemCutter ( )

Destructor.

Definition at line 54 of file elem_cutter.C.

{}

Member Function Documentation

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

cutting algoritm in 1D.

Definition at line 211 of file elem_cutter.C.

Referenced by operator()().

{
  libmesh_not_implemented();
}

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

cutting algoritm in 2D.

Definition at line 219 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::libmesh_assert(), libMesh::Elem::n_vertices(), and libMesh::Elem::point().

Referenced by operator()().

{
#ifndef LIBMESH_HAVE_TRIANGLE

  // current implementation requires triangle!
  libMesh::err << "ERROR: current libMesh ElemCutter 2D implementation requires\n"
               << "       the \"triangle\" library!\n"
               << std::endl;
  libmesh_not_implemented();

#else // OK, LIBMESH_HAVE_TRIANGLE

  std::cout << "Inside cut face element!\n";

  libmesh_assert (_inside_mesh_2D.get()  != NULL);
  libmesh_assert (_outside_mesh_2D.get() != NULL);

  _inside_mesh_2D->clear();
  _outside_mesh_2D->clear();

  for (unsigned int v=0; v<elem.n_vertices(); v++)
    {
      if (vertex_distance_func[v] >= 0.)
        _outside_mesh_2D->add_point (elem.point(v));

      if (vertex_distance_func[v] <= 0.)
        _inside_mesh_2D->add_point (elem.point(v));
    }

  for (std::vector<Point>::const_iterator it=_intersection_pts.begin();
       it != _intersection_pts.end(); ++it)
    {
      _inside_mesh_2D->add_point(*it);
      _outside_mesh_2D->add_point(*it);
    }


  // Customize the variables for the triangulation
  // we will be cutting reference cell, and want as few triangles
  // as possible, so jack this up larger than the area we will be
  // triangulating so we are governed only by accurately defining
  // the boundaries.
  _triangle_inside->desired_area()  = 100.;
  _triangle_outside->desired_area() = 100.;

  // allow for small angles
  _triangle_inside->minimum_angle()  = 5.;
  _triangle_outside->minimum_angle() = 5.;

  // Turn off Laplacian mesh smoothing after generation.
  _triangle_inside->smooth_after_generating()  = false;
  _triangle_outside->smooth_after_generating() = false;

  // Triangulate!
  _triangle_inside->triangulate();
  _triangle_outside->triangulate();

  // std::ostringstream name;

  // name << "cut_face_"
  //  << cut_cntr++
  //  << ".dat";
  // _inside_mesh_2D->write  ("in_"  + name.str());
  // _outside_mesh_2D->write ("out_" + name.str());

  // finally, add the elements to our lists.
  {
    _inside_elem.clear();  _outside_elem.clear();

    MeshBase::const_element_iterator
      it  = _inside_mesh_2D->elements_begin(),
      end = _inside_mesh_2D->elements_end();

    for (; it!=end; ++it)
      _inside_elem.push_back (*it);

    it  = _outside_mesh_2D->elements_begin();
    end = _outside_mesh_2D->elements_end();

    for (; it!=end; ++it)
      _outside_elem.push_back (*it);
  }

#endif
}

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

cutting algoritm in 3D.

Definition at line 308 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::libmesh_assert(), libMesh::Elem::n_vertices(), libMesh::Quality::name(), and libMesh::Elem::point().

Referenced by operator()().

{
#ifndef LIBMESH_HAVE_TETGEN

  // current implementation requires tetgen!
  libMesh::err << "ERROR: current libMesh ElemCutter 3D implementation requires\n"
               << "       the \"tetgen\" library!\n"
               << std::endl;
  libmesh_not_implemented();

#else // OK, LIBMESH_HAVE_TETGEN

  std::cout << "Inside cut cell element!\n";

  libmesh_assert (_inside_mesh_3D.get()  != NULL);
  libmesh_assert (_outside_mesh_3D.get() != NULL);

  _inside_mesh_3D->clear();
  _outside_mesh_3D->clear();

  for (unsigned int v=0; v<elem.n_vertices(); v++)
    {
      if (vertex_distance_func[v] >= 0.)
        _outside_mesh_3D->add_point (elem.point(v));

      if (vertex_distance_func[v] <= 0.)
        _inside_mesh_3D->add_point (elem.point(v));
    }

  for (std::vector<Point>::const_iterator it=_intersection_pts.begin();
       it != _intersection_pts.end(); ++it)
    {
      _inside_mesh_3D->add_point(*it);
      _outside_mesh_3D->add_point(*it);
    }


  // Triangulate!
  _tetgen_inside->triangulate_pointset();
  //_inside_mesh_3D->print_info();
  _tetgen_outside->triangulate_pointset();
  //_outside_mesh_3D->print_info();


  // (below generates some horribly expensive meshes,
  //  but seems immune to the 0 volume problem).
  // _tetgen_inside->pointset_convexhull();
  // _inside_mesh_3D->find_neighbors();
  // _inside_mesh_3D->print_info();
  // _tetgen_inside->triangulate_conformingDelaunayMesh (1.e3, 100.);
  // _inside_mesh_3D->print_info();

  // _tetgen_outside->pointset_convexhull();
  // _outside_mesh_3D->find_neighbors();
  // _outside_mesh_3D->print_info();
  // _tetgen_outside->triangulate_conformingDelaunayMesh (1.e3, 100.);
  // _outside_mesh_3D->print_info();

  std::ostringstream name;

  name << "cut_cell_"
       << cut_cntr++
       << ".dat";
  _inside_mesh_3D->write  ("in_"  + name.str());
  _outside_mesh_3D->write ("out_" + name.str());

  // finally, add the elements to our lists.
  {
    _inside_elem.clear();  _outside_elem.clear();

    MeshBase::const_element_iterator
      it  = _inside_mesh_3D->elements_begin(),
      end = _inside_mesh_3D->elements_end();

    for (; it!=end; ++it)
      if ((*it)->volume() > std::numeric_limits<Real>::epsilon())
        _inside_elem.push_back (*it);

    it  = _outside_mesh_3D->elements_begin();
    end = _outside_mesh_3D->elements_end();

    for (; it!=end; ++it)
      if ((*it)->volume() > std::numeric_limits<Real>::epsilon())
        _outside_elem.push_back (*it);
  }

#endif
}

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 156 of file elem_cutter.C.

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

Referenced by operator()().

{
  _intersection_pts.clear();

  for (unsigned int e=0; e<elem.n_edges(); e++)
    {
      UniquePtr<Elem> edge (elem.build_edge(e));

      // find the element nodes el0, el1 that map
      unsigned int
        el0 = elem.get_node_index(edge->get_node(0)),
        el1 = elem.get_node_index(edge->get_node(1));

      libmesh_assert (elem.is_vertex(el0));
      libmesh_assert (elem.is_vertex(el1));
      libmesh_assert_less (el0, vertex_distance_func.size());
      libmesh_assert_less (el1, vertex_distance_func.size());

      const Real
        d0 = vertex_distance_func[el0],
        d1 = vertex_distance_func[el1];

      // if this egde has a 0 crossing
      if (d0*d1 < 0.)
        {
          libmesh_assert_not_equal_to (d0, d1);

          // then find d_star in [0,1], the
          // distance from el0 to el1 where the 0 lives.
          const Real d_star = d0 / (d0 - d1);


          // Prevent adding nodes trivially close to existing
          // nodes.
          const Real endpoint_tol = 0.01;

          if ( (d_star > endpoint_tol) &&
               (d_star < (1.-endpoint_tol)) )
            {
              const Point x_star = (edge->point(0)*(1.-d_star) +
                                    edge->point(1)*d_star);

              std::cout << "adding cut point (d_star, x_star) = "
                        << d_star << " , " << x_star << std::endl;

              _intersection_pts.push_back (x_star);
            }
        }
    }
}

const std::vector<Elem const*>& libMesh::ElemCutter::inside_elements ( ) const [inline]

Returns a list of in 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 114 of file elem_cutter.h.

References _inside_elem.

  { return _inside_elem; }

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 89 of file elem_cutter.C.

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

Referenced by operator()().

{
  libmesh_assert_equal_to (elem.n_vertices(), vertex_distance_func.size());

  Real
    vmin = vertex_distance_func.front(),
    vmax = vmin;

  for (std::vector<Real>::const_iterator it=vertex_distance_func.begin();
       it!=vertex_distance_func.end(); ++it)
    {
      vmin = std::min (vmin, *it);
      vmax = std::max (vmax, *it);
    }

  // if the distance function changes sign, we're cut.
  return (vmin*vmax < 0.);
}

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 59 of file elem_cutter.C.

References libMesh::Elem::n_vertices().

Referenced by operator()().

{
  libmesh_assert_equal_to (elem.n_vertices(), vertex_distance_func.size());

  for (std::vector<Real>::const_iterator it=vertex_distance_func.begin();
       it!=vertex_distance_func.end(); ++it)
    if (*it > 0.) return false;

  // if the distance function is nonpositive, we are outside
  return true;
}

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 74 of file elem_cutter.C.

References libMesh::Elem::n_vertices().

Referenced by operator()().

{
  libmesh_assert_equal_to (elem.n_vertices(), vertex_distance_func.size());

  for (std::vector<Real>::const_iterator it=vertex_distance_func.begin();
       it!=vertex_distance_func.end(); ++it)
    if (*it < 0.) return false;

  // if the distance function is nonnegative, we are outside
  return true;
}

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 111 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(), libMesh::libmesh_assert(), and libMesh::Elem::n_vertices().

{
  libmesh_assert_equal_to (vertex_distance_func.size(), elem.n_vertices());

  _inside_elem.clear();
  _outside_elem.clear();

  // check for quick return?
  {
    // completely outside?
    if (this->is_outside(elem, vertex_distance_func))
      {
        //std::cout << "element completely outside\n";
        _outside_elem.push_back(&elem);
        return;
      }

    // completely inside?
    else if (this->is_inside(elem, vertex_distance_func))
      {
        //std::cout << "element completely inside\n";
        _inside_elem.push_back(&elem);
        return;
      }

    libmesh_assert (this->is_cut (elem, vertex_distance_func));
  }

  // we now know we are in a cut element, find the intersecting points.
  this->find_intersection_points (elem, vertex_distance_func);

  // and then dispatch the proper method
  switch (elem.dim())
    {
    case 1: this->cut_1D(elem, vertex_distance_func); break;
    case 2: this->cut_2D(elem, vertex_distance_func); break;
    case 3: this->cut_3D(elem, vertex_distance_func); break;
    default: libmesh_error_msg("Invalid element dimension: " << elem.dim());
    }
}

const std::vector<Elem const*>& libMesh::ElemCutter::outside_elements ( ) const [inline]

Returns a list of in 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 122 of file elem_cutter.h.

References _outside_elem.