Commit 1a4564f9 authored by Andreas Dedner's avatar Andreas Dedner

change dune/fem/space/test/test-spaceinterpolation.cc to demonstrate the

issue with the interpolate method on the LagrangeSpace (based on
dune-localfunction) when used with dimRange>1

get the LFESpaces to work correctly with dimR>1
- still need to make sure RT,BDM, etc are doing he right thing
- prolongation and restriction is not yet working (error in dof manager)

remove some debug output

fix an issue with the correct 'block index' in the lfe space. Now RT seems
to work fine again
parent 083e59ab
......@@ -281,6 +281,27 @@ namespace Dune
}
};
///////////////////////////////////////////////////////
template< class Grid >
struct BoundaryIdGetter
{
typedef std::function<int(const typename Grid::template Codim<0>::Geometry::GlobalCoordinate&)> Caller;
BoundaryIdGetter(const Caller &caller)
: caller_(caller) {}
template< class Intersection >
int boundaryId ( const Intersection &intersection ) const
{
if (!intersection.boundary) return 0;
const auto x = intersection.geometry().center();
int val = caller(x);
return (val==0)?BoundaryIdProvider<Grid>::boundaryId(intersection):val;
}
private:
Caller caller_;
};
} // namespace Fem
} // namespace Dune
......
......@@ -4,4 +4,5 @@ dune_install(
interpolation.hh
space.hh
shapefunctionset.hh
localrestrictprolong.hh
)
......@@ -164,7 +164,7 @@ namespace Dune
typedef SingletonList< LFEMapType *, BlockMapperType, BlockMapperSingletonFactory > BlockMapperProviderType;
public:
typedef LocalFiniteElementInterpolation< BasisFunctionSetType, LocalInterpolationType > InterpolationType;
typedef LocalFiniteElementInterpolation< ThisType, LocalInterpolationType, LocalBasisType::dimRange==1 > InterpolationType;
using BaseType::order;
......
......@@ -10,6 +10,8 @@
#include <dune/fem/function/common/localcontribution.hh>
#include <dune/fem/space/basisfunctionset/transformed.hh>
#include <dune/fem/space/basisfunctionset/vectorial.hh>
#include <dune/fem/space/combinedspace/interpolation.hh>
namespace Dune
{
......@@ -19,10 +21,8 @@ namespace Dune
namespace Impl
{
// LocalFunctionWrapper
// --------------------
template< class LocalFunction, class BasisFunctionSet >
struct LocalFunctionWrapper
{
......@@ -45,11 +45,8 @@ namespace Dune
private:
const LocalFunction &lf_;
};
// LocalFunctionWrapper
// --------------------
template< class LocalFunction, class Entity, class ShapeFunctionSet, class Transformation >
struct LocalFunctionWrapper< LocalFunction, TransformedBasisFunctionSet< Entity, ShapeFunctionSet, Transformation > >
{
......@@ -85,13 +82,19 @@ namespace Dune
// LocalFiniteElementInterpolation
// -------------------------------
template< class BasisFunctionSet, class LocalInterpolation >
class LocalFiniteElementInterpolation
template< class Space, class LocalInterpolation, bool scalarSFS >
class LocalFiniteElementInterpolation;
// a vector valued shape basis function set taken from
// dune-localfunction - the vector value 'blow up' is not yet supported
// for this case
template< class Space, class LocalInterpolation >
class LocalFiniteElementInterpolation<Space,LocalInterpolation,false>
{
typedef LocalFiniteElementInterpolation< BasisFunctionSet, LocalInterpolation > ThisType;
typedef LocalFiniteElementInterpolation< Space, LocalInterpolation,false > ThisType;
public:
typedef BasisFunctionSet BasisFunctionSetType;
typedef typename Space::BasisFunctionSetType BasisFunctionSetType;
typedef LocalInterpolation LocalInterpolationType;
private:
......@@ -109,7 +112,7 @@ namespace Dune
public:
explicit LocalFiniteElementInterpolation ( const BasisFunctionSetType &basisFunctionSet,
const LocalInterpolationType &localInterpolation = LocalInterpolationType() )
: basisFunctionSet_( basisFunctionSet ),
: basisFunctionSet_( basisFunctionSet ),
localInterpolation_( localInterpolation )
{}
......@@ -123,16 +126,110 @@ namespace Dune
template< class LocalFunction, class DiscreteFunction, template< class > class Assembly >
void operator() ( const LocalFunction &localFunction, LocalContribution< DiscreteFunction, Assembly > &localContribution ) const
{
LocalFunctionWrapper< LocalFunction > wrapper( localFunction, basisFunctionSet() );
localInterpolation().interpolate( wrapper, localContribution.localDofVector() );
(*this)(localFunction,localContribution.localDofVector());
}
BasisFunctionSetType basisFunctionSet () const { return basisFunctionSet_; }
const LocalInterpolationType &localInterpolation () const { return localInterpolation_; }
private:
BasisFunctionSetType basisFunctionSet_;
LocalInterpolationType localInterpolation_;
};
// a scalar valued shape basis function set taken from
// dune-localfunction - for this the vector value 'blow up' is supported
// as with other spaces
template< class Space, class LocalInterpolation >
class LocalFiniteElementInterpolation<Space,LocalInterpolation,true>
{
typedef LocalFiniteElementInterpolation< Space, LocalInterpolation,true > ThisType;
public:
typedef typename Space::BasisFunctionSetType BasisFunctionSetType;
typedef LocalInterpolation LocalInterpolationType;
private:
typedef typename BasisFunctionSetType::FunctionSpaceType FunctionSpaceType;
typedef typename FunctionSpaceType::RangeType RangeType;
typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
static const int dimRange = FunctionSpaceType::dimRange;
typedef std::size_t size_type;
typedef VerticalDofAlignment< BasisFunctionSetType, RangeType> DofAlignmentType;
struct RangeConverter
{
RangeConverter ( std::size_t range ) : range_( range ) {}
template< class T >
FieldVector< T, 1 > operator() ( const FieldVector< T, dimRange > &in ) const
{
return in[ range_ ];
}
template< class T, int j >
FieldMatrix< T, 1, j > operator() ( const FieldMatrix< T, dimRange, j > &in ) const
{
return in[ range_ ];
}
protected:
std::size_t range_;
};
template <class DofVector, class DofAlignment>
struct SubDofVectorWrapper
: public SubDofVector< DofVector, DofAlignment >
{
typedef SubDofVector< DofVector, DofAlignment > BaseType;
SubDofVectorWrapper( DofVector& dofs, int coordinate, const DofAlignment &dofAlignment )
: BaseType( dofs, coordinate, dofAlignment )
{}
//! do nothing on clear/resize since it's done in apply of this class
void clear() {}
void resize( const unsigned int) {}
};
public:
explicit LocalFiniteElementInterpolation ( const BasisFunctionSetType &basisFunctionSet,
const LocalInterpolationType &localInterpolation = LocalInterpolationType() )
: basisFunctionSet_( basisFunctionSet ),
localInterpolation_( localInterpolation ),
dofAlignment_( basisFunctionSet_ )
{}
template< class LocalFunction, class Dof >
void operator() ( const LocalFunction &localFunction, std::vector< Dof > &localDofVector ) const
{
typedef std::vector< Dof > LocalDofVector;
// clear dofs before something is adedd
localDofVector.clear();
for( std::size_t i = 0; i < dimRange; ++i )
{
SubDofVectorWrapper< LocalDofVector, DofAlignmentType > subLdv( localDofVector, i, dofAlignment_ );
localInterpolation().interpolate(
localFunctionConverter( localFunction, RangeConverter( i ) ), subLdv );
}
}
template< class LocalFunction, class DiscreteFunction, template< class > class Assembly >
void operator() ( const LocalFunction &localFunction, LocalContribution< DiscreteFunction, Assembly > &localContribution ) const
{
(*this)(localFunction,localContribution.localDofVector());
}
BasisFunctionSetType basisFunctionSet () const { return basisFunctionSet_; }
const LocalInterpolationType &localInterpolation () const { return localInterpolation_; }
private:
BasisFunctionSet basisFunctionSet_;
BasisFunctionSetType basisFunctionSet_;
LocalInterpolationType localInterpolation_;
DofAlignmentType dofAlignment_;
};
} // namespace Fem
......
#ifndef DUNE_FEM_LFESPACE_LOCALRESTRICTPROLONG_HH
#define DUNE_FEM_LFESPACE_LOCALRESTRICTPROLONG_HH
#include <dune/common/dynmatrix.hh>
#include <dune/fem/function/localfunction/localfunction.hh>
#include <dune/fem/space/common/localrestrictprolong.hh>
namespace Dune
{
namespace Fem
{
template< class LFEMap, class FunctionSpace, template< class > class Storage >
struct DefaultLocalRestrictProlong< LocalFiniteElementSpace< LFEMap, FunctionSpace, Storage > >
{
typedef LocalFiniteElementSpace< LFEMap, FunctionSpace, Storage > DiscreteFunctionSpaceType;
typedef DefaultLocalRestrictProlong< DiscreteFunctionSpaceType > ThisType;
typedef typename DiscreteFunctionSpaceType::DomainFieldType DomainFieldType;
DefaultLocalRestrictProlong (const DiscreteFunctionSpaceType &space)
: space_( space ), weight_(-1)
{}
/** \brief explicit set volume ratio of son and father
*
* \param[in] weight volume of son / volume of father
*
* \note If this ratio is set, it is assume to be constant.
*/
void setFatherChildWeight ( const DomainFieldType &weight )
{
weight_ = weight;
}
//! restrict data to father
template< class LFFather, class LFSon, class LocalGeometry >
void restrictLocal ( LFFather &lfFather, const LFSon &lfSon,
const LocalGeometry &geometryInFather, bool initialize ) const
{
const DomainFieldType weight = (weight_ < DomainFieldType( 0 ) ? calcWeight( lfFather.entity(), lfSon.entity() ) : weight_);
assert( weight > 0.0 );
const int numDofs = lfFather.numDofs();
assert( lfFather.numDofs() == lfSon.numDofs() );
auto P = calcProlongMatrix(lfFather.entity(), lfSon.entity(), numDofs);
if( initialize )
for( int i = 0; i < numDofs; ++i )
lfFather[ i ] = 0;
// y += alpha A^T x
P.usmtv(weight,lfSon.localDofVector(),lfFather.localDofVector());
}
//! prolong data to children
template< class LFFather, class LFSon, class LocalGeometry >
void prolongLocal ( const LFFather &lfFather, LFSon &lfSon,
const LocalGeometry &geometryInFather, bool initialize ) const
{
const int numDofs = lfFather.numDofs();
assert( lfFather.numDofs() == lfSon.numDofs() );
auto P = calcProlongMatrix(lfFather.entity(), lfSon.entity(), numDofs);
for( int i = 0; i < numDofs; ++i )
P.mv(lfFather.localDofVector(),lfSon.localDofVector());
}
//! do discrete functions need a communication after restriction / prolongation?
bool needCommunication () const { return true; }
protected:
template< class Entity >
static DomainFieldType calcWeight ( const Entity &father, const Entity &son )
{
return son.geometry().volume() / father.geometry().volume();
}
template< class Entity >
DynamicMatrix<double> calcProlongMatrix( const Entity &father, const Entity &son, int numDofs ) const
{
// TODO
DynamicMatrix<double> P(numDofs,numDofs);
return P;
}
const DiscreteFunctionSpaceType &space_;
DomainFieldType weight_;
};
} // namespace Fem
} // namespace Dune
#endif // #ifndef DUNE_FEM_LOCALRESTRICTPROLONG_HH
......@@ -25,6 +25,7 @@
#include <dune/fem/space/localfiniteelement/shapefunctionset.hh>
#include <dune/fem/space/localfiniteelement/capabilities.hh>
#include <dune/fem/space/localfiniteelement/interpolation.hh>
#include <dune/fem/space/localfiniteelement/localrestrictprolong.hh>
namespace Dune
{
......@@ -48,8 +49,12 @@ namespace Dune
typedef GridFunctionSpace< GridPartType, FunctionSpace > FunctionSpaceType;
static constexpr int codimension = 0;
static constexpr bool isScalar = LocalFiniteElementType::Traits::LocalBasisType::Traits::dimRange==1;
typedef Hybrid::IndexRange< int, 1 > LocalBlockIndices;
typedef std::conditional_t<isScalar,
Hybrid::IndexRange< int, FunctionSpace::dimRange >,
Hybrid::IndexRange< int, 1 >
> LocalBlockIndices;
private:
typedef typename GridPartType::template Codim< codimension >::EntityType EntityType;
......@@ -60,8 +65,12 @@ namespace Dune
typedef LocalFunctionsShapeFunctionSet< typename LocalFiniteElementType::Traits::LocalBasisType > LocalFunctionsShapeFunctionSetType;
typedef SelectCachingShapeFunctionSet< LocalFunctionsShapeFunctionSetType, Storage > StoredShapeFunctionSetType;
typedef ShapeFunctionSetProxy< StoredShapeFunctionSetType > ShapeFunctionSetType;
// typedef VectorialShapeFunctionSet< ShapeFunctionSetProxy< StoredShapeFunctionSetType >, typename FunctionSpaceType::RangeType > ShapeFunctionSetType;
typedef ShapeFunctionSetProxy< StoredShapeFunctionSetType > ShapeFunctionSetProxyType;
// only extend to vector valued in case that the original space is scalar
typedef std::conditional_t<isScalar,
VectorialShapeFunctionSet< ShapeFunctionSetProxyType, typename FunctionSpaceType::RangeType >,
ShapeFunctionSetProxyType
> ShapeFunctionSetType;
private:
template< class LFEM >
......@@ -102,9 +111,8 @@ namespace Dune
typedef DiscreteFunctionSpaceDefault< LocalFiniteElementSpaceTraits< LFEMap, FunctionSpace, Storage > >
BaseType;
typedef typename BaseType::Traits Traits;
public:
typedef typename BaseType::Traits Traits;
typedef typename BaseType::FunctionSpaceType FunctionSpaceType;
typedef typename BaseType::GridPartType GridPartType;
......@@ -166,7 +174,7 @@ namespace Dune
typedef SingletonList< LFEMapType *, BlockMapperType, BlockMapperSingletonFactory > BlockMapperProviderType;
public:
typedef LocalFiniteElementInterpolation< BasisFunctionSetType, LocalInterpolationType > InterpolationType;
typedef LocalFiniteElementInterpolation< ThisType, LocalInterpolationType, Traits::isScalar > InterpolationType;
using BaseType::order;
......@@ -181,7 +189,7 @@ namespace Dune
{}
template< class GridPart, std::enable_if_t< std::is_same< GridPart, GridPartType >::value && !std::is_same< KeyType, std::tuple<> >::value, int > = 0 >
explicit LocalFiniteElementSpace ( GridPart &gridPart, const KeyType &key,
explicit LocalFiniteElementSpace ( GridPart &gridPart, const KeyType &key = KeyType(1),
const InterfaceType commInterface = InteriorBorder_All_Interface,
const CommunicationDirection commDirection = ForwardCommunication )
: BaseType( gridPart, commInterface, commDirection ),
......@@ -257,6 +265,13 @@ namespace Dune
std::unique_ptr< BlockMapperType, typename BlockMapperProviderType::Deleter > blockMapper_;
};
template< class LFEMap, class FunctionSpace, template< class > class Storage, int newRange >
struct ToNewDimRangeFunctionSpace<
LocalFiniteElementSpace<LFEMap, FunctionSpace, Storage>, newRange>
{
typedef LocalFiniteElementSpace<LFEMap, typename ToNewDimRangeFunctionSpace<FunctionSpace,newRange>::Type, Storage> Type;
};
} // namespace Fem
} // namespace Dune
......
#ifndef DUNE_FEM_SPACE_RAVIARTTHOMAS_LOCALINTERPOLATION_HH
#define DUNE_FEM_SPACE_RAVIARTTHOMAS_LOCALINTERPOLATION_HH
#error SHOULDN't be needed
// C++ includes
#include <cassert>
#include <vector>
......
......@@ -58,6 +58,7 @@ LINK_LIBRARIES dunefem MPI_RANKS 1 2 3 4 TIMEOUT 9999999 )
dune_add_test( SOURCES test-slavedofs.cc LINK_LIBRARIES dunefem MPI_RANKS 1 2 4 TIMEOUT 300 COMPILE_DEFINITIONS "${DEFAULTFLAGS};USE_COMBINED_SPACE" )
dune_add_test( SOURCES test-raviartthomasinterpolation.cc CMAKE_GUARD dune-localfunctions_FOUND LINK_LIBRARIES dunefem COMPILE_DEFINITIONS "${DEFAULTFLAGS}" )
dune_add_test( SOURCES test-spaceinterpolation.cc LINK_LIBRARIES dunefem COMPILE_DEFINITIONS "${DEFAULTFLAGS}" )
if( ${FEM_TORTURE_TESTSS} )
# test adaptation of combinedspaces
......
......@@ -10,10 +10,12 @@
template< class DiscreteSpace >
struct LocalBasis
{
typedef DiscreteSpace FunctionSpaceType;
typedef typename DiscreteSpace::BasisFunctionSetType BasisFunctionSetType;
typedef typename DiscreteSpace::DomainType DomainType;
typedef typename DiscreteSpace::RangeType RangeType;
typedef typename DiscreteSpace::JacobianRangeType JacobianRangeType;
typedef typename DiscreteSpace::HessianRangeType HessianRangeType;
typedef typename DiscreteSpace::RangeFieldType RangeFieldType;
typedef typename BasisFunctionSetType::EntityType EntityType;
......
......@@ -34,6 +34,8 @@
#include <dune/fem/space/rannacherturek.hh>
#include <dune/fem/space/raviartthomas.hh>
#include <dune/fem/io/file/dataoutput.hh>
#include <dune/fem/space/test/checklocalinterpolation.hh>
......@@ -90,33 +92,31 @@ static const int dimRange = GridPartType::dimensionworld;
typedef Dune::Fem::GridFunctionSpace< GridPartType, Dune::FieldVector< typename GridPartType::ctype, dimRange > > FunctionSpaceType;
typedef std::tuple<
Dune::Fem::FiniteVolumeSpace< FunctionSpaceType, GridPartType >,
Dune::Fem::DiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 0 >,
Dune::Fem::DiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 1 >,
Dune::Fem::DiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 2 >,
Dune::Fem::LagrangeDiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 1 >,
Dune::Fem::LagrangeDiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 2 >,
Dune::Fem::LegendreDiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 1 >,
Dune::Fem::LegendreDiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 2 >,
// Dune::Fem::FiniteVolumeSpace< FunctionSpaceType, GridPartType >,
// Dune::Fem::DiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 0 >,
// Dune::Fem::DiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 1 >,
// Dune::Fem::DiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 2 >,
// Dune::Fem::LagrangeDiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 1 >,
// Dune::Fem::LagrangeDiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 2 >,
// Dune::Fem::LegendreDiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 1 >,
// Dune::Fem::LegendreDiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 2 >,
#if HAVE_DUNE_LOCALFUNCTIONS
Dune::Fem::BrezziDouglasMariniSpace< FunctionSpaceType, GridPartType, 1 >,
Dune::Fem::BrezziDouglasMariniSpace< FunctionSpaceType, GridPartType, 2 >,
Dune::Fem::RaviartThomasSpace< FunctionSpaceType, GridPartType, 0 >,
// Dune::Fem::BrezziDouglasMariniSpace< FunctionSpaceType, GridPartType, 1 >,
// Dune::Fem::BrezziDouglasMariniSpace< FunctionSpaceType, GridPartType, 2 >,
// Dune::Fem::RaviartThomasSpace< FunctionSpaceType, GridPartType, 0 >,
Dune::Fem::RaviartThomasSpace< FunctionSpaceType, GridPartType, 1 >,
// Dune::Fem::LagrangeSpace< FunctionSpaceType, GridPartType >,
Dune::Fem::LagrangeSpace< FunctionSpaceType, GridPartType >,
// Dune::Fem::RannacherTurekSpace< FunctionSpaceType, GridPartType >,
#endif // #if HAVE_DUNE_LOCALFUNCTIONS
Dune::Fem::LagrangeDiscreteFunctionSpace< FunctionSpaceType, GridPartType, 1 >,
Dune::Fem::LagrangeDiscreteFunctionSpace< FunctionSpaceType, GridPartType, 2 >
Dune::Fem::LagrangeDiscreteFunctionSpace< FunctionSpaceType, GridPartType, 1 >
// Dune::Fem::LagrangeDiscreteFunctionSpace< FunctionSpaceType, GridPartType, 2 >
> DiscreteFunctionSpacesType;
typedef ErrorTuple< DiscreteFunctionSpacesType >::Type ErrorTupleType;
// algorithm
// ---------
template< class DiscreteFunctionSpace >
std::pair< Real< DiscreteFunctionSpace >, Real< DiscreteFunctionSpace > >
algorithm ( typename DiscreteFunctionSpace::GridPartType &gridPart )
......@@ -129,6 +129,17 @@ algorithm ( typename DiscreteFunctionSpace::GridPartType &gridPart )
const auto uGridExact = gridFunctionAdapter( "exact solution", uExact, gridPart, 3 );
interpolate( uGridExact, u );
#if 1
{
static int turn = 0;
typedef std::tuple< decltype(u)* > IODataType;
IODataType data( &u );
Dune::Fem::DataOutput< GridType, IODataType > output( gridPart.grid(), data );
output.writeData( turn, "test" );
++turn;
}
#endif
checkLocalInterpolation( space );
Dune::Fem::L2Norm< GridPartType > l2norm( gridPart );
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment