[dune-pdelab] some problems in dune
Robert Jendrny
Robert.Jendrny at math.tu-dortmund.de
Thu Oct 30 15:14:37 CET 2014
Hi Rebecca,
thank you for your answer ...
I am working with the code from pdelab exercise2, example02 (Dune course
this year - 2014)
I hope the attached files are enough
The most important file is example02_Q1.hh, there everything is calculated
I want to solve Poisson (-Laplace(u) = 1) on unit square with AMG
For Q1 I nearly get independent numbers of iterations for mesh size h,
h/2, h/4,.... which is fine
But If I change Q1 to Q2, the numbers of iterations "explode" if the
mesh gets finer and finer
Maybe someone can fix it ;-)
Bye,
Robert
On 10/30/2014 03:00 PM, Rebecca de Cuveland wrote:
> Hey Robert,
> perhaps you can show us your program and output? (verbosity 3 or 4 for AMG and 1 for linear solver)
> AMG is quite sensible for the used options. I only use Q1, but as far as I know Q2 should also work...
> Bye,
> Rebecca
>
--
Robert Jendrny
B.Sc. Technomathematik
TU Dortmund, Fakultät für Mathematik
LS III: Angewandte Mathematik und Numerik
M 1020, Vogelpothsweg 87, 44227 Dortmund
robert.jendrny at math.tu-dortmund.de
-------------- next part --------------
template<typename GV, typename RF>
class ExactSolution
: public Dune::PDELab::GridFunctionBase<
Dune::PDELab::GridFunctionTraits<GV,RF,1,Dune::FieldVector<RF,1> >
,ExactSolution<GV,RF> >
{
public:
typedef Dune::PDELab::GridFunctionTraits<GV,RF,1,Dune::FieldVector<RF,1> > Traits;
//! constructor
ExactSolution
(const GV& g_) : g(g_) {}
//! \copydoc GridFunctionBase::evaluate()
inline void evaluate (const typename Traits::ElementType& e,
const typename Traits::DomainType& x,
typename Traits::RangeType& y) const
{
typename Traits::DomainType xglobal = e.geometry().global(x);
// Für Einheitsquadrat
y = 16.0 * sin(M_PI * xglobal[0]) * sin(M_PI * xglobal[1]);
}
inline const typename Traits::GridViewType& getGridView () const
{
return g;
}
private:
const GV & g;
};
/*! \brief Adapter returning ||f1(x)-f2(x)||^2 for two given grid functions
\tparam T1 a grid function type
\tparam T2 a grid function type
*/
template<typename T1, typename T2>
class DifferenceSquaredAdapter
: public Dune::PDELab::GridFunctionBase<
Dune::PDELab::GridFunctionTraits<typename T1::Traits::GridViewType,
typename T1::Traits::RangeFieldType,
1,Dune::FieldVector<typename T1::Traits::RangeFieldType,1> >
,DifferenceSquaredAdapter<T1,T2> >
{
public:
typedef Dune::PDELab::GridFunctionTraits<typename T1::Traits::GridViewType,
typename T1::Traits::RangeFieldType,
1,Dune::FieldVector<typename T1::Traits::RangeFieldType,1> > Traits;
//! constructor
DifferenceSquaredAdapter (const T1& t1_, const T2& t2_) : t1(t1_), t2(t2_) {}
//! \copydoc GridFunctionBase::evaluate()
inline void evaluate (const typename Traits::ElementType& e,
const typename Traits::DomainType& x,
typename Traits::RangeType& y) const
{
typename Traits::RangeType v1, v2;
t1.evaluate( e, x, v1);
t2.evaluate( e, x, v2);
v1 -= v2; //v1= v1-v2
y = v1*v1;// y = (v1-v2)^2
}
inline const typename Traits::GridViewType& getGridView () const
{
return t1.getGridView();
}
private:
const T1& t1;
const T2& t2;
};
/*! \brief Adapter returning ||f(x)||^2 for the given grid functions
\tparam T a grid function type
*/
template<typename T>
class SquaredAdapter
: public Dune::PDELab::GridFunctionBase<
Dune::PDELab::GridFunctionTraits<typename T::Traits::GridViewType,
typename T::Traits::RangeFieldType,
1,Dune::FieldVector<typename T::Traits::RangeFieldType,1> >
,SquaredAdapter<T> >
{
public:
typedef Dune::PDELab::GridFunctionTraits<typename T::Traits::GridViewType,
typename T::Traits::RangeFieldType,
1,Dune::FieldVector<typename T::Traits::RangeFieldType,1> > Traits;
//! constructor
SquaredAdapter (const T& t_) : t(t_) {}
//! \copydoc GridFunctionBase::evaluate()
inline void evaluate (const typename Traits::ElementType& e,
const typename Traits::DomainType& x,
typename Traits::RangeType& y) const
{
typename T::Traits::RangeType y1;
t.evaluate(e,x,y1);
y = y1.two_norm2();
}
inline const typename Traits::GridViewType& getGridView () const
{
return t.getGridView();
}
private:
const T& t;
};
/*! A function that computes the L2 norm of the given function object
* be careful, this l2norm works only for DGF,
* that means function f should be evaluated in
* !!local coordinates!!
*/
template<class GV,class F>
double l2norm(const GV& gv, const F& f, const unsigned qorder = 3)
{
// first we extract the dimensions of the grid
const int dim = GV::Grid::dimension;
// type used for coordinates in the grid
typedef typename GV::Grid::ctype ctype;
// Extract types
typedef typename GV::Traits::template Codim<0>::Iterator ElementIterator;
typedef typename GV::Traits::template Codim<0>::Entity::Geometry Geometry;
// Initialize norm
double norm = 0.0;
typedef typename F::Traits::RangeType RF;
RF val;
// iterate through all entities of codim 0 of the grid view
for (ElementIterator it = gv.template begin<0>();
it!=gv.template end<0>(); ++it)
{
const Geometry& geo = it->geometry();
Dune::GeometryType gt = geo.type();
// get quadrature rule of order 2
const Dune::QuadratureRule<ctype,dim>&
rule = Dune::QuadratureRules<ctype,dim>::rule(gt, qorder);
// compute approximate integral
for (typename Dune::QuadratureRule<ctype,dim>::const_iterator git=rule.begin();
git!=rule.end(); ++git)
{
// evaluate the given grid functions at integration point
f.evaluate(*it,git->position(),val);
// compute square value of the function
//val*=val; <-- BUG 4 life!
// quadrature element weight
RF weight = git->weight();
// integration element volume
RF detjac = geo.integrationElement(git->position());
// accumulate in a norm
norm += val*weight*detjac;
}
}
// return sqrt of the norm
return std::sqrt(norm);
}
/*! A function that computes the L2 norm of the given function object
* be careful, this l2norm works only function object,
* which evaluates the function in !!global coordinates!!
*/
template<class GV,class F>
double l2norm_global(const GV& gv, const F& f, const unsigned qorder = 3)
{
// first we extract the dimensions of the grid
const int dim = GV::Grid::dimension;
// type used for coordinates in the grid
typedef typename GV::Grid::ctype ctype;
// Extract types
typedef typename GV::Traits::template Codim<0>::Iterator ElementIterator;
typedef typename GV::Traits::template Codim<0>::Entity::Geometry Geometry;
// Initialize norm
double norm = 0.0;
// iterate through all entities of codim 0 of the grid view
for (ElementIterator it = gv.template begin<0>();
it!=gv.template end<0>(); ++it)
{
const Geometry& geo = it->geometry();
Dune::GeometryType gt = geo.type();
// get quadrature rule of order 2
const Dune::QuadratureRule<ctype,dim>&
rule = Dune::QuadratureRules<ctype,dim>::rule(gt, qorder);
// compute approximate integral
for (typename Dune::QuadratureRule<ctype,dim>::const_iterator git=rule.begin();
git!=rule.end(); ++git)
{
const double value = f(geo.global(git->position()));
norm += git->weight()*value*value;
}
}
// return sqrt of the norm
return std::sqrt(norm);
}
template<typename GV, typename RF>
class BETA
: public Dune::PDELab::GridFunctionBase<
Dune::PDELab::GridFunctionTraits<GV,RF,GV::dimension,
Dune::FieldVector<RF,GV::dimension> >,
BETA<GV,RF> >
{
const GV& gv;
public:
typedef Dune::PDELab::GridFunctionTraits<GV,RF,
GV::dimension, Dune::FieldVector<RF,GV::dimension> > Traits;
typedef Dune::PDELab::GridFunctionBase<
Dune::PDELab::GridFunctionTraits<GV, RF,GV::dimension,
Dune::FieldVector<RF,GV::dimension> >,
BETA<GV,RF> > BaseT;
BETA (const GV& gv_ ) : gv(gv_){}
inline void evaluate (const typename Traits::ElementType& e,
const typename Traits::DomainType& x,
typename Traits::RangeType& y) const
{
y = 0;
y[0] = 0.0;
y[1] = 0.0;
// An example for another velocity field maybe...:
//typename Traits::DomainType xg = e.geometry().global(x);
//y[0] = 1.0;
//y[1] = 0.5 * cos(10*xg[0]) ;
}
inline const typename Traits::GridViewType& getGridView () { return gv; }
typename Traits::RangeFieldType
f (const typename Traits::ElementType& e, const typename Traits::DomainType& x, typename Traits::RangeType v) const
{
// (?,?, a Diffusionparameter, v Geschwindigkeitsvektor, c Massparamater)
//xglobal sind Koords. x= xglobal[0], y= xglobal[1]
typename Traits::DomainType xglobal = e.geometry().global(x);
typename Traits::RangeFieldType ret;
ret = //exakt for normal quad
// 32.0 * M_PI * M_PI * sin(M_PI * xglobal[0]) * sin(M_PI * xglobal[1]);
1.0;
return ret;
}
};
/*
##################################################################################################################
##################################################################################################################
LOOP SOLVER AMG SOR
##################################################################################################################
##################################################################################################################
*/
template<typename Real, typename U, typename GV, typename GFS,
typename CC, typename BCTypeParam, typename BETA>
void sdsolver_LS_AMG( U& u, const GV& gv, const GFS& gfs, const CC& cc,
const BCTypeParam& bctype, const BETA& beta, int level, bool damping, const std::string outputname)
{
// <<<3>>> Make grid operator space
typedef Example02LocalOperator<BCTypeParam,BETA> LOP; // operator including boundary
LOP lop( bctype, beta);
typedef Dune::PDELab::ISTLMatrixBackend MBE;
typedef Dune::PDELab::GridOperator<
GFS,GFS, /* ansatz and test space */
LOP, /* local operator */
MBE, /* matrix backend */
Real,Real,Real, /* field types for domain, range and jacobian */
CC,CC /* constraints transformation for ansatz and test space */
> GO;
GO go(gfs,cc,gfs,cc,lop);
//std::cout << "<<<3>>> done" << std::endl;
// <<<5>>> Select a linear solver backend
typedef Dune::PDELab::ISTLBackend_SEQ_LS_AMG_SOR<GO> LS;
LS ls(5000,-1,false,false);
Dune::Amg::Parameters params;
//Cycle: 1 V, 2 W
params.setGamma(2);
//Presmoothing
params.setNoPreSmoothSteps(2);
//Postsmoothing
params.setNoPostSmoothSteps(2);
//Coarsegrid
params.setCoarsenTarget(1);
params.setMaxLevel(level);
ls.setparams(params);
//std::cout << "<<<5>>> done" << std::endl;
// <<<6>>> assemble and solve linear problem
typedef Dune::PDELab::StationaryLinearProblemSolver<GO,LS,U> SLP;
SLP slp(go,ls,1e-6,1e-99,-1);
slp.apply(u);
std::cout<< std::endl;
//std::cout << "<<<6>>> done" << std::endl;
// <<<7>>> graphical output and L2-Error
//Define discrte gridfunction for graphical output and L2-Error
typedef Dune::PDELab::DiscreteGridFunction<GFS,U> DGF;
DGF udgf(gfs,u);
// ??
/*typedef SquaredAdapter<DGF> Squared;
Squared squaredsolution(udgf);
typedef ExactSolution<GV,double> ES;
ES es(gv);
typedef DifferenceSquaredAdapter<DGF,ES> DifferenceSquaredAdapter;
DifferenceSquaredAdapter difference(udgf,es);*/
// #### OUTPUT ####
/* std::cout << "===================Statistics for LS_AMG_SOR===================" <<std::endl;
std::cout << "NVert: "<< gfs.globalSize() << std::endl;
std::cout << "NEL: "<< gv.size(0) << std::endl;
std::cout << "L2-Error: "<< l2norm(gv,difference,2) << std::endl;
std::cout << "Solvertime: "<< slp.result().elapsed << std::endl;
std::cout << "KonvRate: "<< slp.result().conv_rate << std::endl;
std::cout << "iters "<< slp.result().iterations <<std::endl;*/
std::cout << "===================Statistics for LS_AMG_SOR===================" <<std::endl;
std::cout << gfs.globalSize() << std::endl;
std::cout << gv.size(0) << std::endl;
//std::cout << l2norm(gv,difference,2) << std::endl;
//std::cout << slp.result().elapsed << std::endl;
std::cout << slp.result().conv_rate << std::endl;
std::cout << slp.result().iterations <<std::endl;
// graphical output
Dune::VTKWriter<GV> vtkwriter(gv,Dune::VTK::conforming);
vtkwriter.addVertexData(new Dune::PDELab::VTKGridFunctionAdapter<DGF>(udgf,"solution"));
vtkwriter.write(outputname,Dune::VTK::appendedraw);
//std::cout << "<<<7>>> done" << std::endl;
}
/*
##################################################################################################################
##################################################################################################################
BiCGSTAB AMG SSOR
##################################################################################################################
##################################################################################################################
*/
template<typename Real, typename U, typename GV, typename GFS,
typename CC, typename BCTypeParam, typename BETA>
void sdsolver_BICGSTAB_AMG( U& u, const GV& gv, const GFS& gfs, const CC& cc,
const BCTypeParam& bctype, const BETA& beta,int level, bool damping, const std::string outputname)
{
// <<<3>>> Make grid operator space
typedef Example02LocalOperator<BCTypeParam,BETA> LOP; // operator including boundary
LOP lop( bctype, beta);
typedef Dune::PDELab::ISTLMatrixBackend MBE;
typedef Dune::PDELab::GridOperator<
GFS,GFS, /* ansatz and test space */
LOP, /* local operator */
MBE, /* matrix backend */
Real,Real,Real, /* field types for domain, range and jacobian */
CC,CC /* constraints transformation for ansatz and test space */
> GO;
GO go(gfs,cc,gfs,cc,lop);
// <<<5>>> Select a linear solver backend
typedef Dune::PDELab::ISTLBackend_SEQ_BCGS_AMG_SSOR<GO> LS;
LS ls(5000,-1,false,false);
Dune::Amg::Parameters params;
//Cycle: 1 V, 2 W
params.setGamma(2);
//Presmoothing
params.setNoPreSmoothSteps(2);
//Postsmoothing
params.setNoPostSmoothSteps(2);
//Coarsegrid
params.setCoarsenTarget(1);
params.setMaxLevel(level);
ls.setparams(params);
//std::cout << "<<<5>>> done" << std::endl;
// <<<6>>> assemble and solve linear problem
typedef Dune::PDELab::StationaryLinearProblemSolver<GO,LS,U> SLP;
SLP slp(go,ls,1e-6,1e-99,-1);
slp.apply(u);
std::cout<< std::endl;
//std::cout << "<<<6>>> done" << std::endl;
// <<<7>>> graphical output and L2-Error
//Define discrte gridfunction for graphical output and L2-Error
typedef Dune::PDELab::DiscreteGridFunction<GFS,U> DGF;
DGF udgf(gfs,u);
// ??
/*typedef SquaredAdapter<DGF> Squared;
Squared squaredsolution(udgf);
typedef ExactSolution<GV,double> ES;
ES es(gv);
typedef DifferenceSquaredAdapter<DGF,ES> DifferenceSquaredAdapter;
DifferenceSquaredAdapter difference(udgf,es);*/
// #### OUTPUT ####
/* std::cout << "===================Statistics for BiCGSTAB_AMG_SSOR===================" <<std::endl;
std::cout << "NVert: "<< gfs.globalSize() << std::endl;
std::cout << "NEL: "<< gv.size(0) << std::endl;
std::cout << "L2-Error: "<< l2norm(gv,difference,2) << std::endl;
std::cout << "Solvertime: "<< slp.result().elapsed << std::endl;
std::cout << "KonvRate: "<< slp.result().conv_rate << std::endl;
std::cout << "iters "<< slp.result().iterations <<std::endl;*/
std::cout << "===================Statistics for BiCGSTAB_AMG_SSOR===================" <<std::endl;
std::cout << gfs.globalSize() << std::endl;
std::cout << gv.size(0) << std::endl;
//std::cout << l2norm(gv,difference,2) << std::endl;
//std::cout << slp.result().elapsed << std::endl;
std::cout << slp.result().conv_rate << std::endl;
std::cout << slp.result().iterations <<std::endl;
// graphical output
Dune::VTKWriter<GV> vtkwriter(gv,Dune::VTK::conforming);
vtkwriter.addVertexData(new Dune::PDELab::VTKGridFunctionAdapter<DGF>(udgf,"solution"));
vtkwriter.write(outputname,Dune::VTK::appendedraw);
//std::cout << "<<<7>>> done" << std::endl;
}
/*
##################################################################################################################
##################################################################################################################
BiCGSTAB SSOR
##################################################################################################################
##################################################################################################################
*/
template<typename Real, typename U, typename GV, typename GFS,
typename CC, typename BCTypeParam, typename BETA>
void sdsolver_BICGSTAB_SSOR( U& u, const GV& gv, const GFS& gfs, const CC& cc,
const BCTypeParam& bctype, const BETA& beta, bool damping, const std::string outputname)
{
// <<<3>>> Make grid operator space
typedef Example02LocalOperator<BCTypeParam,BETA> LOP; // operator including boundary
LOP lop( bctype, beta);
typedef Dune::PDELab::ISTLMatrixBackend MBE;
typedef Dune::PDELab::GridOperator<
GFS,GFS, /* ansatz and test space */
LOP, /* local operator */
MBE, /* matrix backend */
Real,Real,Real, /* field types for domain, range and jacobian */
CC,CC /* constraints transformation for ansatz and test space */
> GO;
GO go(gfs,cc,gfs,cc,lop);
//std::cout << "<<<3>>> done" << std::endl;
// <<<5>>> Select a linear solver backend
typedef Dune::PDELab::ISTLBackend_SEQ_BCGS_SSOR LS;
LS ls(5000,true);
//std::cout << "<<<5>>> done" << std::endl;
// <<<6>>> assemble and solve linear problem
typedef Dune::PDELab::StationaryLinearProblemSolver<GO,LS,U> SLP;
SLP slp(go,ls,u,1e-6);
slp.apply(u);
std::cout<< std::endl;
//std::cout << "<<<6>>> done" << std::endl;
// <<<7>>> graphical output and L2-Error
//Define discrte gridfunction for graphical output and L2-Error
typedef Dune::PDELab::DiscreteGridFunction<GFS,U> DGF;
DGF udgf(gfs,u);
// ??
/*typedef SquaredAdapter<DGF> Squared;
Squared squaredsolution(udgf);
typedef ExactSolution<GV,double> ES;
ES es(gv);
typedef DifferenceSquaredAdapter<DGF,ES> DifferenceSquaredAdapter;
DifferenceSquaredAdapter difference(udgf,es);*/
// #### OUTPUT ####
/* std::cout << "===================Statistics for LS_AMG_SOR===================" <<std::endl;
std::cout << "NVert: "<< gfs.globalSize() << std::endl;
std::cout << "NEL: "<< gv.size(0) << std::endl;
std::cout << "L2-Error: "<< l2norm(gv,difference,2) << std::endl;
std::cout << "Solvertime: "<< slp.result().elapsed << std::endl;
std::cout << "KonvRate: "<< slp.result().conv_rate << std::endl;
std::cout << "iters "<< slp.result().iterations <<std::endl;*/
std::cout << "===================Statistics for BiCGSTAB_SSOR===================" <<std::endl;
std::cout << gfs.globalSize() << std::endl;
std::cout << gv.size(0) << std::endl;
//std::cout << l2norm(gv,difference,2) << std::endl;
//std::cout << slp.result().elapsed << std::endl;
std::cout << slp.result().conv_rate << std::endl;
std::cout << slp.result().iterations <<std::endl;
// graphical output
Dune::VTKWriter<GV> vtkwriter(gv,Dune::VTK::conforming);
vtkwriter.addVertexData(new Dune::PDELab::VTKGridFunctionAdapter<DGF>(udgf,"solution"));
vtkwriter.write(outputname,Dune::VTK::appendedraw);
//std::cout << "<<<7>>> done" << std::endl;
}
/*
##################################################################################################################
##################################################################################################################
BiCGSTAB JAC
##################################################################################################################
##################################################################################################################
*/
template<typename Real, typename U, typename GV, typename GFS,
typename CC, typename BCTypeParam, typename BETA>
void sdsolver_BICGSTAB_JAC( U& u, const GV& gv, const GFS& gfs, const CC& cc,
const BCTypeParam& bctype, const BETA& beta, bool damping, const std::string outputname)
{
// <<<3>>> Make grid operator space
typedef Example02LocalOperator<BCTypeParam,BETA> LOP; // operator including boundary
LOP lop( bctype, beta);
typedef Dune::PDELab::ISTLMatrixBackend MBE;
typedef Dune::PDELab::GridOperator<
GFS,GFS, /* ansatz and test space */
LOP, /* local operator */
MBE, /* matrix backend */
Real,Real,Real, /* field types for domain, range and jacobian */
CC,CC /* constraints transformation for ansatz and test space */
> GO;
GO go(gfs,cc,gfs,cc,lop);
// <<<5>>> Select a linear solver backend
typedef Dune::PDELab::ISTLBackend_SEQ_BCGS_Jac LS;
LS ls(5000,true);
//std::cout << "<<<5>>> done" << std::endl;
// <<<6>>> assemble and solve linear problem
typedef Dune::PDELab::StationaryLinearProblemSolver<GO,LS,U> SLP;
SLP slp(go,ls,u,1e-6);
slp.apply(u);
std::cout<< std::endl;
//std::cout << "<<<6>>> done" << std::endl;
// <<<7>>> graphical output and L2-Error
//Define discrte gridfunction for graphical output and L2-Error
typedef Dune::PDELab::DiscreteGridFunction<GFS,U> DGF;
DGF udgf(gfs,u);
// ??
/*typedef SquaredAdapter<DGF> Squared;
Squared squaredsolution(udgf);
typedef ExactSolution<GV,double> ES;
ES es(gv);
typedef DifferenceSquaredAdapter<DGF,ES> DifferenceSquaredAdapter;
DifferenceSquaredAdapter difference(udgf,es);*/
// #### OUTPUT ####
/* std::cout << "===================Statistics for BiCGSTAB_AMG_SSOR===================" <<std::endl;
std::cout << "NVert: "<< gfs.globalSize() << std::endl;
std::cout << "NEL: "<< gv.size(0) << std::endl;
std::cout << "L2-Error: "<< l2norm(gv,difference,2) << std::endl;
std::cout << "Solvertime: "<< slp.result().elapsed << std::endl;
std::cout << "KonvRate: "<< slp.result().conv_rate << std::endl;
std::cout << "iters "<< slp.result().iterations <<std::endl;*/
std::cout << "===================Statistics for BiCGSTAB_JAC===================" <<std::endl;
std::cout << gfs.globalSize() << std::endl;
std::cout << gv.size(0) << std::endl;
//std::cout << l2norm(gv,difference,2) << std::endl;
//std::cout << slp.result().elapsed << std::endl;
std::cout << slp.result().conv_rate << std::endl;
std::cout << slp.result().iterations <<std::endl;
// graphical output
Dune::VTKWriter<GV> vtkwriter(gv,Dune::VTK::conforming);
vtkwriter.addVertexData(new Dune::PDELab::VTKGridFunctionAdapter<DGF>(udgf,"solution"));
vtkwriter.write(outputname,Dune::VTK::appendedraw);
//std::cout << "<<<7>>> done" << std::endl;
}
/*
##################################################################################################################
##################################################################################################################
CG JAC
##################################################################################################################
##################################################################################################################
*/
template<typename Real, typename U, typename GV, typename GFS,
typename CC, typename BCTypeParam, typename BETA>
void sdsolver_CG_JAC( U& u, const GV& gv, const GFS& gfs, const CC& cc,
const BCTypeParam& bctype, const BETA& beta, bool damping, const std::string outputname)
{
// <<<3>>> Make grid operator space
typedef Example02LocalOperator<BCTypeParam,BETA> LOP; // operator including boundary
LOP lop( bctype, beta);
typedef Dune::PDELab::ISTLMatrixBackend MBE;
typedef Dune::PDELab::GridOperator<
GFS,GFS, /* ansatz and test space */
LOP, /* local operator */
MBE, /* matrix backend */
Real,Real,Real, /* field types for domain, range and jacobian */
CC,CC /* constraints transformation for ansatz and test space */
> GO;
GO go(gfs,cc,gfs,cc,lop);
// <<<5>>> Select a linear solver backend
typedef Dune::PDELab::ISTLBackend_SEQ_CG_Jac LS;
LS ls(5000,true);
//std::cout << "<<<5>>> done" << std::endl;
// <<<6>>> assemble and solve linear problem
typedef Dune::PDELab::StationaryLinearProblemSolver<GO,LS,U> SLP;
SLP slp(go,ls,u,1e-6);
slp.apply(u);
std::cout<< std::endl;
//std::cout << "<<<6>>> done" << std::endl;
// <<<7>>> graphical output and L2-Error
//Define discrte gridfunction for graphical output and L2-Error
typedef Dune::PDELab::DiscreteGridFunction<GFS,U> DGF;
DGF udgf(gfs,u);
// ??
/*typedef SquaredAdapter<DGF> Squared;
Squared squaredsolution(udgf);
typedef ExactSolution<GV,double> ES;
ES es(gv);
typedef DifferenceSquaredAdapter<DGF,ES> DifferenceSquaredAdapter;
DifferenceSquaredAdapter difference(udgf,es);*/
// #### OUTPUT ####
/* std::cout << "===================Statistics for BiCGSTAB_AMG_SSOR===================" <<std::endl;
std::cout << "NVert: "<< gfs.globalSize() << std::endl;
std::cout << "NEL: "<< gv.size(0) << std::endl;
std::cout << "L2-Error: "<< l2norm(gv,difference,2) << std::endl;
std::cout << "Solvertime: "<< slp.result().elapsed << std::endl;
std::cout << "KonvRate: "<< slp.result().conv_rate << std::endl;
std::cout << "iters "<< slp.result().iterations <<std::endl;*/
std::cout << "===================Statistics for CG_JAC===================" <<std::endl;
std::cout << gfs.globalSize() << std::endl;
std::cout << gv.size(0) << std::endl;
//std::cout << l2norm(gv,difference,2) << std::endl;
//std::cout << slp.result().elapsed << std::endl;
std::cout << slp.result().conv_rate << std::endl;
std::cout << slp.result().iterations <<std::endl;
// graphical output
Dune::VTKWriter<GV> vtkwriter(gv,Dune::VTK::conforming);
vtkwriter.addVertexData(new Dune::PDELab::VTKGridFunctionAdapter<DGF>(udgf,"solution"));
vtkwriter.write(outputname,Dune::VTK::appendedraw);
//std::cout << "<<<7>>> done" << std::endl;
}
/*
##################################################################################################################
##################################################################################################################
LOOP SOLVER AMG SSOR
##################################################################################################################
##################################################################################################################
*/
template<typename Real, typename U, typename GV, typename GFS,
typename CC, typename BCTypeParam, typename BETA>
void sdsolver_LS_AMG_SS( U& u, const GV& gv, const GFS& gfs, const CC& cc,
const BCTypeParam& bctype, const BETA& beta,int level, bool damping, const std::string outputname)
{
// <<<3>>> Make grid operator space
typedef Example02LocalOperator<BCTypeParam,BETA> LOP; // operator including boundary
LOP lop( bctype, beta);
typedef Dune::PDELab::ISTLMatrixBackend MBE;
typedef Dune::PDELab::GridOperator<
GFS,GFS, /* ansatz and test space */
LOP, /* local operator */
MBE, /* matrix backend */
Real,Real,Real, /* field types for domain, range and jacobian */
CC,CC /* constraints transformation for ansatz and test space */
> GO;
GO go(gfs,cc,gfs,cc,lop);
//std::cout << "<<<3>>> done" << std::endl;
// <<<5>>> Select a linear solver backend
typedef Dune::PDELab::ISTLBackend_SEQ_LS_AMG_SSOR<GO> LS;
LS ls(5000,-1,false,false);
Dune::Amg::Parameters params;
//Cycle: 1 V, 2 W
params.setGamma(2);
//Presmoothing
params.setNoPreSmoothSteps(2);
//Postsmoothing
params.setNoPostSmoothSteps(2);
//Coarsegrid
params.setCoarsenTarget(1);
params.setMaxLevel(level);
ls.setparams(params);
//std::cout << "<<<5>>> done" << std::endl;
// <<<6>>> assemble and solve linear problem
typedef Dune::PDELab::StationaryLinearProblemSolver<GO,LS,U> SLP;
SLP slp(go,ls,1e-6,1e-99,-1);
slp.apply(u);
std::cout<< std::endl;
//std::cout << "<<<6>>> done" << std::endl;
// <<<7>>> graphical output and L2-Error
//Define discrte gridfunction for graphical output and L2-Error
typedef Dune::PDELab::DiscreteGridFunction<GFS,U> DGF;
DGF udgf(gfs,u);
// ??
/*typedef SquaredAdapter<DGF> Squared;
Squared squaredsolution(udgf);
typedef ExactSolution<GV,double> ES;
ES es(gv);
typedef DifferenceSquaredAdapter<DGF,ES> DifferenceSquaredAdapter;
DifferenceSquaredAdapter difference(udgf,es);*/
// #### OUTPUT ####
/* std::cout << "===================Statistics for LS_AMG_SOR===================" <<std::endl;
std::cout << "NVert: "<< gfs.globalSize() << std::endl;
std::cout << "NEL: "<< gv.size(0) << std::endl;
std::cout << "L2-Error: "<< l2norm(gv,difference,2) << std::endl;
std::cout << "Solvertime: "<< slp.result().elapsed << std::endl;
std::cout << "KonvRate: "<< slp.result().conv_rate << std::endl;
std::cout << "iters "<< slp.result().iterations <<std::endl;*/
std::cout << "===================Statistics for LS_AMG_SSOR===================" <<std::endl;
std::cout << gfs.globalSize() << std::endl;
std::cout << gv.size(0) << std::endl;
//std::cout << l2norm(gv,difference,2) << std::endl;
//std::cout << slp.result().elapsed << std::endl;
std::cout << slp.result().conv_rate << std::endl;
std::cout << slp.result().iterations <<std::endl;
// graphical output
Dune::VTKWriter<GV> vtkwriter(gv,Dune::VTK::conforming);
vtkwriter.addVertexData(new Dune::PDELab::VTKGridFunctionAdapter<DGF>(udgf,"solution"));
vtkwriter.write(outputname,Dune::VTK::appendedraw);
//std::cout << "<<<7>>> done" << std::endl;
}
/*
##################################################################################################################
##################################################################################################################
BiCGSTAB AMG SOR
##################################################################################################################
##################################################################################################################
*/
template<typename Real, typename U, typename GV, typename GFS,
typename CC, typename BCTypeParam, typename BETA>
void sdsolver_BICGSTAB_AMG_SOR( U& u, const GV& gv, const GFS& gfs, const CC& cc,
const BCTypeParam& bctype, const BETA& beta,int level, bool damping, const std::string outputname)
{
// <<<3>>> Make grid operator space
typedef Example02LocalOperator<BCTypeParam,BETA> LOP; // operator including boundary
LOP lop( bctype, beta);
typedef Dune::PDELab::ISTLMatrixBackend MBE;
typedef Dune::PDELab::GridOperator<
GFS,GFS, /* ansatz and test space */
LOP, /* local operator */
MBE, /* matrix backend */
Real,Real,Real, /* field types for domain, range and jacobian */
CC,CC /* constraints transformation for ansatz and test space */
> GO;
GO go(gfs,cc,gfs,cc,lop);
// <<<5>>> Select a linear solver backend
typedef Dune::PDELab::ISTLBackend_SEQ_BCGS_AMG_SOR<GO> LS;
LS ls(5000,-1,false,false);
Dune::Amg::Parameters params;
//Cycle: 1 V, 2 W
params.setGamma(2);
//Presmoothing
params.setNoPreSmoothSteps(2);
//Postsmoothing
params.setNoPostSmoothSteps(2);
//Coarsegrid
params.setCoarsenTarget(1);
params.setMaxLevel(level);
ls.setparams(params);
//std::cout << "<<<5>>> done" << std::endl;
// <<<6>>> assemble and solve linear problem
typedef Dune::PDELab::StationaryLinearProblemSolver<GO,LS,U> SLP;
SLP slp(go,ls,1e-6,1e-99,-1);
slp.apply(u);
std::cout<< std::endl;
//std::cout << "<<<6>>> done" << std::endl;
// <<<7>>> graphical output and L2-Error
//Define discrte gridfunction for graphical output and L2-Error
typedef Dune::PDELab::DiscreteGridFunction<GFS,U> DGF;
DGF udgf(gfs,u);
// ??
/*typedef SquaredAdapter<DGF> Squared;
Squared squaredsolution(udgf);
typedef ExactSolution<GV,double> ES;
ES es(gv);
typedef DifferenceSquaredAdapter<DGF,ES> DifferenceSquaredAdapter;
DifferenceSquaredAdapter difference(udgf,es);*/
// #### OUTPUT ####
/* std::cout << "===================Statistics for BiCGSTAB_AMG_SSOR===================" <<std::endl;
std::cout << "NVert: "<< gfs.globalSize() << std::endl;
std::cout << "NEL: "<< gv.size(0) << std::endl;
std::cout << "L2-Error: "<< l2norm(gv,difference,2) << std::endl;
std::cout << "Solvertime: "<< slp.result().elapsed << std::endl;
std::cout << "KonvRate: "<< slp.result().conv_rate << std::endl;
std::cout << "iters "<< slp.result().iterations <<std::endl;*/
std::cout << "===================Statistics for BiCGSTAB_AMG_SOR===================" <<std::endl;
std::cout << gfs.globalSize() << std::endl;
std::cout << gv.size(0) << std::endl;
//std::cout << l2norm(gv,difference,2) << std::endl;
//std::cout << slp.result().elapsed << std::endl;
std::cout << slp.result().conv_rate << std::endl;
std::cout << slp.result().iterations <<std::endl;
// graphical output
Dune::VTKWriter<GV> vtkwriter(gv,Dune::VTK::conforming);
vtkwriter.addVertexData(new Dune::PDELab::VTKGridFunctionAdapter<DGF>(udgf,"solution"));
vtkwriter.write(outputname,Dune::VTK::appendedraw);
//std::cout << "<<<7>>> done" << std::endl;
}
/*
##################################################################################################################
##################################################################################################################
CG AMG SSOR
##################################################################################################################
##################################################################################################################
*/
template<typename Real, typename U, typename GV, typename GFS,
typename CC, typename BCTypeParam, typename BETA>
void sdsolver_CG_AMG_SSOR( U& u, const GV& gv, const GFS& gfs, const CC& cc,
const BCTypeParam& bctype, const BETA& beta,int level, bool damping, const std::string outputname)
{
// <<<3>>> Make grid operator space
typedef Example02LocalOperator<BCTypeParam,BETA> LOP; // operator including boundary
LOP lop( bctype, beta);
typedef Dune::PDELab::ISTLMatrixBackend MBE;
typedef Dune::PDELab::GridOperator<
GFS,GFS, /* ansatz and test space */
LOP, /* local operator */
MBE, /* matrix backend */
Real,Real,Real, /* field types for domain, range and jacobian */
CC,CC /* constraints transformation for ansatz and test space */
> GO;
GO go(gfs,cc,gfs,cc,lop);
// <<<5>>> Select a linear solver backend
typedef Dune::PDELab::ISTLBackend_SEQ_CG_AMG_SSOR<GO> LS;
LS ls(5000,-1,false,false);
Dune::Amg::Parameters params;
//Cycle: 1 V, 2 W
params.setGamma(2);
//Presmoothing
params.setNoPreSmoothSteps(2);
//Postsmoothing
params.setNoPostSmoothSteps(2);
//Coarsegrid
params.setCoarsenTarget(1);
params.setMaxLevel(level);
ls.setparams(params);
//std::cout << "<<<5>>> done" << std::endl;
// <<<6>>> assemble and solve linear problem
typedef Dune::PDELab::StationaryLinearProblemSolver<GO,LS,U> SLP;
SLP slp(go,ls,1e-6,1e-99,-1);
slp.apply(u);
std::cout<< std::endl;
//std::cout << "<<<6>>> done" << std::endl;
// <<<7>>> graphical output and L2-Error
//Define discrte gridfunction for graphical output and L2-Error
typedef Dune::PDELab::DiscreteGridFunction<GFS,U> DGF;
DGF udgf(gfs,u);
// ??
/*typedef SquaredAdapter<DGF> Squared;
Squared squaredsolution(udgf);
typedef ExactSolution<GV,double> ES;
ES es(gv);
typedef DifferenceSquaredAdapter<DGF,ES> DifferenceSquaredAdapter;
DifferenceSquaredAdapter difference(udgf,es);*/
// #### OUTPUT ####
/* std::cout << "===================Statistics for BiCGSTAB_AMG_SSOR===================" <<std::endl;
std::cout << "NVert: "<< gfs.globalSize() << std::endl;
std::cout << "NEL: "<< gv.size(0) << std::endl;
std::cout << "L2-Error: "<< l2norm(gv,difference,2) << std::endl;
std::cout << "Solvertime: "<< slp.result().elapsed << std::endl;
std::cout << "KonvRate: "<< slp.result().conv_rate << std::endl;
std::cout << "iters "<< slp.result().iterations <<std::endl;*/
std::cout << "===================Statistics for CG_AMG_SSOR===================" <<std::endl;
std::cout << gfs.globalSize() << std::endl;
std::cout << gv.size(0) << std::endl;
//std::cout << l2norm(gv,difference,2) << std::endl;
//std::cout << slp.result().elapsed << std::endl;
std::cout << slp.result().conv_rate << std::endl;
std::cout << slp.result().iterations <<std::endl;
// graphical output
Dune::VTKWriter<GV> vtkwriter(gv,Dune::VTK::conforming);
vtkwriter.addVertexData(new Dune::PDELab::VTKGridFunctionAdapter<DGF>(udgf,"solution"));
vtkwriter.write(outputname,Dune::VTK::appendedraw);
//std::cout << "<<<7>>> done" << std::endl;
}
/*
##################################################################################################################
##################################################################################################################
CG SSOR
##################################################################################################################
##################################################################################################################
*/
template<typename Real, typename U, typename GV, typename GFS,
typename CC, typename BCTypeParam, typename BETA>
void sdsolver_CG_SSOR( U& u, const GV& gv, const GFS& gfs, const CC& cc,
const BCTypeParam& bctype, const BETA& beta, bool damping, const std::string outputname)
{
// <<<3>>> Make grid operator space
typedef Example02LocalOperator<BCTypeParam,BETA> LOP; // operator including boundary
LOP lop( bctype, beta);
typedef Dune::PDELab::ISTLMatrixBackend MBE;
typedef Dune::PDELab::GridOperator<
GFS,GFS, /* ansatz and test space */
LOP, /* local operator */
MBE, /* matrix backend */
Real,Real,Real, /* field types for domain, range and jacobian */
CC,CC /* constraints transformation for ansatz and test space */
> GO;
GO go(gfs,cc,gfs,cc,lop);
// <<<5>>> Select a linear solver backend
typedef Dune::PDELab::ISTLBackend_SEQ_CG_SSOR LS;
LS ls(5000,true);
//std::cout << "<<<5>>> done" << std::endl;
// <<<6>>> assemble and solve linear problem
typedef Dune::PDELab::StationaryLinearProblemSolver<GO,LS,U> SLP;
SLP slp(go,ls,u,1e-6);
slp.apply(u);
std::cout<< std::endl;
//std::cout << "<<<6>>> done" << std::endl;
// <<<7>>> graphical output and L2-Error
//Define discrte gridfunction for graphical output and L2-Error
typedef Dune::PDELab::DiscreteGridFunction<GFS,U> DGF;
DGF udgf(gfs,u);
// ??
/*typedef SquaredAdapter<DGF> Squared;
Squared squaredsolution(udgf);
typedef ExactSolution<GV,double> ES;
ES es(gv);
typedef DifferenceSquaredAdapter<DGF,ES> DifferenceSquaredAdapter;
DifferenceSquaredAdapter difference(udgf,es);*/
// #### OUTPUT ####
/* std::cout << "===================Statistics for BiCGSTAB_AMG_SSOR===================" <<std::endl;
std::cout << "NVert: "<< gfs.globalSize() << std::endl;
std::cout << "NEL: "<< gv.size(0) << std::endl;
std::cout << "L2-Error: "<< l2norm(gv,difference,2) << std::endl;
std::cout << "Solvertime: "<< slp.result().elapsed << std::endl;
std::cout << "KonvRate: "<< slp.result().conv_rate << std::endl;
std::cout << "iters "<< slp.result().iterations <<std::endl;*/
std::cout << "===================Statistics for CG_SSOR===================" <<std::endl;
std::cout << gfs.globalSize() << std::endl;
std::cout << gv.size(0) << std::endl;
//std::cout << l2norm(gv,difference,2) << std::endl;
//std::cout << slp.result().elapsed << std::endl;
std::cout << slp.result().conv_rate << std::endl;
std::cout << slp.result().iterations <<std::endl;
// graphical output
Dune::VTKWriter<GV> vtkwriter(gv,Dune::VTK::conforming);
vtkwriter.addVertexData(new Dune::PDELab::VTKGridFunctionAdapter<DGF>(udgf,"solution"));
vtkwriter.write(outputname,Dune::VTK::appendedraw);
//std::cout << "<<<7>>> done" << std::endl;
}
template<class GV> void example02_Q1 (const GV& gv)
{
// <<<1>>> Choose domain and range field type
typedef typename GV::Grid::ctype Coord;
typedef double Real;
//std::cout << "<<<1>>> done" << std::endl;
// <<<2>>> Make grid function space
//std::cout << "=================== Q1 - Elements ===================" << std::endl;
//typedef Dune::PDELab::QkLocalFiniteElementMap<GV,Coord,Real,1> FEM; // Quadrilateral Elements
std::cout << "=================== Q2 - Elements ===================" << std::endl;
typedef Dune::PDELab::QkLocalFiniteElementMap<GV,Coord,Real,2> FEM; // Quadrilateral Elements
FEM fem(gv);
typedef Dune::PDELab::ConformingDirichletConstraints CONSTRAINTS; // constraints class
typedef Dune::PDELab::ISTLVectorBackend<> VBE;
typedef Dune::PDELab::GridFunctionSpace<GV,FEM,CONSTRAINTS,VBE> GFS;
GFS gfs(gv,fem);
typedef BETA<GV,Real> BETAType;
BETAType beta(gv);
BCTypeParam<BETAType> bctype(beta); // boundary condition type
typedef typename GFS::template ConstraintsContainer<Real>::Type CC;
CC cc;
Dune::PDELab::constraints( bctype, gfs, cc ); // assemble constraints
//std::cout << "constrained dofs=" << cc.size() << " of " << gfs.globalSize() << std::endl;
//std::cout << "<<<2>>> done" << std::endl;
// <<<3>>> moved into sdsolver()
// <<<4>>> Make FE function extending Dirichlet boundary conditions
// BackendVectorSelector can be used if a vector is needed but no
// GridOperator is available yet
typedef typename Dune::PDELab::BackendVectorSelector<GFS,Real>::Type U;
U u(gfs,0.0);
//Only for use without gmesh, because of dx
//typedef BCExtension_2b<GV,Real,VEC> G; // boundary value + extension
//G g(gv,dx);
//Only for use withgmesh
typedef BCExtension<GV,Real> G; // boundary value + extension
G g(gv);
Dune::PDELab::interpolate(g,gfs,u); // interpolate coefficient vector
//std::cout << "<<<4>>> done" << std::endl;
// <<<5>>> - <<<7>>> moved inside "sdsolver"
int wahl;
printf("LS_AMG_SOR: 1, BICGSTAB_AMG_SSOR: 2, BICGSTAB_SSOR: 3, BICGSTAB_JAC: 4, CG_JAC: 5, LS_AMG_SSOR: 6, BiCGSTAB_AMG_SOR: 7, CG_AMG_SSOR: 8, CG_SSOR: 9 .... choice: ");
scanf("%d",&wahl);
int level;
if (wahl == 1)
{
printf("Level input: ");
scanf("%d",&level);
std::cout << "Used levels:" << level+1 << std::endl;
sdsolver_LS_AMG<Real>( u,gv,gfs,cc,bctype,beta,level, false, "LS_AMG_SOR");
}
if (wahl == 2)
{
printf("Level input: ");
scanf("%d",&level);
std::cout << "Used levels:" << level+1 << std::endl;
sdsolver_BICGSTAB_AMG<Real>( u,gv,gfs,cc,bctype,beta,level, false, "BICGSTAB_AMG_SSOR");
}
if (wahl == 3)
sdsolver_BICGSTAB_SSOR<Real>( u,gv,gfs,cc,bctype,beta, false, "BICGSTAB_SSOR");
if (wahl == 4)
sdsolver_BICGSTAB_JAC<Real>( u,gv,gfs,cc,bctype,beta, false, "BICGSTAB_JAC");
if (wahl == 5)
sdsolver_CG_JAC<Real>( u,gv,gfs,cc,bctype,beta, false, "CG_JAC");
if (wahl == 6)
{
printf("Level input: ");
scanf("%d",&level);
std::cout << "Used levels:" << level+1 << std::endl;
sdsolver_LS_AMG_SS<Real>( u,gv,gfs,cc,bctype,beta,level, false, "LS_AMG_SSOR");
}
if (wahl == 7)
{
printf("Level input: ");
scanf("%d",&level);
std::cout << "Used levels:" << level+1 << std::endl;
sdsolver_BICGSTAB_AMG_SOR<Real>( u,gv,gfs,cc,bctype,beta,level, false, "BICGSTAB_AMG_SOR");
}
if (wahl == 8)
{
printf("Level input: ");
scanf("%d",&level);
std::cout << "Used levels:" << level+1 << std::endl;
sdsolver_CG_AMG_SSOR<Real>( u,gv,gfs,cc,bctype,beta,level, false, "CG_AMG_SSOR");
}
if (wahl == 9)
sdsolver_CG_SSOR<Real>( u,gv,gfs,cc,bctype,beta, false, "CG_SSOR");
}
-------------- next part --------------
/** \brief A function that defines Dirichlet boundary conditions AND
its extension to the interior */
template<typename GV, typename RF>
class BCExtension
: public Dune::PDELab::GridFunctionBase<Dune::PDELab::
GridFunctionTraits<GV,RF,1,Dune::FieldVector<RF,1> >, BCExtension<GV,RF> > {
const GV& gv;
public:
typedef Dune::PDELab::GridFunctionTraits<GV,RF,1,Dune::FieldVector<RF,1> > Traits;
//! construct from grid view
BCExtension (const GV& gv_) : gv(gv_) {}
//! evaluate extended function on element
inline void evaluate (const typename Traits::ElementType& e,
const typename Traits::DomainType& xlocal,
typename Traits::RangeType& y) const
{
const int dim = Traits::GridViewType::Grid::dimension;
typedef typename Traits::GridViewType::Grid::ctype ctype;
Dune::FieldVector<ctype,dim> x = e.geometry().global(xlocal);
y = 0.0;
return;
}
//! get a reference to the grid view
inline const GV& getGridView () {return gv;}
};
-------------- next part --------------
#include<dune/common/fvector.hh>
#include<dune/pdelab/common/function.hh>
#include<dune/pdelab/constraints/constraintsparameters.hh>
//! \brief Parameter class selecting boundary conditions
template<class BETA>
class BCTypeParam
: public Dune::PDELab::DirichletConstraintsParameters
{
const BETA &betaF;
public:
BCTypeParam(const BETA &betaF_) : betaF(betaF_) { }
//! Test whether boundary is Dirichlet-constrained
template<typename I>
bool isDirichlet(const I & intersection,
const Dune::FieldVector<typename I::ctype, I::dimension-1> & coord
) const
{
Dune::FieldVector<typename I::ctype, I::dimension>
xg = intersection.geometry().global( coord );
return true; // Dirichlet b.c. on all other boundaries
}
};
-------------- next part --------------
#include<dune/geometry/quadraturerules.hh>
#include<dune/geometry/referenceelements.hh>
#include<dune/pdelab/localoperator/defaultimp.hh>
#include<dune/pdelab/localoperator/flags.hh>
#include<dune/pdelab/localoperator/pattern.hh>
/** a local operator for solving the equation
*
* - \Delta u = f in \Omega
* u = g on \Gamma_D\subseteq\partial\Omega
*
* with conforming finite elements on all types of grids in any dimension
*
* \tparam BCType parameter class indicating the type of boundary condition
*/
template<class BCType, class BETA>
class Example02LocalOperator :
public Dune::PDELab::NumericalJacobianApplyVolume<Example02LocalOperator<BCType,BETA>>,
public Dune::PDELab::NumericalJacobianVolume<Example02LocalOperator<BCType,BETA> >,
public Dune::PDELab::NumericalJacobianApplyBoundary<Example02LocalOperator<BCType,BETA>>,
public Dune::PDELab::NumericalJacobianBoundary<Example02LocalOperator<BCType,BETA>>,
public Dune::PDELab::FullVolumePattern,
public Dune::PDELab::LocalOperatorDefaultFlags
{
public:
// pattern assembly flags
enum { doPatternVolume = true };
// residual assembly flags
enum { doAlphaVolume = true };
//enum { doAlphaBoundary = true }; // assemble boundary
Example02LocalOperator (const BCType& bctype_, // boundary cond. type
const BETA& betaF_, // convection field
unsigned int intorder_=3) :
bctype(bctype_), betaF(betaF_), intorder(intorder_)
{}
// volume integral depending on test and ansatz functions
template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
{
// pull in some standard functions
using std::pow;
using std::abs;
// assume Galerkin: lfsu == lfsv
// This yields more efficient code since the local functionspace only
// needs to be evaluated once, but would be incorrect for a finite volume
// method
// dimensions
const int dim = EG::Geometry::dimension;
const int dimw = EG::Geometry::dimensionworld;
// extract some types
typedef typename LFSU::Traits::FiniteElementType::
Traits::LocalBasisType::Traits::DomainFieldType DF;
typedef typename LFSU::Traits::FiniteElementType::
Traits::LocalBasisType::Traits::RangeFieldType RF;
typedef typename LFSU::Traits::FiniteElementType::
Traits::LocalBasisType::Traits::JacobianType Jacobian;
typedef typename LFSU::Traits::FiniteElementType::
Traits::LocalBasisType::Traits::RangeType Range;
typedef Dune::FieldVector<RF,dimw> Gradient;
typedef typename LFSU::Traits::SizeType size_type;
// select quadrature rule
Dune::GeometryType gt = eg.geometry().type();
const Dune::QuadratureRule<DF,dim>&
rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
// loop over quadrature points
for (typename Dune::QuadratureRule<DF,dim>::const_iterator
it=rule.begin(); it!=rule.end(); ++it)
{
// evaluate basis functions on reference element
std::vector<Range> phi(lfsu.size());
lfsu.finiteElement().localBasis().evaluateFunction(it->position(),phi);
// compute u at integration point
RF u=0.0;
for (size_type i=0; i<lfsu.size(); ++i)
u += x(lfsu,i)*phi[i];
// evaluate gradient of basis functions on reference element
std::vector<Jacobian> js(lfsu.size());
lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
// transform gradients from reference element to real element
const Dune::FieldMatrix<DF,dimw,dim>
&jac = eg.geometry().jacobianInverseTransposed(it->position());
std::vector<Gradient> gradphi(lfsu.size());
for (size_type i=0; i<lfsu.size(); i++)
jac.mv(js[i][0],gradphi[i]);
// compute gradient of u
Gradient gradu(0.0);
for (size_type i=0; i<lfsu.size(); ++i)
gradu.axpy(x(lfsu,i),gradphi[i]);
// evaluate parameters;
// Dune::FieldVector<RF,dim>
// globalpos = eg.geometry().global(it->position());
// get the vector-field beta:
typename BETA::Traits::RangeType beta;
Dune::FieldVector<DF,dim> localcenter =
Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
betaF.evaluate( eg.entity(), localcenter, beta );
// Get RHS
typename BETA::Traits::RangeFieldType f =betaF.f(eg.entity(),it->position(),beta);
//RF f = 0;
// Standard Galerkin
RF factor = it->weight()*eg.geometry().integrationElement(it->position());
for (size_type i=0; i<lfsu.size(); ++i)
r.accumulate(lfsu, i, factor * ( (gradu*gradphi[i]) - f*phi[i] ));
}
}
private:
const BCType& bctype;
const BETA& betaF;
unsigned int intorder;
};
More information about the dune-pdelab
mailing list