Commit 9b7cf76c authored by Andreas Dedner's avatar Andreas Dedner
Browse files

Merge branch 'feature/Topology-to-GeometryType' into 'master'

Topology-to-GeometryTyp

See merge request !186
parents 2b1c69a0 d8e8a4f9
Pipeline #36312 passed with stage
in 13 minutes and 24 seconds
......@@ -27,6 +27,10 @@
* Imported the Python bindings from the 2.7 branch of dune-python.
* Replaced the combination of function arguments `topologyId` and `dim` with a single `GeometryType` argument.
Tagged the old versions of: `numLagrangePoints`, `equidistantLagrangePoints`, `RTL2InterpolationBuilder::topologyId()`,
`VirtualMonomialBasis(topologyId)`, `VirtualMonomialBasis::topologyId()` as deprecated.
# Release 2.7
* The header `lagrange.hh` now includes all headers of all Lagrange implementations,
......
......@@ -18,52 +18,57 @@ namespace Dune
// numLagrangePoints
// -----------------
inline std::size_t numLagrangePoints ( unsigned int topologyId, int dim, std::size_t order )
inline std::size_t numLagrangePoints ( const GeometryType& gt, std::size_t order )
{
assert( topologyId < Impl::numTopologies( dim ) );
const int dim = gt.dim();
if( dim > 0 )
{
const unsigned int baseId = Impl::baseTopologyId( topologyId, dim );
if( Impl::isPyramid( topologyId, dim ) )
const GeometryType baseGeometryType = Impl::getBase( gt );
if( gt.isConical() )
{
std::size_t size = 0;
for( unsigned int o = 0; o <= order; ++o )
size += numLagrangePoints( baseId, dim-1, o );
size += numLagrangePoints( baseGeometryType, o );
return size;
}
else
return numLagrangePoints( baseId, dim-1, order ) * (order+1);
return numLagrangePoints( baseGeometryType, order ) * (order+1);
}
else
return 1;
}
[[deprecated("Use numLagrangePoints(const GeometryType& gt, std::size_t order ) instead.")]]
inline std::size_t numLagrangePoints ( unsigned int topologyId, unsigned int dim, std::size_t order )
{
return numLagrangePoints ( GeometryType(topologyId, dim), order);
}
// equidistantLagrangePoints
// -------------------------
template< class ct, unsigned int cdim >
inline static unsigned int equidistantLagrangePoints ( unsigned int topologyId, unsigned int dim, unsigned int codim, std::size_t order, unsigned int *count, LagrangePoint< ct, cdim > *points )
inline static unsigned int equidistantLagrangePoints ( const GeometryType& gt, unsigned int codim, std::size_t order, unsigned int *count, LagrangePoint< ct, cdim > *points )
{
const unsigned int dim = gt.dim();
assert( (0 <= codim) && (codim <= dim) && (dim <= cdim) );
assert( topologyId < Impl::numTopologies( dim ) );
if( dim > 0 )
{
const unsigned int baseId = Impl::baseTopologyId( topologyId, dim );
const unsigned int numBaseN = (codim < dim ? Geo::Impl::size( baseId, dim-1, codim ) : 0);
const unsigned int numBaseM = (codim > 0 ? Geo::Impl::size( baseId, dim-1, codim-1 ) : 0);
const GeometryType baseGeometryType = Impl::getBase( gt );
const unsigned int numBaseN = (codim < dim ? Geo::Impl::size( baseGeometryType.id(), baseGeometryType.dim(), codim ) : 0);
const unsigned int numBaseM = (codim > 0 ? Geo::Impl::size( baseGeometryType.id(), baseGeometryType.dim(), codim-1 ) : 0);
if( Impl::isPrism( topologyId, dim ) )
if( gt.isPrismatic() )
{
unsigned int size = 0;
if( codim < dim )
{
for( unsigned int i = 1; i < order; ++i )
{
const unsigned int n = equidistantLagrangePoints( baseId, dim-1, codim, order, count, points );
const unsigned int n = equidistantLagrangePoints( baseGeometryType, codim, order, count, points );
for( unsigned int j = 0; j < n; ++j )
{
LocalKey &key = points->localKey_;
......@@ -77,7 +82,7 @@ namespace Dune
if( codim > 0 )
{
const unsigned int n = equidistantLagrangePoints( baseId, dim-1, codim-1, order, count+numBaseN, points );
const unsigned int n = equidistantLagrangePoints( baseGeometryType, codim-1, order, count+numBaseN, points );
for( unsigned int j = 0; j < n; ++j )
{
LocalKey &key = points[ j ].localKey_;
......@@ -95,7 +100,7 @@ namespace Dune
}
else
{
unsigned int size = (codim > 0 ? equidistantLagrangePoints( baseId, dim-1, codim-1, order, count, points ) : 0);
unsigned int size = (codim > 0 ? equidistantLagrangePoints( baseGeometryType, codim-1, order, count, points ) : 0);
LagrangePoint< ct, cdim > *const end = points + size;
for( ; points != end; ++points )
points->localKey_ = LocalKey( points->localKey_.subEntity(), codim, points->localKey_.index() );
......@@ -104,7 +109,7 @@ namespace Dune
{
for( unsigned int i = order-1; i > 0; --i )
{
const unsigned int n = equidistantLagrangePoints( baseId, dim-1, codim, i, count+numBaseM, points );
const unsigned int n = equidistantLagrangePoints( baseGeometryType, codim, i, count+numBaseM, points );
LagrangePoint< ct, cdim > *const end = points + n;
for( ; points != end; ++points )
{
......@@ -135,6 +140,13 @@ namespace Dune
}
}
template< class ct, unsigned int cdim >
[[deprecated("Use equidistantLagrangePoints ( GeometryType gt, ... ) instead.")]]
inline static unsigned int equidistantLagrangePoints ( unsigned int topologyId, unsigned int dim, unsigned int codim, std::size_t order, unsigned int *count, LagrangePoint< ct, cdim > *points )
{
return equidistantLagrangePoints ( GeometryType(topologyId, dim), codim, order, *count, *points );
}
// EquidistantPointSet
......@@ -156,7 +168,7 @@ namespace Dune
void build ( GeometryType gt )
{
assert( gt.dim() == dimension );
points_.resize( numLagrangePoints( gt.id(), dimension, order() ) );
points_.resize( numLagrangePoints( gt, order() ) );
typename Base::LagrangePoint *p = points_.data();
std::vector< unsigned int > count;
......@@ -164,7 +176,7 @@ namespace Dune
{
count.resize( Geo::Impl::size( gt.id(), dimension, dimension-mydim ) );
std::fill( count.begin(), count.end(), 0u );
p += equidistantLagrangePoints( gt.id(), dimension, dimension-mydim, order(), count.data(), p );
p += equidistantLagrangePoints( gt, dimension-mydim, order(), count.data(), p );
}
const auto &refElement = referenceElement<F,dimension>(gt);
F weight = refElement.volume()/F(double(points_.size()));
......@@ -172,23 +184,22 @@ namespace Dune
p.weight_ = weight;
}
template< class T >
template< GeometryType::Id geometryId >
bool build ()
{
build( GeometryType( T() ) );
build( GeometryType( geometryId ) );
return true;
}
bool buildCube ()
{
using namespace Impl;
return build< typename CubeTopology< dim >::type > ();
return build< GeometryTypes::cube(dim) > ();
}
static bool supports ( GeometryType gt, std::size_t order ) { return true; }
template< class T >
template< GeometryType::Id geometryId>
static bool supports ( std::size_t order ) {
return supports( GeometryType( T() ), order );
return supports( GeometryType( geometryId ), order );
}
private:
......
......@@ -103,17 +103,17 @@ namespace Dune
typedef typename LagrangePointSetFactory::Key Key;
typedef const LocalLagrangeInterpolation< LP,dim,F > Object;
template< class Topology >
template< GeometryType::Id geometryId >
static Object *create ( const Key &key )
{
const LagrangePointSet *lagrangeCoeff
= LagrangePointSetFactory::template create< Topology >( key );
= LagrangePointSetFactory::template create< geometryId >( key );
if ( lagrangeCoeff == 0 )
return 0;
else
return new Object( *lagrangeCoeff );
}
template< class Topology >
template< GeometryType::Id geometryId >
static bool supports ( const Key &key )
{
return true;
......
......@@ -23,14 +23,14 @@ namespace Dune
const typedef LP<F,dim> Object;
typedef std::size_t Key;
template< class T >
template< GeometryType::Id geometryId >
static Object *create ( const Key &order )
{
if (order == 0 || !Object::template supports<T>(order))
if (order == 0 || !Object::template supports<geometryId>(order))
return 0;
typedef typename std::remove_const<Object>::type LagrangeCoefficients;
LagrangeCoefficients *object = new LagrangeCoefficients(order);
if ( !object->template build<T>() )
if ( !object->template build<geometryId>() )
{
delete object;
object = nullptr;
......
......@@ -36,17 +36,17 @@ namespace Dune
typedef unsigned int Key;
typedef const Basis Object;
typedef typename Impl::SimplexTopology< dim >::type SimplexTopology;
static constexpr GeometryType SimplexGeometry = GeometryTypes::simplex(dim);
template< class Topology >
template< GeometryType::Id geometryId >
static Object *create ( const unsigned int order )
{
const MonomialBasisType &monomialBasis = *MonomialBasisProviderType::template create< SimplexTopology >( order );
const MonomialBasisType &monomialBasis = *MonomialBasisProviderType::template create< SimplexGeometry >( order );
static CoefficientMatrix _coeffs;
if( _coeffs.size() <= monomialBasis.size() )
{
ONBCompute::ONBMatrix< Topology, ComputeField > matrix( order );
ONBCompute::ONBMatrix< geometryId, ComputeField > matrix( order );
_coeffs.fill( matrix );
}
......
......@@ -35,64 +35,66 @@ namespace ONBCompute
// Integral
// --------
template< class Topology >
struct Integral;
template< class Base >
struct Integral< Dune::Impl::Pyramid< Base > >
template< Dune::GeometryType::Id geometryId >
struct Integral
{
static constexpr Dune::GeometryType geometry = geometryId;
static constexpr int dimension = geometry.dim();
template< int dim, class scalar_t >
static int compute ( const Dune::MultiIndex< dim, scalar_t > &alpha,
scalar_t &p, scalar_t &q )
{
const int dimension = Base::dimension+1;
int i = alpha.z( Base::dimension );
int ord = Integral< Base >::compute( alpha, p, q );
if constexpr (geometry.isVertex())
{
p = scalar_t( 1 );
q = scalar_t( 1 );
return 0;
}
else
{
if constexpr ( geometry.isPrismatic() )
return computePrismatic(alpha,p,q);
else
return computeConical(alpha,p,q);
}
}
template< int dim, class scalar_t >
static int computeConical ( const Dune::MultiIndex< dim, scalar_t > &alpha,
scalar_t &p, scalar_t &q )
{
int i = alpha.z( dimension-1 );
constexpr Dune::GeometryType::Id baseGeometryId = Dune::Impl::getBase(geometry);
int ord = Integral< baseGeometryId >::compute( alpha, p, q );
p *= factorial< scalar_t >( 1, i );
q *= factorial< scalar_t >( dimension + ord, dimension + ord + i );
return ord + i;
}
};
template< class Base >
struct Integral< Dune::Impl::Prism< Base > >
{
template< int dim, class scalar_t >
static int compute ( const Dune::MultiIndex< dim, scalar_t > &alpha,
scalar_t &p, scalar_t &q )
static int computePrismatic ( const Dune::MultiIndex< dim, scalar_t > &alpha,
scalar_t &p, scalar_t &q )
{
int i = alpha.z( Base::dimension );
int ord = Integral< Base >::compute( alpha, p, q );
//Integral< Base >::compute( alpha, p, q );
int i = alpha.z( dimension-1 );
constexpr Dune::GeometryType::Id baseGeometryId = Dune::Impl::getBase(geometry);
int ord = Integral< baseGeometryId >::compute( alpha, p, q );
//p *= scalar_t( 1 );
q *= scalar_t( i+1 );
return ord + i;
}
};
template<>
struct Integral< Dune::Impl::Point >
{
template< int dim, class scalar_t >
static int compute ( const Dune::MultiIndex< dim, scalar_t > &alpha,
scalar_t &p, scalar_t &q )
{
p = scalar_t( 1 );
q = scalar_t( 1 );
return 0;
}
};
// ONBMatrix
// ---------
template< class Topology, class scalar_t >
template< Dune::GeometryType::Id geometryId, class scalar_t >
class ONBMatrix
: public Dune::LFEMatrix< scalar_t >
{
typedef ONBMatrix< Topology, scalar_t > This;
typedef ONBMatrix< geometryId, scalar_t > This;
typedef Dune::LFEMatrix< scalar_t > Base;
public:
......@@ -102,7 +104,8 @@ namespace ONBCompute
explicit ONBMatrix ( unsigned int order )
{
// get all multiindecies for monomial basis
const unsigned int dim = Topology::dimension;
constexpr Dune::GeometryType geometry = geometryId;
constexpr unsigned int dim = geometry.dim();
typedef Dune::MultiIndex< dim, scalar_t > MI;
Dune::StandardMonomialBasis< dim, MI > basis( order );
const std::size_t size = basis.size();
......@@ -123,7 +126,7 @@ namespace ONBCompute
{
for( std::size_t j = 0; j < size; ++j )
{
Integral< Topology >::compute( y[ i ][ 0 ] * y[ j ][ 0 ], p, q );
Integral< geometryId >::compute( y[ i ][ 0 ] * y[ j ][ 0 ], p, q );
S( i, j ) = p;
S( i, j ) /= q;
}
......
......@@ -68,22 +68,22 @@ namespace Dune
typedef std::size_t Key;
typedef const LocalCoefficientsContainer Object;
template< class Topology >
template< GeometryType::Id geometryId >
static Object *create( const Key &key )
{
typedef RaviartThomasL2InterpolationFactory< dim, double > InterpolationFactory;
if( !supports< Topology >( key ) )
if( !supports< geometryId >( key ) )
return nullptr;
typename InterpolationFactory::Object *interpolation = InterpolationFactory::template create< Topology >( key );
typename InterpolationFactory::Object *interpolation = InterpolationFactory::template create< geometryId >( key );
Object *localKeys = new Object( *interpolation );
InterpolationFactory::release( interpolation );
return localKeys;
}
template< class Topology >
template< GeometryType::Id geometryId >
static bool supports ( const Key &key )
{
return Impl::IsSimplex< Topology >::value;
return GeometryType(geometryId).isSimplex();
}
static void release( Object *object ) { delete object; }
};
......@@ -127,9 +127,10 @@ namespace Dune
TestFaceBasisFactory::release( f.basis_ );
}
unsigned int topologyId () const { return topologyId_; }
[[deprecated("Use type().id() instead.")]]
unsigned int topologyId () const { return type().id(); }
GeometryType type () const { return GeometryType( topologyId(), dimension ); }
GeometryType type () const { return geometry_; }
std::size_t order () const { return order_; }
......@@ -145,13 +146,14 @@ namespace Dune
// normal of face f
const Normal &normal ( unsigned int f ) const { assert( f < faceSize() ); return *(faceStructure_[ f ].normal_); }
template< class Topology >
template< GeometryType::Id geometryId >
void build ( std::size_t order )
{
constexpr GeometryType geometry = geometryId;
geometry_ = geometry;
order_ = order;
topologyId_ = Topology::id;
testBasis_ = (order > 0 ? TestBasisFactory::template create< Topology >( order-1 ) : nullptr);
testBasis_ = (order > 0 ? TestBasisFactory::template create< geometry >( order-1 ) : nullptr);
const auto &refElement = ReferenceElements< Field, dimension >::general( type() );
faceSize_ = refElement.size( 1 );
......@@ -169,7 +171,7 @@ namespace Dune
* And depending on the dynamic face index a different face geometry is needed.
*
*/
TestFaceBasis *faceBasis = Impl::IfTopology< CreateFaceBasis, dimension-1 >::apply( refElement.type( face, 1 ).id(), order );
TestFaceBasis *faceBasis = Impl::IfGeometryType< CreateFaceBasis, dimension-1 >::apply( refElement.type( face, 1 ), order );
faceStructure_.emplace_back( faceBasis, refElement.integrationOuterNormal( face ) );
}
assert( faceStructure_.size() == faceSize_ );
......@@ -186,15 +188,16 @@ namespace Dune
const Dune::FieldVector< Field, dimension > *normal_;
};
template< class FaceTopology >
template< GeometryType::Id faceGeometryId >
struct CreateFaceBasis
{
static TestFaceBasis *apply ( std::size_t order ) { return TestFaceBasisFactory::template create< FaceTopology >( order ); }
static TestFaceBasis *apply ( std::size_t order ) { return TestFaceBasisFactory::template create< faceGeometryId >( order ); }
};
std::vector< FaceStructure > faceStructure_;
TestBasis *testBasis_ = nullptr;
unsigned int topologyId_, faceSize_;
GeometryType geometry_;
unsigned int faceSize_;
std::size_t order_;
};
......@@ -250,12 +253,12 @@ namespace Dune
{
return size_;
}
template <class Topology>
template <GeometryType::Id geometryId>
void build( std::size_t order )
{
size_ = 0;
order_ = order;
builder_.template build<Topology>(order_);
builder_.template build<geometryId>(order_);
if (builder_.testBasis())
size_ += dimension*builder_.testBasis()->size();
for ( unsigned int f=0; f<builder_.faceSize(); ++f )
......@@ -283,7 +286,7 @@ namespace Dune
template< class Func, class Container, bool type >
void interpolate ( typename Base::template Helper<Func,Container,type> &func ) const
{
const Dune::GeometryType geoType( builder_.topologyId(), dimension );
const Dune::GeometryType geoType = builder_.type();
std::vector< Field > testBasisVal;
......@@ -424,19 +427,20 @@ namespace Dune
typedef const RaviartThomasL2Interpolation<dim,Field> Object;
typedef std::size_t Key;
typedef typename std::remove_const<Object>::type NonConstObject;
template <class Topology>
template <GeometryType::Id geometryId>
static Object *create( const Key &key )
{
if ( !supports<Topology>(key) )
if ( !supports<geometryId>(key) )
return 0;
NonConstObject *interpol = new NonConstObject();
interpol->template build<Topology>(key);
interpol->template build<geometryId>(key);
return interpol;
}
template< class Topology >
template< GeometryType::Id geometryId >
static bool supports ( const Key &key )
{
return Impl::IsSimplex<Topology>::value;
return GeometryType(geometryId).isSimplex();
}
static void release( Object *object ) { delete object; }
};
......
......@@ -12,7 +12,7 @@
namespace Dune
{
template < class Topology, class Field >
template < GeometryType::Id geometryId, class Field >
struct RTVecMatrix;
template <unsigned int dim, class Field>
......@@ -31,11 +31,11 @@ namespace Dune
{
typedef MonomialBasisProvider<dd,FF> Type;
};
template< class Topology >
template< GeometryType::Id geometryId >
static Object *create ( const Key &order )
{
RTVecMatrix<Topology,Field> vecMatrix(order);
MBasis *mbasis = MBasisFactory::template create<Topology>(order+1);
RTVecMatrix<geometryId,Field> vecMatrix(order);
MBasis *mbasis = MBasisFactory::template create<geometryId>(order+1);
typename std::remove_const<Object>::type *tmBasis = new typename std::remove_const<Object>::type(*mbasis);
tmBasis->fill(vecMatrix);
return tmBasis;
......@@ -43,12 +43,13 @@ namespace Dune
static void release( Object *object ) { delete object; }
};
template <class Topology, class Field>
template <GeometryType::Id geometryId, class Field>
struct RTVecMatrix
{
static const unsigned int dim = Topology::dimension;
static constexpr GeometryType geometry = geometryId;
static const unsigned int dim = geometry.dim();
typedef MultiIndex<dim,Field> MI;
typedef MonomialBasis<Topology,MI> MIBasis;
typedef MonomialBasis<geometryId,MI> MIBasis;
RTVecMatrix(std::size_t order)
{
/*
......
......@@ -2,6 +2,8 @@
// vi: set et ts=4 sw=2 sts=2:
#include <config.h>
#include <dune/geometry/type.hh>
#include <dune/localfunctions/utility/field.hh>
#include <dune/localfunctions/utility/basisprint.hh>
......@@ -14,9 +16,9 @@
* shape functions on simplices.
*
* The topology can be chosen at compile time by setting TOPOLOGY
* to a string like
* to a Dune::GeometryType like
* \code
* Pyramid<Pyramid<Point> > >
* GeometryTypes::simplex(2)
* \endcode
* which generates a 2d simplex. If TOPOLOGY is not set, all
* topologies up to 4d are tested. Note, this may lead to prolonged
......@@ -72,24 +74,25 @@ bool test(const Basis &basis, const Points &points, bool verbose)
return ret;
}
template <class Topology>
template <Dune::GeometryType::Id geometryId>
bool test(unsigned int order, bool verbose = false)
{
typedef Dune::LagrangeBasisFactory<Dune::EquidistantPointSet,Topology::dimension,StorageField,ComputeField> BasisFactory;
typedef Dune::LagrangeCoefficientsFactory< Dune::EquidistantPointSet, Topology::dimension,double > LagrangeCoefficientsFactory;
constexpr Dune::GeometryType geometry = geometryId;
typedef Dune::LagrangeBasisFactory<Dune::EquidistantPointSet, geometry.dim(), StorageField, ComputeField> BasisFactory;
typedef Dune::LagrangeCoefficientsFactory< Dune::EquidistantPointSet, geometry.dim(), double > LagrangeCoefficientsFactory;
bool ret = true;