[Dune] dune-grid-howto broken
Martin Drohmann
mdrohmann at uni-muenster.de
Fri Mar 6 18:08:43 CET 2009
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.
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
> > +
>
>
>
--
Martin Drohmann
Institut für Numerische und Angewandte Mathematik
Westfälische Wilhelms-Universität Münster
Einsteinstr. 62, D-48149 Münster
Tel. +49 (0) 251 83-32753
Fax. +49 (0) 251 83-32729
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.dune-project.org/pipermail/dune/attachments/20090306/c0721232/attachment.htm>
More information about the Dune
mailing list