Commit df6c54ac authored by Martin Nolte's avatar Martin Nolte

adapt to geomtry being returned as object

[[Imported from SVN: r370]]
parent aaf53171
......@@ -38,7 +38,7 @@ void adaptiveintegration (Grid& grid, const Functor& f)
double value=0; /*@\label{aic:int0}@*/
for (ElementLeafIterator it = gridView.template begin<0>();
it!=gridView.template end<0>(); ++it)
value += integrateentity(it,f,highorder); /*@\label{aic:int1}@*/
value += integrateEntity(*it,f,highorder); /*@\label{aic:int1}@*/
// print result
double estimated_error = std::abs(value-oldvalue);
......@@ -71,16 +71,16 @@ void adaptiveintegration (Grid& grid, const Functor& f)
it!=grid.template leafend<0>(); ++it)
{
// error on this entity
double lowresult=integrateentity(it,f,loworder);
double highresult=integrateentity(it,f,highorder);
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 fatherlowresult=integrateentity(it->father(),f,loworder);
double fatherhighresult=integrateentity(it->father(),f,highorder);
double fatherlowresult=integrateEntity(*(it->father()),f,loworder);
double fatherhighresult=integrateEntity(*(it->father()),f,highorder);
double fathererror = std::abs(fatherlowresult-fatherhighresult);
// local extrapolation
......@@ -93,8 +93,8 @@ void adaptiveintegration (Grid& grid, const Functor& f)
for (ElementLeafIterator it = gridView.template begin<0>(); /*@\label{aic:mark0}@*/
it!=gridView.template end<0>(); ++it)
{
double lowresult=integrateentity(it,f,loworder);
double highresult=integrateentity(it,f,highorder);
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);
} /*@\label{aic:mark1}@*/
......
......@@ -48,11 +48,11 @@ void elementdata (const G& grid, const F& f)
for (ElementLeafIterator it = gridView.template begin<0>(); /*@\label{edh:loop0}@*/
it!=gridView.template end<0>(); ++it)
{
// cell geometry type
const LeafGeometry & gt = it->geometry();
// cell geometry
const LeafGeometry geo = it->geometry();
// get global coordinate of cell center
Dune::FieldVector<ct,dimworld> global = gt.center();
Dune::FieldVector<ct,dimworld> global = geo.center();
// evaluate functor and store value
c[mapper.map(*it)] = f(global); /*@\label{edh:feval}@*/
......
......@@ -47,11 +47,11 @@ void evolve (const G& grid, const M& mapper, V& c, double t, double& dt)
for (LeafIterator it = gridView.template begin<0>(); it!=endit; ++it)
{
// cell geometry
const LeafGeometry &gt = it->geometry();
const LeafGeometry geo = it->geometry();
// cell volume, assume linear map here
double volume = gt.volume();
double volume = geo.volume();
// cell index
int indexi = mapper.map(*it);
......@@ -64,15 +64,15 @@ void evolve (const G& grid, const M& mapper, V& c, double t, double& dt)
for (IntersectionIterator is = gridView.ibegin(*it); is!=isend; ++is)
{
// get geometry type of face
const IntersectionGeometry &gtf = is->geometry();
const IntersectionGeometry igeo = is->geometry();
// get normal vector scaled with volume
Dune::FieldVector<ct,dimworld> integrationOuterNormal
= is->centerUnitOuterNormal();
integrationOuterNormal *= gtf.volume();
integrationOuterNormal *= igeo.volume();
// center of face in global coordinates
Dune::FieldVector<ct,dimworld> faceglobal = gtf.center();
Dune::FieldVector<ct,dimworld> faceglobal = igeo.center();
// evaluate velocity at face center
Dune::FieldVector<double,dimworld> velocity = u(faceglobal,t);
......@@ -96,8 +96,8 @@ void evolve (const G& grid, const M& mapper, V& c, double t, double& dt)
(it->level()==outside->level() && indexi<indexj) )
{
// compute factor in neighbor
const LeafGeometry &nbgt = outside->geometry();
double nbvolume = nbgt.volume();
const LeafGeometry nbgeo = outside->geometry();
double nbvolume = nbgeo.volume();
double nbfactor = velocity*integrationOuterNormal/nbvolume;
if (factor<0) // inflow
......
......@@ -33,10 +33,10 @@ void initialize (const G& grid, const M& mapper, V& c)
for (LeafIterator it = gridView.template begin<0>(); it!=endit; ++it)
{
// get geometry
const Geometry &gt = it->geometry();
const Geometry geo = it->geometry();
// get global coordinate of cell center
Dune::FieldVector<ct,dimworld> global = gt.center();
Dune::FieldVector<ct,dimworld> global = geo.center();
// initialize cell concentration
c[mapper.map(*it)] = c0(global);
......
......@@ -7,21 +7,24 @@
#include <dune/geometry/quadraturerules.hh>
//! compute integral of function over entity with given order
template<class Iterator, class Functor>
double integrateentity (const Iterator& it, const Functor& f, int p)
template<class Entity, class Function>
double integrateEntity (const Entity &entity, const Function &f, int p)
{
// dimension of the entity
const int dim = Iterator::Entity::dimension;
const int dim = Entity::dimension;
// type used for coordinates in the grid
typedef typename Iterator::Entity::ctype ct;
typedef typename Entity::ctype ctype;
// get geometry
const typename Entity::Geometry geometry = entity.geometry();
// get geometry type
Dune::GeometryType gt = it->type();
const Dune::GeometryType gt = geometry.type();
// get quadrature rule of order p
const Dune::QuadratureRule<ct,dim>&
rule = Dune::QuadratureRules<ct,dim>::rule(gt,p); /*@\label{ieh:qr}@*/
const Dune::QuadratureRule<ctype,dim>&
rule = Dune::QuadratureRules<ctype,dim>::rule(gt,p); /*@\label{ieh:qr}@*/
// ensure that rule has at least the requested order
if (rule.order()<p)
......@@ -29,16 +32,17 @@ double integrateentity (const Iterator& it, const Functor& f, int p)
// compute approximate integral
double result=0;
for (typename Dune::QuadratureRule<ct,dim>::const_iterator i=rule.begin(); /*@\label{ieh:for}@*/
for (typename Dune::QuadratureRule<ctype,dim>::const_iterator i=rule.begin(); /*@\label{ieh:for}@*/
i!=rule.end(); ++i)
{
double fval = f(it->geometry().global(i->position())); /*@\label{ieh:fval}@*/
double fval = f(geometry.global(i->position())); /*@\label{ieh:fval}@*/
double weight = i->weight(); /*@\label{ieh:weight}@*/
double detjac = it->geometry().integrationElement(i->position()); /*@\label{ieh:detjac}@*/
double detjac = geometry.integrationElement(i->position()); /*@\label{ieh:detjac}@*/
result += fval * weight * detjac; /*@\label{ieh:result}@*/
}
// return result
return result;
}
#endif
......@@ -34,7 +34,7 @@ void uniformintegration (Grid& grid)
double value = 0.0;
LeafIterator eendit = gridView.template end<0>();
for (LeafIterator it = gridView.template begin<0>(); it!=eendit; ++it)
value += integrateentity(it,f,1); /*@\label{ic:call}@*/
value += integrateEntity(*it,f,1); /*@\label{ic:call}@*/
// print result and error estimate
std::cout << "elements="
......
......@@ -59,10 +59,10 @@ void parevolve (const G& grid, const M& mapper, V& c, double t, double& dt)
for (LeafIterator it = gridView.template begin<0,Dune::All_Partition>(); it!=endit; ++it) /*@\label{peh:begin}@*/
{
// cell geometry
const LeafGeometry & gt = it->geometry();
const LeafGeometry geo = it->geometry();
// cell volume
double volume = gt.volume();
double volume = geo.volume();
// cell index
int indexi = mapper.map(*it);
......@@ -77,15 +77,15 @@ void parevolve (const G& grid, const M& mapper, V& c, double t, double& dt)
const Intersection &intersection = *is;
// get geometry type of face
const IntersectionGeometry & gtf = intersection.geometry();
const IntersectionGeometry igeo = intersection.geometry();
// get normal vector scaled with volume
Dune::FieldVector< ct, dimworld > integrationOuterNormal
= intersection.centerUnitOuterNormal();
integrationOuterNormal *= gtf.volume();
integrationOuterNormal *= igeo.volume();
// center of face in global coordinates
Dune::FieldVector< ct, dimworld > faceglobal = gtf.center();
Dune::FieldVector< ct, dimworld > faceglobal = igeo.center();
// evaluate velocity at face center
Dune::FieldVector<double,dim> velocity = u(faceglobal,t);
......@@ -111,8 +111,8 @@ void parevolve (const G& grid, const M& mapper, V& c, double t, double& dt)
|| ((insideLevel == outsideLevel) && (indexi < indexj)) )
{
// compute factor in neighbor
const LeafGeometry & nbgt = outside->geometry();
double nbvolume = nbgt.volume();
const LeafGeometry nbgeo = outside->geometry();
double nbvolume = nbgeo.volume();
double nbfactor = velocity*integrationOuterNormal/nbvolume;
if( factor < 0 ) // inflow
......
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