referenceelements.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_GEOMETRY_REFERENCEELEMENTS_HH
#define DUNE_GEOMETRY_REFERENCEELEMENTS_HH
 
#include <algorithm>
#include <limits>
 
#include <dune/common/array.hh>
#include <dune/common/forloop.hh>
#include <dune/common/nullptr.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/visibility.hh>
#include <dune/common/unused.hh>
 
#include <dune/geometry/affinegeometry.hh>
#include <dune/geometry/genericgeometry/codimtable.hh>
#include <dune/geometry/genericgeometry/subtopologies.hh>
#include <dune/geometry/genericgeometry/referencedomain.hh>
 
namespace Dune
{
 
  // Internal Forward Declarations
  // -----------------------------
 
  template< class ctype, int dim >
  class ReferenceElementContainer;
 
  template< class ctype, int dim >
  struct ReferenceElements;
 
 
 
  template< class ctype, int dim >
  class ReferenceElement
  {
    typedef ReferenceElement< ctype, dim > This;
 
    friend class ReferenceElementContainer< ctype, dim >;
 
    struct SubEntityInfo;
 
    // make copy constructor private
    ReferenceElement ( const This & );
 
    ReferenceElement () {}
 
    template< int codim > struct CreateGeometries;
 
  public:
    template< int codim >
    struct Codim
    {
      typedef AffineGeometry< ctype, dim-codim, dim > Geometry;
    };
 
    int size ( int c ) const
    {
      assert( (c >= 0) && (c <= dim) );
      return info_[ c ].size();
    }
 
    int size ( int i, int c, int cc ) const
    {
      assert( (i >= 0) && (i < size( c )) );
      return info_[ c ][ i ].size( cc );
    }
 
    int subEntity ( int i, int c, int ii, int cc ) const
    {
      assert( (i >= 0) && (i < size( c )) );
      return info_[ c ][ i ].number( ii, cc );
    }
 
    const GeometryType &type ( int i, int c ) const
    {
      assert( (i >= 0) && (i < size( c )) );
      return info_[ c ][ i ].type();
    }
 
    const GeometryType &type () const { return type( 0, 0 ); }
 
    const FieldVector< ctype, dim > &position( int i, int c ) const
    {
      assert( (c >= 0) && (c <= dim) );
      return baryCenters_[ c ][ i ];
    }
 
    bool checkInside ( const FieldVector< ctype, dim > &local ) const
    {
      const ctype tolerance = ctype( 64 ) * std::numeric_limits< ctype >::epsilon();
      return GenericGeometry::template checkInside< ctype, dim >( type().id(), dim, local, tolerance );
    }
 
    template< int codim >
    bool
    DUNE_DEPRECATED_MSG("Method will be removed after dune-geometry-2.4")
    checkInside ( const FieldVector< ctype, dim-codim > &local, int i ) const
    {
      const ctype tolerance = ctype( 64 ) * std::numeric_limits< ctype >::epsilon();
      return GenericGeometry::template checkInside< ctype, dim-codim >( type( i, codim ).id(), dim-codim, local, tolerance );
    }
 
    template< int codim >
    typename Codim< codim >::Geometry geometry ( int i ) const
    {
      integral_constant< int, codim > codimVariable;
      return geometries_[ codimVariable ][ i ];
    }
 
    ctype volume () const
    {
      return volume_;
    }
 
    const FieldVector< ctype, dim > &integrationOuterNormal ( int face ) const
    {
      assert( (face >= 0) && (face < int( integrationNormals_.size() )) );
      return integrationNormals_[ face ];
    }
 
