genericreferenceelements.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_GEOMETRY_GENERICREFERENCEELEMENTS_HH
#define DUNE_GEOMETRY_GENERICREFERENCEELEMENTS_HH
 
#include <dune/common/deprecated.hh>
#include <dune/common/forloop.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/visibility.hh>
 
#include <dune/geometry/genericgeometry/codimtable.hh>
#include <dune/geometry/genericgeometry/subtopologies.hh>
#include <dune/geometry/genericgeometry/referencedomain.hh>
#include <dune/geometry/genericgeometry/hybridmapping.hh>
#include <dune/geometry/genericgeometry/mappingprovider.hh>
 
namespace Dune
{
 
  // Internal Forward Declarations
  // -----------------------------
 
  template< class ctype, int dim >
  class GenericReferenceElementContainer;
 
 
 
  // GenericReferenceElement
  // -----------------------
 
  template< class ctype, int dim >
  class GenericReferenceElement
  {
    typedef GenericReferenceElement< ctype, dim > This;
 
    friend class GenericReferenceElementContainer< ctype, dim >;
 
    // make copy constructor private
    GenericReferenceElement(const GenericReferenceElement &);
 
    GenericReferenceElement () {}
 
    ~GenericReferenceElement ()
    {
      ForLoop< Destroy, 0, dim >::apply( mappings_ );
      integral_constant< int, 0 > codim0Variable;
      if(mappings_[  codim0Variable ].size())
        delete mappings_[ codim0Variable ][ 0 ];
    }
 
    class SubEntityInfo;
    template< class Topology > class CornerStorage;
    template< class Topology > struct Initialize;
    template< int codim > struct Destroy;
 
    struct GeometryTraits
      : public GenericGeometry::DefaultGeometryTraits< ctype, dim, dim >
    {
      typedef GenericGeometry::DefaultGeometryTraits< ctype, dim, dim > Base;
 
      typedef typename Base::CoordTraits CoordTraits;
 
      template< class Topology >
      struct Mapping
      {
        typedef GenericGeometry::CornerMapping< CoordTraits, Topology, dim, CornerStorage< Topology >, true > type;
      };
 
      struct Caching
      {
        static const GenericGeometry::EvaluationType evaluateJacobianTransposed = GenericGeometry::PreCompute;
        static const GenericGeometry::EvaluationType evaluateJacobianInverseTransposed = GenericGeometry::PreCompute;
        static const GenericGeometry::EvaluationType evaluateIntegrationElement = GenericGeometry::PreCompute;
        static const GenericGeometry::EvaluationType evaluateNormal = GenericGeometry::PreCompute;
      };
 
    };
 
  public:
    template< int codim >
    struct Codim
    {
      typedef GenericGeometry::HybridMapping< dim-codim, GeometryTraits > Mapping;
    };
 
  private:
    template< int codim >
    struct MappingArray
      : public std::vector< typename Codim< codim >::Mapping * >
    {};
 
    typedef GenericGeometry::CodimTable< MappingArray, dim > MappingsTable;
 
    std::vector< SubEntityInfo > info_[ dim+1 ];
 
    ctype volume_;
    std::vector< FieldVector< ctype, dim > > volumeNormals_;
 
    MappingsTable mappings_;
 
  public:
    int size ( int c ) const
    {
      assert( (c >= 0) && (c <= dim) );
      return info_[ c ].size();
    }
 
    int size ( int i, int c, int cc ) const
    {
      assert( (c >= 0) && (c <= dim) );
      return info_[ c ][ i ].size( cc );
    }
 
    int subEntity ( int i, int c, int ii, int cc ) const
    {
      assert( (c >= 0) && (c <= dim) );
      return info_[ c ][ i ].number( ii, cc );
    }
 
    const FieldVector< ctype, dim > &position( int i, int c ) const
    {
      assert( (c >= 0) && (c <= dim) );
      return info_[ c ][ i ].position();
    }
 
    bool checkInside ( const FieldVector< ctype, dim > &local ) const
    {
      return checkInside< 0 >( local, 0 );
    }
 
    template< int codim >
    bool checkInside ( const FieldVector< ctype, dim-codim > &local, int i ) const
    {
      return mapping< codim >( i ).checkInside( local );
    }
 
    template< int codim >
    FieldVector< ctype, dim >
    global( const FieldVector< ctype, dim-codim > &local, int i, int c ) const
    {
      if( c != codim )
        DUNE_THROW( Exception, "Local Coordinate Type does not correspond to codimension c." );
      assert( c == codim );
      return mapping< codim >( i ).global( local );
    }
 
    template< int codim >
    FieldVector< ctype, dim >
    global( const FieldVector< ctype, dim-codim > &local, int i ) const
    {
      return mapping< codim >( i ).global( local );
    }
 
    template< int codim >
    typename Codim< codim >::Mapping &mapping( int i ) const
    {
      integral_constant< int, codim > codimVariable;
      return *(mappings_[ codimVariable ][ i ]);
    }
 
