Commit c2d08489 authored by Henrik Stolzmann's avatar Henrik Stolzmann Committed by Andreas Dedner
Browse files

Cleanup dune-localfunctions

Replaced old `Topology` template parameter based on template recursion
with a `GeometryType::Id` parameter.
parent 2b1c69a0
......@@ -172,23 +172,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 ( Dune::Impl::isTensor(geometry) )
return computeTensor(alpha,p,q);
else
return computeCone(alpha,p,q);
}
}
template< int dim, class scalar_t >
static int computeCone ( 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 computeTensor ( 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; }
};
......@@ -145,13 +145,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;
order_ = order;
topologyId_ = Topology::id;
topologyId_ = geometry.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 +170,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 ).id(), order );
faceStructure_.emplace_back( faceBasis, refElement.integrationOuterNormal( face ) );
}
assert( faceStructure_.size() == faceSize_ );
......@@ -186,10 +187,10 @@ 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_;
......@@ -250,12 +251,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 )
......@@ -424,19 +425,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);
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;
for (unsigned int o = 0; o <= order; ++o)
{
const typename LagrangeCoefficientsFactory::Object *pointsPtr = LagrangeCoefficientsFactory::template create< Topology >( o );
const typename LagrangeCoefficientsFactory::Object *pointsPtr = LagrangeCoefficientsFactory::template create< geometry >( o );
if ( pointsPtr == 0)
continue;
std::cout << "# Testing " << Topology::name() << " in dimension " << Topology::dimension << " with order " << o << std::endl;
std::cout << "Testing " << geometry << " with order " << o << std::endl;
typename BasisFactory::Object &basis = *BasisFactory::template create<Topology>(o);
typename BasisFactory::Object &basis = *BasisFactory::template create<geometry>(o);
ret |= test(basis,*pointsPtr,verbose);
......@@ -97,10 +100,10 @@ bool test(unsigned int order, bool verbose = false)
// derivatives in a human readabible form (aka LaTeX source)
#ifdef TEST_OUTPUT_FUNCTIONS
std::stringstream name;
name << "lagrange_" << Topology::name() << "_p" << o << ".basis";
name << "lagrange_" << geometry << "_p" << o << ".basis";
std::ofstream out(name.str().c_str());
Dune::basisPrint<0,BasisFactory,typename BasisFactory::StorageField,Topology>(out,basis);
Dune::basisPrint<1,BasisFactory,typename BasisFactory::StorageField,Topology>(out,basis);
Dune::basisPrint<0,BasisFactory,typename BasisFactory::StorageField,geometry>(out,basis);
Dune::basisPrint<1,BasisFactory,typename BasisFactory::StorageField,geometry>(out,basis);
#endif // TEST_OUTPUT_FUNCTIONS
LagrangeCoefficientsFactory::release( pointsPtr );
......@@ -152,31 +155,29 @@ int main ( int argc, char **argv )
bool tests = true;
#ifdef CHECKDIM1
tests &= test<Prism<Point> > (order);
tests &= test<Pyramid<Point> > (order);
tests &= test<GeometryTypes::cube(1)> (order);
tests &= test<GeometryTypes::simplex(1)> (order);
#endif
#ifdef CHECKDIM2
tests &= test<Prism<Prism<Point> > > (order);
tests &= test<Pyramid<Pyramid<Point> > >(order);
tests &= test<GeometryTypes::cube(2)> (order);
tests &= test<GeometryTypes::simplex(2)> (order);
#endif
#ifdef CHECKDIM3
tests &= test<Prism<Prism<Prism<Point> > > >(order);
tests &= test<Prism<Pyramid<Pyramid<Point> > > >(order);
tests &= test<Pyramid<Pyramid<Pyramid<Point> > > >(order);
tests &= test<GeometryTypes::cube(3)> (order);
tests &= test<GeometryTypes::prism> (order);
tests &= test<GeometryTypes::pyramid> (order);
tests &= test<GeometryTypes::simplex(3)> (order);
#endif
// tests &= test<Pyramid<Prism<Prism<Point> > > >(order);
std::cout << "NOT CHECKING PYRAMID!" << std::endl;
// reduce tested order to 4 in 4d unless explicitly asked for more
if (argc < 2)
order = 4;
#ifdef CHECKDIM4
tests &= test<Prism<Prism<Prism<Prism<Point> > > > >(order);
tests &= test<Pyramid<Pyramid<Pyramid<Pyramid<Point> > > > >(order);
tests &= test<GeometryTypes::cube(4)> (order);
tests &= test<GeometryTypes::simplex(4)> (order);
#endif
return (tests ? 0 : 1);
......
......@@ -15,9 +15,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
......@@ -36,16 +36,16 @@ typedef double StorageField;
typedef double ComputeField;
#endif
template <class Topology>
template< Dune::GeometryType::Id geometryId >
bool test(unsigned int order)
{
bool ret = true;
Dune::GeometryType gt(Topology::id, Topology::dimension);
constexpr Dune::GeometryType geometry = geometryId;
for (unsigned int o = 0; o <= order; ++o)
{
std::cout << "Testing " << Topology::name() << " in dimension " << Topology::dimension << " with order " << o << std::endl;
typedef Dune::OrthonormalBasisFactory<Topology::dimension,StorageField,ComputeField> BasisFactory;
const typename BasisFactory::Object &basis = *BasisFactory::template create<Topology>(o);
std::cout << "Testing " << geometry << " with order " << o << std::endl;
typedef Dune::OrthonormalBasisFactory<geometry.dim(),StorageField,ComputeField> BasisFactory;
const typename BasisFactory::Object &basis = *BasisFactory::template create<geometry>(o);
const unsigned int size = basis.size( );
......@@ -55,8 +55,8 @@ bool test(unsigned int order)
for( unsigned int i = 0; i < size * size; ++i )
m[ i ] = 0;
const Dune::QuadratureRule<double,Topology::dimension> &quadrature =
Dune::QuadratureRules<double,Topology::dimension>::rule(gt,2*order+1);
const Dune::QuadratureRule<double,geometry.dim()> &quadrature =
Dune::QuadratureRules<double,geometry.dim()>::rule(geometry,2*order+1);
const unsigned int quadratureSize = quadrature.size();
for( unsigned int qi = 0; qi < quadratureSize; ++qi )
{
......@@ -84,10 +84,10 @@ bool test(unsigned int order)
// derivatives in a human readabible form (aka LaTeX source)
#ifdef TEST_OUTPUT_FUNCTIONS
std::stringstream name;
name << "orthonormal_" << Topology::name() << "_p" << o << ".basis";
name << "orthonormal_" << geometry << "_p" << o << ".basis";
std::ofstream out(name.str().c_str());
Dune::basisPrint<0,BasisFactory,typename BasisFactory::StorageField,Topology>(out,basis);
Dune::basisPrint<1,BasisFactory,typename BasisFactory::StorageField,Topology>(out,basis);
Dune::basisPrint<0,BasisFactory,typename BasisFactory::StorageField,geometry>(out,basis);
Dune::basisPrint<1,BasisFactory,typename BasisFactory::StorageField,geometry>(out,basis);
#endif // TEST_OUTPUT_FUNCTIONS
BasisFactory::release(&basis);
......@@ -136,25 +136,25 @@ int main ( int argc, char **argv )
bool tests = true;
#ifdef CHECKDIM1
tests &= test<Prism<Point> > (order);
tests &= test<Pyramid<Point> > (order);
tests &= test<GeometryTypes::cube(1)> (order);
tests &= test<GeometryTypes::simplex(1)> (order);
#endif
#ifdef CHECKDIM2
tests &= test<Prism<Prism<Point> > > (order);
tests &= test<Pyramid<Pyramid<Point> > >(order);
tests &= test<GeometryTypes::cube(2)> (order);
tests &= test<GeometryTypes::simplex(2)> (order);
#endif
#ifdef CHECKDIM3
tests &= test<Prism<Prism<Prism<Point> > > >(order);
tests &= test<Prism<Pyramid<Pyramid<Point> > > >(order);
tests &= test<Pyramid<Prism<Prism<Point> > > >(order);
tests &= test<Pyramid<Pyramid<Pyramid<Point> > > >(order);
tests &= test<GeometryTypes::cube(3)> (order);
tests &= test<GeometryTypes::prism> (order);
tests &= test<GeometryTypes::pyramid> (order);
tests &= test<GeometryTypes::simplex(3)> (order);
#endif
#ifdef CHECKDIM4
tests &= test<Prism<Prism<Prism<Prism<Point> > > > >(order);
tests &= test<Pyramid<Pyramid<Pyramid<Pyramid<Point> > > > >(order);
tests &= test<GeometryTypes::cube(4)> (order);
tests &= test<GeometryTypes::simplex(4)> (order);
#endif