Commit 02b2a11f authored by Robert K's avatar Robert K

[!304] Bugfix/inverse operator if

Merge branch 'bugfix/inverse-operator-if' into 'master'

See merge request [!304]

  [!304]: Nonedune-fem/dune-fem/merge_requests/304
parents 93467cf8 669f7428
Pipeline #17668 passed with stage
in 44 minutes and 42 seconds
......@@ -17,6 +17,8 @@
#include <dune/fem/function/localfunction/mutable.hh>
#include <dune/fem/misc/petsc/petscvector.hh>
#include <dune/fem/space/common/restrictprolonginterface.hh>
namespace Dune
{
namespace Fem
......@@ -204,6 +206,123 @@ namespace Dune
DofVectorType& dofVector_;
};
template <class DiscreteFunctionSpace>
class RestrictProlongDefault< PetscDiscreteFunction< DiscreteFunctionSpace > >
: public RestrictProlongInterfaceDefault< RestrictProlongTraits< RestrictProlongDefault< PetscDiscreteFunction< DiscreteFunctionSpace > >,
typename DiscreteFunctionSpace::DomainFieldType > >
{
typedef PetscDiscreteFunction< DiscreteFunctionSpace > DiscreteFunction;
typedef RestrictProlongDefault< DiscreteFunction > ThisType;
typedef RestrictProlongInterfaceDefault< RestrictProlongTraits< ThisType, typename DiscreteFunction::DomainFieldType > > BaseType;
typedef AdaptiveDiscreteFunction< DiscreteFunctionSpace > AdaptiveDiscreteFunctionType;
typedef RestrictProlongDefault < AdaptiveDiscreteFunctionType > AdaptiveRestrictProlongType;
public:
typedef DiscreteFunction DiscreteFunctionType;
typedef typename BaseType::DomainFieldType DomainFieldType;
typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
typedef typename DiscreteFunctionType::LocalFunctionType LocalFunctionType;
typedef typename DiscreteFunctionType::GridPartType GridPartType;
typedef DefaultLocalRestrictProlong< DiscreteFunctionSpaceType > LocalRestrictProlongType;
explicit RestrictProlongDefault ( DiscreteFunctionType &discreteFunction )
: discreteFunction_( discreteFunction ),
adaptiveFunction_( discreteFunction_.name()+"-adaptive", discreteFunction_.space() ),
rpOp_( adaptiveFunction_ ),
initialized_( false )
{
}
RestrictProlongDefault ( const RestrictProlongDefault& other )
: discreteFunction_( const_cast< DiscreteFunction& > (other.discreteFunction_) ),
adaptiveFunction_( other.adaptiveFunction_.name()+"-copy", discreteFunction_.space() ),
rpOp_( adaptiveFunction_ ),
initialized_( false )
{
}
void setFatherChildWeight ( const DomainFieldType &weight ) const
{
rpOp_.setFatherChildWeight( weight );
}
//! restrict data to father
template< class Entity >
void restrictLocal ( const Entity &father, const Entity &son, bool initialize ) const
{
assert( initialized_ );
rpOp_.restrictLocal( father, son, initialize );
}
//! prolong data to children
template< class Entity >
void prolongLocal ( const Entity &father, const Entity &son, bool initialize ) const
{
assert( initialized_ );
rpOp_.prolongLocal( father, son, initialize );
}
//! add discrete function to communicator with given unpack operation
template< class Communicator, class Operation >
void addToList ( Communicator &comm, const Operation& op)
{
rpOp_.addToList( comm, op );
}
//! add discrete function to communicator
template< class Communicator >
void addToList ( Communicator &comm )
{
rpOp_.addToList( comm );
}
//! remove discrete function from communicator
template< class Communicator >
void removeFromList ( Communicator &comm )
{
rpOp_.removeFromList( comm );
}
//! add discrete function to load balancer
template< class LoadBalancer >
void addToLoadBalancer ( LoadBalancer& lb )
{
rpOp_.addToLoadBalancer( lb );
}
void initialize ()
{
adaptiveFunction_.assign( discreteFunction_ );
rpOp_.initialize();
initialized_ = true ;
}
void finalize ()
{
// only finalize if previously initialized
assert( initialized_ );
rpOp_.finalize();
discreteFunction_.assign( adaptiveFunction_ );
initialized_ = false ;
}
protected:
using BaseType::calcWeight;
using BaseType::entitiesAreCopies;
protected:
DiscreteFunctionType& discreteFunction_;
AdaptiveDiscreteFunctionType adaptiveFunction_;
AdaptiveRestrictProlongType rpOp_;
bool initialized_;
};
} // namespace Fem
} // namespace Dune
......
......@@ -65,14 +65,36 @@ namespace Dune
const MapperType& mapper )
: BaseType( space.grid(), mapper, myArray_ ),
myArray_( space )
{}
{
}
void resize( const bool )
{ // do nothing here, only in compress
}
void reserve( int )
{
// do nothing here, only in compress
}
void dofCompress( const bool clearResizedArrays )
{
myArray_.resize();
if( clearResizedArrays )
{
myArray_.clear();
}
}
//! enable dof compression for dof storage (default is empty)
void enableDofCompression()
{
std::cerr << "WARNING: PetscVector cannot handle dof compression!" << std::endl;
}
protected:
DofArrayType myArray_;
};
// PetscVector
// -----------
......@@ -155,7 +177,7 @@ namespace Dune
std::size_t size () const { return mappers().ghostMapper().size(); }
void resize( const std::size_t newsize )
void resize( const std::size_t newsize = 0 )
{
// TODO: keep old data stored in current vector
// remove old vectors
......@@ -370,14 +392,17 @@ namespace Dune
{
Vec& vec = *getGhostedVector();
const PetscInt blocks = other.size();
const PetscInt bs = blockSize ;
const PetscInt b = 0;
const PetscScalar* vecData = static_cast< const PetscScalar* > (other.data());
::Dune::Petsc::VecSetBlockSize( vec, bs * blocks );
::Dune::Petsc::VecSetValuesBlocked( vec, 1, &b, vecData, INSERT_VALUES );
::Dune::Petsc::VecSetBlockSize( vec, blockSize );
// use mapping from the ghost mapper which removes duplicate dofs
const auto& mapping = mappers_.ghostMapper().mapping();
const size_t blocks = mapping.size();
assert( blocks == other.size() );
for( size_t b=0, bs = 0; b<blocks; ++b, bs += blockSize)
{
PetscInt block = mapping[ b ];
const PetscScalar* values = static_cast< const PetscScalar* > (other.data()+bs);
::Dune::Petsc::VecSetValuesBlocked( vec, 1, &block, values, INSERT_VALUES );
}
::Dune::Petsc::VecGhostGetLocalForm( vec_, &ghostedVec_ );
updateGhostRegions();
}
......@@ -389,12 +414,17 @@ namespace Dune
assert( DofBlock :: dimension == blockSize );
Vec& vec = *getGhostedVector();
const PetscInt blocks = other.size();
for( PetscInt b=0; b<blocks; ++b )
// use mapping from the ghost mapper which removes duplicate dofs
const auto& mapping = mappers_.ghostMapper().mapping();
const size_t blocks = mapping.size();
assert( blocks == other.size() );
for( size_t b=0; b<blocks; ++b )
{
PetscInt block = mapping[ b ];
const PetscScalar* values = static_cast< const PetscScalar* > (&(other[ b ][ 0 ])) ;
::Dune::Petsc::VecSetValuesBlocked( vec, 1, &b, values, INSERT_VALUES );
::Dune::Petsc::VecSetValuesBlocked( vec, 1, &block, values, INSERT_VALUES );
}
::Dune::Petsc::VecGhostGetLocalForm( vec_, &ghostedVec_ );
updateGhostRegions();
}
......@@ -403,9 +433,24 @@ namespace Dune
template <class Container>
void copyTo( SimpleBlockVector< Container, blockSize >& other ) const
{
PetscScalar *array = nullptr;
VecGetArray( ghostedVec_, &array );
std::copy_n( array, blockSize * other.size(), other.data() );
typedef typename Container::FieldType FieldType;
const PetscScalar *array = nullptr;
VecGetArrayRead( ghostedVec_, &array );
// use mapping from the ghost mapper which removes duplicate dofs
const auto& mapping = mappers_.ghostMapper().mapping();
const size_t blocks = mapping.size();
assert( blocks == other.size() );
for( size_t b=0; b<blocks; ++b )
{
const PetscScalar* petscBlock = array + (blockSize * mapping[ b ]);
FieldType* otherBlock = other.data() + (b * blockSize);
for( int i=0; i<blockSize; ++i )
{
otherBlock[ i ] = petscBlock[ i ];
}
}
VecRestoreArrayRead( ghostedVec_, &array );
}
// assign from other given ISTLBlockVector with same block size
......@@ -413,19 +458,23 @@ namespace Dune
void copyTo ( ISTLBlockVector< DofBlock >& other ) const
{
assert( DofBlock :: dimension == blockSize );
PetscScalar *array = nullptr;
VecGetArray( ghostedVec_, &array );
const PetscInt blocks = other.size();
for( PetscInt b=0, id = 0; b<blocks; ++b )
const PetscScalar *array = nullptr;
VecGetArrayRead( ghostedVec_, &array );
// use mapping from the ghost mapper which removes duplicate dofs
const auto& mapping = mappers_.ghostMapper().mapping();
const size_t blocks = mapping.size();
assert( blocks == other.size() );
for( size_t b=0; b<blocks; ++b )
{
auto& block = other[ b ];
for( int d=0; d<blockSize; ++d, ++id )
const PetscScalar* petscBlock = array + (blockSize * mapping[ b ]);
DofBlock& otherBlock = other[ b ];
for( int i=0; i<blockSize; ++i )
{
block[ d ] = array[ id ];
otherBlock[ i ] = petscBlock[ i ];
}
}
VecRestoreArrayRead( ghostedVec_, &array );
}
PetscVector& operator= ( const ThisType& other )
......
......@@ -252,7 +252,7 @@ namespace Dune
virtual int preconditionerIteration () const
{
return parameter_.getValue< int >( keyPrefix_ + "preconditioning.iterations", 0 );
return parameter_.getValue< int >( keyPrefix_ + "preconditioning.iterations", 1 );
}
virtual int preconditionerLevel () const
......
......@@ -434,6 +434,7 @@ int main(int argc, char** argv)
pass &= Algorithm< InverseOperator, LinearOperator >::apply( grid, designation, verboseSolver, &param);
}
/*
// PetscInverseOperator + PetscLinearOperator + AdaptiveDiscreteFunction
// this case is appearing for adaptive simulations with PETSc as solver backend
{
......@@ -459,6 +460,7 @@ int main(int argc, char** argv)
));
pass &= Algorithm< InverseOperator, LinearOperator >::apply( grid, designation, verboseSolver, &param);
}
*/
#endif // HAVE_PETSC
/*
......
......@@ -366,6 +366,9 @@ namespace Dune
template <PartitionIteratorType pitype>
void genericAdapt () const
{
// initialize restrict prolong operator (e.g. PetscRestrictProlong... )
rpOp_.initialize();
// call pre-adapt, returns true if at least
// one element is marked for coarsening
bool restr = grid_.preAdapt();
......@@ -380,6 +383,7 @@ namespace Dune
if(restr)
{
// get macro grid view
MacroGridView macroView = grid_.levelGridView( 0 );
......@@ -453,6 +457,9 @@ namespace Dune
// do cleanup
grid_.postAdapt();
// finalize restrict prolong operator (e.g. PetscRestrictProlong... )
rpOp_.finalize();
}
private:
......@@ -606,6 +613,8 @@ namespace Dune
typedef AdaptationManagerBase<GridType,RestProlOperatorImp> BaseType;
typedef LoadBalancer<GridType> Base2Type;
using BaseType :: rpOp_;
// reference counter to ensure only one instance per grid exists
ObjectType& referenceCounter_;
......@@ -682,8 +691,15 @@ namespace Dune
/** @copydoc LoadBalancerInterface::loadBalance */
virtual bool loadBalance ()
{
// same as for the adapt method
rpOp_.initialize () ;
// call load balance
return Base2Type :: loadBalance();
const bool result = Base2Type :: loadBalance( );
// finalize rp object (mostly RestrictProlongDefault for PetscDF)
rpOp_.finalize () ;
return result ;
}
/** @copydoc LoadBalancerInterface::loadBalanceTime */
......
......@@ -76,6 +76,9 @@ namespace Dune
// if preAdapt was already called just return
if( initializeCalled_ ) return ;
// initialize restrict-prolong object
rpOp_.initialize();
// unset was changed
wasChanged_ = false;
// reserve memory
......@@ -105,6 +108,9 @@ namespace Dune
wasChanged_ = false;
}
// initialize restrict-prolong object
rpOp_.finalize();
// set postAdaptCalled flag
finalizeCalled_ = true ;
......
......@@ -46,6 +46,20 @@ namespace Dune
//! \brief field type of domain vector space
typedef typename Traits::DomainFieldType DomainFieldType;
/** \brief initialize restrict prolong object (if necessary) before adaptation takes place
*/
void initialize ()
{
CHECK_AND_CALL_INTERFACE_IMPLEMENTATION( asImp().initialize() );
}
/** \brief finalize restrict prolong object (if necessary) after adaptation and dof compression was finished
*/
void finalize ()
{
CHECK_AND_CALL_INTERFACE_IMPLEMENTATION( asImp().finalize() );
}
/** \brief explicit set volume ratio of son and father
*
* \param[in] weight volume of son / volume of father
......@@ -157,9 +171,13 @@ namespace Dune
return (indexSet.index( father ) == indexSet.index( son ));
}
public:
/** \copydoc RestrictProlongInterface::setFatherChildWeight(const DomainFieldType &weight) const*/
void setFatherChildWeight ( const DomainFieldType &weight ) const {}
void initialize () {}
void finalize () {}
};
......
......@@ -52,6 +52,8 @@ namespace Dune
template< int i > struct ProlongLocal;
template< int i > struct RestrictLocal;
template< int i > struct SetFatherChildWeight;
template< int i > struct Initialize;
template< int i > struct Finalize;
public:
/** \copydoc Dune::Fem::RestrictProlongInterface::DomainFieldType */
......@@ -75,6 +77,18 @@ namespace Dune
* \{
*/
/** \copydoc Dune::Fem::RestrictProlongInterface::initialize */
void initialize ()
{
Dune::Fem::ForLoop< Initialize, 0, sizeof...( Tail ) >::apply( tuple_ );
}
/** \copydoc Dune::Fem::RestrictProlongInterface::setFatherChildWeight */
void finalize ()
{
Dune::Fem::ForLoop< Finalize, 0, sizeof...( Tail ) >::apply( tuple_ );
}
/** \copydoc Dune::Fem::RestrictProlongInterface::setFatherChildWeight */
void setFatherChildWeight ( const DomainFieldType &weight ) const
{
......@@ -240,6 +254,34 @@ namespace Dune
};
// RestrictProlongTuple< Head, Tail... >::Initialize
// -------------------------------------------------
template< class Head, class... Tail >
template< int i >
struct RestrictProlongTuple< Head, Tail... >::Initialize
{
static void apply ( std::tuple< Head, Tail... > &tuple )
{
std::get< i >( tuple ).initialize();
}
};
// RestrictProlongTuple< Head, Tail... >::Finalize
// -----------------------------------------------
template< class Head, class... Tail >
template< int i >
struct RestrictProlongTuple< Head, Tail... >::Finalize
{
static void apply ( std::tuple< Head, Tail... > &tuple )
{
std::get< i >( tuple ).finalize();
}
};
// RestrictProlongDefaultTuple
// ---------------------------
......
......@@ -47,7 +47,6 @@ foreach( test ${TESTS} )
LINK_LIBRARIES dunefem )
endforeach()
# dgcomm tests should run on GRIDDIM=3 meshes
dune_add_test( NAME dgcomm SOURCES dgcomm.cc
COMPILE_DEFINITIONS "GRIDDIM=3;WORLDDIM=3;${GRIDTYPE};DIMRANGE=3;POLORDER=2;WANT_GRAPE=0;WANT_CACHED_COMM_MANAGER=0"
......
......@@ -5,6 +5,7 @@
using namespace Dune;
#include <dune/fem/function/adaptivefunction.hh>
#include <dune/fem/function/petscdiscretefunction.hh>
#include <dune/fem/function/vectorfunction.hh>
#include <dune/fem/function/blockvectorfunction.hh>
#include <dune/fem/function/blockvectordiscretefunction.hh>
......@@ -82,7 +83,16 @@ typedef DiscontinuousGalerkinSpace< FunctionSpace < double , double, MyGridType:
typedef typename DiscreteFunctionSpaceType :: FunctionSpaceType FunctionSpaceType;
//! define the type of discrete function we are using , see
#if HAVE_PETSC
// PetscDiscreteFunction uses AdaptiveDiscreteFunction for dof prolongation and
// resttriction
typedef PetscDiscreteFunction< DiscreteFunctionSpaceType > DiscreteFunctionType;
#else
typedef AdaptiveDiscreteFunction< DiscreteFunctionSpaceType > DiscreteFunctionType;
#endif
//typedef Dune::Fem::ManagedDiscreteFunction< VectorDiscreteFunction< DiscreteFunctionSpaceType, std::vector< double > > > DiscreteFunctionType;
//typedef ISTLBlockVectorDiscreteFunction< DiscreteFunctionSpaceType > DiscreteFunctionType;
//typedef Dune::Fem::ReferenceBlockVector< FunctionSpaceType::RangeFieldType,
......@@ -149,7 +159,8 @@ void adapt( MyGridType &grid, DiscreteFunctionType &solution, int step )
for( const auto& entity : space)
grid.mark( mark, entity );
adop.adapt();
std::cout << message << std::endl;
if( Parameter::verbose() )
std::cout << message << std::endl;
}
if( Parameter :: verbose() )
......@@ -166,12 +177,14 @@ double algorithm ( MyGridType &grid, DiscreteFunctionType &solution, int step, i
DiscreteFunctionType tmp ( solution );
ExactSolution f;
auto gridFunc = gridFunctionAdapter(f, solution.space().gridPart(), 2);
{
ExactSolution f;
Dune::Fem::interpolate( f, solution );
Dune::Fem::interpolate( gridFunc, solution );
Dune :: Fem :: L2Norm< GridPartType > l2norm ( solution.space().gridPart(), 2*order+2 ) ;
double new_error = l2norm.distance( f ,solution );
std::cout << "before ref." << new_error << "\n\n";
if( Parameter::verbose() )
std::cout << "before ref." << new_error << "\n\n";
}
adapt(grid,solution,step);
......@@ -194,12 +207,11 @@ double algorithm ( MyGridType &grid, DiscreteFunctionType &solution, int step, i
}
#endif
ExactSolution f;
// calculation L2 error on refined grid
// pol ord for calculation the error chould by higher than
// pol for evaluation the basefunctions
Dune :: Fem :: L2Norm< GridPartType > l2norm ( solution.space().gridPart(), 2*order+2 ) ;
double error = l2norm.distance( f, solution );
double error = l2norm.distance( gridFunc, solution );
#if USE_GRAPE
// if Grape was found, then display last solution
......@@ -212,9 +224,12 @@ double algorithm ( MyGridType &grid, DiscreteFunctionType &solution, int step, i
#endif
//! perform l2-projection to refined grid
Dune::Fem::interpolate( f, solution );
double new_error = l2norm.distance( f, solution );
std::cout << "\nL2 Error : " << error << " on new grid " << new_error << "\n\n";
Dune::Fem::interpolate( gridFunc, solution );
double new_error = l2norm.distance( gridFunc, solution );
if( Parameter::verbose() )
{
std::cout << "\nL2 Error : " << error << " on new grid " << new_error << "\n\n";
}
#if USE_GRAPE
// if Grape was found, then display last solution
if(turn > 0)
......@@ -271,7 +286,8 @@ try {
DiscreteFunctionType solution ( "sol", space );
solution.clear();
std::cout << "------------ Refining:" << std::endl;
if( Parameter::verbose() )
std::cout << "------------ Refining:" << std::endl;
for(int i=0; i<ml; i+=1)
{
error[i] = algorithm ( grid , solution, step, (i==ml-1));
......@@ -280,7 +296,8 @@ try {
if ( isLocallyAdaptive )
{
double eoc = log( error[i-1]/error[i]) / M_LN2;
std::cout << "EOC = " << eoc << std::endl;
if( Parameter::verbose() )
std::cout << "EOC = " << eoc << std::endl;
if( std::abs( eoc - (space.order()+eocThreshold) ) < 0 )
{
DUNE_THROW(InvalidStateException,"EOC check of refinement failed");
......@@ -290,7 +307,8 @@ try {
std::cout << "no EOC for non-adaptive grid" << std::endl;
}
}
std::cout << "------------ Coarsening:" << std::endl;
if( Parameter::verbose() )
std::cout << "------------ Coarsening:" << std::endl;
for(int i=ml-1; i>=0; i-=1)
{
error[i] = algorithm ( grid , solution,-step, 1);
......@@ -299,7 +317,8 @@ try {
if( isLocallyAdaptive )
{
double eoc = log( error[i+1]/error[i]) / M_LN2;
std::cout << "EOC = " << eoc << std::endl;
if( Parameter::verbose() )
std::cout << "EOC = " << eoc << std::endl;
if( std::abs( eoc + (space.order()+eocThreshold) ) < 0 )
{
DUNE_THROW(InvalidStateException,"EOC check of coarsening failed");
......
......@@ -25,6 +25,7 @@ const int polOrder = POLORDER;
#include <dune/fem/space/lagrange.hh>
// #include <dune/fem/space/padaptivespace.hh>
#include <dune/fem/function/adaptivefunction.hh>
#include <dune/fem/function/petscdiscretefunction.hh>
#include <dune/fem/quadrature/cachingquadrature.hh>
#include <dune/fem/space/common/interpolate.hh>
#include <dune/fem/misc/l2norm.hh>
......@@ -178,7 +179,11 @@ typedef LagrangeDiscreteFunctionSpace< FunctionSpaceType, GridPartType, polOrder
// DiscreteFunctionSpaceType;
//! type of the discrete function we are using
#if HAVE_PETSC
typedef AdaptiveDiscreteFunction< DiscreteFunctionSpaceType > DiscreteFunctionType;
#else
typedef AdaptiveDiscreteFunction< DiscreteFunctionSpaceType > DiscreteFunctionType;
#endif
typedef ExactSolution< FunctionSpaceType > ExactSolutionType;
......@@ -256,9 +261,13 @@ void algorithm ( GridPartType &gridPart,
double preL2error = l2norm.distance( f, solution );
double preH1error = h1norm.distance( f, solution );
std::cout << "Unknowns before adaptation: " << solution.space().size() << std::endl;
std::cout << "L2 error before adaptation: " << preL2error << std::endl;
std::cout << "H1 error before adaptation: " << preH1error << std::endl;
size_t unknowns = gridPart.comm().sum( solution.space().size() );
if( Parameter::verbose() )
{
std::cout << "Unknowns before adaptation: " << unknowns << std::endl;
std::cout << "L2 error before adaptation: " << preL2error << std::endl;
std::cout << "H1 error before adaptation: " << preH1error << std::endl;
}
adapt( gridPart.grid(), solution, step, locallyAdaptive );
......@@ -269,16 +278,20 @@ void algorithm ( GridPartType &gridPart,
double postL2error = l2norm.distance( solution, fexact );
double postH1error = h1norm.distance( f, solution );
std::cout << "Unknowns after "
<< (step < 0 ? "restriction" : "prolongation")
<< ": " << solution.space().size() << std::endl;
std::cout << "L2 error after "
<< (step < 0 ? "restriction" : "prolongation")
<< ": " << postL2error << std::endl;
std::cout << "H1 error after "
unknowns = gridPart.comm().sum( solution.space().size() );
if( Parameter::verbose() )
{
std::cout << "Unknowns after "
<< (step < 0 ? "restriction" : "prolongation")
<< ": " << postH1error << std::endl;
<< ": " << unknowns << std::endl;
std::cout << "L2 error after "
<< (step < 0 ? "restriction" : "prolongation")
<< ": " << postL2error << std::endl;
std::cout << "H1 error after "
<< (step < 0 ? "restriction" : "prolongation")
<< ": " << postH1error << std::endl;
}
#if USE_GRAPE && SHOW_RESTRICT_PROLONG
if( turn > 0 ) {
GrapeDataDisplay< MyGridType > grape( gridPart.grid() );
......@@ -286,12 +299,15 @@ void algorithm ( GridPartType &gridPart,
}
#endif
interpolate( f, solution );
interpolate( f, solution );
double newL2error = l2norm.distance( f, solution );
double newH1error = h1norm.distance( f, solution );
std :: cout << "L2 error for interpolation after adaption: " << newL2error << std :: endl;
std :: cout << "H1 error for interpolation after adaption: " << newH1error << std :: endl;
if( Parameter::verbose() )
{
std :: cout << "L2 error for interpolation after adaption: " << newL2error << std :: endl;
std :: cout << "H1 error for interpolation after adaption: " << newH1error << std :: endl;
}
#if USE_GRAPE && SHOW_INTERPOLATION
if( turn > 0 ) {
......@@ -303,8 +319,11 @@ void algorithm ( GridPartType &gridPart,
double l2eoc = -log( newL2error / preL2error) / M_LN2;
double h1eoc = -log( newH1error / preH1error) / M_LN2;
std :: cout << "L2 EOC: " << l2eoc << std :: endl;
std :: cout << "H1 EOC: " << h1eoc << std :: endl;
if( Parameter::verbose() )
{
std :: cout << "L2 EOC: " << l2eoc << std :: endl;
std :: cout << "H1 EOC: " << h1eoc << std :: endl;
}
const bool isLocallyAdaptive = Dune::Fem::Capabilities::isLocallyAdaptive< GridPartType :: GridType > :: v ;
// threshold for EOC difference to predicted value
......@@ -329,7 +348,9 @@ void algorithm ( GridPartType &gridPart,
dataio.writeData( solution, xdr, "sol", turn );
#endif
std :: cout << std :: endl;