equidistantpoints.hh

Go to the documentation of this file.

#ifndef DUNE_EQUIDISTANTPOINTS_HH
#define DUNE_EQUIDISTANTPOINTS_HH

#include <vector>

#include <dune/common/fvector.hh>
#include <dune/common/forloop.hh>
#include <dune/geometry/topologyfactory.hh>
#include <dune/geometry/genericgeometry/topologytypes.hh>
#include <dune/geometry/genericgeometry/subtopologies.hh>

#include <dune/localfunctions/lagrange/emptypoints.hh>
#include <dune/localfunctions/utility/field.hh>

namespace Dune
{

  // Internal Forward Declarations
  // -----------------------------

  template< class F, unsigned int dim >
  class EquidistantPointSet;

  // EquidistantPointSetImpl
  // ----------------------------

  template< class Topology, class F >
  class EquidistantPointSetImpl;

  template< class F >
  class EquidistantPointSetImpl< GenericGeometry::Point, F >
  {
    typedef EquidistantPointSetImpl< GenericGeometry::Point, F > This;

    typedef GenericGeometry::Point Topology;

    friend class EquidistantPointSet< F, Topology::dimension >;
    friend class EquidistantPointSetImpl< GenericGeometry::Prism< Topology >, F >;
    friend class EquidistantPointSetImpl< GenericGeometry::Pyramid< Topology >, F >;

  public:
    typedef F Field;

    static const unsigned int dimension = Topology::dimension;

    static unsigned int size ( const unsigned int order )
    {
      return 1;
    }

  private:
    template< unsigned int codim, unsigned int dim >
    static unsigned int setup ( const unsigned int order,
                                unsigned int *count,
                                LagrangePoint< Field, dim > *points )
    {
      assert( codim == 0 );
      points->localKey_ = LocalKey( 0, 0, count[ 0 ]++ );
      points->point_ = Field( 0 );
      return 1;
    }
  };

  template< class BaseTopology, class F >
  class EquidistantPointSetImpl< GenericGeometry::Prism< BaseTopology >, F >
  {
    typedef EquidistantPointSetImpl< GenericGeometry::Prism< BaseTopology >, F > This;

    typedef GenericGeometry::Prism< BaseTopology > Topology;

    friend class EquidistantPointSet< F, Topology::dimension >;
    friend class EquidistantPointSetImpl< GenericGeometry::Prism< Topology >, F >;
    friend class EquidistantPointSetImpl< GenericGeometry::Pyramid< Topology >, F >;

    typedef EquidistantPointSetImpl< BaseTopology, F > BaseImpl;

  public:
    typedef F Field;

    static const unsigned int dimension = Topology::dimension;

    static unsigned int size ( const unsigned int order )
    {
      return BaseImpl::size( order ) * (order+1);
    }

  // private:
    template< unsigned int codim, unsigned int dim >
    static unsigned int setup ( const unsigned int order,
                                unsigned int *count,
                                LagrangePoint< Field, dim > *points )
    {
      unsigned int size = 0;
      unsigned int numBaseN = 0;

      if( codim < dimension )
      {
        const unsigned int vcodim = (codim < dimension ? codim : dimension-1);
        numBaseN = GenericGeometry::Size< BaseTopology, vcodim >::value;
        for( unsigned int i = 1; i < order; ++i )
        {
          const unsigned int n = BaseImpl::template setup< vcodim, dim >( order, count, points );
          for( unsigned int j = 0; j < n; ++j )
          {
            LocalKey &key = points->localKey_;
            key = LocalKey( key.subEntity(), codim, key.index() );
            points->point_[ dimension-1 ] = Field( i ) / Field( order );
            ++points;
          }
          size += n;
        }
      }

      if( codim > 0 )
      {
        const unsigned int vcodim = (codim > 0 ? codim : 1);
        const unsigned int numBaseM = GenericGeometry::Size< BaseTopology, vcodim-1 >::value;
        const unsigned int n = BaseImpl::template setup< vcodim-1, dim >( order, count+numBaseN, points );
        for( unsigned int j = 0; j < n; ++j )
        {
          LocalKey &key = points[ j ].localKey_;
          key = LocalKey( key.subEntity() + numBaseN, codim, key.index() );
          points[ j + n ].point_ = points[ j ].point_;
          points[ j + n ].point_[ dimension-1 ] = Field( 1 );
          points[ j + n ].localKey_ = LocalKey( key.subEntity() + numBaseM, codim, key.index() );
          ++count[ key.subEntity() + numBaseM ];
        }
        size += 2*n;
      }
      return size;
    }
  };

  template< class BaseTopology, class F >
  class EquidistantPointSetImpl< GenericGeometry::Pyramid< BaseTopology >, F >
  {
  align="bottom">Dune:;EquidistantPointSetImpl< GenericGeometry::Point, F >
new Fieldspan>{