    const GeometryType &type ( int i, int c ) const
    {
      assert( (c >= 0) && (c <= dim) );
      return info_[ c ][ i ].type();
    }
 
    const GeometryType &type () const { return type( 0, 0 ); }
 
    ctype volume () const
    {
      return volume_;
    }
 
    const FieldVector< ctype, dim > &volumeOuterNormal ( int face ) const
    {
      assert( (face >= 0) && (face < int( volumeNormals_.size())) );
      return volumeNormals_[ face ];
    }
 
    template< class Topology >
    void initializeTopology ()
    {
      dune_static_assert( (Topology::dimension == dim),
                          "Cannot initialize reference element for different dimension." );
      typedef Initialize< Topology > Init;
      typedef GenericGeometry::VirtualMapping< Topology, GeometryTraits > VirtualMapping;
 
      // set up subentities
      integral_constant< int, 0 > codim0Variable;
      mappings_[ codim0Variable ].resize( 1 );
      mappings_[ codim0Variable ][ 0 ]  = new VirtualMapping( codim0Variable );
 
      Dune::ForLoop< Init::template Codim, 0, dim >::apply( info_, mappings_ );
 
      // compute reference element volume
      typedef GenericGeometry::ReferenceDomain< Topology > ReferenceDomain;
      volume_ = ReferenceDomain::template volume< ctype >();
 
      // compute normals
      volumeNormals_.resize( ReferenceDomain::numNormals );
      for( unsigned int i = 0; i < ReferenceDomain::numNormals; ++i )
        ReferenceDomain::integrationOuterNormal( i ,volumeNormals_[ i ] );
    }
  };
 
 
  template< class ctype, int dim >
  class GenericReferenceElement< ctype, dim >::SubEntityInfo
  {
    template< class Topology, int codim > struct Initialize
    {
      template< int subcodim > struct SubCodim;
    };
 
    std::vector< int > numbering_[ dim+1 ];
    FieldVector< ctype, dim > baryCenter_;
    GeometryType type_;
 
  public:
    int size ( int cc ) const
    {
      assert( cc <= dim );
      return numbering_[ cc ].size();
    }
 
    int number ( int ii, int cc ) const
    {
      assert( cc <= dim );
      return numbering_[ cc ][ ii ];
    }
 
    const FieldVector< ctype, dim > &position () const
    {
      return baryCenter_;
    }
 
    const GeometryType &type () const
    {
      return type_;
    }
 
    template< class Topology, unsigned int codim, unsigned int i >
    DUNE_EXPORT void initialize ()
    {
      typedef Initialize< Topology, codim > Init;
      typedef GenericGeometry::ReferenceDomain< Topology > RefDomain;
 
      const unsigned int iVariable = i;
      Dune::ForLoop< Init::template SubCodim, 0, dim-codim >::apply( iVariable, numbering_ );
 
      baryCenter_ = ctype( 0 );
      static const unsigned int numCorners = size( dim );
      for( unsigned int j = 0; j < numCorners; ++j )
      {
        FieldVector< ctype, dim > corner;
        RefDomain::corner( number( j, dim ), corner );
        baryCenter_ += corner;
      }
      baryCenter_ *= ctype( 1 ) / ctype( numCorners );
 
      typedef typename GenericGeometry::SubTopology< Topology, codim, i >::type SubTopology;
      type_ = GeometryType( SubTopology::id, SubTopology::dimension );
    }
  };
 
 
  template< class ctype, int dim >
  template< class Topology >
  class GenericReferenceElement< ctype, dim >::CornerStorage
  {
    typedef GenericGeometry::ReferenceDomain< Topology > RefDomain;
 
  public:
    static const unsigned int size = Topology::numCorners;
 
    template< class SubTopology >
    struct SubStorage
    {
      typedef CornerStorage< SubTopology > type;
    };
 
    explicit CornerStorage ( const integral_constant< int, 0 > & )
    {
      for( unsigned int i = 0; i < size; ++i )
        RefDomain::corner( i, coords_[ i ] );
    }
 
    template< class Mapping, unsigned int codim >
    explicit
    CornerStorage ( const GenericGeometry::SubMappingCoords< Mapping, codim > &coords )
    {
      for( unsigned int i = 0; i < size; ++i )
        coords_[ i ] = coords[ i ];
    }
 
    const FieldVector< ctype, dim > &operator[] ( unsigned int i ) const
    {
      return coords_[ i ];
    }
 
  private:
    FieldVector< ctype, dim > coords_[ size ];
  };
 
 
  template< class ctype, int dim >
  template< class Topology, int codim >
  template< int subcodim >
  struct GenericReferenceElement< ctype, dim >::SubEntityInfo::Initialize< Topology, codim >::SubCodim
  {
    typedef GenericGeometry::SubTopologySize< Topology, codim, subcodim > SubSize;
    typedef GenericGeometry::GenericSubTopologyNumbering< Topology, codim, subcodim > SubNumbering;
 
