Commit 6979c466 authored by Andreas Dedner's avatar Andreas Dedner

completed implementation of lfe restrict prolong - not working yet

tried an approach to fix the generic restriction.
restrictFinalize is not called after the data on the father has been
computed. Still doesn't quite work....

LocalRestrictProlong does not have an implemented default class so the new
   prolongFinalize( LocalFunctionOnFather )
needs to be implemented (empty) in each specialization. I probably missed some
parent 900ae801
......@@ -37,6 +37,7 @@ namespace Dune
// helper structs
template< int > struct RestrictLocal;
template< int > struct RestrictFinalize;
template< int > struct ProlongLocal;
template< std::size_t ... i >
......@@ -68,7 +69,11 @@ namespace Dune
{
Fem::ForLoop< RestrictLocal, 0, setSize >::apply( lfFather, lfSon, geometryInFather, initialize, localRestrictProlongTuple_ );
}
template< class LFFather >
void restrictFinalize ( LFFather &lfFather ) const
{
Fem::ForLoop< RestrictFinalize, 0, setSize >::apply( lfFather, lfSon, localRestrictProlongTuple_ );
}
template< class LFFather, class LFSon, class LocalGeometry >
void prolongLocal ( const LFFather &lfFather, LFSon &lfSon,
......@@ -132,6 +137,35 @@ namespace Dune
};
// RestrictFinalize
// ----------------
template< class ... DiscreteFunctionSpaces >
template< int i >
struct TupleLocalRestrictProlong< DiscreteFunctionSpaces ... >::
RestrictFinalize
{
template< class LFFather, class Tuple >
static void apply ( LFFather &lfFather,
const Tuple &tuple )
{
typedef SubVector< typename LFFather::LocalDofVectorType, OffsetSubMapper > SubDofVectorTypeFather;
typedef typename LFFather::BasisFunctionSetType::template SubBasisFunctionSet< i >::type SubFatherBasisFunctionSetType;
SubFatherBasisFunctionSetType subFatherBasisFunctionSet = lfFather.basisFunctionSet().template subBasisFunctionSet< i >();
std::size_t fatherBasisSetOffset = lfFather.basisFunctionSet().offset(i);
SubDofVectorTypeFather fatherSubDofVector( lfFather.localDofVector(), OffsetSubMapper( subFatherBasisFunctionSet.size(),
fatherBasisSetOffset ) );
LocalFunction< SubFatherBasisFunctionSetType, SubDofVectorTypeFather > subLFFather( subFatherBasisFunctionSet, fatherSubDofVector );
std::get< i >( tuple ).restrictFinalize( subLFFather );
}
};
// RestrictLocal
// -------------
......
......@@ -509,6 +509,7 @@ namespace Dune
// reset initialize flag
initialize = false;
}
restop.restrictFinalize(entity);
}
}
}
......
......@@ -133,6 +133,7 @@ namespace Dune
restrictLocal( father, *it, initialize );
initialize = false;
}
rpOp_.restrictFinalize( father );
}
}
......
......@@ -697,6 +697,10 @@ namespace Dune
}
}
template <class EntityType>
inline void restrictFinalize( const EntityType &father ) const
{}
//! prolong data to children and resize memory if doResize is true
template <class EntityType>
inline void prolongLocal ( const EntityType & father, const EntityType & son , bool initialize ) const
......
......@@ -68,6 +68,9 @@ namespace Dune
lfFather[ i ] += weight * lfSon[ i ];
}
}
template< class LFFather >
void restrictFinalize ( LFFather &lfFather ) const
{}
//! prolong data to children
template< class LFFather, class LFSon, class LocalGeometry >
......@@ -122,6 +125,9 @@ namespace Dune
void restrictLocal ( LFFather &lfFather, const LFSon &lfSon,
const LocalGeometry &geometryInFather, bool initialize ) const
{}
template< class LFFather >
void restrictFinalize ( LFFather &lfFather ) const
{}
//! prolong data to children
template< class LFFather, class LFSon, class LocalGeometry >
......
......@@ -166,6 +166,7 @@ namespace Dune
DUNE_THROW( GridError, "Cannot restrict over more than one level." );
initialize = false;
}
localRestrictProlong_.restrictFinalize(father);
}
template< class Function >
......
......@@ -87,6 +87,13 @@ namespace Dune
CHECK_AND_CALL_INTERFACE_IMPLEMENTATION( asImp().restrictLocal( father, son, geometryInFather, initialize ) );
}
//! finalize restriction on father
template <class Entity>
void restrictFinalize(const Entity &father) const
{
CHECK_AND_CALL_INTERFACE_IMPLEMENTATION( asImp().restrictFinalize( father ) );
}
//! prolong data to children
template< class Entity >
void prolongLocal ( const Entity &father, const Entity &son, bool initialize ) const
......@@ -257,6 +264,18 @@ namespace Dune
localRP_.restrictLocal( lfFather, constLf_, geometryInFather, initialize );
}
//! finalize restriction on father
template <class Entity>
void restrictFinalize(const Entity &father) const
{
typedef typename GridPartType::template Codim< Entity::codimension >::EntityType GridPartEntityType;
const GridPartType &gridPart = discreteFunction_.gridPart();
const GridPartEntityType &gpFather = gridPart.convert( father );
LocalFunctionType lfFather = discreteFunction_.localFunction( father );
localRP_.restrictFinalize( lfFather );
}
//! prolong data to children
template< class Entity >
void prolongLocal ( const Entity &father, const Entity &son, bool initialize ) const
......
......@@ -110,6 +110,11 @@ namespace Dune
Dune::Fem::ForLoop< RestrictLocal, 0, sizeof...( Tail ) >::apply( father, child, geometryInFather, initialize, tuple_ );
}
template< class Entity >
void restrictFinalize ( const Entity &father ) const
{
}
/** \copydoc Dune::Fem::RestrictProlongInterface::prolongLocal */
template< class Entity >
void prolongLocal ( const Entity &father, const Entity &child, bool initialize ) const
......
......@@ -111,6 +111,9 @@ namespace Dune
lfFather += temp_;
}
}
template< class LFFather >
void restrictFinalize ( LFFather &lfFather ) const
{}
template< class LFFather, class LFSon, class LocalGeometry >
void prolongLocal ( const LFFather &lfFather, LFSon &lfSon,
......
......@@ -78,6 +78,9 @@ namespace Dune
}
}
}
template< class LFFather >
void restrictFinalize ( LFFather &lfFather ) const
{}
template< class LFFather, class LFSon, class LocalGeometry >
void prolongLocal ( const LFFather &lfFather, LFSon &lfSon,
......
......@@ -203,12 +203,13 @@ namespace Dune
dofAlignment_( basisFunctionSet_ )
{}
template< class LocalFunction, class Dof >
void operator() ( const LocalFunction &localFunction, std::vector< Dof > &localDofVector ) const
template< class LocalFunction, class Vector >
void operator() ( const LocalFunction &localFunction, Vector &localDofVector ) const
{
typedef std::vector< Dof > LocalDofVector;
typedef Vector LocalDofVector;
// clear dofs before something is adedd
localDofVector.clear();
// localDofVector.clear(); // does not exist on DynVector so use 'fill' instead
std::fill(localDofVector.begin(),localDofVector.end(),0);
for( std::size_t i = 0; i < dimRange; ++i )
{
SubDofVectorWrapper< LocalDofVector, DofAlignmentType > subLdv( localDofVector, i, dofAlignment_ );
......
......@@ -11,6 +11,45 @@ namespace Dune
namespace Fem
{
namespace Impl
{
// BasisFunctionWrapper
// --------------------
template< class BasisFunctionSet, class LocalGeometry >
struct BasisFunctionWrapper
{
typedef std::remove_reference_t< BasisFunctionSet > BasisFunctionSetType;
typedef std::remove_reference_t< LocalGeometry > LocalGeometryType;
struct Traits
{
typedef typename BasisFunctionSetType::DomainType DomainType;
typedef typename BasisFunctionSetType::RangeType RangeType;
};
typedef typename BasisFunctionSetType::EntityType EntityType;
typedef typename BasisFunctionSetType::FunctionSpaceType FunctionSpaceType;
typedef typename BasisFunctionSetType::DomainType DomainType;
typedef typename BasisFunctionSetType::RangeType RangeType;
typedef typename BasisFunctionSetType::JacobianRangeType JacobianRangeType;
typedef typename BasisFunctionSetType::HessianRangeType HessianRangeType;
BasisFunctionWrapper ( const LocalGeometryType &localGeo, const BasisFunctionSetType &bset, int i )
: localGeo_( localGeo ), bset_( bset ), i_( i ), values_( bset.size() ) {}
template< class Arg >
void evaluate ( const Arg &x, typename Traits::RangeType &y ) const
{
bset_.evaluateAll( localGeo_.global(x), values_);
y = values_[i_];
}
private:
const LocalGeometryType &localGeo_;
const BasisFunctionSetType &bset_;
const int i_;
mutable std::vector<RangeType> values_;
};
} // namespce Impl
template< class LFEMap, class FunctionSpace, template< class > class Storage >
struct DefaultLocalRestrictProlong< LocalFiniteElementSpace< LFEMap, FunctionSpace, Storage > >
{
......@@ -20,7 +59,8 @@ namespace Dune
typedef typename DiscreteFunctionSpaceType::DomainFieldType DomainFieldType;
DefaultLocalRestrictProlong (const DiscreteFunctionSpaceType &space)
: space_( space ), weight_(-1)
: space_( space ), weight_(-1),
lfFather_(0)
{}
/** \brief explicit set volume ratio of son and father
......@@ -45,13 +85,35 @@ namespace Dune
const int numDofs = lfFather.numDofs();
assert( lfFather.numDofs() == lfSon.numDofs() );
if (initialize)
lfFather_.resize(numDofs,0);
std::cout << "initialize=" << initialize;
for (int r=0;r<numDofs;++r)
std::cout << " " << lfFather_[r] << " , " << lfSon[r];
auto P = calcProlongMatrix(lfFather.entity(), lfSon.entity(), numDofs);
P.usmv(weight, lfSon.localDofVector(), lfFather_);
if( initialize )
for( int i = 0; i < numDofs; ++i )
lfFather[ i ] = 0;
// y += alpha A^T x
P.usmtv(weight,lfSon.localDofVector(),lfFather.localDofVector());
std::cout << " -> ";
for (int r=0;r<numDofs;++r)
std::cout << " " << lfFather_[r];
std::cout << std::endl;
}
template <class LFFather>
void restrictFinalize( LFFather &lfFather ) const
{
const int numDofs = lfFather.numDofs();
assert( numDofs == lfFather_.size() );
for (int r=0;r<numDofs;++r)
lfFather[r] = lfFather_[r];
std::cout << "finalize";
for (int r=0;r<numDofs;++r)
std::cout << " " << lfFather[r];
std::cout << std::endl << std::endl;
lfFather_.clear();
}
//! prolong data to children
......@@ -62,8 +124,17 @@ namespace Dune
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());
// Note: we want lfSon = P^T lfFather
// but lfSon and lfFather might share dofs (e.g. Lagrange space)
// and then the necessary lfSon.clear() before computing the matrix
// vector multiplication would lead to errors so we first need to
// copy the lfFather dofs:
DynamicVector<double> lfCopy(numDofs);
for (int r=0;r<numDofs;++r)
lfCopy[r] = lfFather[r];
// y = P^T x
P.mtv(lfCopy,lfSon.localDofVector());
}
//! do discrete functions need a communication after restriction / prolongation?
......@@ -75,16 +146,35 @@ namespace Dune
{
return son.geometry().volume() / father.geometry().volume();
}
// note that this method returns
// P = ( l(phi_r)_c )_rc so it is the transposed of the matrix needed
template< class Entity >
DynamicMatrix<double> calcProlongMatrix( const Entity &father, const Entity &son, int numDofs ) const
{
// TODO
DynamicMatrix<double> P(numDofs,numDofs);
const auto &localBasis = space_.basisFunctionSet(father);
const auto &localInterpolation = space_.interpolation(son);
const auto &localGeo = son.geometryInFather();
for (int r=0;r<numDofs;++r)
{
Impl::BasisFunctionWrapper< decltype(localBasis), decltype(localGeo) >
basisWrapper( localGeo, localBasis, r );
localInterpolation( basisWrapper, P[r] );
}
/*
for (int r=0;r<numDofs;++r)
{
for (int c=0;c<numDofs;++c)
std::cout << P[c][r] << " ";
std::cout << std::endl;
}
std::cout << std::endl << std::endl;
*/
return P;
}
const DiscreteFunctionSpaceType &space_;
DomainFieldType weight_;
mutable std::vector<double> lfFather_;
};
} // namespace Fem
......
......@@ -71,6 +71,9 @@ namespace Dune
}
}
}
template< class LFFather >
void restrictFinalize ( LFFather &lfFather ) const
{}
template< class LFFather, class LFSon, class LocalGeometry >
......
......@@ -16,7 +16,7 @@
#endif
// polynomial order of base functions
const int polOrder = 2; // POLORDER;
const int polOrder = 1; // POLORDER;
#include <iostream>
#include <sstream>
......@@ -400,7 +400,7 @@ try
if( Parameter::verbose() )
std :: cout << std :: endl << "Coarsening:" << std::endl;
for( int i = ml - 1; i >= 0; --i )
// for( int i = ml - 1; i >= 0; --i )
algorithm( gridPart, solution, -step, 1, locallyAdaptive );
return 0;
......
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