  private:
    void initialize ( unsigned int topologyId )
    {
      assert( topologyId < GenericGeometry::numTopologies( dim ) );
 
      // set up subentities
      for( int codim = 0; codim <= dim; ++codim )
      {
        const unsigned int size = GenericGeometry::size( topologyId, dim, codim );
        info_[ codim ].resize( size );
        for( unsigned int i = 0; i < size; ++i )
          info_[ codim ][ i ].initialize( topologyId, codim, i );
      }
 
      // compute corners
      const unsigned int numVertices = size( dim );
      baryCenters_[ dim ].resize( numVertices );
      GenericGeometry::referenceCorners( topologyId, dim, &(baryCenters_[ dim ][ 0 ]) );
 
      // compute barycenters
      for( int codim = 0; codim < dim; ++codim )
      {
        baryCenters_[ codim ].resize( size(codim) );
        for( int i = 0; i < size( codim ); ++i )
        {
          baryCenters_[ codim ][ i ] = FieldVector< ctype, dim >( ctype( 0 ) );
          const unsigned int numCorners = size( i, codim, dim );
          for( unsigned int j = 0; j < numCorners; ++j )
            baryCenters_[ codim ][ i ] += baryCenters_[ dim ][ subEntity( i, codim, j, dim ) ];
          baryCenters_[ codim ][ i ] *= ctype( 1 ) / ctype( numCorners );
        }
      }
 
      // compute reference element volume
      volume_ = GenericGeometry::template referenceVolume< ctype >( topologyId, dim );
 
      // compute integration outer normals
      if( dim > 0 )
      {
        integrationNormals_.resize( size( 1 ) );
        GenericGeometry::referenceIntegrationOuterNormals( topologyId, dim, &(integrationNormals_[ 0 ]) );
      }
 
      // set up geometries
      Dune::ForLoop< CreateGeometries, 0, dim >::apply( *this, geometries_ );
    }
 
    template< int codim >
    struct GeometryArray
    : public std::vector< typename Codim< codim >::Geometry >
    {};
 
    typedef GenericGeometry::CodimTable< GeometryArray, dim > GeometryTable;
 
    ctype volume_;
 
    std::vector< FieldVector< ctype, dim > > baryCenters_[ dim+1 ];
    std::vector< FieldVector< ctype, dim > > integrationNormals_;
 
    GeometryTable geometries_;
 
    std::vector< SubEntityInfo > info_[ dim+1 ];
  };
 
  template< class ctype, int dim >
  struct ReferenceElement< ctype, dim >::SubEntityInfo
  {
    SubEntityInfo ()
      : numbering_( nullptr )
    {
      std::fill( offset_.begin(), offset_.end(), 0 );
    }
 
    SubEntityInfo ( const SubEntityInfo &other )
      : offset_( other.offset_ ),
        type_( other.type_ )
    {
      numbering_ = allocate();
      std::copy( other.numbering_, other.numbering_ + capacity(), numbering_ );
    }
 
    ~SubEntityInfo () { deallocate( numbering_ ); }
 
    const SubEntityInfo &operator= ( const SubEntityInfo &other )
    {
      type_ = other.type_;
      offset_ = other.offset_;
 
      deallocate( numbering_ );
      numbering_ = allocate();
      std::copy( other.numbering_, other.numbering_ + capacity(), numbering_ );
 
      return *this;
    }
 
    int size ( int cc ) const
    {
      assert( (cc >= codim()) && (cc <= dim) );
      return (offset_[ cc+1 ] - offset_[ cc ]);
    }
 
    int number ( int ii, int cc ) const
    {
      assert( (ii >= 0) && (ii < size( cc )) );
      return numbering_[ offset_[ cc ] + ii ];
    }
 
    const GeometryType &type () const { return type_; }
 
    void initialize ( unsigned int topologyId, int codim, unsigned int i )
    {
      const unsigned int subId = GenericGeometry::subTopologyId( topologyId, dim, codim, i );
      type_ = GeometryType( subId, dim-codim );
 
      // compute offsets
      for( int cc = 0; cc <= codim; ++cc )
        offset_[ cc ] = 0;
      for( int cc = codim; cc <= dim; ++cc )
        offset_[ cc+1 ] = offset_[ cc ] + GenericGeometry::size( subId, dim-codim, cc-codim );
 
      // compute subnumbering
      deallocate( numbering_ );
      numbering_ = allocate();
      for( int cc = codim; cc <= dim; ++cc )
        GenericGeometry::subTopologyNumbering( topologyId, dim, codim, i, cc-codim, numbering_+offset_[ cc ], numbering_+offset_[ cc+1 ] );
    }
 
  protected:
    int codim () const { return dim - type().dim(); }
 
