Commit a5fe4687 authored by Andreas Dedner's avatar Andreas Dedner

added hessian method to LFE shape function set

bug fix

bug fix

fix the lfe space restrict prolong

need Traits class in LFConverter so that it can be correctly passed into
dune-localfunction's interpolate methods

final clean up
parent 2460b27f
......@@ -99,6 +99,16 @@ namespace Dune
static const int dimDomain = FunctionSpaceType::dimDomain;
struct Traits
{
typedef typename FunctionSpaceType::DomainType DomainType;
typedef typename FunctionSpaceType::RangeType RangeType;
typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;
typedef typename FunctionSpaceType::DomainFieldType DomainFieldType;
typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
};
LocalFunctionConverter ( const HostLocalFunction &hostLocalFunction, const Converter &converter = Converter() )
: BaseType( hostLocalFunction ), converter_( converter )
{}
......
......@@ -281,27 +281,6 @@ namespace Dune
}
};
///////////////////////////////////////////////////////
template< class Grid >
struct BoundaryIdGetter
{
typedef std::function<int(const typename Grid::template Codim<0>::Geometry::GlobalCoordinate&)> Caller;
BoundaryIdGetter(const Caller &caller)
: caller_(caller) {}
template< class Intersection >
int boundaryId ( const Intersection &intersection ) const
{
if (!intersection.boundary) return 0;
const auto x = intersection.geometry().center();
int val = caller(x);
return (val==0)?BoundaryIdProvider<Grid>::boundaryId(intersection):val;
}
private:
Caller caller_;
};
} // namespace Fem
} // namespace Dune
......
......@@ -72,7 +72,7 @@ namespace Dune
template< class LFFather >
void restrictFinalize ( LFFather &lfFather ) const
{
Fem::ForLoop< RestrictFinalize, 0, setSize >::apply( lfFather, lfSon, localRestrictProlongTuple_ );
Fem::ForLoop< RestrictFinalize, 0, setSize >::apply( lfFather, localRestrictProlongTuple_ );
}
template< class LFFather, class LFSon, class LocalGeometry >
......
......@@ -271,7 +271,7 @@ namespace Dune
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 );
LocalFunctionType lfFather = discreteFunction_.localFunction( gpFather );
localRP_.restrictFinalize( lfFather );
}
......
......@@ -203,18 +203,19 @@ namespace Dune
dofAlignment_( basisFunctionSet_ )
{}
template< class LocalFunction, class Vector >
template< class LocalFunction, class Vector>
void operator() ( const LocalFunction &localFunction, Vector &localDofVector ) const
{
typedef Vector LocalDofVector;
// clear dofs before something is adedd
// 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_ );
localInterpolation().interpolate(
localFunctionConverter( localFunction, RangeConverter( i ) ), subLdv );
SubDofVectorWrapper< Vector, DofAlignmentType > subLdv( localDofVector, i, dofAlignment_ );
(*this)(
localFunctionConverter( localFunction, RangeConverter(i) ),
subLdv, PriorityTag<42>()
);
}
}
......@@ -228,6 +229,22 @@ namespace Dune
const LocalInterpolationType &localInterpolation () const { return localInterpolation_; }
private:
template< class LocalFunction, class Vector>
auto operator() ( const LocalFunction &localFunction, Vector &localDofVector, PriorityTag<1> ) const
-> void_t< decltype(
std::declval<LocalInterpolationType>().interpolate(
localFunction, localDofVector)) >
{
localInterpolation().interpolate( localFunction, localDofVector );
}
template< class LocalFunction, class Vector>
void operator() ( const LocalFunction &localFunction, Vector &localDofVector, PriorityTag<0> ) const
{
std::vector<typename Vector::value_type> tmp(basisFunctionSet_.size());
localInterpolation().interpolate( localFunction, tmp);
for (unsigned int i=0;i<tmp.size();++i)
localDofVector[i] = tmp[i];
}
BasisFunctionSetType basisFunctionSet_;
LocalInterpolationType localInterpolation_;
DofAlignmentType dofAlignment_;
......
......@@ -64,6 +64,9 @@ namespace Dune
{
values_.reserve( size() );
jacobians_.reserve( size() );
hessians_.resize( DomainType::dimension*(DomainType::dimension+1)/2. );
for (unsigned int i=0;i<hessians_.size();++i)
hessians_[i].reserve( size() );
}
int order () const { return localBasis_.order(); }
......@@ -89,7 +92,22 @@ namespace Dune
template< class Point, class Functor >
void hessianEach ( const Point &x, Functor f ) const
{
DUNE_THROW( NotImplemented, "Method hessianEach not implemented" );
std::array<unsigned int, DomainType::dimension> multiIndex;
std::fill(multiIndex.begin(),multiIndex.end(),0);
unsigned int k = 0;
for (unsigned int i=0;i<DomainType::dimension;++i)
{
multiIndex[i] = 1;
for (unsigned int j=i;j<DomainType::dimension;++j)
{
multiIndex[j] += 1;
localBasis_.partial(multiIndex,coordinate(x),hessians_[k]);
multiIndex[j] -= 1;
++k;
}
multiIndex[i] -= 1;
}
callFunctor( hessians_, f );
}
private:
......@@ -101,10 +119,33 @@ namespace Dune
for( Iterator it = v.begin(); it != v.end(); ++it )
f( i++, *it );
}
template< class T, class Functor >
static void callFunctor ( const std::vector< std::vector<T> > &v, Functor f )
{
HessianRangeType h;
for (unsigned int b=0;b<v[0].size();++b)
{
unsigned int k = 0;
for (unsigned int i=0;i<DomainType::dimension;++i)
{
for (unsigned int j=i;j<DomainType::dimension;++j)
{
for (unsigned int r=0;r<RangeType::dimension;++r)
{
h[r][i][j] = v[b][k][r];
h[r][j][i] = v[b][k][r];
}
++k;
}
}
f( b, h );
}
}
const LocalBasis& localBasis_;
mutable std::vector< RangeType > values_;
mutable std::vector< JacobianRangeType > jacobians_;
mutable std::vector< std::vector< RangeType > > hessians_;
};
} // namespace Fem
......
#ifndef DUNE_FEM_SPACE_RAVIARTTHOMAS_LOCALINTERPOLATION_HH
#define DUNE_FEM_SPACE_RAVIARTTHOMAS_LOCALINTERPOLATION_HH
#error SHOULDN't be needed
// alternative interpolation used for testing
// C++ includes
#include <cassert>
......
if( NOT GRIDTYPE )
set( GRIDTYPE ALUGRID_CONFORM )
set( GRIDTYPE YASPGRID )
endif()
if( NOT GRIDDIM )
......@@ -58,7 +58,6 @@ LINK_LIBRARIES dunefem MPI_RANKS 1 2 3 4 TIMEOUT 9999999 )
dune_add_test( SOURCES test-slavedofs.cc LINK_LIBRARIES dunefem MPI_RANKS 1 2 4 TIMEOUT 300 COMPILE_DEFINITIONS "${DEFAULTFLAGS};USE_COMBINED_SPACE" )
dune_add_test( SOURCES test-raviartthomasinterpolation.cc CMAKE_GUARD dune-localfunctions_FOUND LINK_LIBRARIES dunefem COMPILE_DEFINITIONS "${DEFAULTFLAGS}" )
dune_add_test( SOURCES test-spaceinterpolation.cc LINK_LIBRARIES dunefem COMPILE_DEFINITIONS "${DEFAULTFLAGS}" )
if( ${FEM_TORTURE_TESTSS} )
# test adaptation of combinedspaces
......
#include <config.h>
#ifdef YASPGRID
#undef YASPGRID
#endif
#define ALUGRID_CONFORM
#define SHOW_INTERPOLATION 0
#define SHOW_RESTRICT_PROLONG 1
#include <config.h>
// to write out the data, set WRITE_DATA to 1
#define WRITE_DATA 0
// to use grape, set to WANT_GRAPE to 1
#ifndef WANT_GRAPE
#define WANT_GRAPE 0
#endif
#if WANT_GRAPE
#define USE_GRAPE 1
#endif
#define USE_LFE 0
#define USE_PSPACE 0
// polynomial order of base functions
const int polOrder = 1; // POLORDER;
const int polOrder = 2; // POLORDER;
#include <iostream>
#include <sstream>
......@@ -27,7 +24,7 @@ const int polOrder = 1; // POLORDER;
#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
#include <dune/fem/space/common/adaptationmanager.hh>
#include <dune/fem/space/lagrange.hh>
// #include <dune/fem/space/padaptivespace.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>
......@@ -36,9 +33,6 @@ const int polOrder = 1; // POLORDER;
#include <dune/fem/misc/h1norm.hh>
#include <dune/fem/io/file/dataoutput.hh>
#if USE_GRAPE
#include <dune/grid/io/visual/grapedatadisplay.hh>
#endif
#include <dune/fem/io/parameter.hh>
#include <dune/fem/test/testgrid.hh>
......@@ -131,8 +125,10 @@ public:
{
phi = 1;
for( int i = 0; i < DomainType :: dimension; ++i )
// phi[ 0 ] += x[ i ] * x[ i ];
{
phi[ 0 ] *= sin( M_PI * x[ i ] );
phi[ 1 ] += x[ i ] * x[ i ];
}
}
void evaluate ( const DomainType &x, RangeFieldType t, RangeType &phi ) const
......@@ -145,8 +141,10 @@ public:
Dphi = 1;
for( int i = 0; i < DomainType :: dimension; ++i )
for( int j = 0; j < DomainType :: dimension; ++j )
// Dphi[ 0 ][ j ] *= ((i != j) ? 1. : 2.*x[i]);
{
Dphi[ 0 ][ j ] *= ((i != j) ? sin( M_PI * x[ i ]) : M_PI * cos( M_PI * x[ i ] ));
Dphi[ 1 ][ j ] *= ((i != j) ? 1. : 2.*x[i]);
}
}
void jacobian( const DomainType &x, RangeFieldType t, JacobianRangeType &Dphi ) const
......@@ -165,15 +163,19 @@ typedef Dune::GridSelector::GridType MyGridType;
typedef CheckGridEnabled< MyGridType >::GridPartType GridPartType;
//! type of the function space
typedef FunctionSpace< double, double, MyGridType::dimensionworld, 1 > FunctionSpaceType;
typedef FunctionSpace< double, double, MyGridType::dimensionworld, 2 > FunctionSpaceType;
//! type of the discrete function space our unkown belongs to
#if USE_LFE
typedef LagrangeSpace< FunctionSpaceType, GridPartType >
// typedef LagrangeDiscreteFunctionSpace< FunctionSpaceType, GridPartType, polOrder >
DiscreteFunctionSpaceType;
//! type of the discrete function space our unkown belongs to
//typedef Fem :: PAdaptiveLagrangeSpace< FunctionSpaceType, GridPartType, polOrder >
// DiscreteFunctionSpaceType;
#elif USE_PSPACE
typedef Fem :: PAdaptiveLagrangeSpace< FunctionSpaceType, GridPartType, polOrder >
DiscreteFunctionSpaceType;
#else
typedef LagrangeDiscreteFunctionSpace< FunctionSpaceType, GridPartType, polOrder >
DiscreteFunctionSpaceType;
#endif
//! type of the discrete function we are using
#if HAVE_PETSC
......@@ -288,13 +290,7 @@ void algorithm ( GridPartType &gridPart,
<< ": " << postH1error << std::endl;
}
#if USE_GRAPE && SHOW_RESTRICT_PROLONG
if( turn > 0 ) {
GrapeDataDisplay< MyGridType > grape( gridPart.grid() );
grape.dataDisplay( solution );
}
#endif
#if 1
#if WRITE_DATA
{
static int turn = 0;
typedef std::tuple< DiscreteFunctionType* > IODataType;
......@@ -315,13 +311,6 @@ void algorithm ( GridPartType &gridPart,
std :: cout << "H1 error for interpolation after adaption: " << newH1error << std :: endl;
}
#if USE_GRAPE && SHOW_INTERPOLATION
if( turn > 0 ) {
GrapeDataDisplay< MyGridType > grape( gridPart.grid );
grape.dataDisplay( solution );
}
#endif
double l2eoc = -log( newL2error / preL2error) / M_LN2;
double h1eoc = -log( newH1error / preH1error) / M_LN2;
......@@ -335,7 +324,7 @@ void algorithm ( GridPartType &gridPart,
// threshold for EOC difference to predicted value
const double eocThreshold = Parameter :: getValue("adapt.eocthreshold", double(0.25) );
if( 0 && isLocallyAdaptive )
if( isLocallyAdaptive )
{
const double sign = step / std::abs( step );
if( std::abs( l2eoc - h1eoc - sign ) > eocThreshold )
......@@ -348,12 +337,6 @@ void algorithm ( GridPartType &gridPart,
DUNE_THROW( InvalidStateException,"Wrong H1-EOC for " << refcrs );
}
#if WRITE_DATA
GrapeDataIO< MyGridType > dataio;
dataio.writeGrid( gridPart.grid(), xdr, "gridout", 0, turn );
dataio.writeData( solution, xdr, "sol", turn );
#endif
if( Parameter::verbose() )
std :: cout << std :: endl;
......@@ -387,7 +370,11 @@ try
const int step = Dune::Fem::TestGrid::refineStepsForHalf();
GridPartType gridPart( grid );
#if USE_LFE
DiscreteFunctionSpaceType discreteFunctionSpace( gridPart, polOrder );
#else
DiscreteFunctionSpaceType discreteFunctionSpace( gridPart );
#endif
DiscreteFunctionType solution( "solution", discreteFunctionSpace );
solution.clear();
......@@ -400,7 +387,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;
......
......@@ -92,24 +92,24 @@ static const int dimRange = GridPartType::dimensionworld;
typedef Dune::Fem::GridFunctionSpace< GridPartType, Dune::FieldVector< typename GridPartType::ctype, dimRange > > FunctionSpaceType;
typedef std::tuple<
// Dune::Fem::FiniteVolumeSpace< FunctionSpaceType, GridPartType >,
// Dune::Fem::DiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 0 >,
// Dune::Fem::DiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 1 >,
// Dune::Fem::DiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 2 >,
// Dune::Fem::LagrangeDiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 1 >,
// Dune::Fem::LagrangeDiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 2 >,
// Dune::Fem::LegendreDiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 1 >,
// Dune::Fem::LegendreDiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 2 >,
Dune::Fem::FiniteVolumeSpace< FunctionSpaceType, GridPartType >,
Dune::Fem::DiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 0 >,
Dune::Fem::DiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 1 >,
Dune::Fem::DiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 2 >,
Dune::Fem::LagrangeDiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 1 >,
Dune::Fem::LagrangeDiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 2 >,
Dune::Fem::LegendreDiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 1 >,
Dune::Fem::LegendreDiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 2 >,
#if HAVE_DUNE_LOCALFUNCTIONS
// Dune::Fem::BrezziDouglasMariniSpace< FunctionSpaceType, GridPartType, 1 >,
// Dune::Fem::BrezziDouglasMariniSpace< FunctionSpaceType, GridPartType, 2 >,
// Dune::Fem::RaviartThomasSpace< FunctionSpaceType, GridPartType, 0 >,
// Dune::Fem::RaviartThomasSpace< FunctionSpaceType, GridPartType, 1 >,
Dune::Fem::LagrangeSpace< FunctionSpaceType, GridPartType >
// Dune::Fem::RannacherTurekSpace< FunctionSpaceType, GridPartType >,
Dune::Fem::BrezziDouglasMariniSpace< FunctionSpaceType, GridPartType, 1 >,
Dune::Fem::BrezziDouglasMariniSpace< FunctionSpaceType, GridPartType, 2 >,
Dune::Fem::RaviartThomasSpace< FunctionSpaceType, GridPartType, 0 >,
Dune::Fem::RaviartThomasSpace< FunctionSpaceType, GridPartType, 1 >,
Dune::Fem::LagrangeSpace< FunctionSpaceType, GridPartType >,
Dune::Fem::RannacherTurekSpace< FunctionSpaceType, GridPartType >,
#endif // #if HAVE_DUNE_LOCALFUNCTIONS
// Dune::Fem::LagrangeDiscreteFunctionSpace< FunctionSpaceType, GridPartType, 1 >
// Dune::Fem::LagrangeDiscreteFunctionSpace< FunctionSpaceType, GridPartType, 2 >
Dune::Fem::LagrangeDiscreteFunctionSpace< FunctionSpaceType, GridPartType, 1 >,
Dune::Fem::LagrangeDiscreteFunctionSpace< FunctionSpaceType, GridPartType, 2 >
> DiscreteFunctionSpacesType;
typedef ErrorTuple< DiscreteFunctionSpacesType >::Type ErrorTupleType;
......@@ -121,7 +121,7 @@ template< class DiscreteFunctionSpace >
std::pair< Real< DiscreteFunctionSpace >, Real< DiscreteFunctionSpace > >
algorithm ( typename DiscreteFunctionSpace::GridPartType &gridPart )
{
DiscreteFunctionSpace space( gridPart,2 );
DiscreteFunctionSpace space( gridPart );
Dune::Fem::AdaptiveDiscreteFunction< DiscreteFunctionSpace > u( "solution", space );
// interpolate a function
......@@ -129,7 +129,7 @@ algorithm ( typename DiscreteFunctionSpace::GridPartType &gridPart )
const auto uGridExact = gridFunctionAdapter( "exact solution", uExact, gridPart, 3 );
interpolate( uGridExact, u );
#if 1
#if 0
{
static int turn = 0;
typedef std::tuple< decltype(u)* > IODataType;
......
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