triangulationwrapper.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_MMESH_CGAL_TRIANGULATIONWRAPPER_HH
#define DUNE_MMESH_CGAL_TRIANGULATIONWRAPPER_HH
 
#include <dune/grid/common/exceptions.hh>
#include <dune/mmesh/grid/pointfieldvector.hh>
 
namespace Dune
{
 
  template< class Coordinate >
  class Plane
  {
  public:
    Plane( const Coordinate& a, const Coordinate& b )
    {
      n_ = a;
      n_ -= b;
      std::swap( n_[0], n_[1] );
      n_[0] *= -1.0;
      n_ /= n_.two_norm();
 
      p_ = a;
    }
 
    template< class Vertex >
    Plane( const Vertex& a, const Vertex& b )
     : Plane( makeFieldVector( a->point() ), makeFieldVector( b->point() ) )
    {}
 
    auto signedDistance( const Coordinate& x ) const
    {
      return n_ * (x - p_);
    }
 
    template< class Vertex >
    auto signedDistance( const Vertex& vh ) const
    {
      return signedDistance( makeFieldVector( vh->point() ) );
    }
 
  private:
    Coordinate p_, n_;
  };
 
  template<>
  class TriangulationWrapper<2> :
    public MMeshDefaults::Triangulation<2>::type
  {
      using ThisType = TriangulationWrapper<2>;
      using BaseType = typename MMeshDefaults::Triangulation<2>::type;
 
      using CellHandle = typename BaseType::Face_handle;
      using FacetHandle = std::pair<typename BaseType::Face_handle, int>;
      using VertexHandle = typename BaseType::Vertex_handle;
 
      using Vertex_pair = std::pair< VertexHandle, VertexHandle >;
      using Vertex_pair_Facet_map = std::map< Vertex_pair, FacetHandle >;
      using Vertex_handle_unique_hash_map = std::map< VertexHandle, VertexHandle >;
 
  public:
    template< class VertexHandle, class Elements >
    void removeAndGiveNewElements (const VertexHandle& vh, Elements& elements)
    {
      std::list<FacetHandle> hole;
      this->make_hole(vh, hole);
      this->fill_hole_delaunay(hole, std::back_inserter(elements));
      this->delete_vertex(vh);
    }
 
    template< class VertexHandle, class Elements >
    void remeshHoleConstrained(const VertexHandle& vh, std::list<FacetHandle>& hole, Elements& elements, const std::vector<VertexHandle> constraint)
    {
      assert( constraint.size() == 2 );
 
      std::list<FacetHandle> hole1;
      std::list<FacetHandle> hole2;
 
      Plane< FieldVector<double, 2> > plane ( constraint[0], constraint[1] );
 
      VertexHandle right;
      VertexHandle left;
 
      auto helperface = create_face();
      FacetHandle helperfacet (helperface, 0);
 
      // collect facets corresponding to their side of the constraint
      for( const FacetHandle& facet : hole )
      {
        const auto& vh = facet.first->vertex( (facet.second+1)%3 );
 
        if ( vh == constraint[0] || vh == constraint[1] )
        {
          const auto& vh2 = facet.first->vertex( (facet.second+2)%3 );
          if ( plane.signedDistance( vh2 ) > 0 )
          {
            hole1.push_back( facet );
            left = vh;
            hole1.push_back( helperfacet );
          }
          else
          {
            hole2.push_back( facet );
            right = vh;
            hole2.push_back( helperfacet );
          }
        }
        else
        {
          if ( plane.signedDistance( vh ) > 0 )
            hole1.push_back( facet );
          else
            hole2.push_back( facet );
        }
      }
 
      // fill hole1
      helperface->set_vertices( VertexHandle(), right, left );
      this->fill_hole(vh, hole1, std::back_inserter(elements));
 
      // fill hole2
      helperface->set_vertices( VertexHandle(), left, right );
      this->fill_hole(vh, hole2, std::back_inserter(elements));
 
      // glue the facets at helperface together
      Face_handle f1; int i1;
      for ( auto f : elements )
        for ( int i = 0; i < 3; ++i )
          if (f->neighbor(i) == helperface)
          {
            if (f1 == Face_handle())
            {
              f1 = f;
              i1 = i;
            }
            else
            {
              f->set_neighbor(i, f1);
              f1->set_neighbor(i1, f);
 
              delete_vertex(vh);
              delete_face(helperface);
              return;
            }
          }
 
      DUNE_THROW( GridError, "Could not remesh hole with interface." );
    }
 
  private:
    Vertex_pair make_vertex_pair ( FacetHandle f ) const
    {
      Vertex_pair vp;
      vp.first = f.first->vertex( (f.second+1)%3 );
      vp.second = f.first->vertex( (f.second+2)%3 );
 
      if ( vp.first > vp.second )
        std::swap( vp.first, vp.second );
 
      return vp;
    }
  };
 
 
  template<>
  class TriangulationWrapper<3> :
    public MMeshDefaults::Triangulation<3>::type
  {
    using ThisType = TriangulationWrapper<3>;
    using BaseType = typename MMeshDefaults::Triangulation<3>::type;
 
  public:
    template< class VertexHandle, class Elements >
    void removeAndGiveNewElements (const VertexHandle& v, Elements fit)
    {
      DUNE_THROW( GridError, "Removal is not supported in 3D!" );
    }
  }; // end of TriangulationWrapper<3>
 
 
  template<int dim>
  class DelaunayTriangulationWrapper :
    public MMeshDefaults::Delaunay<dim>::type
  {
      using ThisType = DelaunayTriangulationWrapper<dim>;
      using BaseType = typename MMeshDefaults::Delaunay<dim>::type;
  };
 
} // namespace Dune
 
#endif