    static void apply ( unsigned int i, std::vector< int > (&numbering)[ dim+1 ] )
    {
      const unsigned int size = SubSize::size( i );
      numbering[ codim+subcodim ].resize( size );
      for( unsigned int j = 0; j < size; ++j )
        numbering[ codim+subcodim ][ j ] = SubNumbering::number( i, j );
    }
  };
 
 
  template< class ctype, int dim >
  template< class Topology >
  struct GenericReferenceElement< ctype, dim >::Initialize
  {
    typedef Dune::GenericReferenceElement< ctype, dim > GenericReferenceElement;
 
    typedef typename GenericReferenceElement::template Codim< 0 >::Mapping ReferenceMapping;
 
    template< int codim >
    struct Codim
    {
      template< int i >
      struct SubTopology
      {
        static void apply ( std::vector< SubEntityInfo > &info )
        {
          info[ i ].template initialize< Topology, codim, i >();
        }
      };
 
      static void
      apply ( std::vector< SubEntityInfo > (&info)[ dim+1 ],
              MappingsTable &mappings )
      {
        const unsigned int size = GenericGeometry::Size< Topology, codim >::value;
        info[ codim ].resize( size );
        Dune::ForLoop< SubTopology, 0, size-1 >::apply( info[ codim ] );
 
        if( codim > 0 )
        {
          integral_constant< int, 0 > codim0Variable;
          const ReferenceMapping &refMapping = *(mappings[ codim0Variable ][ 0 ]);
 
          typedef typename GenericGeometry::MappingProvider< ReferenceMapping, codim > MappingProvider;
 
          integral_constant< int, codim > codimVariable;
          mappings[ codimVariable ].resize( size );
          for( unsigned int i = 0; i < size; ++i ) {
            char* storage = new char[MappingProvider::maxMappingSize];
            mappings[ codimVariable ][ i ] = refMapping.template trace< codim >( i, storage );
          }
        }
      }
    };
  };
 
 
 
  template< class ctype, int dim >
  template< int codim >
  struct GenericReferenceElement< ctype, dim >::Destroy
  {
    static void apply ( MappingsTable &mappings )
    {
      if (codim > 0 )
      {
        integral_constant< int, codim > codimVariable;
        for( size_t i = 0; i < mappings[ codimVariable ].size(); ++i ) {
          typedef typename Codim<codim>::Mapping Mapping;
          mappings[ codimVariable ][ i ]->~Mapping();
          char* storage = (char*)mappings[ codimVariable ][ i ];
          delete[](storage);
        }
      }
    }
  };
 
 
  // GenericReferenceElementContainer
  // --------------------------------
 
  template< class ctype, int dim >
  class GenericReferenceElementContainer
  {
    static const unsigned int numTopologies = (1u << dim);
 
  public:
    typedef GenericReferenceElement< ctype, dim > value_type;
    typedef const value_type *const_iterator;
 
    DUNE_DEPRECATED_MSG("GenericReferenceElementContainer was renamed to ReferenceElementContainer.")
    GenericReferenceElementContainer ()
    {
      ForLoop< Builder, 0, numTopologies-1 >::apply( values_ );
    }
 
    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:
    template< int topologyId >
    struct Builder
    {
      static void apply ( value_type (&values)[ numTopologies ] )
      {
        typedef typename GenericGeometry::Topology< topologyId, dim >::type Topology;
        values[ topologyId ].template initializeTopology< Topology >();
      }
    };
 
    value_type values_[ numTopologies ];
  };
 
 
  // GenericReferenceElements
  // ------------------------
 
  template< class ctype, int dim >
  struct GenericReferenceElements
  {
    typedef typename GenericReferenceElementContainer< ctype, dim >::const_iterator Iterator;
 
    DUNE_DEPRECATED_MSG("GenericReferenceElements was renamed to ReferenceElements.")
    GenericReferenceElements() {}
 
    static const GenericReferenceElement< ctype, dim > &
    DUNE_DEPRECATED_MSG("GenericReferenceElements was renamed to ReferenceElements.")
    general ( const GeometryType &type )
    {
      return container() ( type );
    }
 
    static const GenericReferenceElement< ctype, dim > &
    DUNE_DEPRECATED_MSG("GenericReferenceElements was renamed to ReferenceElements.")
    simplex ()
    {
      return container().simplex();
    }
 
    static const GenericReferenceElement< ctype, dim > &
    DUNE_DEPRECATED_MSG("GenericReferenceElements was renamed to ReferenceElements.")
    cube ()
    {
      return container().cube();
    }
 
    static Iterator begin () { return container().begin(); }
    static Iterator end () { return container().end(); }
 
  private:
    DUNE_EXPORT static const GenericReferenceElementContainer< ctype, dim > &container ()
    {
      static GenericReferenceElementContainer< ctype, dim > container;
      return container;
    }
  };
 
} // namespace Dune
 
#endif // #ifndef DUNE_GEOMETRY_GENERICREFERENCEELEMENTS_HH