Commit 4de2d4db authored by Christian Engwer's avatar Christian Engwer

The IntersectionIterators are now extracted from

 GridViews instead of the grid directly.

creadits to Martin Drohmann

[[Imported from SVN: r238]]
parent a2baaad5
......@@ -21,8 +21,13 @@
template<class Grid, class Functor>
void adaptiveintegration (Grid& grid, const Functor& f)
{
// get grid view type for leaf grid part
typedef typename Grid::LeafGridView GridView;
// get iterator type
typedef typename Grid::template Codim<0>::LeafIterator ElementLeafIterator;
typedef typename GridView::template Codim<0>::Iterator ElementLeafIterator;
// get grid view on leaf part
GridView gridView = grid.leafView();
// algorithm parameters
const double tol=1E-8;
......@@ -35,8 +40,8 @@ void adaptiveintegration (Grid& grid, const Functor& f)
{
// compute integral on current mesh
double value=0; /*@\label{aic:int0}@*/
for (ElementLeafIterator it = grid.template leafbegin<0>();
it!=grid.template leafend<0>(); ++it)
for (ElementLeafIterator it = gridView.template begin<0>();
it!=gridView.template end<0>(); ++it)
value += integrateentity(it,f,highorder); /*@\label{aic:int1}@*/
// print result
......@@ -89,8 +94,8 @@ void adaptiveintegration (Grid& grid, const Functor& f)
double kappa = std::min(maxextrapolatederror,0.5*maxerror); /*@\label{aic:kappa1}@*/
// mark elements for refinement
for (ElementLeafIterator it = grid.template leafbegin<0>(); /*@\label{aic:mark0}@*/
it!=grid.template leafend<0>(); ++it)
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);
......@@ -105,7 +110,7 @@ void adaptiveintegration (Grid& grid, const Functor& f)
}
// write grid in VTK format
Dune::VTKWriter<typename Grid::LeafGridView> vtkwriter(grid.leafView());
Dune::VTKWriter<typename Grid::LeafGridView> vtkwriter(gridView);
vtkwriter.write("adaptivegrid",Dune::VTKOptions::binaryappended);
}
......
......@@ -26,7 +26,11 @@ void elementdata (const G& grid, const F& f)
const int dim = G::dimension;
const int dimworld = G::dimensionworld;
typedef typename G::ctype ct;
typedef typename G::template Codim<0>::LeafIterator ElementLeafIterator;
typedef typename G::LeafGridView GridView;
typedef typename GridView::template Codim<0>::Iterator ElementLeafIterator;
// get grid view on leaf part
GridView gridView = grid.leafView();
// make a mapper for codim 0 entities in the leaf grid
Dune::LeafMultipleCodimMultipleGeomTypeMapper<G,P0Layout>
......@@ -36,8 +40,8 @@ void elementdata (const G& grid, const F& f)
std::vector<double> c(mapper.size()); /*@\label{edh:c}@*/
// iterate through all entities of codim 0 at the leafs
for (ElementLeafIterator it = grid.template leafbegin<0>(); /*@\label{edh:loop0}@*/
it!=grid.template leafend<0>(); ++it)
for (ElementLeafIterator it = gridView.template begin<0>(); /*@\label{edh:loop0}@*/
it!=gridView.template end<0>(); ++it)
{
// cell geometry type
Dune::GeometryType gt = it->type();
......@@ -55,7 +59,7 @@ void elementdata (const G& grid, const F& f)
// generate a VTK file
// Dune::LeafP0Function<G,double> cc(grid,c);
Dune::VTKWriter<typename G::LeafGridView> vtkwriter(grid.leafView()); /*@\label{edh:vtk0}@*/
Dune::VTKWriter<typename G::LeafGridView> vtkwriter(gridView); /*@\label{edh:vtk0}@*/
vtkwriter.addCellData(c,"data");
vtkwriter.write("elementdata",Dune::VTKOptions::binaryappended); /*@\label{edh:vtk1}@*/
......@@ -69,7 +73,7 @@ void elementdata (const G& grid, const F& f)
// display data
grape.displayVector("concentration", // name of data that appears in grape
c, // data vector
grid.leafIndexSet(), // used index set
gridView.indexSet(), // used index set
polynomialOrder, // polynomial order of data
dimRange); // dimRange of data
}
......
......@@ -12,15 +12,21 @@ void evolve (const G& grid, const M& mapper, V& c, double t, double& dt)
// type used for coordinates in the grid
typedef typename G::ctype ct;
// iterator type
typedef typename G::template Codim<0>::LeafIterator LeafIterator;
// type of grid view on leaf part
typedef typename G::LeafGridView GridView;
// element iterator type
typedef typename GridView::template Codim<0>::Iterator LeafIterator;
// intersection iterator type
typedef typename G::template Codim<0>::LeafIntersectionIterator IntersectionIterator;
typedef typename GridView::IntersectionIterator IntersectionIterator;
// entity pointer type
typedef typename G::template Codim<0>::EntityPointer EntityPointer;
// get grid view on leaf part
GridView gridView = grid.leafView();
// allocate a temporary vector for the update
V update(c.size()); /*@\label{evh:update}@*/
for (typename V::size_type i=0; i<c.size(); i++) update[i] = 0;
......@@ -29,8 +35,8 @@ void evolve (const G& grid, const M& mapper, V& c, double t, double& dt)
dt = 1E100;
// compute update vector and optimum dt in one grid traversal
LeafIterator endit = grid.template leafend<0>(); /*@\label{evh:loop0}@*/
for (LeafIterator it = grid.template leafbegin<0>(); it!=endit; ++it)
LeafIterator endit = gridView.template end<0>(); /*@\label{evh:loop0}@*/
for (LeafIterator it = gridView.template begin<0>(); it!=endit; ++it)
{
// cell geometry type
Dune::GeometryType gt = it->type();
......@@ -54,8 +60,8 @@ void evolve (const G& grid, const M& mapper, V& c, double t, double& dt)
double sumfactor = 0.0;
// run through all intersections with neighbors and boundary
IntersectionIterator isend = it->ileafend(); /*@\label{evh:flux0}@*/
for (IntersectionIterator is = it->ileafbegin(); is!=isend; ++is)
IntersectionIterator isend = gridView.iend(*it); /*@\label{evh:flux0}@*/
for (IntersectionIterator is = gridView.ibegin(*it); is!=isend; ++is)
{
// get geometry type of face
Dune::GeometryType gtf = is->intersectionSelfLocal().type();
......
......@@ -23,27 +23,35 @@ bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax, int k)
// type used for coordinates in the grid
typedef typename G::ctype ct;
// grid view types
typedef typename G::LeafGridView LeafGridView;
typedef typename G::LevelGridView LevelGridView;
// iterator types
typedef typename G::template Codim<0>::LeafIterator LeafIterator;
typedef typename G::template Codim<0>::LevelIterator LevelIterator;
typedef typename LeafGridView::template Codim<0>::Iterator LeafIterator;
typedef typename LevelGridView::template Codim<0>::Iterator LevelIterator;
// entity and entity pointer
typedef typename G::template Codim<0>::Entity Entity;
typedef typename G::template Codim<0>::EntityPointer EntityPointer;
// intersection iterator type
typedef typename G::template Codim<0>::LeafIntersectionIterator IntersectionIterator;
typedef typename LeafGridView::IntersectionIterator LeafIntersectionIterator;
// global id set types
// global id set types, local means that the numbering is unique in a single process only.
typedef typename G::template Codim<0>::LocalIdSet IdSet;
// type for the index set, note that this is _not_ an integer
typedef typename IdSet::IdType IdType;
// get grid view on leaf grid
LeafGridView leafView = grid.leafView();
// compute cell indicators
V indicator(c.size(),-1E100);
double globalmax = -1E100;
double globalmin = 1E100;
for (LeafIterator it = grid.template leafbegin<0>(); /*@\label{fah:loop0}@*/
it!=grid.template leafend<0>(); ++it)
for (LeafIterator it = leafView.template begin<0>(); /*@\label{fah:loop0}@*/
it!=leafView.template end<0>(); ++it)
{
// my index
int indexi = mapper.map(*it);
......@@ -52,10 +60,10 @@ bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax, int k)
globalmax = std::max(globalmax,c[indexi]);
globalmin = std::min(globalmin,c[indexi]);
IntersectionIterator isend = it->ileafend();
for (IntersectionIterator is = it->ileafbegin(); is!=isend; ++is)
LeafIntersectionIterator isend = leafView.iend(*it);
for (LeafIntersectionIterator is = leafView.ibegin(*it); is!=isend; ++is)
{
const typename IntersectionIterator::Intersection &intersection = *is;
const typename LeafIntersectionIterator::Intersection &intersection = *is;
if( !intersection.neighbor() )
continue;
......@@ -78,8 +86,8 @@ bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax, int k)
// mark cells for refinement/coarsening
double globaldelta = globalmax-globalmin;
int marked=0;
for (LeafIterator it = grid.template leafbegin<0>(); /*@\label{fah:loop2}@*/
it!=grid.template leafend<0>(); ++it)
for (LeafIterator it = leafView.template begin<0>(); /*@\label{fah:loop2}@*/
it!=leafView.template end<0>(); ++it)
{
if (indicator[mapper.map(*it)]>refinetol*globaldelta
&& (it.level()<lmax || !it->isRegular()))
......@@ -87,10 +95,10 @@ bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax, int k)
const Entity &entity = *it;
grid.mark( 1, entity );
++marked;
IntersectionIterator isend = entity.ileafend();
for( IntersectionIterator is = entity.ileafbegin(); is != isend; ++is )
LeafIntersectionIterator isend = leafView.iend(entity);
for( LeafIntersectionIterator is = leafView.ibegin(entity); is != isend; ++is )
{
const typename IntersectionIterator::Intersection intersection = *is;
const typename LeafIntersectionIterator::Intersection intersection = *is;
if( !intersection.neighbor() )
continue;
......@@ -113,8 +121,11 @@ bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax, int k)
std::map<IdType,RestrictedValue> restrictionmap; // restricted concentration /*@\label{fah:loop4}@*/
const IdSet& idset = grid.localIdSet();
for (int level=grid.maxLevel(); level>=0; level--)
for (LevelIterator it = grid.template lbegin<0>(level);
it!=grid.template lend<0>(level); ++it)
{
// get grid view on level grid
LevelGridView levelView = grid.levelView(level);
for (LevelIterator it = levelView.template begin<0>();
it!=levelView.template end<0>(); ++it)
{
// get your map entry
IdType idi = idset.id(*it);
......@@ -138,6 +149,7 @@ bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax, int k)
rvf.count += 1;
}
} /*@\label{fah:loop5}@*/
}
grid.preAdapt();
// adapt mesh and mapper
......@@ -147,14 +159,17 @@ bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax, int k)
// interpolate new cells, restrict coarsened cells
for (int level=0; level<=grid.maxLevel(); level++) /*@\label{fah:loop6}@*/
for (LevelIterator it = grid.template lbegin<0>(level);
it!=grid.template lend<0>(level); ++it)
{
LevelGridView levelView = grid.levelView(level);
for (LevelIterator it = levelView.template begin<0>();
it!=levelView.template end<0>(); ++it)
{
// get your id
IdType idi = idset.id(*it);
// check map entry
typename std::map<IdType,RestrictedValue>::iterator rit = restrictionmap.find(idi);
typename std::map<IdType,RestrictedValue>::iterator rit
= restrictionmap.find(idi);
if (rit!=restrictionmap.end())
{
// entry is in map, write in leaf
......@@ -166,7 +181,7 @@ bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax, int k)
}
else
{
// value is not in map, interpolate
// value is not in map, interpolate from father element
if (it.level()>0)
{
EntityPointer ep = it->father();
......@@ -187,6 +202,7 @@ bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax, int k)
}
}
} /*@\label{fah:loop7}@*/
}
grid.postAdapt();
return rv;
......
......@@ -13,12 +13,18 @@ void initialize (const G& grid, const M& mapper, V& c)
// type used for coordinates in the grid
typedef typename G::ctype ct;
// type of grid view on leaf part
typedef typename G::LeafGridView GridView;
// leaf iterator type
typedef typename G::template Codim<0>::LeafIterator LeafIterator;
typedef typename GridView::template Codim<0>::Iterator LeafIterator;
// get grid view on leaf part
GridView gridView = grid.leafView();
// iterate through leaf grid an evaluate c0 at cell center
LeafIterator endit = grid.template leafend<0>();
for (LeafIterator it = grid.template leafbegin<0>(); it!=endit; ++it)
LeafIterator endit = gridView.template end<0>();
for (LeafIterator it = gridView.template begin<0>(); it!=endit; ++it)
{
// get geometry type
Dune::GeometryType gt = it->type();
......
......@@ -20,8 +20,14 @@ void uniformintegration (Grid& grid)
// function to integrate
Exp<typename Grid::ctype,Grid::dimension> f;
// get GridView on leaf grid - type
typedef typename Grid :: LeafGridView GridView;
// get GridView instance
GridView gridView = grid.leafView();
// get iterator type
typedef typename Grid::template Codim<0>::LeafIterator LeafIterator;
typedef typename GridView :: template Codim<0> :: Iterator LeafIterator;
// loop over grid sequence
double oldvalue=1E100;
......@@ -29,8 +35,8 @@ void uniformintegration (Grid& grid)
{
// compute integral with some order
double value = 0.0;
LeafIterator eendit = grid.template leafend<0>();
for (LeafIterator it = grid.template leafbegin<0>(); it!=eendit; ++it)
LeafIterator eendit = gridView.template end<0>();
for (LeafIterator it = gridView.template begin<0>(); it!=eendit; ++it)
value += integrateentity(it,f,1); /*@\label{ic:call}@*/
// print result and error estimate
......
......@@ -15,12 +15,15 @@ void parevolve (const G& grid, const M& mapper, V& c, double t, double& dt)
// type used for coordinates in the grid
typedef typename G::ctype ct;
// type for grid view on leaf part
typedef typename G::LeafGridView GridView;
// iterator type
typedef typename G::template Codim<0>::
template Partition<Dune::All_Partition>::LeafIterator LeafIterator; /*@\label{peh:pit}@*/
typedef typename GridView::template Codim<0>::
template Partition<Dune::All_Partition>::Iterator LeafIterator; /*@\label{peh:pit}@*/
// intersection iterator type
typedef typename G::template Codim<0>::LeafIntersectionIterator IntersectionIterator;
typedef typename GridView::IntersectionIterator IntersectionIterator;
// type of intersection
typedef typename IntersectionIterator::Intersection Intersection;
......@@ -35,10 +38,13 @@ void parevolve (const G& grid, const M& mapper, V& c, double t, double& dt)
// initialize dt very large
dt = 1E100;
// get grid view instance on leaf grid
GridView gridView = grid.leafView();
// compute update vector and optimum dt in one grid traversal
// iterate over all entities, but update is only used on interior entities
LeafIterator endit = grid.template leafend<0,Dune::All_Partition>(); /*@\label{peh:end}@*/
for (LeafIterator it = grid.template leafbegin<0,Dune::All_Partition>(); it!=endit; ++it) /*@\label{peh:begin}@*/
LeafIterator endit = gridView.template end<0,Dune::All_Partition>(); /*@\label{peh:end}@*/
for (LeafIterator it = gridView.template begin<0,Dune::All_Partition>(); it!=endit; ++it) /*@\label{peh:begin}@*/
{
// cell geometry type
Dune::GeometryType gt = it->type();
......@@ -62,8 +68,8 @@ void parevolve (const G& grid, const M& mapper, V& c, double t, double& dt)
double sumfactor = 0.0;
// run through all intersections with neighbors and boundary
const IntersectionIterator isend = it->ileafend();
for( IntersectionIterator is = it->ileafbegin(); is != isend; ++is )
const IntersectionIterator isend = gridView.iend(*it);
for( IntersectionIterator is = gridView.ibegin(*it); is != isend; ++is )
{
const Intersection &intersection = *is;
......
......@@ -26,20 +26,25 @@ void traversal (G& grid)
// Leaf Traversal
std::cout << "*** Traverse codim 0 leaves" << std::endl;
// the grid has an iterator providing the access to
// all elements (better codim 0 entities) which are leafs
// of the refinement tree.
// Note the use of the typename keyword and the traits class
typedef typename G::template Codim<0>::LeafIterator ElementLeafIterator; /*@\label{tc:ittype}@*/
// type of the GridView used for traversal
// every grid exports a LeafGridView and a LevelGridView
typedef typename G :: LeafGridView LeafGridView; /*@\label{tc:lfgv}@*/
// get the instance of the LeafGridView
LeafGridView leafView = grid.leafView(); /*@\label{tc:lfv}@*/
// Get the iterator type
// Note the use of the typename and template keywords
typedef typename LeafGridView::template Codim<0>::Iterator ElementLeafIterator; /*@\label{tc:ittype}@*/
// iterate through all entities of codim 0 at the leafs
int count = 0;
for (ElementLeafIterator it = grid.template leafbegin<0>(); /*@\label{tc:forel}@*/
it!=grid.template leafend<0>(); ++it)
for (ElementLeafIterator it = leafView.template begin<0>(); /*@\label{tc:forel}@*/
it!=leafView.template end<0>(); ++it)
{ /*@\label{tc:forel0}@*/
Dune::GeometryType gt = it->type(); /*@\label{tc:reftype}@*/
std::cout << "visiting leaf " << gt /*@\label{tc:print}@*/
<< " with first vertex at " << it->geometry()[0]
<< " with first vertex at " << it->geometry().corner(0)
<< std::endl;
count++; /*@\label{tc:count}@*/
} /*@\label{tc:forel1}@*/
......@@ -52,16 +57,17 @@ void traversal (G& grid)
// Get the iterator type
// Note the use of the typename and template keywords
typedef typename G::template Codim<dim>::LeafIterator VertexLeafIterator; /*@\label{tc:vertit}@*/
typedef typename LeafGridView :: template Codim<dim>
:: Iterator VertexLeafIterator; /*@\label{tc:vertit}@*/
// iterate through all entities of codim 0 on the given level
count = 0;
for (VertexLeafIterator it = grid.template leafbegin<dim>(); /*@\label{tc:forve}@*/
it!=grid.template leafend<dim>(); ++it)
for (VertexLeafIterator it = leafView.template begin<dim>(); /*@\label{tc:forve}@*/
it!=leafView.template end<dim>(); ++it)
{
Dune::GeometryType gt = it->type();
std::cout << "visiting " << gt
<< " at " << it->geometry()[0]
<< " at " << it->geometry().corner(0)
<< std::endl;
count++;
}
......@@ -72,20 +78,28 @@ void traversal (G& grid)
std::cout << std::endl;
std::cout << "*** Traverse codim 0 level-wise" << std::endl;
// type of the GridView used for traversal
// every grid exports a LeafGridView and a LevelGridView
typedef typename G :: LevelGridView LevelGridView; /*@\label{tc:level0}@*/
// Get the iterator type
// Note the use of the typename and template keywords
typedef typename G::template Codim<0>::LevelIterator ElementLevelIterator; /*@\label{tc:level0}@*/
typedef typename LevelGridView :: template Codim<0>
:: Iterator ElementLevelIterator;
// iterate through all entities of codim 0 on the given level
for (int level=0; level<=grid.maxLevel(); level++)
{
// get the instance of the LeafGridView
LevelGridView levelView = grid.levelView(level);
count = 0;
for (ElementLevelIterator it = grid.template lbegin<0>(level);
it!=grid.template lend<0>(level); ++it)
for (ElementLevelIterator it = levelView.template begin<0>();
it!=levelView.template end<0>(); ++it)
{
Dune::GeometryType gt = it->type();
std::cout << "visiting " << gt
<< " with first vertex at " << it->geometry()[0]
<< " with first vertex at " << it->geometry().corner(0)
<< std::endl;
count++;
}
......
......@@ -25,8 +25,12 @@ void vertexdata (const G& grid, const F& f)
// get dimension and coordinate type from Grid
const int dim = G::dimension;
typedef typename G::ctype ct;
typedef typename G::LeafGridView GridView;
// dertermine type of LeafIterator for codimension = dimension
typedef typename G::template Codim<dim>::LeafIterator VertexLeafIterator;
typedef typename GridView::template Codim<dim>::Iterator VertexLeafIterator;
// get grid view on the leaf part
GridView gridView = grid.leafView();
// make a mapper for codim 0 entities in the leaf grid
Dune::LeafMultipleCodimMultipleGeomTypeMapper<G,P1Layout>
......@@ -36,8 +40,8 @@ void vertexdata (const G& grid, const F& f)
std::vector<double> c(mapper.size());
// iterate through all entities of codim 0 at the leafs
for (VertexLeafIterator it = grid.template leafbegin<dim>();
it!=grid.template leafend<dim>(); ++it)
for (VertexLeafIterator it = gridView.template begin<dim>();
it!=gridView.template end<dim>(); ++it)
{
// evaluate functor and store value
c[mapper.map(*it)] = f(it->geometry().corner(0));
......@@ -59,7 +63,7 @@ void vertexdata (const G& grid, const F& f)
// display data
grape.displayVector("concentration", // name of data that appears in grape
c, // data vector
grid.leafIndexSet(), // used index set
gridView.indexSet(), // used index set
polynomialOrder, // polynomial order of data
dimRange); // dimRange of data
}
......
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