    unsigned int *allocate () { return (capacity() != 0 ? new unsigned int[ capacity() ] : nullptr); }
    void deallocate ( unsigned int *ptr ) { delete[] ptr; }
    unsigned int capacity () const { return offset_[ dim+1 ]; }
 
  private:
    unsigned int *numbering_;
    array< unsigned int, dim+2 > offset_;
    GeometryType type_;
  };
 
 
  template< class ctype, int dim >
  template< int codim >
  struct ReferenceElement< ctype, dim >::CreateGeometries
  {
    template< int cc >
    static const ReferenceElement< ctype, dim-cc > &
    subRefElement( const ReferenceElement< ctype, dim > &refElement, int i, integral_constant< int, cc > )
    {
      return ReferenceElements< ctype, dim-cc >::general( refElement.type( i, cc ) );
    }
 
    static const ReferenceElement< ctype, dim > &
    subRefElement( const ReferenceElement< ctype, dim > &refElement, int i, integral_constant< int, 0 > )
    {
      DUNE_UNUSED_PARAMETER(i);
      return refElement;
    }
 
    static void apply ( const ReferenceElement< ctype, dim > &refElement, GeometryTable &geometries )
    {
      const int size = refElement.size( codim );
      std::vector< FieldVector< ctype, dim > > origins( size );
      std::vector< FieldMatrix< ctype, dim - codim, dim > > jacobianTransposeds( size );
      GenericGeometry::referenceEmbeddings( refElement.type().id(), dim, codim, &(origins[ 0 ]), &(jacobianTransposeds[ 0 ]) );
 
      integral_constant< int, codim > codimVariable;
      geometries[ codimVariable ].reserve( size );
      for( int i = 0; i < size; ++i )
      {
        typename Codim< codim >::Geometry geometry( subRefElement( refElement, i, codimVariable ), origins[ i ], jacobianTransposeds[ i ] );
        geometries[ codimVariable ].push_back( geometry );
      }
    }
  };
 
 
 
  // ReferenceElementContainer
  // -------------------------
 
  template< class ctype, int dim >
  class ReferenceElementContainer
  {
    static const unsigned int numTopologies = (1u << dim);
 
  public:
    typedef ReferenceElement< ctype, dim > value_type;
    typedef const value_type *const_iterator;
 
    ReferenceElementContainer ()
    {
      for( unsigned int topologyId = 0; topologyId < numTopologies; ++topologyId )
        values_[ topologyId ].initialize( topologyId );
    }
 
    const value_type &operator() ( const GeometryType &type ) const
    {
      assert( type.dim() == dim );
      return values_[ type.id() ];
    }
 
    const value_type &simplex () const
    {
      return values_[ GenericGeometry::SimplexTopology< dim >::type::id ];
    }
 
    const value_type &cube () const
    {
      return values_[ GenericGeometry::CubeTopology< dim >::type::id ];
    }
 
    const value_type &pyramid () const
    {
      return values_[ GenericGeometry::PyramidTopology< dim >::type::id ];
    }
 
    const value_type &prism () const
    {
      return values_[ GenericGeometry::PrismTopology< dim >::type::id ];
    }
 
    const_iterator begin () const { return values_; }
    const_iterator end () const { return values_ + numTopologies; }
 
  private:
    value_type values_[ numTopologies ];
  };
 
 
 
  // ReferenceElements
  // ------------------------
 
  template< class ctype, int dim >
  struct ReferenceElements
  {
    typedef typename ReferenceElementContainer< ctype, dim >::const_iterator Iterator;
 
    static const ReferenceElement< ctype, dim > &
    general ( const GeometryType &type )
    {
      return container() ( type );
    }
 
    static const ReferenceElement< ctype, dim > &simplex ()
    {
      return container().simplex();
    }
 
    static const ReferenceElement< ctype, dim > &cube ()
    {
      return container().cube();
    }
 
    static Iterator begin () { return container().begin(); }
    static Iterator end () { return container().end(); }
 
  private:
    DUNE_EXPORT static const ReferenceElementContainer< ctype, dim > &container ()
    {
      static ReferenceElementContainer< ctype, dim > container;
      return container;
    }
  };
 
} // namespace Dune
 
#endif // #ifndef DUNE_GEOMETRY_REFERENCEELEMENTS_HH