Skip to content
Snippets Groups Projects
Commit dd77f4c7 authored by Robert K's avatar Robert K
Browse files

Merged from release/2.4 to make stokes compile.

parent e9ce4269
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 564 additions and 296 deletions
dune/fem-dg/algorithm/sub/elliptic.hh
+ 127
− 136
View file @ dd77f4c7
...@@ -150,30 +150,30 @@ namespace Fem ...@@ -150,30 +150,30 @@ namespace Fem
std::shared_ptr< Matrix > matrix_; std::shared_ptr< Matrix > matrix_;
}; };
template< class ErrorEstimatorImp >
template< class DiscreteFunctionImp, template<class> class SigmaDiscreteFunctionChooserImp, class AssemblerImp, int polOrder>
class PoissonSigmaEstimator class PoissonSigmaEstimator
{ {
public: public:
typedef ErrorEstimatorImp ErrorEstimatorType;
typedef typename ErrorEstimatorType::DiscreteFunctionType DiscreteFunctionType;
typedef typename ErrorEstimatorType::SigmaDiscreteFunctionType SigmaDiscreteFunctionType;
typedef typename ErrorEstimatorType::SigmaDiscreteFunctionSpaceType SigmaDiscreteFunctionSpaceType;
typedef typename ErrorEstimatorType::DGOperatorType DGOperatorType;
static const int polynomialOrder = ErrorEstimatorType::polynomialOrder;
typedef typename DGOperatorType::ContainerType ContainerType;
typedef DiscreteFunctionImp DiscreteFunctionType; typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType; typedef typename DiscreteFunctionType::GridPartType GridPartType;
typedef typename DiscreteFunctionType::GridPartType GridPartType; typedef typename GridPartType::GridType GridType;
typedef typename GridPartType::GridType GridType;
typedef AssemblerImp AssemblerType; PoissonSigmaEstimator( ContainerType& container,
typedef AssemblerImp DGOperatorType; const DGOperatorType& assembler,
static const int polynomialOrder = polOrder;
typedef typename DiscreteFunctionSpaceType ::
template ToNewDimRange< GridType::dimension * DiscreteFunctionSpaceType::FunctionSpaceType::dimRange >::NewFunctionSpaceType
SigmaFunctionSpaceType;
PoissonSigmaEstimator( GridPartType& gridPart,
const DiscreteFunctionType& solution,
const AssemblerType& assembler,
const std::string name = "" ) const std::string name = "" )
: gridPart_( gridPart ), : container_( container ),
solution_( solution ), gridPart_( container_.gridPart() ),
solution_( *container_.solution() ),
assembler_( assembler ), assembler_( assembler ),
sigmaSpace_( gridPart_ ), sigmaSpace_( gridPart_ ),
sigmaDiscreteFunction_( "sigma-"+name, sigmaSpace_ ), sigmaDiscreteFunction_( "sigma-"+name, sigmaSpace_ ),
...@@ -183,44 +183,6 @@ namespace Fem ...@@ -183,44 +183,6 @@ namespace Fem
sigmaEstimateFunction_( "function 4 estimate-"+name, sigmaLocalEstimate_, gridPart_, solution_.space().order() ) sigmaEstimateFunction_( "function 4 estimate-"+name, sigmaLocalEstimate_, gridPart_, solution_.space().order() )
{} {}
//- hybrid spaces use PAdaptiveDGSpace
template <class Grid, int topoId>
struct SigmaSpaceChooser
{
typedef Fem::PAdaptiveDGSpace< SigmaFunctionSpaceType, GridPartType, polOrder > Type;
};
//- cubes use LegendreDGSpace
template <class Grid>
struct SigmaSpaceChooser< Grid, 1 >
{
typedef Dune::Fem::LegendreDiscontinuousGalerkinSpace< SigmaFunctionSpaceType, GridPartType, polOrder > Type;
};
//- cubes use OrthonormalDGSpace
template <class Grid>
struct SigmaSpaceChooser< Grid, 0 >
{
typedef Dune::Fem::DiscontinuousGalerkinSpace< SigmaFunctionSpaceType, GridPartType, polOrder > Type;
};
// work arround internal compiler error
enum { simplexTopoId = GenericGeometry::SimplexTopology< GridType::dimension >::type::id,
cubeTopoId = GenericGeometry::CubeTopology< GridType::dimension >::type::id,
myTopo = Dune::Capabilities::hasSingleGeometryType< GridType >::topologyId
};
enum { topoId = (simplexTopoId == myTopo) ? 0 : // 0 = simplex, 1 = cube, -1 = hybrid
(myTopo == cubeTopoId) ? 1 : -1 };
typedef typename SigmaSpaceChooser< GridType, topoId >::Type SigmaDiscreteFunctionSpaceType;
typedef typename SigmaDiscreteFunctionChooserImp<SigmaDiscreteFunctionSpaceType>::type
SigmaDiscreteFunctionType;
// compute the function sigma = grad u + sum_e r_e // compute the function sigma = grad u + sum_e r_e
template <class DF, class Operator> template <class DF, class Operator>
struct SigmaLocal : public Fem::LocalFunctionAdapterHasInitialize struct SigmaLocal : public Fem::LocalFunctionAdapterHasInitialize
...@@ -367,9 +329,9 @@ namespace Fem ...@@ -367,9 +329,9 @@ namespace Fem
SigmaLocalType sigmaLocal_; SigmaLocalType sigmaLocal_;
}; };
typedef Dune::Fem::LocalFunctionAdapter< SigmaLocal<DiscreteFunctionType, AssemblerType> > typedef Dune::Fem::LocalFunctionAdapter< SigmaLocal<DiscreteFunctionType, DGOperatorType> >
SigmaEstimateFunctionType; SigmaEstimateFunctionType;
typedef SigmaLocal<DiscreteFunctionType, AssemblerType> SigmaLocalType; typedef SigmaLocal<DiscreteFunctionType, DGOperatorType> SigmaLocalType;
typedef SigmaLocalFunction<SigmaLocalType > SigmaLocalFunctionType; typedef SigmaLocalFunction<SigmaLocalType > SigmaLocalFunctionType;
typedef Dune::Fem::LocalFunctionAdapter<SigmaLocalFunctionType> SigmaLocalFunctionAdapterType; typedef Dune::Fem::LocalFunctionAdapter<SigmaLocalFunctionType> SigmaLocalFunctionAdapterType;
...@@ -391,13 +353,14 @@ namespace Fem ...@@ -391,13 +353,14 @@ namespace Fem
public: public:
ContainerType& container_;
GridPartType& gridPart_; GridPartType& gridPart_;
const DiscreteFunctionType& solution_; const DiscreteFunctionType& solution_;
const AssemblerType& assembler_; const DGOperatorType& assembler_;
SigmaDiscreteFunctionSpaceType sigmaSpace_; SigmaDiscreteFunctionSpaceType sigmaSpace_;
SigmaDiscreteFunctionType sigmaDiscreteFunction_; SigmaDiscreteFunctionType sigmaDiscreteFunction_;
SigmaLocal<DiscreteFunctionType, AssemblerType> SigmaLocal<DiscreteFunctionType, DGOperatorType>
sigmaLocalEstimate_; sigmaLocalEstimate_;
SigmaLocalFunctionType sigmaLocalFunction_; SigmaLocalFunctionType sigmaLocalFunction_;
SigmaLocalFunctionAdapterType sigma_; SigmaLocalFunctionAdapterType sigma_;
...@@ -407,41 +370,7 @@ namespace Fem ...@@ -407,41 +370,7 @@ namespace Fem
template< int polOrder > template< class DiscreteFunctionSpaceImp, int polOrder, class SigmaEstimatorImp >
class NoEstimator
{
public:
NoEstimator(){}
int minimalOrder()
{
return polOrder;
}
bool isPadaptive()
{
return false;
}
template< class ProblemImp>
double estimate( ProblemImp problem )
{ return 0.0; }
template< class EntityImp >
int newOrder( double tolerance, EntityImp& entity )
{
return minimalOrder();
}
bool mark( double tolerance ){ return false; }
void closure(){}
};
template< int polOrder, class DiscreteFunctionSpaceImp, class EstimatorImp = NoEstimator<polOrder> >
class PAdaptivity class PAdaptivity
{ {
struct PolOrderStructure struct PolOrderStructure
...@@ -461,21 +390,27 @@ namespace Fem ...@@ -461,21 +390,27 @@ namespace Fem
public: public:
typedef DiscreteFunctionSpaceImp DiscreteFunctionSpaceType; typedef DiscreteFunctionSpaceImp DiscreteFunctionSpaceType;
typedef typename DiscreteFunctionSpaceType::GridType GridType; typedef typename DiscreteFunctionSpaceType::GridType GridType;
typedef EstimatorImp EstimatorType; typedef SigmaEstimatorImp SigmaEstimatorType;
typedef typename SigmaEstimatorType::ErrorEstimatorType ErrorEstimatorType;
typedef typename ErrorEstimatorType::DGOperatorType DGOperatorType;
typedef typename SigmaEstimatorType::ContainerType ContainerType;
typedef PersistentContainer<GridType,PolOrderStructure> PolOrderContainer; typedef PersistentContainer<GridType,PolOrderStructure> PolOrderContainer;
PAdaptivity( GridType& grid, const DiscreteFunctionSpaceType& space, EstimatorType& estimator, const std::string name = "" ) PAdaptivity( ContainerType& container, DGOperatorType& assembler, const std::string name = "" )
: polOrderContainer_( grid, 0 ), : polOrderContainer_( container.grid(), 0 ),
space_( space ), space_( container.space() ),
estimator_( estimator ), sigmaEstimator_( container, assembler, name ),
errorEstimator_( *container.solution(), assembler, sigmaEstimator_.sigma() ),
param_( AdaptationParameters() ) param_( AdaptationParameters() )
{ {
#ifdef PADAPTSPACE #ifdef PADAPTSPACE
// we start with max order // we start with max order
typedef typename PolOrderContainer::Iterator Iterator ; typedef typename PolOrderContainer::Iterator Iterator ;
const Iterator end = polOrderContainer_.end(); const Iterator end = polOrderContainer_.end();
const int minimalOrder = estimator_.minimalOrder() ; const int minimalOrder = errorEstimator_.minimalOrder() ;
for( Iterator it = polOrderContainer_.begin(); it != end; ++it ) for( Iterator it = polOrderContainer_.begin(); it != end; ++it )
{ {
(*it).value() = minimalOrder ; (*it).value() = minimalOrder ;
...@@ -485,19 +420,19 @@ namespace Fem ...@@ -485,19 +420,19 @@ namespace Fem
bool adaptive() const bool adaptive() const
{ {
return estimator_.isPadaptive(); return errorEstimator_.isPadaptive();
} }
void prepare() void prepare()
{ {
#ifdef PADAPTSPACE #ifdef PADAPTSPACE
const int minimalOrder = estimator_.minimalOrder(); const int minimalOrder = errorEstimator_.minimalOrder();
// only implemented for PAdaptiveSpace // only implemented for PAdaptiveSpace
std::vector<int> polOrderVec( space_.gridPart().indexSet().size(0) ); std::vector<int> polOrderVec( space_.gridPart().indexSet().size(0) );
std::fill( polOrderVec.begin(), polOrderVec.end(), polOrder ); std::fill( polOrderVec.begin(), polOrderVec.end(), polOrder );
polOrderContainer_.resize(); polOrderContainer_.resize();
if ( estimator_.isPadaptive() ) if ( errorEstimator_.isPadaptive() )
{ {
typedef typename DiscreteFunctionSpaceType::IteratorType IteratorType; typedef typename DiscreteFunctionSpaceType::IteratorType IteratorType;
typedef typename IteratorType::Entity EntityType ; typedef typename IteratorType::Entity EntityType ;
...@@ -528,12 +463,12 @@ namespace Fem ...@@ -528,12 +463,12 @@ namespace Fem
template< class ProblemImp > template< class ProblemImp >
bool estimateMark( const ProblemImp& problem ) bool estimateMark( const ProblemImp& problem )
{ {
double tolerance = param_.refinementTolerance();
#ifdef PADAPTSPACE #ifdef PADAPTSPACE
double tolerance = param_.refinementTolerance();
// resize container // resize container
polOrderContainer_.resize(); polOrderContainer_.resize();
const double error = estimator_.estimate( problem ); const double error = errorEstimator_.estimate( problem );
std::cout << "ESTIMATE: " << error << std::endl; std::cout << "ESTIMATE: " << error << std::endl;
typedef typename DiscreteFunctionSpaceType::IteratorType IteratorType; typedef typename DiscreteFunctionSpaceType::IteratorType IteratorType;
const IteratorType end = space_.end(); const IteratorType end = space_.end();
...@@ -541,9 +476,9 @@ namespace Fem ...@@ -541,9 +476,9 @@ namespace Fem
{ {
const typename IteratorType::Entity &entity = *it; const typename IteratorType::Entity &entity = *it;
polOrderContainer_[entity].value() = polOrderContainer_[entity].value() =
estimator_.newOrder( 0.98*tolerance, entity ); errorEstimator_.newOrder( 0.98*tolerance, entity );
} }
return (error < std::abs(tolerance) ? false : estimator_.mark( 0.98 * tolerance)); return (error < std::abs(tolerance) ? false : errorEstimator_.mark( 0.98 * tolerance));
#else #else
return false; return false;
#endif #endif
...@@ -551,52 +486,109 @@ namespace Fem ...@@ -551,52 +486,109 @@ namespace Fem
void closure() void closure()
{ {
#ifdef PADAPTSPACE #ifdef PADAPTSPACE
estimator_.closure(); errorEstimator_.closure();
#endif #endif
} }
ErrorEstimatorType& errorEstimator()
{
return errorEstimator_;
}
SigmaEstimatorType& sigmaEstimator()
{
return sigmaEstimator_;
}
private: private:
PolOrderContainer polOrderContainer_; PolOrderContainer polOrderContainer_;
const DiscreteFunctionSpaceType& space_; const DiscreteFunctionSpaceType& space_;
EstimatorType& estimator_; SigmaEstimatorType sigmaEstimator_;
ErrorEstimatorType errorEstimator_;
AdaptationParameters param_; AdaptationParameters param_;
}; };
/**
* \brief Adaptation indicator doing no indication and marking of the entities.
*/
class NoPAdaptIndicator
{
public:
typedef uint64_t UInt64Type;
template< class... Args >
NoPAdaptIndicator( Args&& ...a )
{}
bool adaptive() const
{
return false;
}
size_t numberOfElements() const
{
return 0;
}
UInt64Type globalNumberOfElements() const { return 0; }
template< class TimeProviderImp >
void setAdaptation( TimeProviderImp& tp ) {}
void preAdapt() {}
void estimateMark( const bool initialAdapt = false ) {}
void postAdapt(){}
void finalize(){}
int minNumberOfElements() const { return 0; }
int maxNumberOfElements() const { return 0; }
const int finestLevel() const { return 0; }
// // return some info
// const typename SigmaEstimatorType::SigmaDiscreteFunctionType& sigma()
// {
// return pAdapt_.sigmaEstimator().sigma();
// }
};
/** /**
* \brief Adaptation indicator doing no indication and marking of the entities. * \brief Adaptation indicator doing no indication and marking of the entities.
*/ */
template< class EstimatorImp, class SigmaEstimatorImp, class ProblemImp > template< class PAdaptivityImp, class ProblemImp >
class PAdaptIndicator class PAdaptIndicator
{ {
protected:
typedef PAdaptivityImp PAdaptivityType;
typedef ProblemImp ProblemType; typedef ProblemImp ProblemType;
typedef SigmaEstimatorImp SigmaEstimatorType; typedef typename PAdaptivityType::SigmaEstimatorType SigmaEstimatorType;
typedef EstimatorImp EstimatorType; typedef typename PAdaptivityType::ErrorEstimatorType ErrorEstimatorType;
typedef typename SigmaEstimatorType::DiscreteFunctionType DiscreteFunctionType; typedef typename SigmaEstimatorType::DiscreteFunctionType DiscreteFunctionType;
typedef typename SigmaEstimatorType::AssemblerType AssemblerType; typedef typename SigmaEstimatorType::DGOperatorType DGOperatorType;
typedef typename SigmaEstimatorType::SigmaFunctionSpaceType SigmaFunctionSpaceType; typedef typename DGOperatorType::ContainerType ContainerType;
typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType; typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
typedef typename DiscreteFunctionType::GridPartType GridPartType; typedef typename DiscreteFunctionType::GridPartType GridPartType;
typedef typename GridPartType::GridType GridType; typedef typename GridPartType::GridType GridType;
static const int polOrder = SigmaEstimatorType::polynomialOrder; static const int polOrder = SigmaEstimatorType::polynomialOrder;
typedef PAdaptivity< polOrder, DiscreteFunctionSpaceType, EstimatorImp > PAdaptivityType;
enum { dimension = GridType::dimension }; enum { dimension = GridType::dimension };
public: public:
typedef uint64_t UInt64Type; typedef uint64_t UInt64Type;
PAdaptIndicator( GridPartType& gridPart, PAdaptIndicator( ContainerType& container,
DiscreteFunctionType& solution, DGOperatorType& assembler,
const ProblemType& problem, const ProblemType& problem,
AssemblerType& assembler,
const std::string name = "" ) const std::string name = "" )
: sigmaEstimator_( gridPart, solution, assembler, name ), : pAdapt_( container, assembler, name ),
estimator_( solution, sigmaEstimator_.sigma(), assembler, gridPart.grid(), name /*, AdaptationParameters( param ) */),
pAdapt_( gridPart.grid(), solution.space(), estimator_ ),
problem_( problem ) problem_( problem )
{} {}
...@@ -625,20 +617,16 @@ namespace Fem ...@@ -625,20 +617,16 @@ namespace Fem
void estimateMark( const bool initialAdapt = false ) void estimateMark( const bool initialAdapt = false )
{ {
// calculate new sigma // calculate new sigma
sigmaEstimator_.update(); pAdapt_.sigmaEstimator().update();
const bool marked = pAdapt_.estimateMark( problem_ ); const bool marked = pAdapt_.estimateMark( problem_ );
if( marked ) if( marked )
pAdapt_.closure(); pAdapt_.closure();
} }
void postAdapt() void postAdapt(){}
{
}
void finalize() void finalize(){}
{
}
int minNumberOfElements() const { return 0; } int minNumberOfElements() const { return 0; }
...@@ -649,13 +637,11 @@ namespace Fem ...@@ -649,13 +637,11 @@ namespace Fem
// return some info // return some info
const typename SigmaEstimatorType::SigmaDiscreteFunctionType& sigma() const typename SigmaEstimatorType::SigmaDiscreteFunctionType& sigma()
{ {
return sigmaEstimator_.sigma(); return pAdapt_.sigmaEstimator().sigma();
} }
private: private:
SigmaEstimatorType sigmaEstimator_;
EstimatorType estimator_;
PAdaptivityType pAdapt_; PAdaptivityType pAdapt_;
const ProblemType& problem_; const ProblemType& problem_;
}; };
...@@ -738,7 +724,7 @@ namespace Fem ...@@ -738,7 +724,7 @@ namespace Fem
space_( container_.space() ), space_( container_.space() ),
assembler_( container_, model() ), assembler_( container_, model() ),
matrix_( container_.matrix() ), matrix_( container_.matrix() ),
adaptIndicator_( std::make_unique<AdaptIndicatorType>( gridPart_, *container_.solution().get(), problem(), assembler_, name() ) ), adaptIndicator_( Std::make_unique<AdaptIndicatorType>( container_, assembler_, problem(), name() ) ),
step_( 0 ), step_( 0 ),
time_( 0 ) time_( 0 )
{ {
...@@ -765,6 +751,11 @@ namespace Fem ...@@ -765,6 +751,11 @@ namespace Fem
return assembler_; return assembler_;
} }
AssemblerType& assembler ()
{
return assembler_;
}
//ADAPTATION //ADAPTATION
virtual AdaptIndicatorType* adaptIndicator() virtual AdaptIndicatorType* adaptIndicator()
{ {
...@@ -814,7 +805,7 @@ namespace Fem ...@@ -814,7 +805,7 @@ namespace Fem
//! finalize computation by calculating errors and EOCs //! finalize computation by calculating errors and EOCs
virtual void doFinalize ( const int loop ) override virtual void doFinalize ( const int loop ) override
{ {
AnalyticalTraits::addEOCErrors( solution(), model(), problem().exactSolution(), adaptIndicator()->sigma() ); //AnalyticalTraits::addEOCErrors( solution(), model(), problem().exactSolution(), adaptIndicator()->sigma() );
} }
......
dune/fem-dg/examples/poisson/algorithmcreator.hh
+ 11
− 4
View file @ dd77f4c7
...@@ -128,14 +128,21 @@ namespace Fem ...@@ -128,14 +128,21 @@ namespace Fem
template< class SigmaDFSpaceType > struct SigmaFunctionChooser template< class SigmaDFSpaceType > struct SigmaFunctionChooser
{ typedef typename AC::template DiscreteFunctions< SigmaDFSpaceType > type; }; { typedef typename AC::template DiscreteFunctions< SigmaDFSpaceType > type; };
typedef PoissonSigmaEstimator< DiscreteFunctionType, SigmaFunctionChooser, typename Operator::AssemblerType, polOrd > //TODO improve -> algorithm configurator
SigmaEstimatorType;
typedef ErrorEstimator< SigmaEstimatorType > EstimatorType;
typedef typename SigmaDiscreteFunctionSelector< DiscreteFunctionType, SigmaFunctionChooser >::type SigmaDiscreteFunctionType;
typedef ErrorEstimator< DiscreteFunctionType, SigmaDiscreteFunctionType, typename Operator::AssemblerType >
ErrorEstimatorType;
typedef PoissonSigmaEstimator< ErrorEstimatorType > SigmaEstimatorType;
typedef PAdaptivity<DFSpaceType, polOrd, SigmaEstimatorType > PAdaptivityType;
public: public:
typedef SubSolverMonitor< SolverMonitor > SolverMonitorType; typedef SubSolverMonitor< SolverMonitor > SolverMonitorType;
typedef SubDiagnostics< Diagnostics > DiagnosticsType; typedef SubDiagnostics< Diagnostics > DiagnosticsType;
typedef PAdaptIndicator< EstimatorType, SigmaEstimatorType, ProblemInterfaceType > AdaptIndicatorType; typedef PAdaptIndicator< PAdaptivityType, ProblemInterfaceType > AdaptIndicatorType;
}; };
template <int polOrd> template <int polOrd>
......
dune/fem-dg/examples/stokes/algorithmcreator.hh
+ 30
− 2
View file @ dd77f4c7
...@@ -146,6 +146,21 @@ namespace Fem ...@@ -146,6 +146,21 @@ namespace Fem
typedef typename AC::template LinearSolvers< DFSpaceType, true> type; typedef typename AC::template LinearSolvers< DFSpaceType, true> type;
}; };
private:
//small helper class
template< class SigmaDFSpaceType > struct SigmaFunctionChooser
{ typedef typename AC::template DiscreteFunctions< SigmaDFSpaceType > type; };
public:
typedef typename SigmaDiscreteFunctionSelector< DiscreteFunctionType, SigmaFunctionChooser >::type SigmaDiscreteFunctionType;
typedef ErrorEstimator< DiscreteFunctionType, SigmaDiscreteFunctionType, typename Operator::AssemblerType >
ErrorEstimatorType;
typedef PoissonSigmaEstimator< ErrorEstimatorType > SigmaEstimatorType;
typedef PAdaptivity<DFSpaceType, polOrd, SigmaEstimatorType > PAdaptivityType;
typedef NoPAdaptIndicator AdaptIndicatorType;
typedef SubSolverMonitor< SolverMonitor > SolverMonitorType; typedef SubSolverMonitor< SolverMonitor > SolverMonitorType;
typedef SubDiagnostics< Diagnostics > DiagnosticsType; typedef SubDiagnostics< Diagnostics > DiagnosticsType;
}; };
...@@ -225,8 +240,21 @@ namespace Fem ...@@ -225,8 +240,21 @@ namespace Fem
static_assert( (int)DFSpaceType::FunctionSpaceType::dimRange == 1 , "pressure dimrange does not fit"); static_assert( (int)DFSpaceType::FunctionSpaceType::dimRange == 1 , "pressure dimrange does not fit");
typedef SubSolverMonitor< SolverMonitor > SolverMonitorType; private:
typedef SubDiagnostics< Diagnostics > DiagnosticsType; typedef typename SubPoissonAlgorithmCreator::template DiscreteTraits<polOrd>::ErrorEstimatorType PoissonErrorEstimatorType;
typedef typename SubPoissonAlgorithmCreator::template DiscreteTraits<polOrd>::SigmaEstimatorType PoissonSigmaEstimatorType;
typedef typename SubPoissonAlgorithmCreator::template DiscreteTraits<polOrd>::PAdaptivityType PoissonPAdaptivityType;
typedef StokesSigmaEstimator< PoissonErrorEstimatorType, typename Operator::AssemblerType > StokesSigmaEstimatorType;
typedef typename StokesSigmaEstimatorType::StokesErrorEstimatorType StokesErrorEstimatorType;
typedef StokesPAdaptivity< PoissonPAdaptivityType,
DFSpaceType, polOrd, StokesSigmaEstimatorType > StokesPAdaptivityType;
public:
typedef SubSolverMonitor< SolverMonitor > SolverMonitorType;
typedef SubDiagnostics< Diagnostics > DiagnosticsType;
typedef StokesPAdaptIndicator< StokesPAdaptivityType, ProblemInterfaceType > AdaptIndicatorType;
}; };
template <int polOrd> template <int polOrd>
......
dune/fem-dg/examples/stokes/stokesalgorithm.hh
+ 266
− 81
View file @ dd77f4c7
...@@ -33,19 +33,27 @@ namespace Dune ...@@ -33,19 +33,27 @@ namespace Dune
namespace Fem namespace Fem
{ {
template< class GridPartImp, class DiscreteVelocityFunctionImp, class DiscretePressureFunctionImp, class SigmaFunctionSpaceImp, class AssemblerImp, class ModelImp, int polOrder>
class StokesSigmaEstimator: public PoissonSigmaEstimator<GridPartImp, DiscreteVelocityFunctionImp, SigmaFunctionSpaceImp, AssemblerImp, polOrder > template< class PoissonErrorEstimatorImp, class AssemblerImp >
class StokesSigmaEstimator
: public PoissonSigmaEstimator< PoissonErrorEstimatorImp >
{ {
typedef PoissonSigmaEstimator<GridPartImp, DiscreteVelocityFunctionImp, SigmaFunctionSpaceImp, AssemblerImp, polOrder > BaseType; typedef PoissonSigmaEstimator< PoissonErrorEstimatorImp > BaseType;
public: public:
typedef DiscreteVelocityFunctionImp DiscreteVelocityFunctionType; typedef typename BaseType::DiscreteFunctionType DiscreteVelocityFunctionType;
typedef DiscretePressureFunctionImp DiscretePressureFunctionType; typedef typename BaseType::DGOperatorType BaseAssemblerType;
typedef SigmaFunctionSpaceImp SigmaFunctionSpaceType; typedef AssemblerImp AssemblerType;
typedef GridPartImp GridPartType; typedef typename AssemblerType::DiscretePressureFunctionType DiscretePressureFunctionType;
typedef typename GridPartType::GridType GridType; typedef typename BaseType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
typedef AssemblerImp AssemblerType; typedef typename BaseType::GridPartType GridPartType;
typedef ModelImp ModelType; typedef typename BaseType::GridType GridType;
typedef typename BaseType::DGOperatorType DGOperatorType;
//typedef typename BaseType::AssemblerType AssemblerType;
static const int polynomialOrder = BaseType::polynomialOrder;
typedef typename AssemblerType::ModelType ModelType;
using BaseType::sigmaSpace_; using BaseType::sigmaSpace_;
using BaseType::sigmaDiscreteFunction_; using BaseType::sigmaDiscreteFunction_;
...@@ -95,16 +103,20 @@ namespace Fem ...@@ -95,16 +103,20 @@ namespace Fem
} }
}; };
struct StokesFlux struct StokesFlux
{ {
typedef typename GridPartType::GridType::template Codim<0>::Entity EntityType; typedef typename GridPartType::GridType::template Codim<0>::Entity EntityType;
typedef typename GridPartType::IntersectionIteratorType IntersectionIteratorType; typedef typename GridPartType::IntersectionIteratorType IntersectionIteratorType;
typedef typename GridPartType::IntersectionType IntersectionType; typedef typename GridPartType::IntersectionType IntersectionType;
typedef typename AssemblerType::FaceQuadratureType FaceQuadratureType; typedef typename AssemblerType::FaceQuadratureType FaceQuadratureType;
typedef typename AssemblerType::VolumeQuadratureType VolumeQuadratureType; typedef typename AssemblerType::VolumeQuadratureType VolumeQuadratureType;
StokesFlux(const DiscretePressureFunctionType &p, const AssemblerType &oper) typedef typename AssemblerType::ContainerType ContainerType;
StokesFlux(const DiscretePressureFunctionType &p, const BaseAssemblerType &oper)
: p_(p), : p_(p),
oper_(oper) oper_(oper)
{} {}
...@@ -195,31 +207,33 @@ namespace Fem ...@@ -195,31 +207,33 @@ namespace Fem
} }
private: private:
const DiscretePressureFunctionType &p_; const DiscretePressureFunctionType &p_;
const AssemblerType &oper_; const BaseAssemblerType &oper_;
}; };
typedef StokesFlux StokesFluxType;
typedef Dune::Fem::LocalFunctionAdapter< SigmaEval<DiscreteVelocityFunctionType,DiscretePressureFunctionType,BaseAssemblerType> > StokesEstimateFunction;
typedef StokesErrorEstimator< DiscreteVelocityFunctionType, StokesEstimateFunction, StokesFluxType > StokesErrorEstimatorType;
typedef Dune::Fem::LocalFunctionAdapter< SigmaEval<DiscreteVelocityFunctionType,DiscretePressureFunctionType,AssemblerType> > StokesEstimateFunction; typedef Dune::Fem::LocalFunctionAdapter< StokesErrorEstimatorType > StokesEstimateDataType;
typedef StokesErrorEstimator< DiscreteVelocityFunctionType, StokesEstimateFunction, StokesFlux > StokesEstimatorType;
typedef Dune::Fem::LocalFunctionAdapter< StokesEstimatorType > StokesEstimateDataType;
typedef typename AssemblerType::ContainerType ContainerType;
StokesSigmaEstimator( GridPartType& gridPart, StokesSigmaEstimator( ContainerType& container,
const DiscreteVelocityFunctionType& solution, const BaseAssemblerType& ass,
const DiscretePressureFunctionType& pressureSolution,
const AssemblerType& assembler, const AssemblerType& assembler,
const std::string keyPrefix = "" ) const std::string keyPrefix = "" )
: BaseType( gridPart, solution, assembler, keyPrefix ), : BaseType( container.template adapter<0>(), ass, keyPrefix ),
stkFlux_( pressureSolution, assembler_), container_( container ),
stkLocalEstimate_( solution_, pressureSolution, assembler_ ), stkFlux_( *container_.template solution<1>(), ass ),
stkEstimateFunction_("stokes estimate", stkLocalEstimate_, gridPart_, solution_.space().order() ), stkLocalEstimate_( *container_.template solution<0>(), *container_.template solution<1>(), ass ),
stkEstimator_( solution_, stkEstimateFunction_, stkFlux_, gridPart_.grid() ), stkEstimateFunction_("stokes estimate", stkLocalEstimate_, container_.template solution<0>()->gridPart(), container_.template space<0>().order() ),
stkEstimateData_("stokesEstimator", stkEstimator_, gridPart_, solution_.space().order() ) stkEstimator_( *container_.template solution<0>(), stkFlux_, stkEstimateFunction_ ),
stkEstimateData_("stokesEstimator", stkEstimator_, container_.template solution<0>()->gridPart(), container_.template space<0>().order() )
{} {}
StokesEstimatorType& estimator() const StokesErrorEstimatorType& estimator() const
{ {
return stkEstimator_; return stkEstimator_;
}; };
...@@ -231,16 +245,228 @@ namespace Fem ...@@ -231,16 +245,228 @@ namespace Fem
private: private:
StokesFlux stkFlux_; ContainerType& container_;
SigmaEval<DiscreteVelocityFunctionType,DiscretePressureFunctionType,AssemblerType> stkLocalEstimate_; StokesFluxType stkFlux_;
SigmaEval<DiscreteVelocityFunctionType,DiscretePressureFunctionType,BaseAssemblerType> stkLocalEstimate_;
StokesEstimateFunction stkEstimateFunction_; StokesEstimateFunction stkEstimateFunction_;
StokesEstimatorType stkEstimator_; StokesErrorEstimatorType stkEstimator_;
StokesEstimateDataType stkEstimateData_; StokesEstimateDataType stkEstimateData_;
}; };
template< class PAdaptivityImp, class DiscreteFunctionSpaceImp, int polOrder, class SigmaEstimatorImp >
class StokesPAdaptivity
{
public:
typedef PAdaptivityImp PAdaptivityType;
typedef DiscreteFunctionSpaceImp DiscreteFunctionSpaceType;
typedef SigmaEstimatorImp SigmaEstimatorType;
typedef typename SigmaEstimatorType::ErrorEstimatorType ErrorEstimatorType;
typedef typename DiscreteFunctionSpaceType::GridType GridType;
typedef typename SigmaEstimatorType::DGOperatorType DGOperatorType;
typedef typename SigmaEstimatorType::AssemblerType AssemblerType;
typedef typename SigmaEstimatorType::BaseAssemblerType BaseAssemblerType;
typedef typename SigmaEstimatorType::ContainerType ContainerType;
typedef typename PAdaptivityType::PolOrderContainer PolOrderContainer;
StokesPAdaptivity( ContainerType& container, BaseAssemblerType& ass, AssemblerType& assembler, const std::string name = "" )
: pAdapt_( container.template adapter<0>(), ass, name ),
polOrderContainer_( container.template solution<0>()->gridPart().grid(), 0 ),
space_( container.template space<1>() ),
sigmaEstimator_( container, ass, assembler, name ),
errorEstimator_( *container.template solution<0>(), ass, sigmaEstimator_.sigma() ),
param_( AdaptationParameters() )
{
#ifdef PADAPTSPACE
// we start with max order
typedef typename PolOrderContainer::Iterator Iterator ;
const Iterator end = polOrderContainer_.end();
const int minimalOrder = errorEstimator_.minimalOrder() ;
for( Iterator it = polOrderContainer_.begin(); it != end; ++it )
{
(*it).value() = minimalOrder ;
}
#endif
}
bool adaptive() const
{
return errorEstimator_.isPadaptive() && pAdapt_.adaptive();
}
void prepare()
{
pAdapt_.prepare();
#ifdef PADAPTSPACE
const int minimalOrder = errorEstimator_.minimalOrder();
// only implemented for PAdaptiveSpace
std::vector<int> polOrderVec( space_.gridPart().indexSet().size(0) );
std::fill( polOrderVec.begin(), polOrderVec.end(), polOrder );
polOrderContainer_.resize();
if ( errorEstimator_.isPadaptive() )
{
typedef typename DiscreteFunctionSpaceType::IteratorType IteratorType;
typedef typename IteratorType::Entity EntityType ;
const IteratorType end = space_.end();
for( IteratorType it = space_.begin(); it != end; ++it )
{
const EntityType& entity = *it;
int order = polOrderContainer_[ entity ].value();
while (order == -1) // is a new element
{
if ( entity.level() == 0)
order = minimalOrder;
else
{
// don't call father twice
order = polOrderContainer_[ entity.father() ].value();
assert(order > 0);
}
}
polOrderVec[space_.gridPart().indexSet().index(entity)] = order;
}
}
space_.adapt( polOrderVec );
#endif
}
template< class ProblemImp >
bool estimateMark( const ProblemImp& problem )
{
return pAdapt_.estimateMark( problem );
//note: no h-Adaptation regarding pressure space!?
}
void closure()
{
pAdapt_.closure();
//note: no closure regarding pressure space!?
}
ErrorEstimatorType& errorEstimator()
{
return errorEstimator_;
}
SigmaEstimatorType& sigmaEstimator()
{
return sigmaEstimator_;
}
private:
PAdaptivityType pAdapt_;
PolOrderContainer polOrderContainer_;
const DiscreteFunctionSpaceType& space_;
SigmaEstimatorType sigmaEstimator_;
ErrorEstimatorType errorEstimator_;
AdaptationParameters param_;
};
/**
* \brief Adaptation indicator doing no indication and marking of the entities.
*/
template< class PAdaptivityImp, class ProblemImp >
class StokesPAdaptIndicator
{
protected:
typedef PAdaptivityImp PAdaptivityType;
typedef ProblemImp ProblemType;
typedef typename PAdaptivityType::SigmaEstimatorType SigmaEstimatorType;
typedef typename PAdaptivityType::ErrorEstimatorType ErrorEstimatorType;
typedef typename SigmaEstimatorType::DiscreteFunctionType DiscreteFunctionType;
typedef typename SigmaEstimatorType::AssemblerType AssemblerType;
typedef typename SigmaEstimatorType::BaseAssemblerType BaseAssemblerType;
typedef typename SigmaEstimatorType::DGOperatorType DGOperatorType;
typedef typename AssemblerType::ContainerType ContainerType;
typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
typedef typename DiscreteFunctionType::GridPartType GridPartType;
typedef typename GridPartType::GridType GridType;
static const int polOrder = SigmaEstimatorType::polynomialOrder;
enum { dimension = GridType::dimension };
public:
typedef uint64_t UInt64Type;
StokesPAdaptIndicator( ContainerType& container,
BaseAssemblerType& ass,
AssemblerType& assembler,
const ProblemType& problem,
const std::string name = "" )
: pAdapt_( container, ass, assembler, name ),
problem_( problem )
{}
bool adaptive() const
{
return pAdapt_.adaptive();
}
size_t numberOfElements() const
{
return 0;
}
UInt64Type globalNumberOfElements() const { return 0; }
template< class TimeProviderImp >
void setAdaptation( TimeProviderImp& tp )
{
}
void preAdapt()
{
pAdapt_.prepare();
}
void estimateMark( const bool initialAdapt = false )
{
// calculate new sigma
pAdapt_.sigmaEstimator().update();
const bool marked = pAdapt_.estimateMark( problem_ );
if( marked )
pAdapt_.closure();
}
void postAdapt(){}
void finalize(){}
int minNumberOfElements() const { return 0; }
int maxNumberOfElements() const { return 0; }
const int finestLevel() const { return 0; }
// return some info
const typename SigmaEstimatorType::SigmaDiscreteFunctionType& sigma()
{
return pAdapt_.sigmaEstimator().sigma();
}
private:
PAdaptivityType pAdapt_;
const ProblemType& problem_;
};
/** /**
* \brief Algorithm for solving the Stokes equation. * \brief Algorithm for solving the Stokes equation.
* *
...@@ -284,14 +510,9 @@ namespace Fem ...@@ -284,14 +510,9 @@ namespace Fem
enum { dimension = GridType::dimension }; enum { dimension = GridType::dimension };
typedef typename DiscreteVelocityFunctionSpaceType :: typedef typename BaseType::AdaptIndicatorType AdaptIndicatorType;
template ToNewDimRange< dimension * ElliptProblemTraits::AnalyticalTraits::ModelType::dimRange >::NewFunctionSpaceType
SigmaFunctionSpaceType;
typedef StokesSigmaEstimator< GridPartType, DiscreteVelocityFunctionType, DiscreteFunctionType, typedef typename BaseType::AdaptationDiscreteFunctionType AdaptationDiscreteFunctionType;
SigmaFunctionSpaceType, typename EllipticalAlgorithmType::AssemblerType,
typename ElliptProblemTraits::AnalyticalTraits::ModelType, polOrd >
StokesSigmaEstimatorType;
typedef CovariantTuple< typename BaseType::IOTupleType::type, typename EllipticalAlgorithmType::IOTupleType::type > typedef CovariantTuple< typename BaseType::IOTupleType::type, typename EllipticalAlgorithmType::IOTupleType::type >
IOTupleType; IOTupleType;
...@@ -321,40 +542,19 @@ namespace Fem ...@@ -321,40 +542,19 @@ namespace Fem
explicit SubStokesAlgorithm( GridType& grid, ContainerType& container ) : explicit SubStokesAlgorithm( GridType& grid, ContainerType& container ) :
BaseType( grid, container.template adapter<1>() ), BaseType( grid, container.template adapter<1>() ),
container_( container ), container_( container ),
space_( container_.template adapter<1>().space() ), space_( container_.template space<1>() ),
assembler_( container_, model() ), assembler_( container_, model() ),
ellAlg_( grid, container_.template adapter<0>() ), ellAlg_( grid, container_.template adapter<0>() ),
ioTuple_( new IOTupleType( *BaseType::dataTuple(), *ellAlg_.dataTuple() ) ), adaptIndicator_( Std::make_unique<AdaptIndicatorType>( container_, ellAlg_.assembler(), assembler_, problem(), name() ) ),
stokesSigmaEstimator_( new StokesSigmaEstimatorType( container_.template adapter<1>().gridPart(), ellAlg_.solution(), solution(), ellAlg_.assembler(), name() ) ) ioTuple_( new IOTupleType( *BaseType::dataTuple(), *ellAlg_.dataTuple() ) )
{ {
} }
virtual IOTupleType& dataTuple () virtual IOTupleType& dataTuple () override
{ {
return *ioTuple_; return *ioTuple_;
} }
bool adaptation(const double tolerance)
{
#if 0
#if PADAPTSPACE
// update to current
polOrderContainer_.resize();
const double error = stokesSigmaEstimator_.estimator().estimate( problem() );
typedef typename DiscreteVelocityFunctionSpaceType::IteratorType IteratorType;
const IteratorType end = ellAlg_.space().end();
for( IteratorType it = ellAlg_.space().begin(); it != end; ++it )
{
const typename IteratorType::Entity &entity = *it;
polOrderContainer_[entity].value() =
stokesSigmaEstimator_.estimator().newOrder( 0.98*tolerance, entity );
}
#endif
#endif
return (error < std::abs(tolerance) ? false : stokesSigmaEstimator_->estimator().mark( 0.98 * tolerance));
}
private: private:
virtual std::shared_ptr< SolverType > doCreateSolver() override virtual std::shared_ptr< SolverType > doCreateSolver() override
{ {
...@@ -369,18 +569,6 @@ namespace Fem ...@@ -369,18 +569,6 @@ namespace Fem
std::cout << "Solver (Stokes) assemble time: " << timer.elapsed() << std::endl; std::cout << "Solver (Stokes) assemble time: " << timer.elapsed() << std::endl;
double absLimit = Dune::Fem::Parameter::getValue<double>("istl.absLimit",1.e-10); double absLimit = Dune::Fem::Parameter::getValue<double>("istl.absLimit",1.e-10);
#if 0
#ifdef PADAPTSPACE
int polOrder = Dune::Fem::Parameter::getValue<double>("femdg.polOrd",1);
// only implemented for PAdaptiveSpace
std::vector<int> polOrderVec( ellAlg_.space().gridPart().indexSet().size(0) );
std::vector<int> polOrderVecPressure( space_.gridPart().indexSet().size(0) );
std::fill( polOrderVec.begin(), polOrderVec.end(), polOrder );
std::fill( polOrderVecPressure.begin(), polOrderVecPressure.end(), polOrder-1 );
ellAlg_.space().adapt( polOrderVec );
space_.adapt( polOrderVecPressure);
#endif
#endif
return std::make_shared< SolverType >( container_, *ellAlg_.solver(), absLimit, 3*ellAlg_.solution().space().size() ); return std::make_shared< SolverType >( container_, *ellAlg_.solver(), absLimit, 3*ellAlg_.solution().space().size() );
} }
...@@ -399,9 +587,6 @@ namespace Fem ...@@ -399,9 +587,6 @@ namespace Fem
virtual void doSolve( const int loop ) override virtual void doSolve( const int loop ) override
{ {
BaseType::doSolve( loop ); BaseType::doSolve( loop );
// TODO check wheather we need the following line
stokesSigmaEstimator_->update();
} }
virtual void doPostSolve( const int loop ) override virtual void doPostSolve( const int loop ) override
...@@ -414,7 +599,7 @@ namespace Fem ...@@ -414,7 +599,7 @@ namespace Fem
virtual void doFinalize( const int loop ) override virtual void doFinalize( const int loop ) override
{ {
ellAlg_.finalize( loop ); ellAlg_.finalize( loop );
AnalyticalTraits::addEOCErrors( solution(), ellAlg_.model(), problem().get<1>().exactSolution() ); AnalyticalTraits::addEOCErrors( solution(), ellAlg_.model(), problem().template get<1>().exactSolution() );
} }
protected: protected:
...@@ -423,8 +608,8 @@ namespace Fem ...@@ -423,8 +608,8 @@ namespace Fem
AssemblerType assembler_; AssemblerType assembler_;
EllipticalAlgorithmType ellAlg_; EllipticalAlgorithmType ellAlg_;
std::unique_ptr< AdaptIndicatorType > adaptIndicator_;
std::unique_ptr< IOTupleType > ioTuple_; std::unique_ptr< IOTupleType > ioTuple_;
std::unique_ptr< StokesSigmaEstimatorType > stokesSigmaEstimator_;
}; };
......
dune/fem-dg/examples/stokes/stokesassembler.hh
+ 4
− 2
View file @ dd77f4c7
...@@ -44,7 +44,7 @@ namespace Fem ...@@ -44,7 +44,7 @@ namespace Fem
virtual double relaxation () const virtual double relaxation () const
{ {
return Fem::Parameter::getValue< int >( keyPrefix_ + "preconditioning.relaxation", 1.1 ); return Fem::Parameter::getValue< double >( keyPrefix_ + "preconditioning.relaxation", 1.1 );
} }
virtual int method () const virtual int method () const
...@@ -291,8 +291,10 @@ namespace Fem ...@@ -291,8 +291,10 @@ namespace Fem
typedef typename GridType::template Codim<0>::Geometry GeometryType; typedef typename GridType::template Codim<0>::Geometry GeometryType;
//! type of quadrature to be used //! type of quadrature to be used
public:
typedef typename OperatorTraits::VolumeQuadratureType VolumeQuadratureType; typedef typename OperatorTraits::VolumeQuadratureType VolumeQuadratureType;
typedef typename OperatorTraits::FaceQuadratureType FaceQuadratureType; typedef typename OperatorTraits::FaceQuadratureType FaceQuadratureType;
protected:
typedef typename ContainerType::template Matrix<0,1> PressureGradMatType; typedef typename ContainerType::template Matrix<0,1> PressureGradMatType;
typedef typename ContainerType::template Matrix<1,0> PressureDivMatType; typedef typename ContainerType::template Matrix<1,0> PressureDivMatType;
...@@ -738,7 +740,7 @@ namespace Fem ...@@ -738,7 +740,7 @@ namespace Fem
double intWeight=faceQuadInner.weight(l); double intWeight=faceQuadInner.weight(l);
DomainType quadInEn=edge.geometry().global(faceQuadInner.localPoint(l)); DomainType quadInEn=edge.geometry().global(faceQuadInner.localPoint(l));
problem_.get<0>().g(quadInEn,dirichletValue); problem_.template get<0>().g(quadInEn,dirichletValue);
double pressureDirichlet; double pressureDirichlet;
for(unsigned int n = 0; n <(unsigned int)localPressure.numDofs() ; ++n) for(unsigned int n = 0; n <(unsigned int)localPressure.numDofs() ; ++n)
......
dune/fem-dg/operator/adaptation/poissonestimator.hh
+ 86
− 33
View file @ dd77f4c7
...@@ -91,6 +91,59 @@ namespace Fem ...@@ -91,6 +91,59 @@ namespace Fem
} }
}; };
template< class DiscreteFunctionImp, template<class> class SigmaDiscreteFunctionChooserImp >
class SigmaDiscreteFunctionSelector
{
typedef DiscreteFunctionImp DiscreteFunctionType;
typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
typedef typename DiscreteFunctionType::GridPartType GridPartType;
typedef typename GridPartType::GridType GridType;
static const int polOrder = DiscreteFunctionSpaceType::polynomialOrder;
typedef typename DiscreteFunctionSpaceType::
template ToNewDimRange< GridType::dimension * DiscreteFunctionSpaceType::FunctionSpaceType::dimRange >::NewFunctionSpaceType
SigmaFunctionSpaceType;
//- hybrid spaces use PAdaptiveDGSpace
template <class Grid, int topoId>
struct SigmaSpaceChooser
{
typedef Fem::PAdaptiveDGSpace< SigmaFunctionSpaceType, GridPartType, polOrder > type;
};
//- cubes use LegendreDGSpace
template <class Grid>
struct SigmaSpaceChooser< Grid, 1 >
{
typedef Dune::Fem::LegendreDiscontinuousGalerkinSpace< SigmaFunctionSpaceType, GridPartType, polOrder > type;
};
//- cubes use OrthonormalDGSpace
template <class Grid>
struct SigmaSpaceChooser< Grid, 0 >
{
typedef Dune::Fem::DiscontinuousGalerkinSpace< SigmaFunctionSpaceType, GridPartType, polOrder > type;
};
// work arround internal compiler error
enum { simplexTopoId = GenericGeometry::SimplexTopology< GridType::dimension >::type::id,
cubeTopoId = GenericGeometry::CubeTopology< GridType::dimension >::type::id,
myTopo = Dune::Capabilities::hasSingleGeometryType< GridType >::topologyId
};
enum { topoId = (simplexTopoId == myTopo) ? 0 : // 0 = simplex, 1 = cube, -1 = hybrid
(myTopo == cubeTopoId) ? 1 : -1 };
typedef typename SigmaSpaceChooser< GridType, topoId >::type SigmaDiscreteFunctionSpaceType;
public:
typedef typename SigmaDiscreteFunctionChooserImp<SigmaDiscreteFunctionSpaceType>::type
type;
};
// Estimator // Estimator
// --------- // ---------
// Template Arguments: // Template Arguments:
...@@ -106,18 +159,21 @@ namespace Fem ...@@ -106,18 +159,21 @@ namespace Fem
// fluxEn, dfluxEn, fluxNb, dfluxNb): // fluxEn, dfluxEn, fluxNb, dfluxNb):
// method to compute -hatK, only fluxEn and fluxNb is used // method to compute -hatK, only fluxEn and fluxNb is used
template<class SigmaEstimatorImp> template< class DiscreteFunctionImp, class SigmaDiscreteFunctionImp, class DGOperatorImp >
class ErrorEstimator class ErrorEstimator
{ {
typedef ErrorEstimator< SigmaEstimatorImp> ThisType; typedef ErrorEstimator< DiscreteFunctionImp, SigmaDiscreteFunctionImp, DGOperatorImp > ThisType;
public: public:
typedef SigmaEstimatorImp SigmaEstimatorType; typedef DiscreteFunctionImp DiscreteFunctionType;
typedef typename SigmaEstimatorType::DiscreteFunctionType DiscreteFunctionType;
typedef typename SigmaEstimatorType::SigmaDiscreteFunctionType SigmaDiscreteFunctionType;
typedef typename SigmaEstimatorType::DGOperatorType DGOperatorType;
typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType; typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
typedef typename DiscreteFunctionType::GridPartType GridPartType;
typedef typename GridPartType::GridType GridType;
typedef DGOperatorImp DGOperatorType;
typedef SigmaDiscreteFunctionImp SigmaDiscreteFunctionType;
typedef typename SigmaDiscreteFunctionType::DiscreteFunctionSpaceType SigmaDiscreteFunctionSpaceType;
typedef typename DiscreteFunctionType::LocalFunctionType LocalFunctionType; typedef typename DiscreteFunctionType::LocalFunctionType LocalFunctionType;
typedef typename SigmaDiscreteFunctionType::LocalFunctionType SigmaLocalFunctionType; typedef typename SigmaDiscreteFunctionType::LocalFunctionType SigmaLocalFunctionType;
typedef typename SigmaLocalFunctionType::RangeType GradientRangeType; typedef typename SigmaLocalFunctionType::RangeType GradientRangeType;
...@@ -127,13 +183,11 @@ namespace Fem ...@@ -127,13 +183,11 @@ namespace Fem
typedef typename DiscreteFunctionSpaceType::DomainType DomainType; typedef typename DiscreteFunctionSpaceType::DomainType DomainType;
typedef typename DiscreteFunctionSpaceType::RangeType RangeType; typedef typename DiscreteFunctionSpaceType::RangeType RangeType;
typedef typename DiscreteFunctionSpaceType::JacobianRangeType JacobianRangeType; typedef typename DiscreteFunctionSpaceType::JacobianRangeType JacobianRangeType;
typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
typedef typename DiscreteFunctionSpaceType::IteratorType IteratorType; typedef typename DiscreteFunctionSpaceType::IteratorType IteratorType;
typedef Dune::Fem::FieldMatrixConverter< GradientRangeType, JacobianRangeType > typedef Dune::Fem::FieldMatrixConverter< GradientRangeType, JacobianRangeType >
SigmaConverterType; SigmaConverterType;
typedef typename GridPartType::GridType GridType;
typedef typename GridPartType::IndexSetType IndexSetType; typedef typename GridPartType::IndexSetType IndexSetType;
typedef typename GridPartType::IntersectionIteratorType IntersectionIteratorType; typedef typename GridPartType::IntersectionIteratorType IntersectionIteratorType;
...@@ -186,19 +240,18 @@ namespace Fem ...@@ -186,19 +240,18 @@ namespace Fem
public: public:
ErrorEstimator (const DiscreteFunctionType &uh, ErrorEstimator (DiscreteFunctionType& uh,
const SigmaDiscreteFunctionType &sigma,
const DGOperatorType &oper, const DGOperatorType &oper,
GridType &grid, const SigmaDiscreteFunctionType &sigma,
const AdaptationParameters& param = AdaptationParameters() ) const AdaptationParameters& param = AdaptationParameters() )
: oper_(oper), : oper_(oper),
uh_( uh ), uh_( uh ),
sigma_( sigma ), sigma_( sigma ),
beta_(1.), beta_(1.),
dfSpace_( uh.space() ), dfSpace_( uh_.space() ),
gridPart_( dfSpace_.gridPart() ), gridPart_( dfSpace_.gridPart() ),
indexSet_( gridPart_.indexSet() ), indexSet_( gridPart_.indexSet() ),
grid_( grid ), grid_( gridPart_.grid() ),
dgSimplexSpace_(), dgSimplexSpace_(),
dgCubeSpace_(), dgCubeSpace_(),
indicator_( indexSet_.size( 0 )), indicator_( indexSet_.size( 0 )),
...@@ -548,18 +601,18 @@ namespace Fem ...@@ -548,18 +601,18 @@ namespace Fem
} }
void calcSmoothness( int polOrder, const ElementType &entity, void calcSmoothness( int pOrder, const ElementType &entity,
const typename ElementType::Geometry &geometry, const typename ElementType::Geometry &geometry,
const LocalFunctionType &uLocal, double &smoothness ) const LocalFunctionType &uLocal, double &smoothness )
{ {
if (geometry.type().isSimplex()) if (geometry.type().isSimplex())
calcSmoothness<HierarchicSimplexDGSpaceType>(dgSimplexSpace(),polOrder,entity,geometry,uLocal,smoothness); calcSmoothness<HierarchicSimplexDGSpaceType>(dgSimplexSpace(),pOrder,entity,geometry,uLocal,smoothness);
else else
calcSmoothness<HierarchicCubeDGSpaceType>(dgCubeSpace(),polOrder,entity,geometry,uLocal,smoothness); calcSmoothness<HierarchicCubeDGSpaceType>(dgCubeSpace(),pOrder,entity,geometry,uLocal,smoothness);
} }
template <class HDGSpaceType> template <class HDGSpaceType>
void calcSmoothness( const typename HDGSpaceType::Type &dgSpace, void calcSmoothness( const typename HDGSpaceType::Type &dgSpace,
int polOrder, const ElementType &entity, int pOrder, const ElementType &entity,
const typename ElementType::Geometry &geometry, const typename ElementType::Geometry &geometry,
const LocalFunctionType &uLocal, double &smoothness ) const LocalFunctionType &uLocal, double &smoothness )
{ {
...@@ -585,14 +638,14 @@ namespace Fem ...@@ -585,14 +638,14 @@ namespace Fem
localMassMatrix.applyInverse( udg ); localMassMatrix.applyInverse( udg );
assert(RangeType::dimension == 1); assert(RangeType::dimension == 1);
// only dofs up to bstart[polOrder] should be used // only dofs up to bstart[pOrder] should be used
// bool correctProjection = true; // bool correctProjection = true;
for (int i=HDGSpaceType::size(polOrder); i<udg.numDofs(); ++i) for (int i=HDGSpaceType::size(pOrder); i<udg.numDofs(); ++i)
{ {
if ( std::abs(udg[i]) > 1e-10 ) if ( std::abs(udg[i]) > 1e-10 )
{ {
std::cout << "for i=" << HDGSpaceType::size(polOrder) << " to " << udg.numDofs() ; std::cout << "for i=" << HDGSpaceType::size(pOrder) << " to " << udg.numDofs() ;
std::cout << " projection wrong for polorder = " << polOrder std::cout << " projection wrong for polorder = " << pOrder
<< " dof[" << i << "] = " << udg[i] << std::endl; << " dof[" << i << "] = " << udg[i] << std::endl;
//correctProjection = false; //correctProjection = false;
} }
...@@ -603,35 +656,35 @@ namespace Fem ...@@ -603,35 +656,35 @@ namespace Fem
{ {
double a = 0; double a = 0;
double vol = geometry.volume(); double vol = geometry.volume();
for (int i=HDGSpaceType::size(polOrder-1) ; i<udg.numDofs(); ++i) for (int i=HDGSpaceType::size(pOrder-1) ; i<udg.numDofs(); ++i)
{ {
assert( udg.numDofs() > i ); assert( udg.numDofs() > i );
a += std::abs(udg[i]); a += std::abs(udg[i]);
} }
a *= sqrt(vol); a *= sqrt(vol);
if (std::abs(a)<1e-8 || polOrder==1) if (std::abs(a)<1e-8 || pOrder==1)
smoothness = -1; // polOrder+1; smoothness = -1; // pOrder+1;
else else
smoothness = log( (2.*polOrder+1)/(2.*a*a) ) / (2.*log(polOrder))-1; smoothness = log( (2.*pOrder+1)/(2.*a*a) ) / (2.*log(pOrder))-1;
} }
else if (padaptiveMethod_ == prior2p) else if (padaptiveMethod_ == prior2p)
{ {
if (polOrder<3) if (pOrder<3)
smoothness = polOrder+1; smoothness = pOrder+1;
else else
{ {
typedef Dune::Fem::H1Norm<GridPartType> H1NormType ; typedef Dune::Fem::H1Norm<GridPartType> H1NormType ;
typedef typename H1NormType :: IntegratorType IntegratorType ; typedef typename H1NormType :: IntegratorType IntegratorType ;
IntegratorType integrator(2*polOrder); IntegratorType integrator(2*pOrder);
typedef typename H1NormType :: template FunctionJacobianSquare< DGLFType > FctJacType; typedef typename H1NormType :: template FunctionJacobianSquare< DGLFType > FctJacType;
FctJacType udg2( udg ); FctJacType udg2( udg );
typename FctJacType::RangeType eta1(0),eta2(0); typename FctJacType::RangeType eta1(0),eta2(0);
for (int i=0 ; i<HDGSpaceType::size(polOrder-2); ++i) for (int i=0 ; i<HDGSpaceType::size(pOrder-2); ++i)
udg[i] = 0; udg[i] = 0;
integrator.integrate( entity, udg2, eta2 ); // square of H1 norm integrator.integrate( entity, udg2, eta2 ); // square of H1 norm
for (int i=HDGSpaceType::size(polOrder-2) ; i<HDGSpaceType::size(polOrder-1); ++i) for (int i=HDGSpaceType::size(pOrder-2) ; i<HDGSpaceType::size(pOrder-1); ++i)
udg[i] = 0; udg[i] = 0;
integrator.integrate( entity, udg2, eta1 ); // square of H1 norm integrator.integrate( entity, udg2, eta1 ); // square of H1 norm
if (std::abs(eta1) < 1e-12 || std::abs(eta2) < 1e-12) if (std::abs(eta1) < 1e-12 || std::abs(eta2) < 1e-12)
...@@ -641,11 +694,11 @@ namespace Fem ...@@ -641,11 +694,11 @@ namespace Fem
// approximated with a much lower polynomial (2 order // approximated with a much lower polynomial (2 order
// are small) - so coarsen // are small) - so coarsen
else else
smoothness = polOrder+1; smoothness = pOrder+1;
} }
else else
{ {
smoothness = -0.5*std::log(eta1[0]/eta2[0]) / std::log( double(polOrder-1)/double(polOrder-2) ); smoothness = -0.5*std::log(eta1[0]/eta2[0]) / std::log( double(pOrder-1)/double(pOrder-2) );
assert( smoothness == smoothness && smoothness < 100.); assert( smoothness == smoothness && smoothness < 100.);
} }
} }
......
dune/fem-dg/operator/adaptation/stokesestimator.hh
+ 40
− 38
View file @ dd77f4c7
...@@ -26,42 +26,45 @@ namespace Fem ...@@ -26,42 +26,45 @@ namespace Fem
// fluxEn, dfluxEn, fluxNb, dfluxNb): // fluxEn, dfluxEn, fluxNb, dfluxNb):
// method to compute -hatK, only fluxEn and fluxNb is used // method to compute -hatK, only fluxEn and fluxNb is used
template<class UFunction, class SigmaFunction, class DGOperator> template< class DiscreteFunctionImp, class SigmaDiscreteFunctionImp, class DGOperatorImp >
class StokesErrorEstimator : public ErrorEstimator<UFunction,SigmaFunction,DGOperator> class StokesErrorEstimator
: public ErrorEstimator< DiscreteFunctionImp, SigmaDiscreteFunctionImp, DGOperatorImp >
{ {
typedef ErrorEstimator< UFunction, SigmaFunction, DGOperator> BaseType; typedef ErrorEstimator< DiscreteFunctionImp, SigmaDiscreteFunctionImp, DGOperatorImp >
typedef StokesErrorEstimator< UFunction, SigmaFunction, DGOperator> ThisType; BaseType;
public: public:
typedef UFunction DiscreteFunctionType; typedef typename BaseType::DGOperatorType DGOperatorType;
typedef typename BaseType::DiscreteFunctionType DiscreteFunctionType;
typedef typename DiscreteFunctionType :: DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
typedef typename DiscreteFunctionType :: LocalFunctionType LocalFunctionType; typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
typedef typename SigmaFunction :: LocalFunctionType SigmaLocalFunctionType; typedef typename DiscreteFunctionType::LocalFunctionType LocalFunctionType;
typedef typename BaseType::SigmaDiscreteFunctionType SigmaDiscreteFunctionType;
typedef typename DiscreteFunctionSpaceType :: DomainFieldType DomainFieldType; typedef typename BaseType::LocalFunctionType SigmaLocalFunctionType;
typedef typename DiscreteFunctionSpaceType :: RangeFieldType RangeFieldType;
typedef typename DiscreteFunctionSpaceType :: DomainType DomainType; typedef typename DiscreteFunctionSpaceType::DomainFieldType DomainFieldType;
typedef typename DiscreteFunctionSpaceType :: RangeType RangeType; typedef typename DiscreteFunctionSpaceType::RangeFieldType RangeFieldType;
typedef typename DiscreteFunctionSpaceType :: JacobianRangeType JacobianRangeType; typedef typename DiscreteFunctionSpaceType::DomainType DomainType;
typedef typename DiscreteFunctionSpaceType :: GridPartType GridPartType; typedef typename DiscreteFunctionSpaceType::RangeType RangeType;
typedef typename DiscreteFunctionSpaceType :: IteratorType IteratorType; typedef typename DiscreteFunctionSpaceType::JacobianRangeType JacobianRangeType;
typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
typedef typename GridPartType :: GridType GridType; typedef typename DiscreteFunctionSpaceType::IteratorType IteratorType;
typedef typename GridPartType :: IndexSetType IndexSetType;
typedef typename GridPartType :: IntersectionIteratorType IntersectionIteratorType; typedef typename GridPartType::GridType GridType;
typedef typename GridPartType::IndexSetType IndexSetType;
typedef typename IntersectionIteratorType :: Intersection IntersectionType; typedef typename GridPartType::IntersectionIteratorType IntersectionIteratorType;
typedef typename GridType :: template Codim< 0 > :: Entity ElementType; typedef typename IntersectionIteratorType::Intersection IntersectionType;
typedef typename GridType :: template Codim< 0 > :: EntityPointer ElementPointerType;
typedef typename ElementType::Geometry GeometryType; typedef typename GridType::template Codim< 0 >::Entity ElementType;
typedef typename GridType::template Codim< 0 >::EntityPointer ElementPointerType;
static const int dimension = GridType :: dimension; typedef typename ElementType::Geometry GeometryType;
typedef FieldMatrix<double,dimension,dimension> JacobianInverseType;
typedef Dune::Fem::CachingQuadrature< GridPartType, 0 > ElementQuadratureType; static const int dimension = GridType::dimension;
typedef Dune::Fem::CachingQuadrature< GridPartType, 1 > FaceQuadratureType; typedef FieldMatrix<double,dimension,dimension> JacobianInverseType;
typedef std :: vector< double > ErrorIndicatorType; typedef Dune::Fem::CachingQuadrature< GridPartType, 0 > ElementQuadratureType;
typedef Dune::Fem::CachingQuadrature< GridPartType, 1 > FaceQuadratureType;
typedef std::vector< double > ErrorIndicatorType;
private: private:
using BaseType::uh_; using BaseType::uh_;
...@@ -76,12 +79,11 @@ namespace Fem ...@@ -76,12 +79,11 @@ namespace Fem
using BaseType::theta_; using BaseType::theta_;
typename BaseType::ErrorIndicatorType Rdiv_; typename BaseType::ErrorIndicatorType Rdiv_;
public: public:
StokesErrorEstimator (const DiscreteFunctionType &uh, StokesErrorEstimator (DiscreteFunctionType& df,
const SigmaFunction &sigma, const DGOperatorType &oper,
const DGOperator &oper, const SigmaDiscreteFunctionType &sigma,
GridType &grid,
const AdaptationParameters& param = AdaptationParameters() ) const AdaptationParameters& param = AdaptationParameters() )
: BaseType(uh,sigma,oper,grid,param), : BaseType(df,oper,sigma,param),
Rdiv_( this->indexSet_.size( 0 )) Rdiv_( this->indexSet_.size( 0 ))
{ {
} }
......
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