Commit ec65a618 authored by Peter Bastian's avatar Peter Bastian

reorganization of integration completed

[[Imported from SVN: r21]]
parent bea625a8
......@@ -12,13 +12,15 @@
//! adaptive refinement test
template<class Grid, class Functor>
void adaptive_integration (Grid& grid, const Functor& f)
void adaptiveintegration (Grid& grid, const Functor& f)
{
// get iterator type
typedef typename Grid::template Codim<0>::LeafIterator ElementLeafIterator;
// algorithm parameters
const double tol=1E-8; // relative accuracy
const double tol=1E-8;
const int loworder=1;
const int highorder=3;
// loop over grid sequence
double oldvalue=1E100;
......@@ -28,10 +30,7 @@ void adaptive_integration (Grid& grid, const Functor& f)
double value=0;
for (ElementLeafIterator it = grid.template leafbegin<0>();
it!=grid.template leafend<0>(); ++it)
{
double error;
value += integrate(it,f,error);
}
value += integrateentity(it,f,highorder);
// print result
double estimated_error = std::abs(value-oldvalue);
......@@ -64,13 +63,17 @@ void adaptive_integration (Grid& grid, const Functor& f)
it!=grid.template leafend<0>(); ++it)
{
// error on this entity
double error;
integrate(it,f,error);
double lowresult=integrateentity(it,f,loworder);
double highresult=integrateentity(it,f,highorder);
double error = std::abs(lowresult-highresult);
// max over whole grid
maxerror = std::max(maxerror,error);
// error on father entity
double fathererror;
integrate(it->father(),f,fathererror);
double fatherlowresult=integrateentity(it->father(),f,loworder);
double fatherhighresult=integrateentity(it->father(),f,highorder);
double fathererror = std::abs(lowresult-highresult);
// local extrapolation
double extrapolatederror = error*error/(fathererror+1E-30);
......@@ -82,8 +85,9 @@ void adaptive_integration (Grid& grid, const Functor& f)
for (ElementLeafIterator it = grid.template leafbegin<0>();
it!=grid.template leafend<0>(); ++it)
{
double error;
integrate(it,f,error);
double lowresult=integrateentity(it,f,loworder);
double highresult=integrateentity(it,f,highorder);
double error = std::abs(lowresult-highresult);
if (error>kappa) grid.mark(1,it);
}
......@@ -102,7 +106,7 @@ void adaptive_integration (Grid& grid, const Functor& f)
template<class Grid>
void dowork (Grid& grid)
{
adaptive_integration(grid,Needle<typename Grid::ctype,Grid::dimension>());
adaptiveintegration(grid,Needle<typename Grid::ctype,Grid::dimension>());
}
int main(int argc, char **argv)
......@@ -121,7 +125,7 @@ int main(int argc, char **argv)
UnitCube<Dune::UGGrid<2,2>,2> uc6;
#endif
dowork(uc1.grid());
dowork(uc6.grid());
#if HAVE_MPI
MPI_Finalize();
......
......@@ -1032,12 +1032,6 @@ numberstyle=\tiny, numbersep=5pt]{../integration.cc}
\section{Adaptive integration}
\begin{lst}[File dune-grid-howto/integrate\_entity.hh] \mbox{}
\nopagebreak
\lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
numberstyle=\tiny, numbersep=5pt]{../integrate_entity.hh}
\end{lst}
\begin{lst}[File dune-grid-howto/adaptiveintegration.cc] \mbox{}
\nopagebreak
\lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
......
......@@ -78,13 +78,15 @@ bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax, int k)
for (LeafIterator it = grid.template leafbegin<0>();
it!=grid.template leafend<0>(); ++it)
{
if (indicator[mapper.map(*it)]>refinetol*globaldelta && (it.level()<lmax || !it->isRegular()))
if (indicator[mapper.map(*it)]>refinetol*globaldelta
&& (it.level()<lmax || !it->isRegular()))
{
grid.mark(1,it);
IntersectionIterator isend = it->iend();
for (IntersectionIterator is = it->ibegin(); is!=isend; ++is)
if (is.neighbor())
if (is.outside()->isLeaf() && (is.outside().level()<lmax || !is.outside()->isRegular()))
if (is.outside()->isLeaf() &&
(is.outside().level()<lmax || !is.outside()->isRegular()))
grid.mark(1,is.outside());
}
if (indicator[mapper.map(*it)]<coarsentol*globaldelta && it.level()>lmin)
......
......@@ -6,61 +6,6 @@
#include "dune/common/exceptions.hh"
#include "dune/quadrature/quadraturerules.hh"
//! compute integral of function over entity with error estimate
template<class Iterator, class Functor>
double integrate (const Iterator& it, const Functor& f, double& localerror)
{
// dimension of the entity
const int dim = Iterator::Entity::dimension;
// type used for coordinates in the grid
typedef typename Iterator::Entity::ctype ct;
// get geometry type
Dune::GeometryType gt = it->geometry().type();
// get low order quadrature
const Dune::QuadratureRule<ct,dim>&
lowrule = Dune::QuadratureRules<ct,dim>::rule(gt,1);
// get higher order quadrature
const Dune::QuadratureRule<ct,dim>&
highrule = Dune::QuadratureRules<ct,dim>::rule(gt,3);
// ensure that highrule has higher order
if (!(highrule.order()>lowrule.order()))
DUNE_THROW(Dune::Exception,"high order rule not better than low order rule");
// compute integral with low order rule
double lowresult=0;
for (typename Dune::QuadratureRule<ct,dim>::const_iterator i=lowrule.begin();
i!=lowrule.end(); ++i)
{
double fval = f(it->geometry().global(i->position()));
double weight = i->weight();
double detjac = it->geometry().integrationElement(i->position());
lowresult += fval * weight * detjac;
}
// compute integral with high order rule
double highresult=0;
for (typename Dune::QuadratureRule<ct,dim>::const_iterator i=highrule.begin();
i!=highrule.end(); ++i)
{
double fval = f(it->geometry().global(i->position()));
double weight = i->weight();
double detjac = it->geometry().integrationElement(i->position());
highresult += fval * weight * detjac;
}
// compute estimate for local error
localerror = std::abs(lowresult-highresult);
// return the more accurate integral
return highresult;
}
//! compute integral of function over entity with given order
template<class Iterator, class Functor>
double integrateentity (const Iterator& it, const Functor& f, int p)
......@@ -96,5 +41,4 @@ double integrateentity (const Iterator& it, const Functor& f, int p)
// return the more accurate integral
return result;
}
#endif
......@@ -138,7 +138,9 @@ void parevolve (const G& grid, const M& mapper, V& c, double t, double& dt)
// exchange update
VectorExchange<M,V> dh(mapper,update);
grid.template communicate<VectorExchange<M,V> >(dh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
grid.template
communicate<VectorExchange<M,V> >(dh,Dune::InteriorBorder_All_Interface,
Dune::ForwardCommunication);
// update the concentration vector
for (unsigned int i=0; i<c.size(); ++i)
......
......@@ -53,7 +53,8 @@ void partimeloop (const G& grid, double tend)
k++;
parevolve(grid,mapper,c,t,dt);
t += dt;
if (grid.comm().rank()==0) std::cout << "k=" << k << " t=" << t << " dt=" << dt << std::endl;
if (grid.comm().rank()==0)
std::cout << "k=" << k << " t=" << t << " dt=" << dt << std::endl;
if (k%20==0) vtkout(grid,c,"pconc",k/20);
}
vtkout(grid,c,"pconc",k/20);
......
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