[Dune] dune-grid-howto broken
Christian Engwer
christi at uni-hd.de
Fri Mar 6 18:12:17 CET 2009
On Fri, Mar 06, 2009 at 06:08:43PM +0100, Martin Drohmann wrote:
> Hi Chrisitan,
>
> as far as i know, parallelization on a yaspgrid only works when the initial
> grid is coarse enough. With the UnitCube it is quite easy to change the
> coarseness for a YaspGrid cube through the second template parameter. That
> is why I left it in the code.
>
> But of course, you could also change the dgf file. So if you think this part
> is superfluous, just delete it.
Well, in this case it sounds more like the DGFparser isn't really
usefull in this situtation.
Christian
>
> Martin
>
> 2009/3/6 Christian Engwer <christi at uni-hd.de>
>
> > Hi Martin,
> >
> > On Fri, Mar 06, 2009 at 05:42:53PM +0100, Martin Drohmann wrote:
> > > Hello Christian,
> > >
> > > sorry for the whitespace changes. I built a new patch and excluded all
> > lines
> > > with "whitespace-only" changes.
> > >
> > > Additionally I resorted the big patch, made small ones out of it and
> > > documented them individually. These patches are in the attached
> > > small_patches.tgz. In the header of each patch you will find a
> > description
> > > of what I changed. I hope this helps.
> >
> > On mall question is left concerning
> > 0005-updated-the-parfinitevolume-example-now-using-dgf-f.patch
> >
> > Is there a reason for having both version, UnitCube and DGFParser in
> > the example? Or can I just remove the UnitCube part?
> >
> > Christian
> >
> > >
> > > Bye,
> > > Martin
> > >
> > > 2009/3/6 Christian Engwer <christi at uni-hd.de>
> > >
> > > > Dear Martin,
> > > >
> > > > thank you for your work on the dune-grid-howto. The changes in the tex
> > > > file include some references to a basicunitcube.hh. The code you sent
> > > > include a lot of white space changes and it is not obvious what all
> > > > the changes are.
> > > >
> > > > Please update your patches to avoid whitespace changes and explain a
> > > > little bit more about the source-code changes.
> > > >
> > > > In the meantime I'll remove the paragraph quoting basicunitcube, in
> > > > order to make the howto work again.
> > > >
> > > > Thank you
> > > > Christian
> > > >
> > > >
> > >
> >
> > > diff --git a/adaptivefinitevolume.cc b/adaptivefinitevolume.cc
> > > index a9f2241..e1425d4 100644
> > > --- a/adaptivefinitevolume.cc
> > > +++ b/adaptivefinitevolume.cc
> > > @@ -32,43 +32,6 @@ struct P0Layout
> > > };
> > >
> > > template<class G>
> > > -void gnuplot (G& grid, std::vector<double>& c)
> > > -{
> > > - // first we extract the dimensions of the grid
> > > - const int dim = G::dimension;
> > > - const int dimworld = G::dimensionworld;
> > > -
> > > - // type used for coordinates in the grid
> > > - // such a type is exported by every grid implementation
> > > - typedef typename G::ctype ct;
> > > -
> > > - // 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;
> > > -
> > > - // make a mapper for codim 0 entities in the leaf grid
> > > - Dune::LeafMultipleCodimMultipleGeomTypeMapper<G,P0Layout>
> > > - mapper(grid);
> > > -
> > > - // iterate through all entities of codim 0 at the leafs
> > > - int count = 0;
> > > - for (ElementLeafIterator it = grid.template leafbegin<0>();
> > > - it!=grid.template leafend<0>(); ++it)
> > > - {
> > > - Dune::GeometryType gt = it->type();
> > > - const Dune::FieldVector<ct,dim>&
> > > - local =
> > Dune::ReferenceElements<ct,dim>::general(gt).position(0,0);
> > > - Dune::FieldVector<ct,dimworld>
> > > - global = it->geometry().global(local);
> > > - std::cout << global[0] << " " << c[mapper.map(*it)] << std::endl;
> > > - count++;
> > > - }
> > > -}
> > > -
> > > -
> > > -template<class G>
> > > void timeloop (G& grid, double tend, int lmin, int lmax)
> > > {
> > > // make a mapper for codim 0 entities in the leaf grid
> > > @@ -83,12 +46,12 @@ void timeloop (G& grid, double tend, int lmin, int
> > lmax)
> > > for (int i=grid.maxLevel(); i<lmax; i++)
> > > {
> > > if (grid.maxLevel()>=lmax) break;
> > > - finitevolumeadapt(grid,mapper,c,lmin,lmax,0);
> > > + finitevolumeadapt(grid,mapper,c,lmin,lmax,0);
> > /*@\label{afv:in}@*/
> > > initialize(grid,mapper,c);
> > > }
> > >
> > > // write initial data
> > > - vtkout(grid,c,"concentration",0);
> > > + vtkout(grid,c,"concentration",0,0);
> > >
> > > // variables for time, timestep etc.
> > > double dt, t=0;
> > > @@ -113,7 +76,7 @@ void timeloop (G& grid, double tend, int lmin, int
> > lmax)
> > > if (t >= saveStep)
> > > {
> > > // write data
> > > - vtkout(grid,c,"concentration",counter);
> > > + vtkout(grid,c,"concentration",counter,t);
> > >
> > > // increase counter and saveStep for next interval
> > > saveStep += saveInterval;
> > > @@ -121,17 +84,20 @@ void timeloop (G& grid, double tend, int lmin, int
> > lmax)
> > > }
> > >
> > > // print info about time, timestep size and counter
> > > - std::cout << "s=" << grid.size(0) << " k=" << k << " t=" << t << "
> > dt=" << dt << std::endl;
> > > + std::cout << "s=" << grid.size(0)
> > > + << " k=" << k << " t=" << t << " dt=" << dt << std::endl;
> > >
> > > // for unstructured grids call adaptation algorithm
> > > if( Dune :: Capabilities :: IsUnstructured<G> :: v )
> > > {
> > > - finitevolumeadapt(grid,mapper,c,lmin,lmax,k);
> > > + finitevolumeadapt(grid,mapper,c,lmin,lmax,k);
> > /*@\label{afv:ad}@*/
> > > }
> > > }
> > >
> > > // write last time step
> > > - vtkout(grid,c,"concentration",counter);
> > > + vtkout(grid,c,"concentration",counter,tend);
> > > +
> > > + // write
> > > }
> > >
> > > //===============================================================
> > > diff --git a/adaptiveintegration.cc b/adaptiveintegration.cc
> > > index ded8d9c..5396c44 100644
> > > --- a/adaptiveintegration.cc
> > > +++ b/adaptiveintegration.cc
> > > @@ -19,8 +19,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;
> > > @@ -33,8 +38,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
> > > @@ -87,8 +92,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);
> > > @@ -103,7 +108,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);
> > > }
> > >
> > > diff --git a/basicunitcube.hh b/basicunitcube.hh
> > > new file mode 100644
> > > index 0000000..65a30dd
> > > --- /dev/null
> > > +++ b/basicunitcube.hh
> > > @@ -0,0 +1,111 @@
> > > +#ifndef BASICUNITCUBE_HH
> > > +#define BASICUNITCUBE_HH
> > > +
> > > +#include <dune/grid/common/gridfactory.hh>
> > > +
> > > +// declaration of a basic unit cube that uses the GridFactory
> > > +template< int dim >
> > > +class BasicUnitCube;
> > > +
> > > +// unit cube in two dimensions with 2 variants: triangle and rectangle
> > elements
> > > +template<>
> > > +class BasicUnitCube< 2 >
> > > +{
> > > +protected:
> > > + template< class Grid >
> > > + static void insertVertices ( Dune::GridFactory< Grid > &factory )
> > > + {
> > > + Dune::FieldVector<double,2> pos;
> > > +
> > > + pos[0] = 0; pos[1] = 0;
> > > + factory.insertVertex(pos); /*@\label{uc:iv}@
> > */
> > > +
> > > + pos[0] = 1; pos[1] = 0;
> > > + factory.insertVertex(pos);
> > > +
> > > + pos[0] = 0; pos[1] = 1;
> > > + factory.insertVertex(pos);
> > > +
> > > + pos[0] = 1; pos[1] = 1;
> > > + factory.insertVertex(pos);
> > > + }
> > > +
> > > + template< class Grid >
> > > + static void insertSimplices ( Dune::GridFactory< Grid > &factory )
> > > + {
> > > + const Dune::GeometryType type( Dune::GeometryType::simplex, 2 );
> > > + std::vector< unsigned int > cornerIDs( 3 );
> > > +
> > > + cornerIDs[0] = 0; cornerIDs[1] = 1; cornerIDs[2] = 2;
> > > + factory.insertElement( type, cornerIDs ); /*@\label{uc:ie}@
> > */
> > > +
> > > + cornerIDs[0] = 2; cornerIDs[1] = 1; cornerIDs[2] = 3;
> > > + factory.insertElement( type, cornerIDs );
> > > + }
> > > +
> > > + template< class Grid >
> > > + static void insertCubes ( Dune::GridFactory< Grid > &factory )
> > > + {
> > > + const Dune::GeometryType type( Dune::GeometryType::cube, 2 );
> > > + std::vector< unsigned int > cornerIDs( 4 );
> > > + for( int i = 0; i < 4; ++i )
> > > + cornerIDs[ i ] = i;
> > > + factory.insertElement( type, cornerIDs );
> > > + }
> > > +};
> > > +
> > > +// unit cube in 3 dimensions with two variants: tetraheda and hexahedra
> > > +template<>
> > > +class BasicUnitCube< 3 >
> > > +{
> > > +protected:
> > > + template< class Grid >
> > > + static void insertVertices ( Dune::GridFactory< Grid > &factory )
> > > + {
> > > + Dune::FieldVector< double, 3 > pos;
> > > +
> > > + pos[0] = 0; pos[1] = 0; pos[2] = 0; factory.insertVertex(pos);
> > > + pos[0] = 1; pos[1] = 0; pos[2] = 0; factory.insertVertex(pos);
> > > + pos[0] = 0; pos[1] = 1; pos[2] = 0; factory.insertVertex(pos);
> > > + pos[0] = 1; pos[1] = 1; pos[2] = 0; factory.insertVertex(pos);
> > > + pos[0] = 0; pos[1] = 0; pos[2] = 1; factory.insertVertex(pos);
> > > + pos[0] = 1; pos[1] = 0; pos[2] = 1; factory.insertVertex(pos);
> > > + pos[0] = 0; pos[1] = 1; pos[2] = 1; factory.insertVertex(pos);
> > > + pos[0] = 1; pos[1] = 1; pos[2] = 1; factory.insertVertex(pos);
> > > + }
> > > +
> > > + template< class Grid >
> > > + static void insertSimplices ( Dune::GridFactory< Grid > &factory )
> > > + {
> > > + const Dune::GeometryType type( Dune::GeometryType::simplex, 3 );
> > > + std::vector< unsigned int > cornerIDs( 4 );
> > > +
> > > + cornerIDs[0] = 0; cornerIDs[1] = 1; cornerIDs[2] = 2;
> > cornerIDs[3] = 4;
> > > + factory.insertElement( type, cornerIDs );
> > > +
> > > + cornerIDs[0] = 1; cornerIDs[1] = 3; cornerIDs[2] = 2;
> > cornerIDs[3] = 7;
> > > + factory.insertElement( type, cornerIDs );
> > > +
> > > + cornerIDs[0] = 1; cornerIDs[1] = 7; cornerIDs[2] = 2;
> > cornerIDs[3] = 4;
> > > + factory.insertElement( type, cornerIDs );
> > > +
> > > + cornerIDs[0] = 1; cornerIDs[1] = 7; cornerIDs[2] = 4;
> > cornerIDs[3] = 5;
> > > + factory.insertElement( type, cornerIDs );
> > > +
> > > + cornerIDs[0] = 4; cornerIDs[1] = 7; cornerIDs[2] = 2;
> > cornerIDs[3] = 6;
> > > + factory.insertElement( type, cornerIDs );
> > > + }
> > > +
> > > + template< class Grid >
> > > + static void insertCubes ( Dune::GridFactory< Grid > &factory )
> > > + {
> > > + const Dune::GeometryType type( Dune::GeometryType::cube, 3 );
> > > + std::vector< unsigned int > cornerIDs( 8 );
> > > + for( int i = 0; i < 8; ++i )
> > > + cornerIDs[ i ] = i;
> > > + factory.insertElement( type, cornerIDs );
> > > + }
> > > +};
> > > +
> > > +#endif /*BASICUNITCUBE_HH*/
> > > +
> > > diff --git a/dune.module b/dune.module
> > > index ef05236..c5e2876 100644
> > > --- a/dune.module
> > > +++ b/dune.module
> > > @@ -1,5 +1,5 @@
> > > # parameters for dune control
> > > Module: dune-grid-howto
> > > Maintainer: dune at dune-project.org
> > > -Version: 1.2svn
> > > +Version: 1.2
> > > Depends: dune-common dune-grid
> > > diff --git a/elementdata.hh b/elementdata.hh
> > > index 1696aa9..a4c8b88 100644
> > > --- a/elementdata.hh
> > > +++ b/elementdata.hh
> > > @@ -24,7 +24,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>
> > > @@ -34,8 +38,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();
> > > @@ -53,7 +57,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}@*/
> > >
> > > @@ -67,7 +71,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
> > > }
> > > diff --git a/evolve.hh b/evolve.hh
> > > index 1b4262a..af0aac4 100644
> > > --- a/evolve.hh
> > > +++ b/evolve.hh
> > > @@ -10,15 +10,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;
> > > @@ -27,8 +33,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();
> > > @@ -52,8 +58,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();
> > > diff --git a/finitevolume.cc b/finitevolume.cc
> > > index 315c72e..be9ab43 100644
> > > --- a/finitevolume.cc
> > > +++ b/finitevolume.cc
> > > @@ -41,14 +41,14 @@ void timeloop (const G& grid, double tend)
> > >
> > > // initialize concentration with initial values
> > > initialize(grid,mapper,c);
> > /*@\label{fvc:init}@*/
> > > - vtkout(grid,c,"concentration",0);
> > > + vtkout(grid,c,"concentration",0,0.0);
> > >
> > > // now do the time steps
> > > double t=0,dt;
> > > int k=0;
> > > const double saveInterval = 0.1;
> > > double saveStep = 0.1;
> > > - int counter = 0;
> > > + int counter = 1;
> > >
> > > while (t<tend)
> > /*@\label{fvc:loop0}@*/
> > > {
> > > @@ -65,7 +65,7 @@ void timeloop (const G& grid, double tend)
> > > if (t >= saveStep)
> > > {
> > > // write data
> > > - vtkout(grid,c,"concentration",counter);
> > > + vtkout(grid,c,"concentration",counter,t);
> > >
> > > // increase counter and saveStep for next interval
> > > saveStep += saveInterval;
> > > @@ -73,11 +73,12 @@ void timeloop (const G& grid, double tend)
> > > }
> > >
> > > // print info about time, timestep size and counter
> > > - std::cout << "s=" << grid.size(0) << " k=" << k << " t=" << t << "
> > dt=" << dt << std::endl;
> > > + std::cout << "s=" << grid.size(0)
> > > + << " k=" << k << " t=" << t << " dt=" << dt << std::endl;
> > > }
> > /*@\label{fvc:loop1}@*/
> > >
> > > // output results
> > > - vtkout(grid,c,"concentration",counter);
> > /*@\label{fvc:file}@*/
> > > + vtkout(grid,c,"concentration",counter,tend); /*@\label{fvc:file}@
> > */
> > > }
> > >
> > > //===============================================================
> > > @@ -93,7 +94,7 @@ int main (int argc , char ** argv)
> > > try {
> > > using namespace Dune;
> > >
> > > - // use unitcube from grids
> > > + // use unitcube from dgf grids
> > > std::stringstream dgfFileName;
> > > dgfFileName << "grids/unitcube" << GridType :: dimension << ".dgf";
> > >
> > > diff --git a/finitevolumeadapt.hh b/finitevolumeadapt.hh
> > > index 42b0f35..dd79b49 100644
> > > --- a/finitevolumeadapt.hh
> > > +++ b/finitevolumeadapt.hh
> > > @@ -21,27 +21,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);
> > > @@ -50,10 +58,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;
> > >
> > > @@ -76,8 +84,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()))
> > > @@ -85,10 +93,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;
> > >
> > > @@ -111,8 +119,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);
> > > @@ -136,6 +147,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
> > > @@ -145,14 +157,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
> > > @@ -164,7 +179,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();
> > > @@ -185,6 +200,7 @@ bool finitevolumeadapt (G& grid, M& mapper, V& c, int
> > lmin, int lmax, int k)
> > > }
> > > }
> > > }
> > /*@\label{fah:loop7}@*/
> > > + }
> > > grid.postAdapt();
> > >
> > > return rv;
> > > diff --git a/grids/3dgrid.al b/grids/3dgrid.al
> > > deleted file mode 100644
> > > index 2bde72f..0000000
> > > --- a/grids/3dgrid.al
> > > +++ /dev/null
> > > @@ -1,39 +0,0 @@
> > > -DIM: 3
> > > -DIM_OF_WORLD: 3
> > > -
> > > -number of vertices: 8
> > > -number of elements: 6
> > > -
> > > -vertex coordinates:
> > > - 0.0 0.0 0.0
> > > - 1.0 0.0 0.0
> > > - 0.0 0.0 1.0
> > > - 1.0 0.0 1.0
> > > - 1.0 1.0 0.0
> > > - 1.0 1.0 1.0
> > > - 0.0 1.0 0.0
> > > - 0.0 1.0 1.0
> > > -
> > > -element vertices:
> > > - 0 5 4 1
> > > - 0 5 3 1
> > > - 0 5 3 2
> > > - 0 5 4 6
> > > - 0 5 7 6
> > > - 0 5 7 2
> > > -
> > > -element boundaries:
> > > - 1 1 0 0
> > > - 1 1 0 0
> > > - 1 1 0 0
> > > - 1 1 0 0
> > > - 1 1 0 0
> > > - 1 1 0 0
> > > -
> > > -element neighbours:
> > > - -1 -1 1 3
> > > - -1 -1 0 2
> > > - -1 -1 5 1
> > > - -1 -1 4 0
> > > - -1 -1 3 5
> > > - -1 -1 2 4
> > > diff --git a/grids/Makefile.am b/grids/Makefile.am
> > > index 6527044..8e7c6a0 100644
> > > --- a/grids/Makefile.am
> > > +++ b/grids/Makefile.am
> > > @@ -1,4 +1,4 @@
> > > -EXTRA_DIST = 2dgrid.al 2dsimplex.alu 3dgrid.al cube.hexa cube.tetra \
> > > +EXTRA_DIST = 2dgrid.al 2dsimplex.alu \
> > > unitcube1.dgf unitcube2.dgf unitcube3.dgf
> > >
> > > griddir = $(datadir)/doc/dune-grid-howto/examples/grids
> > > diff --git a/grids/cube.hexa b/grids/cube.hexa
> > > deleted file mode 100644
> > > index 2085520..0000000
> > > --- a/grids/cube.hexa
> > > +++ /dev/null
> > > @@ -1,31 +0,0 @@
> > > -!Hexahedra
> > > -
> > > -8
> > > -0.000000 0.000000 0.000000
> > > -1.000000 0.000000 0.000000
> > > -1.000000 1.000000 0.000000
> > > -0.000000 1.000000 0.000000
> > > -0.000000 0.000000 1.000000
> > > -1.000000 0.000000 1.000000
> > > -1.000000 1.000000 1.000000
> > > -0.000000 1.000000 1.000000
> > > -
> > > -1
> > > -0 1 2 3 4 5 6 7
> > > -
> > > -6
> > > --2 4 0 3 7 4
> > > --3 4 1 5 6 2
> > > --1 4 0 4 5 1
> > > --1 4 3 2 6 7
> > > --1 4 0 1 2 3
> > > --1 4 5 4 7 6
> > > -
> > > -0 -1
> > > -1 -1
> > > -2 -1
> > > -3 -1
> > > -4 -1
> > > -5 -1
> > > -6 -1
> > > -7 -1
> > > diff --git a/grids/cube.tetra b/grids/cube.tetra
> > > deleted file mode 100644
> > > index c9eeb98..0000000
> > > --- a/grids/cube.tetra
> > > +++ /dev/null
> > > @@ -1,41 +0,0 @@
> > > -!Tetrahedra
> > > -8
> > > -0.000000 0.000000 0.000000
> > > -1.000000 0.000000 0.000000
> > > -0.000000 1.000000 1.000000
> > > -0.000000 0.000000 1.000000
> > > -0.000000 1.000000 0.000000
> > > -1.000000 0.000000 1.000000
> > > -1.000000 1.000000 0.000000
> > > -1.000000 1.000000 1.000000
> > > -
> > > -6
> > > -0 1 2 3
> > > -0 2 1 4
> > > -1 5 2 3
> > > -1 6 4 2
> > > -1 6 2 7
> > > -1 7 2 5
> > > -
> > > -12
> > > --2 3 3 2 0
> > > --2 3 1 3 0
> > > --1 3 4 1 0
> > > --2 3 2 4 0
> > > --1 3 2 3 5
> > > --2 3 5 3 1
> > > --2 3 4 2 6
> > > --1 3 4 6 1
> > > --2 3 2 7 6
> > > --2 3 6 7 1
> > > --1 3 2 5 7
> > > --2 3 7 5 1
> > > -
> > > -0 -1
> > > -1 -1
> > > -2 -1
> > > -3 -1
> > > -4 -1
> > > -5 -1
> > > -6 -1
> > > -7 -1
> > > diff --git a/initialize.hh b/initialize.hh
> > > index cc4872a..96bbbd2 100644
> > > --- a/initialize.hh
> > > +++ b/initialize.hh
> > > @@ -11,12 +11,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();
> > > diff --git a/integration.cc b/integration.cc
> > > index 4bbab9b..8de5237 100644
> > > --- a/integration.cc
> > > +++ b/integration.cc
> > > @@ -18,8 +18,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;
> > > @@ -27,8 +33,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
> > > diff --git a/parevolve.hh b/parevolve.hh
> > > index c2c0f96..6d724c9 100644
> > > --- a/parevolve.hh
> > > +++ b/parevolve.hh
> > > @@ -13,12 +13,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;
> > > @@ -33,10 +36,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();
> > > @@ -60,8 +66,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;
> > >
> > > diff --git a/parfinitevolume.cc b/parfinitevolume.cc
> > > index 9f1f2f2..cc801de 100644
> > > --- a/parfinitevolume.cc
> > > +++ b/parfinitevolume.cc
> > > @@ -5,9 +5,13 @@
> > > #include<dune/grid/common/mcmgmapper.hh> // mapper class
> > > #include <dune/common/mpihelper.hh> // include mpi helper class
> > >
> > > +
> > > +// checks for defined gridtype and inlcudes appropriate dgfparser
> > implementation
> > > +#include <dune/grid/io/file/dgfparser/dgfgridtype.hh>
> > > +
> > > #include"vtkout.hh"
> > > #include"unitcube.hh"
> > > -#include"transportproblem.hh"
> > > +#include "transportproblem2.hh"
> > > #include"initialize.hh"
> > > #include"parfvdatahandle.hh"
> > > #include"parevolve.hh"
> > > @@ -40,21 +44,41 @@ void partimeloop (const G& grid, double tend)
> > >
> > > // initialize concentration with initial values
> > > initialize(grid,mapper,c);
> > > - vtkout(grid,c,"pconc",0);
> > > + vtkout(grid,c,"pconc",0,0.0,grid.comm().rank());
> > >
> > > // now do the time steps
> > > double t=0,dt;
> > > int k=0;
> > > + const double saveInterval = 0.1;
> > > + double saveStep = 0.1;
> > > + int counter = 1;
> > > while (t<tend)
> > > {
> > > + // augment time step counter
> > > k++;
> > > +
> > > + // apply finite volume scheme
> > > parevolve(grid,mapper,c,t,dt);
> > > +
> > > + // augment time
> > > t += dt;
> > > +
> > > + // check if data should be written
> > > + if (t >= saveStep)
> > > + {
> > > + // write data
> > > + vtkout(grid,c,"pconc",counter,t,grid.comm().rank());
> > > +
> > > + //increase counter and saveStep for next interval
> > > + saveStep += saveInterval;
> > > + ++counter;
> > > + }
> > > +
> > > + // print info about time, timestep size and counter
> > > if (grid.comm().rank()==0)
> > /*@\label{pfc:rank0}@*/
> > > 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);
> > > + vtkout(grid,c,"pconc",counter,tend,grid.comm().rank());
> > > }
> > >
> > > //===============================================================
> > > @@ -68,9 +92,36 @@ int main (int argc , char ** argv)
> > >
> > > // start try/catch block to get error messages from dune
> > > try {
> > > - UnitCube<Dune::YaspGrid<2>,64> uc;
> > > + using namespace Dune;
> > > +
> > > + /*
> > > + UnitCube<YaspGrid<2>,64> uc;
> > > uc.grid().globalRefine(2);
> > > partimeloop(uc.grid(),0.5);
> > > + */
> > > +
> > > + // use unitcube from dgf grids
> > > + std::stringstream dgfFileName;
> > > + dgfFileName << "grids/unitcube" << GridType :: dimension << ".dgf";
> > > +
> > > + // create grid pointer, GridType is defined by gridtype.hh
> > > + GridPtr<GridType> gridPtr( dgfFileName.str() );
> > > +
> > > + // grid reference
> > > + GridType& grid = *gridPtr;
> > > +
> > > + // half grid width 4 times
> > > + int level = 6 * DGFGridInfo<GridType>::refineStepsForHalf();
> > > +
> > > + // refine grid until upper limit of level
> > > + grid.globalRefine(level);
> > > +
> > > + // re-partition grid for better load balancing
> > > + grid.loadBalance();
> > /*@\label{pfv:lb}@*/
> > > +
> > > + // do time loop until end time 0.5
> > > + partimeloop(grid, 0.5);
> > > +
> > > }
> > > catch (std::exception & e) {
> > > std::cout << "STL ERROR: " << e.what() << std::endl;
> > > diff --git a/traversal.cc b/traversal.cc
> > > index de444bb..f061529 100644
> > > --- a/traversal.cc
> > > +++ b/traversal.cc
> > > @@ -24,20 +24,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}@*/
> > > @@ -50,16 +55,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++;
> > > }
> > > @@ -70,20 +76,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++;
> > > }
> > > diff --git a/unitcube.hh b/unitcube.hh
> > > index 6683610..c1ef749 100644
> > > --- a/unitcube.hh
> > > +++ b/unitcube.hh
> > > @@ -4,116 +4,8 @@
> > > #include <dune/common/exceptions.hh>
> > > #include <dune/common/fvector.hh>
> > >
> > > -#include <dune/grid/common/gridfactory.hh>
> > > -
> > > -// UGGrid 3d, variant 2 (tetrahedra) specialization
> > > -template< int dim >
> > > -class BasicUnitCube;
> > > -
> > > -
> > > -template<>
> > > -class BasicUnitCube< 2 >
> > > -{
> > > -protected:
> > > - template< class Grid >
> > > - static void insertVertices ( Dune::GridFactory< Grid > &factory )
> > > - {
> > > - Dune::FieldVector<double,2> pos;
> > > -
> > > - pos[0] = 0; pos[1] = 0;
> > > - factory.insertVertex(pos);
> > > -
> > > - pos[0] = 1; pos[1] = 0;
> > > - factory.insertVertex(pos);
> > > -
> > > - pos[0] = 0; pos[1] = 1;
> > > - factory.insertVertex(pos);
> > > -
> > > - pos[0] = 1; pos[1] = 1;
> > > - factory.insertVertex(pos);
> > > - }
> > > -
> > > - template< class Grid >
> > > - static void insertSimplices ( Dune::GridFactory< Grid > &factory )
> > > - {
> > > - const Dune::GeometryType type( Dune::GeometryType::simplex, 2 );
> > > - std::vector< unsigned int > cornerIDs( 3 );
> > > -
> > > - cornerIDs[0] = 0; cornerIDs[1] = 1; cornerIDs[2] = 2;
> > > - factory.insertElement( type, cornerIDs );
> > > -
> > > - cornerIDs[0] = 2; cornerIDs[1] = 1; cornerIDs[2] = 3;
> > > - factory.insertElement( type, cornerIDs );
> > > - }
> > > -
> > > - template< class Grid >
> > > - static void insertCubes ( Dune::GridFactory< Grid > &factory )
> > > - {
> > > - const Dune::GeometryType type( Dune::GeometryType::cube, 2 );
> > > - std::vector< unsigned int > cornerIDs( 4 );
> > > - for( int i = 0; i < 4; ++i )
> > > - cornerIDs[ i ] = i;
> > > - factory.insertElement( type, cornerIDs );
> > > - }
> > > -};
> > > -
> > > -
> > > -template<>
> > > -class BasicUnitCube< 3 >
> > > -{
> > > -protected:
> > > - template< class Grid >
> > > - static void insertVertices ( Dune::GridFactory< Grid > &factory )
> > > - {
> > > - Dune::FieldVector< double, 3 > pos;
> > > -
> > > - pos[0] = 0; pos[1] = 0; pos[2] = 0; factory.insertVertex(pos);
> > > - pos[0] = 1; pos[1] = 0; pos[2] = 0; factory.insertVertex(pos);
> > > - pos[0] = 0; pos[1] = 1; pos[2] = 0; factory.insertVertex(pos);
> > > - pos[0] = 1; pos[1] = 1; pos[2] = 0; factory.insertVertex(pos);
> > > - pos[0] = 0; pos[1] = 0; pos[2] = 1; factory.insertVertex(pos);
> > > - pos[0] = 1; pos[1] = 0; pos[2] = 1; factory.insertVertex(pos);
> > > - pos[0] = 0; pos[1] = 1; pos[2] = 1; factory.insertVertex(pos);
> > > - pos[0] = 1; pos[1] = 1; pos[2] = 1; factory.insertVertex(pos);
> > > - }
> > > -
> > > - template< class Grid >
> > > - static void insertSimplices ( Dune::GridFactory< Grid > &factory )
> > > - {
> > > - const Dune::GeometryType type( Dune::GeometryType::simplex, 3 );
> > > - std::vector< unsigned int > cornerIDs( 4 );
> > > -
> > > - cornerIDs[0] = 0; cornerIDs[1] = 1; cornerIDs[2] = 2;
> > cornerIDs[3] = 4;
> > > - factory.insertElement( type, cornerIDs );
> > > -
> > > - cornerIDs[0] = 1; cornerIDs[1] = 3; cornerIDs[2] = 2;
> > cornerIDs[3] = 7;
> > > - factory.insertElement( type, cornerIDs );
> > > -
> > > - cornerIDs[0] = 1; cornerIDs[1] = 7; cornerIDs[2] = 2;
> > cornerIDs[3] = 4;
> > > - factory.insertElement( type, cornerIDs );
> > > -
> > > - cornerIDs[0] = 1; cornerIDs[1] = 7; cornerIDs[2] = 4;
> > cornerIDs[3] = 5;
> > > - factory.insertElement( type, cornerIDs );
> > > -
> > > - cornerIDs[0] = 4; cornerIDs[1] = 7; cornerIDs[2] = 2;
> > cornerIDs[3] = 6;
> > > - factory.insertElement( type, cornerIDs );
> > > - }
> > > -
> > > - template< class Grid >
> > > - static void insertCubes ( Dune::GridFactory< Grid > &factory )
> > > - {
> > > - const Dune::GeometryType type( Dune::GeometryType::cube, 3 );
> > > - std::vector< unsigned int > cornerIDs( 8 );
> > > - for( int i = 0; i < 8; ++i )
> > > - cornerIDs[ i ] = i;
> > > - factory.insertElement( type, cornerIDs );
> > > - }
> > > -};
> > > -
> > > -
> > > -
> > > // default implementation for any template parameter
> > > -template<typename T, int variant>
> > > +template<typename T, int variant>
> > /*@\label{uc:uc0}@*/
> > > class UnitCube
> > > {
> > > public:
> > > @@ -133,7 +25,10 @@ public:
> > > private:
> > > // the constructed grid object
> > > T grid_;
> > > -};
> > > +};
> > /*@\label{uc:uc1}@*/
> > > +
> > > +// include basic unitcube using GridFactory concept
> > > +#include "basicunitcube.hh"
> > >
> > > // include specializations
> > > #include"unitcube_onedgrid.hh"
> > > diff --git a/unitcube_alugrid.hh b/unitcube_alugrid.hh
> > > index 86f64b6..9e2efb1 100644
> > > --- a/unitcube_alugrid.hh
> > > +++ b/unitcube_alugrid.hh
> > > @@ -3,25 +3,38 @@
> > >
> > > #if HAVE_ALUGRID
> > > #include <dune/grid/alugrid.hh>
> > > +#include <dune/grid/alugrid/3d/alu3dgridfactory.hh>
> > >
> > > -// ALU3dGrid tetrahedra specialization. Note: element type determined by
> > type
> > > +// ALU3dGrid and ALU2dGrid simplex specialization.
> > > +// Note: element type determined by type
> > > template<>
> > > class UnitCube<Dune::ALUSimplexGrid<3,3>,1>
> > > +: public BasicUnitCube< 3 >
> > > {
> > > public:
> > > typedef Dune::ALUSimplexGrid<3,3> GridType;
> > >
> > > - UnitCube () : filename("grids/cube.tetra"), grid_(filename.c_str())
> > > - {}
> > > +private:
> > > + GridType * grid_;
> > >
> > > - GridType& grid ()
> > > +public:
> > > + UnitCube ()
> > > {
> > > - return grid_;
> > > + Dune::GridFactory< GridType > factory;
> > > + BasicUnitCube< 3 >::insertVertices( factory );
> > > + BasicUnitCube< 3 >::insertSimplices( factory );
> > > + grid_ = factory.createGrid( );
> > > }
> > >
> > > -private:
> > > - std::string filename;
> > > - GridType grid_;
> > > + ~UnitCube()
> > > + {
> > > + delete grid_;
> > > + }
> > > +
> > > + GridType &grid ()
> > > + {
> > > + return *grid_;
> > > + }
> > > };
> > >
> > > // ALU2SimplexGrid 2d specialization. Note: element type determined by
> > type
> > > @@ -47,21 +60,32 @@ private:
> > > // ALU3dGrid hexahedra specialization. Note: element type determined by
> > type
> > > template<>
> > > class UnitCube<Dune::ALUCubeGrid<3,3>,1>
> > > +: public BasicUnitCube< 3 >
> > > {
> > > public:
> > > typedef Dune::ALUCubeGrid<3,3> GridType;
> > >
> > > - UnitCube () : filename("grids/cube.hexa"), grid_(filename.c_str())
> > > - {}
> > > +private:
> > > + GridType * grid_;
> > >
> > > - GridType& grid ()
> > > +public:
> > > + UnitCube ()
> > > {
> > > - return grid_;
> > > + Dune::GridFactory< GridType > factory;
> > > + BasicUnitCube< 3 >::insertVertices( factory );
> > > + BasicUnitCube< 3 >::insertCubes( factory );
> > > + grid_ = factory.createGrid( );
> > > }
> > >
> > > -private:
> > > - std::string filename;
> > > - GridType grid_;
> > > + ~UnitCube()
> > > + {
> > > + delete grid_;
> > > + }
> > > +
> > > + GridType &grid ()
> > > + {
> > > + return *grid_;
> > > + }
> > > };
> > > #endif
> > >
> > > diff --git a/unitcube_uggrid.hh b/unitcube_uggrid.hh
> > > index 79c6757..bc7c02b 100644
> > > --- a/unitcube_uggrid.hh
> > > +++ b/unitcube_uggrid.hh
> > > @@ -24,7 +24,8 @@ public:
> > > else if( variant == 2 )
> > > BasicUnitCube< dim >::insertSimplices( factory );
> > > else
> > > - DUNE_THROW( Dune::NotImplemented, "Variant " << variant << " of UG
> > unit cube not implemented." );
> > > + DUNE_THROW( Dune::NotImplemented, "Variant "
> > > + << variant << " of UG unit cube not implemented." );
> > > grid_ = factory.createGrid();
> > > }
> > >
> > > diff --git a/vertexdata.hh b/vertexdata.hh
> > > index 8fecd47..c4bb8bb 100644
> > > --- a/vertexdata.hh
> > > +++ b/vertexdata.hh
> > > @@ -23,8 +23,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>
> > > @@ -34,8 +38,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));
> > > @@ -57,7 +61,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
> > > }
> > > diff --git a/visualization.cc b/visualization.cc
> > > index 3ccf52b..c3af5d6 100644
> > > --- a/visualization.cc
> > > +++ b/visualization.cc
> > > @@ -43,28 +43,33 @@ int main(int argc, char **argv)
> > > {
> > > /*
> > > UnitCube<Dune::OneDGrid,1> uc0;
> > > - UnitCube<Dune::YaspGrid<3>,1> uc1;
> > > - UnitCube<Dune::YaspGrid<2>,1> uc2;
> > > - */
> > > + UnitCube<Dune::YaspGrid<dimGrid>,1> uc1;
> > > +
> > > #if HAVE_UG
> > > - UnitCube< Dune::UGGrid< dimGrid >, 2 > uc6;
> > > - dowork( uc6.grid(), 3 );
> > > + UnitCube< Dune::UGGrid< dimGrid >, 2 > uc2;
> > > + dowork( uc2.grid(), 3 );
> > > #endif
> > >
> > > #ifdef GRIDDIM
> > > #if HAVE_ALBERTA
> > > - UnitCube< Dune::AlbertaGrid< dimGrid, dimGrid >, 1 > uc7;
> > > + UnitCube< Dune::AlbertaGrid< dimGrid, dimGrid >, 1 > uc3;
> > > // note: The 3d cube cannot be bisected recursively
> > > - dowork( uc7.grid(), (dimGrid < 3 ? 6 : 0 ) );
> > > + dowork( uc3.grid(), (dimGrid < 3 ? 6 : 0 ) );
> > > #endif
> > > #endif
> > > + */
> > >
> > > UnitCube< Dune::SGrid< dimGrid, dimGrid >, 1 > uc4;
> > > dowork( uc4.grid(), 3 );
> > >
> > > #if HAVE_ALUGRID
> > > - UnitCube< Dune::ALUSimplexGrid< dimGrid, dimGrid > , 1 > uc8;
> > > - dowork( uc8.grid(), 3 );
> > > + UnitCube< Dune::ALUSimplexGrid< dimGrid, dimGrid > , 1 > uc5;
> > > + dowork( uc5.grid(), 3 );
> > > +
> > > +#if GRIDDIM == 3
> > > + UnitCube< Dune::ALUCubeGrid< dimGrid, dimGrid > , 1 > uc6;
> > > + dowork( uc6.grid(), 3 );
> > > +#endif
> > > #endif
> > > }
> > > catch (std::exception & e) {
> > > diff --git a/vtkout.hh b/vtkout.hh
> > > index 562cb71..a810ec9 100644
> > > --- a/vtkout.hh
> > > +++ b/vtkout.hh
> > > @@ -2,11 +2,20 @@
> > > #include <stdio.h>
> > >
> > > template<class G, class V>
> > > -void vtkout (const G& grid, const V& c, const char* name, int k)
> > > +void vtkout (const G& grid, const V& c, const char* name, int k, double
> > time=0.0, int rank=0)
> > > {
> > > Dune::VTKWriter<typename G::LeafGridView> vtkwriter(grid.leafView());
> > > char fname[128];
> > > + char sername[128];
> > > sprintf(fname,"%s-%05d",name,k);
> > > + sprintf(sername,"%s.series",name);
> > > vtkwriter.addCellData(c,"celldata");
> > > vtkwriter.write(fname,Dune::VTKOptions::ascii);
> > > +
> > > + if ( rank == 0)
> > > + {
> > > + std::ofstream serstream(sername, (k==0 ? std::ios_base::out :
> > std::ios_base::app));
> > > + serstream << k << " " << fname << ".vtu " << time << std::endl;
> > > + serstream.close();
> > > + }
> > > }
> > > diff --git a/writePVD b/writePVD
> > > new file mode 100644
> > > index 0000000..14525da
> > > --- /dev/null
> > > +++ b/writePVD
> > > @@ -0,0 +1,64 @@
> > > +#! /bin/bash
> > > +
> > > +# read parameters
> > > +if [ -z "$1" ]; then
> > > + echo "Usage: writePVD [file_prefix]"
> > > + exit
> > > +fi
> > > +
> > > +param=$1
> > > +SERIESFILE="${param##*/}"
> > > +DIR="${param%$SERIESFILE}"
> > > +DIR="${DIR:-.}"
> > > +FILEPREFIX="${SERIESFILE%.*}"
> > > +SERIESFILE=$1
> > > +
> > > +echo $SERIESFILE $DIR $FILEPREFIX
> > > +
> > > +THISDIR=$PWD
> > > +# cd $DIR
> > > +
> > > +# output file name
> > > +OUTPUTFILE=$FILEPREFIX.pvd
> > > +
> > > +#write header
> > > +echo "Writing file: $OUTPUTFILE"
> > > +echo "<?xml version=\"1.0\"?>" > $OUTPUTFILE
> > > +echo "<VTKFile type=\"Collection\" version=\"0.1\"
> > byte_order=\"LittleEndian\">" >> $OUTPUTFILE
> > > +echo " <Collection>" >> $OUTPUTFILE
> > > +echo "" >> $OUTPUTFILE
> > > +
> > > +# Set loop separator to end of line
> > > +BAKIFS=$IFS
> > > +IFS=$(echo -en "\n\b")
> > > +# the following applied to the series file
> > > +exec 3<&0
> > > +exec 0<$SERIESFILE
> > > +while read line
> > > +do
> > > + # file name and time stamp
> > > + NAME=$(echo $line | cut -d " " -f 2)
> > > + TIME=$(echo $line | cut -d " " -f 3)
> > > + # strip unnecessary path from file name
> > > + while ! test -e $NAME
> > > + do
> > > + NAME="${NAME#*/}"
> > > + if [ $NAME == "" ];
> > > + then # file does not exists
> > > + NAME=$(echo $line | cut -d " " -f 2)
> > > + echo $NAME " not found!"
> > > + exit
> > > + fi
> > > + done
> > > + # write line into output file
> > > + echo "<DataSet timestep=\"$TIME\" group=\"\" part=\"0\"
> > file=\"$NAME\"/>" >> $OUTPUTFILE
> > > +done
> > > +exec 0<&3
> > > +# restore $IFS which was used to determine what the field separators are
> > > +IFS=$BAKIFS
> > > +
> > > +# write the end of file
> > > +echo "" >> $OUTPUTFILE
> > > +echo " </Collection>" >> $OUTPUTFILE
> > > +echo "</VTKFile>" >> $OUTPUTFILE
> > > +
> >
> >
> >
>
>
> _______________________________________________
> Dune mailing list
> Dune at dune-project.org
> http://lists.dune-project.org/mailman/listinfo/dune
More information about the Dune
mailing list