[Dune] dune-grid-howto broken
Christian Engwer
christi at uni-hd.de
Fri Mar 6 16:11:22 CET 2009
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
On Fri, Mar 06, 2009 at 04:03:27PM +0100, Oliver Sander wrote:
> Hi Christian!
> Martin Drohmann hatte mir ein großes Patch geschickt, wo Text
> und Code geändert wird. Die Code-Änderungen erschienen zweifel-
> haft (viele WhiteSpace-Änderungen, ich hatte aber auch nicht so
> viel Zeit, mir das anzuschauen). Was dieses basisccube ding soll
> weiß ich aber auch nicht.
>
> Ich bin schon auf dem Weg in dem Urlaub. Anbei findest Du das
> Patch von Martin. Bitte kümmer Du Dich darum.
>
> Vielen Dank, viele Grüße,
> Oliver
>
> Christian Engwer schrieb:
>> The last changes of rund-grid-howto broke the howto.
>>
>> I'm not the author of these changes, so I'm not sure what was the
>> rational for basicunitcube.hh.
>>
>> LaTeX Warning: Reference `tc:lfgv' on page 12 undefined on input line
>> 441.
>>
>> [12] [13] [14]
>> Chapter 3.
>> [15] [16] [17] [18]
>>
>> LaTeX Warning: Reference `ch:mappers' on page 19 undefined on input
>> line 849.
>>
>>
>> LaTeX Warning: Reference `ch:mappers' on page 19 undefined on input
>> line 859.
>>
>> [19]
>> Chapter 4.
>> (../unitcube.hh [20] [21]) (../unitcube_onedgrid.hh [22])
>> (../unitcube_sgrid.hh
>> ) (../unitcube_yaspgrid.hh [23])
>>
>> ! Package Listings Error: File `../basicunitcube(.hh)' not found.
>>
>> There is no file basicunitcube.hh. Is this just some kind of wildcard
>> for any of the unicube files? Or is this a real file, but wasn't
>> commited?
>>
>> Christian
>>
>> _______________________________________________
>> Dune mailing list
>> Dune at dune-project.org
>> http://lists.dune-project.org/mailman/listinfo/dune
>>
>
>
>
> diff --git a/adaptivefinitevolume.cc b/adaptivefinitevolume.cc
> index a9f2241..cd1a607 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;
> @@ -99,39 +62,42 @@ void timeloop (G& grid, double tend, int lmin, int lmax)
>
> std::cout << "s=" << grid.size(0) << " k=" << k << " t=" << t << std::endl;
> while (t<tend)
> - {
> - // augment time step counter
> - ++k;
> -
> - // apply finite volume scheme
> - evolve(grid,mapper,c,t,dt);
> -
> - // augment time
> - t += dt;
> -
> - // check if data should be written
> - if (t >= saveStep)
> {
> - // write data
> - vtkout(grid,c,"concentration",counter);
> -
> - // increase counter and saveStep for next interval
> - saveStep += saveInterval;
> - ++counter;
> + // augment time step counter
> + ++k;
> +
> + // apply finite volume scheme
> + evolve(grid,mapper,c,t,dt);
> +
> + // augment time
> + t += dt;
> +
> + // check if data should be written
> + if (t >= saveStep)
> + {
> + // write data
> + vtkout(grid,c,"concentration",counter,t);
> +
> + // increase counter and saveStep for next interval
> + saveStep += saveInterval;
> + ++counter;
> + }
> +
> + // print info about time, timestep size and counter
> + 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); /*@\label{afv:ad}@*/
> + }
> }
>
> - // print info about time, timestep size and counter
> - 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);
> - }
> - }
> -
> // write last time step
> - vtkout(grid,c,"concentration",counter);
> + vtkout(grid,c,"concentration",counter,tend);
> +
> + // write
> }
>
> //===============================================================
> @@ -140,14 +106,14 @@ void timeloop (G& grid, double tend, int lmin, int lmax)
>
> int main (int argc , char ** argv)
> {
> - // initialize MPI, finalize is done automatically on exit
> + // initialize MPI, finalize is done automatically on exit
> Dune::MPIHelper::instance(argc,argv);
>
> // start try/catch block to get error messages from dune
> try {
> using namespace Dune;
>
> - // use unitcube from grids
> + // use unitcube from grids
> std::stringstream dgfFileName;
> dgfFileName << "grids/unitcube" << GridType :: dimension << ".dgf";
>
> @@ -161,12 +127,12 @@ int main (int argc , char ** argv)
> int minLevel = 2 * DGFGridInfo<GridType>::refineStepsForHalf();
>
> // refine grid until upper limit of level
> - grid.globalRefine(minLevel);
> + grid.globalRefine(minLevel);
>
> - // maximal allowed level during refinement
> + // maximal allowed level during refinement
> int maxLevel = minLevel + 3 * DGFGridInfo<GridType>::refineStepsForHalf();
>
> - // do time loop until end time 0.5
> + // do time loop until end time 0.5
> timeloop(grid, 0.5, minLevel, maxLevel);
> }
> catch (std::exception & e) {
> diff --git a/adaptiveintegration.cc b/adaptiveintegration.cc
> index ded8d9c..0fa792c 100644
> --- a/adaptiveintegration.cc
> +++ b/adaptiveintegration.cc
> @@ -1,26 +1,31 @@
> // $Id$
>
> -#include"config.h"
> -#include<iostream>
> -#include<iomanip>
> -#include<dune/grid/io/file/vtk/vtkwriter.hh> // VTK output routines
> +#include "config.h"
> +#include <iostream>
> +#include <iomanip>
> +#include <dune/grid/io/file/vtk/vtkwriter.hh> // VTK output routines
> #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"unitcube.hh"
> -#include"functors.hh"
> -#include"integrateentity.hh"
> +#include "unitcube.hh"
> +#include "functors.hh"
> +#include "integrateentity.hh"
>
>
> //! adaptive refinement test
> 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;
> @@ -32,9 +37,9 @@ void adaptiveintegration (Grid& grid, const Functor& f)
> for (int k=0; k<100; k++)
> {
> // compute integral on current mesh
> - double value=0; /*@\label{aic:int0}@*/
> - for (ElementLeafIterator it = grid.template leafbegin<0>();
> - it!=grid.template leafend<0>(); ++it)
> + double value=0; /*@\label{aic:int0}@*/
> + for (ElementLeafIterator it = gridView.template begin<0>();
> + it!=gridView.template end<0>(); ++it)
> value += integrateentity(it,f,highorder); /*@\label{aic:int1}@*/
>
> // 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..5a15cb2
> --- /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/doc/grid-howto.bib b/doc/grid-howto.bib
> index e37fb33..964b930 100644
> --- a/doc/grid-howto.bib
> +++ b/doc/grid-howto.bib
> @@ -104,7 +104,7 @@
>
>
> @Book{Deuflhard93,
> - author = {P Deuflhard and A. Hohmann},
> + author = {P. Deuflhard and A. Hohmann},
> title = {Numerische Mathematik I},
> publisher = {Walter de Gruyter},
> year = {1993},
> diff --git a/doc/grid-howto.tex b/doc/grid-howto.tex
> index fb7e2f0..68fc9ae 100644
> --- a/doc/grid-howto.tex
> +++ b/doc/grid-howto.tex
> @@ -249,7 +249,7 @@ Documentation of the files and classes in \Dune{} is provided in code and
> can be extracted using the
> doxygen\footnote{\href{http://www.stack.nl/~dimitri/doxygen/}{\texttt{http://www.stack.nl/$\sim$dimitri/doxygen/}}}
> software available elsewhere. The code documentation can either be built
> -locally on your machine (in html and other formats, e.~g.~\LaTeX)
> +locally on your machine (in html and other formats, e.g.~\LaTeX)
> or its latest version is available at
> \begin{center}
> \href{http://www.dune-project.org/doc/}%
> @@ -385,6 +385,14 @@ this case all entities are of type ``cube''.
> modest. Also try dimensions 4 and 5!
> \end{exc}
>
> +\begin{exc}
> + The static methods \lstinline!Dune::gridlevellist! and
> + \lstinline!Dune::gridleaflist! produce a very detailed output of the grid's
> + elements on a specified grid level. Change the code and print out this
> + information for the leaf grid or a grid on lower level. Try to understand the
> + output.
> +\end{exc}
> +
> \section{Traversing a grid --- A first look at the grid interface}
>
> After looking at very first simple example we are now ready to go on
> @@ -397,7 +405,7 @@ numberstyle=\tiny, numbersep=5pt]{../traversal.cc}
> \end{lst}
>
> The \lstinline!main! function near the end of the listing is pretty
> -similar to previous one except that we use a 2d grid for the unit
> +similar to the previous one except that we use a 2d grid for the unit
> square that just consists of one cell. In line \ref{tc:refine} this
> cell is refined once using the standard method of grid refinement of
> the implementation. Here, the cell is refined into four smaller cells.
> @@ -426,22 +434,28 @@ several. As in the STL the elements of any container can be accessed
> with iterators which are generalized pointers. Again, a \Dune{} grid
> knows several different iterators which provide access to the
> different kinds of entities and which also provide different patterns
> -of access.
> +of access.
> +
> +As we usually do not want to use the entire hierarchy of the grid, we first
> +define a view on that part of the grid we are interested in. This can be a
> +level or the leaf part of the grid. In line \ref{tc:lfgv} a type for a
> +\lstinline!GridView! on the leaf grid is defined.
>
> -Line \ref{tc:ittype} extracts the type of an iterator from the grid
> +Line \ref{tc:ittype} extracts the type of an iterator from this view
> class. \lstinline!Codim! is a \lstinline!struct! within the grid class
> that takes an integer template parameter specifying the codimension
> over which to iterate. Within the \lstinline!Codim! structure the type
> -\lstinline!LeafIterator! is defined. Since we specified codimension 0
> +\lstinline!Iterator! is defined. Since we specified codimension 0
> this iterator is used to iterate
> -over the elements which are not refined any further, i.~e.~which are
> +over the elements which are not refined any further, i.e.~which are
> the leaves of the refinement trees.
>
> The \lstinline!for!-loop in line \ref{tc:forel} now visits every such
> -element. The \lstinline!leafbegin! and \lstinline!leafend! on the grid
> -class deliver the first leaf element and one past the last leaf
> -element. Note that the \lstinline!template! keyword must be used and
> -template parameters are passed explicitely. Within the loop body in
> +element. The \lstinline!begin! and \lstinline!end! on the
> +\lstinline!LeafGridView!
> +class deliver the first leaf element and one past the last leaf element. Note
> +that the \lstinline!template! keyword must be used and template parameters are
> +passed explicitely. Within the loop body in
> lines \ref{tc:forel0}-\ref{tc:forel1} the iterator \lstinline!it! acts
> like a pointer to an entity of dimension \lstinline!dim! and
> codimension 0. The exact type would be
> @@ -457,8 +471,8 @@ mapping is also called the ``local to global'' mapping. The
> corresponding reference element has a certain type which is extracted
> in line \ref{tc:reftype}. Since the reference elements are polyhedra
> they consist of a finite number of corners. The images of the corners
> -under the local to global map can be accessed via an
> -\lstinline!operator[]!. Line \ref{tc:print} prints the geometry type
> +under the local to global map can be accessed via the
> +\lstinline!corner(int n)! method. Line \ref{tc:print} prints the geometry type
> and the position of the first corner of the element. Then line
> \ref{tc:count} just counts the number of elements visited.
>
> @@ -480,12 +494,13 @@ degrees of freedom associated to, say, faces).
>
> Finally, we show in lines \ref{tc:level0}-\ref{tc:level1} how the
> hierarchic structure of the mesh can be accessed. To that end a
> -\lstinline!LevelIterator! is used. It provides access to all entities
> -of a given codimension (here 0) on a given grid level. The coarsest
> +\lstinline!LevelGridView! is used. It provides via an \lstinline!Iterator!
> +access to all entities of a given codimension (here 0) on a given grid level.
> +The coarsest
> grid level (the initial macro grid) has number zero and the number of
> the finest grid level is returned by the \lstinline!maxLevel()! method
> -of the grid. The methods \lstinline!lbegin()! and \lstinline!lend()!
> -on the grid deliver iterators to the first and one-past-the-last
> +of the grid. The methods \lstinline!begin()! and \lstinline!end()!
> +on the view deliver iterators to the first and one-past-the-last
> entity of a given grid level supplied as an integer argument to these
> methods.
>
> @@ -674,7 +689,7 @@ definition it is clear that the concept does not necessarily exist as
> program text.
>
> A class that implements a concept is called a
> -\textit{model} of the concept. E.~g.~in the standard template library (STL)
> +\textit{model} of the concept. E.g.~in the standard template library (STL)
> the class \lstinline!std::vector<int>! is a model of the concept
> ``container''. If all instances of a class template are a model of
> a given concept we can also say that the class template is a model of
> @@ -685,7 +700,7 @@ In standard OO language a concept would be formulated as
> an abstract base class and all the models would be implemented as
> derived classes. However, for reasons of efficiency we do not want to
> use dynamic polymorphism. Moreover, concepts are more powerful because
> -the models of a concept can use different types, e.~g.~as return types of
> +the models of a concept can use different types, e.g.~as return types of
> methods. As an example consider the STL where the begin method on a
> vector of int returns \lstinline!std::vector<int>::iterator! and on a
> list of int it returns \lstinline!std::list<int>::iterator! which may
> @@ -772,6 +787,14 @@ you to the documentation of the corresponding wrapper class template
> classes are typically aggregated in the grid class and accessed via
> iterators.
>
> +\minisec{\href{TODO}{GridView}}
> +
> + The GridView gives a view on a level or the leaf part of a grid. It
> + provides iterators for access to the elements of this view and a
> + communication method for parallel computations. Alternatively, a
> + LevelIterator of a LeafIterator can be directly accessed from a grid.
> + These iterator types are described below.
> +
> \minisec{\href{http://www.dune-project.org/doc/doxygen/html/classDune_1_1Entity.html}{Entity}}
>
> The entity class encapsulates
> @@ -793,61 +816,47 @@ you to the documentation of the corresponding wrapper class template
> entity. Moreover it is immutable, i.e. the referenced entity can
> not be modified.
>
> -\minisec{\href{http://www.dune-project.org/doc/doxygen/html/classDune_1_1LevelIterator.html}{LevelIterator}}
> +\minisec{\href{http://www.dune-project.org/doc/doxygen/html/classDune_1_1Iterator.html}{Iterator}}
>
> - LevelIterator is an
> - immutable iterator that provides access to an entity. It can be
> - incremented to visit all entities of a given
> - codimension and level of the grid. An EntityPointer is
> - assignable from a LevelIterator.
> -
> -\minisec{\href{http://www.dune-project.org/doc/doxygen/html/classDune_1_1LeafIterator.html}{LeafIterator}}
> -
> - LeafIterator is an
> + Iterator is an
> immutable iterator that provides access to an entity. It can by
> incremented to visit all entities of a given
> - codimension of the leaf grid. An EntityPointer is assignable
> - from a LeafIterator.
> -
> -\minisec{\href{http://www.dune-project.org/doc/doxygen/html/classDune_1_1HierarchicIterator.html}{HierarchicIterator}}
> + codimension of a GridView. An EntityPointer is assignable
> + from an Iterator.
>
> - HierarchicIterator
> - is an immutable iterator that provides access to an entity. It can
> - be incremented to visit all entities of
> - codimension 0 that resulted from subdivision of a given entity of
> - codimension 0. An EntityPointer is assignable from a
> - HierarchicIterator.
> -
> -\minisec{\href{http://www.dune-project.org/doc/doxygen/dune-grid-html/group__GIIntersectionIterator.html}{LevelIntersectionIterator}}
> +\minisec{\href{http://www.dune-project.org/doc/doxygen/dune-grid-html/group__GIIntersectionIterator.html}{IntersectionIterator}}
>
> IntersectionIterator provides access to all entities of
> codimension 0 that have an intersection of codimension 1 with a
> given entity of codimension 0. In a conforming mesh these are the
> face neighbors of an element. For two entities with a common
> - intersection the IntersectionIterator also provides information
> - about the geometric location of the intersection. Furthermore it
> - also provides information about intersections of an entity with
> + intersection the IntersectionIterator can be dereferenced as an Intersection
> + object which in turn provides information about the geometric location of
> + the intersection. Furthermore this Intersection class also provides
> + information about intersections of an entity
> + with
> the internal or external boundaries.
> - The LevelIntersectionIterator provides intersections between
> - codimension 0 entities having the same level.
> + The IntersectionIterator provides intersections between
> + codimension 0 entities among the same GridView.
>
> -\minisec{\href{http://www.dune-project.org/doc/doxygen/dune-grid-html/group__GIIntersectionIterator.html}{LeafIntersectionIterator}}
> -
> - This iterator has the same properties as the LevelIntersectionIterator
> - but provides intersections between leaf entities of the grid.
> -
> \minisec{\href{http://www.dune-project.org/doc/doxygen/html/classDune_1_1IndexSet.html}{LevelIndexSet}, \href{http://www.dune-project.org/doc/doxygen/html/classDune_1_1IndexSet.html}{LeafIndexSet}}
>
> LevelIndexSet and LeafIndexSet which are both models of
> Dune::IndexSet are used to attach any kind of user-defined data to
> (subsets of) entities of the grid. This data is supposed to be
> - stored in one-dimensional arrays for reasons of efficiency.
> + stored in one-dimensional arrays for reasons of efficiency. An IndexSet is
> + usually not used directly but through a Mapper (c.f.~chapter
> + \ref{ch:mappers}).
>
> \minisec{\href{http://www.dune-project.org/doc/doxygen/html/classDune_1_1IdSet.html}{LocalIdSet}, \href{http://www.dune-project.org/doc/doxygen/html/classDune_1_1IdSet.html}{GlobalIdSet}}
>
> LocalIdSet and GlobalIdSet which are both models of Dune::IdSet
> are used to save user data during a grid refinement phase and
> - during dynamic load balancing in the parallel case.
> + during dynamic load balancing in the parallel case. The LocalIdSet is
> + unique for all entities on the current partition, whereas the GlobalIdSet
> + gives a unique mapping over all grid partitions. An IdSet is usually not used
> + directly but through a Mapper (c.f.~chapter
> + \ref{ch:mappers}).
>
>
> \section{Propagation of type information}
> @@ -865,8 +874,8 @@ compiler. As in the Standard Template Library strange error messages
> may occur if you try to mix these types.
>
> In order to avoid these problems the related types of an
> -implementation are as public types from most classes of an
> -implementation. E.~g., in order to extract the
> +implementation are provided as public types by most classes of an
> +implementation. E.g., in order to extract the
> \lstinline!XLevelIterator! (for codimension 0) from the
> \lstinline!XGrid! class you would write
> \begin{lstlisting}
> @@ -951,11 +960,28 @@ used.
> numberstyle=\tiny, numbersep=5pt]{../unitcube_yaspgrid.hh}
> \end{lst}
>
> +\minisec{GridFactory}
> +
> +The file \lstinline!basicunitcube.hh!\ provides an interface class for the
> +set-up of a unit cube in two or three dimensions with simplicial or cubic
> +elements with help of the dune-grid class \lstinline!GridFactory!. This class
> +hides grid specific methods for insertion of vertices and elements in the macro
> +grid. Specializations of the \lstinline!GridFactory!\ exist for
> +\lstinline!UGGrid!, \lstinline!AlbertaGrid!\ and the three dimensional
> +\lstinline!ALUGrid!\ objects.
> +
> +\begin{lst}[File dune-grid-howto/basicunitcube.hh] \mbox{}
> +\nopagebreak
> +\lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
> +numberstyle=\tiny, numbersep=5pt]{../basicunitcube.hh}
> +\end{lst}
> +
> \minisec{UGGrid}
>
> -The following listing shows how to create \lstinline!UGGrid!
> -objects. Two and three-dimensional versions are available. The
> -\lstinline!variant! parameter can take on two values: 1 for
> +The following listing shows how to create \lstinline!UGGrid!\
> +objects based on the \lstinline!BasicUnitCube!. Two and three-dimensional
> +versions are available. The \lstinline!variant! parameter can take on two
> +values: 1 for
> quadrilateral/hexahedral grids and 2 for triangular/tetrahedral
> grids. The initial grids are read in AmiraMesh format.
>
> @@ -969,14 +995,17 @@ numberstyle=\tiny, numbersep=5pt]{../unitcube_uggrid.hh}
>
> The following listing contains specializations of the
> \lstinline!UnitCube! template for Alberta in two and three
> -dimensions. When using Alberta the \Dune{} framework has to be
> -configured with a dimension (\lstinline!--with-alberta-dim=2!,
> +dimensions. When using Alberta versions less than $2.0$ the \Dune{} framework
> +has to be configured with a dimension (\lstinline!--with-alberta-dim=2!,
> \lstinline!--with-alberta-world-dim=2!) and only this dimension can then be
> used.
> The dimension from the configure run is available in the macro
> \lstinline!ALBERTA_DIM! and \lstinline!ALBERTA_WORLD_DIM!
> in the file \lstinline!config.h! (see
> -next section). The \lstinline!variant! parameter must be 1.
> +next section). The \lstinline!variant! parameter must be 1. Like in the
> +\lstinline!UGGrid!\ implementation the grid factory concept is used as the
> +\lstinline!UnitCube!\ specializations are based on the
> +\lstinline!BasicUnitCube!.
>
> \begin{lst}[File dune-grid-howto/unitcube\_albertagrid.hh] \mbox{}
> \nopagebreak
> @@ -986,13 +1015,15 @@ numberstyle=\tiny, numbersep=5pt]{../unitcube_albertagrid.hh}
>
> \minisec{ALUGrid}
>
> -The next listing shows the instantiation of \lstinline!ALUSimplexGrid! or
> -\lstinline!ALUCubeGrid! objects.
> -The ALUGrid implementation supports either simplicial grids ,i.e.
> +The next listing shows the instantiation of \lstinline!ALUSimplexGrid!\ or
> +\lstinline!ALUCubeGrid!\ objects.
> +The ALUGrid implementation supports either simplicial grids, i.e.\
> tetrahedral or triangular grids, and hexahedral grids and the
> element type has to be chosen at compile-time. This is done by choosing
> either \lstinline!ALUSimplexGrid! or \lstinline!ALUCubeGrid!.
> -The \lstinline!variant! parameter must be 1.
> +The \lstinline!variant!\ parameter must be 1.
> +For three dimensional versions the grid objects are set up with help of the
> +\lstinline!BasicUnitCube!\ class.
>
> \begin{lst}[File dune-grid-howto/unitcube\_alugrid.hh] \mbox{}
> \nopagebreak
> @@ -1003,20 +1034,18 @@ numberstyle=\tiny, numbersep=5pt]{../unitcube_alugrid.hh}
> \section{Using configuration information provided by configure}
>
> The \lstinline!./configure! script in the application
> -(\lstinline!dune-grid-howto! here) produces a file
> -\lstinline!config.h! that contains information about the configuration
> -parameters. E.~g.~which of the optional grid implementations is
> +(\lstinline!dune-grid-howto!\ here) produces a file
> +\lstinline!config.h!\ that contains information about the configuration
> +parameters. E.g.~which of the optional grid implementations is
> available and which dimension has been selected (if applicable). This
> information can then be used at compile-time to include header files
> or code that depend on optional packages.
>
> -As an example, the macro \lstinline!HAVE_UG! can be used to compile
> +As an example, the macro \lstinline!HAVE_UG!\ can be used to compile
> UG-specific code as in
> \begin{lstlisting}[basicstyle=\ttfamily\scriptsize]
> #if HAVE_UG
> #include"dune/grid/uggrid.hh"
> -#include"dune/io/file/amirameshreader.hh"
> -#endif
> \end{lstlisting}
>
> It is important that the file \lstinline!config.h! is the first
> @@ -1028,28 +1057,28 @@ Dune has its own macro grid format, the \underline{D}une \underline{G}rid \under
> A detailed description of the DGF and how to use it can be found on the
> homepage of Dune under the documentation section (see
> \href{http://www.dune-project.org/doc/doxygen/dune-grid-html/group__DuneGridFormatParser.html}%
> -{{\small\texttt{http://www.dune-project.org/doc/doxygen/dune-grid-html/group\_\_DuneGridFormatParser.html}}}).
> -
> -Here we only give a short introduction. To use the DGF parser the configuration
> -option \\
> -\lstinline!--with-grid-dim={1,2,3}! must be provided
> -during configuration run. Optional \\
> -\lstinline!--with-grid-type=ALBERTAGRID!.
> -Furthermore, \lstinline!ALUGRID_CUBE,ALUGRID_SIMPLEX,ALUGRID_CONFORM!,
> -\lstinline!ONEDGRID,SGRID,UGGRID!, and \lstinline!YASPGRID! can be chosen as grid types.
> -can be chosen. Note that both values will also be changeable later.
> +{DuneGridFormatParser}\footnote{\href{http://www.dune-project.org/doc/doxygen/dune-grid-html/group__DuneGridFormatParser.html}%
> +{{\small\texttt{http://www.dune-project.org/doc/doxygen/dune-grid-html/group\_\_DuneGridFormatParser.html}}}}).
> +
> +Here we only give a short introduction. The configuration
> +\lstinline!--with-grid-dim={1,2,3}! must be provided during configuration run
> +in order to use the DGF parser. A default grid type can be chosen with the
> +configuration option \lstinline!--with-grid-type=ALBERTAGRID!.
> +Further possible grid types are \lstinline!ALUGRID_CUBE,ALUGRID_SIMPLEX,ALUGRID_CONFORM!,
> +\lstinline!ONEDGRID,SGRID,UGGRID!, and \lstinline!YASPGRID!.
> +Note that both values will also be changeable later.
> If the \lstinline!--with-grid-dim! option was not provided during configuration the
> DFG grid type definition will not work. Nevertheless, the grid parser will work
> but the grid type has to be defined by the user and the appropriate DGF parser
> specialization has to be included.
> -Assuming the \lstinline!--with-grid-dim! was provided the DGF grid type
> +Assuming the \lstinline!--with-grid-dim!\ was provided the DGF grid type
> definition works by first including \lstinline!gridtype.hh!.
> \begin{lstlisting}[basicstyle=\ttfamily\scriptsize]
> -#include <dune/grid/io/file/dgfparser/gridtype.hh>
> +#include <dune/grid/io/file/dgfparser/dgfgridtype.hh>
> \end{lstlisting}
> -Depending on the pre-configured values of \lstinline!GRIDDIM! and
> +Depending on the pre-configured values of \lstinline!GRIDDIM!\ and
> \lstinline!GRIDTYPE! a typedef for the grid to use will be provided by
> -including \lstinline!gridtype.hh!. The follwoing example show how an
> +including \lstinline!dgfgridtype.hh!. The follwoing example show how an
> instance of the defined grid is generated. Given a DGF file, for example
> \lstinline!unitcube2.dgf!, a grid pointer is created as follows.
> \begin{lstlisting}[basicstyle=\ttfamily\scriptsize]
> @@ -1184,7 +1213,7 @@ refined. Then
> E \approx |I_f-I_c|
> \end{equation}
> is an estimate for the error. If
> -the refinement is such that every element is halfened in every
> +the refinement is such that every element is bisected in every
> coordinate direction, the function to be integrated is sufficiently
> smooth and the order of the quadrature rule is $p+1$,
> then the error should be reduced by a factor of $(1/2)^p$ after
> @@ -1242,7 +1271,7 @@ whether an entity has already been visited by some algorithm or
> not. In this chapter we will show with some examples how arbitrary
> user data can be attached to a grid.
>
> -\section{Mappers}
> +\section{Mappers}\label{ch:mappers}
>
> The general situation is that a user wants to store some arbitrary
> data with a subset of the entities of a grid. Remember that entities
> @@ -1269,7 +1298,7 @@ index. The set of indices starts at zero and is consecutive.
>
> Let us assume that the set of all entities in the grid is $E$ and
> that $E^\prime\subseteq E$ is the subset of entities for which data is
> -to be stored. E.~g.~this could be all the vertices in the leaf grid in
> +to be stored. E.g.~this could be all the vertices in the leaf grid in
> the case of $P_1$ finite elements. Then the access from grid entities
> to user data is a two stage process: A so-called \textit{mapper}
> provides a map
> @@ -1306,7 +1335,7 @@ shown below.
>
> \minisec{Id based mappers}
>
> -Id-based mapper can also be used while a grid changes, i.~e.~it is
> +Id-based mappers can also be used while a grid changes, i.e.~it is
> ensured that the map $m$ can still be evaluated for all entities $e$
> that are still in the grid after modification. For that it
> has to be implemented on the basis of a \lstinline!Dune::IdSet!. This may be
> @@ -1323,7 +1352,7 @@ id-based mapper would be used to transfer the required data
>
> \section{Visualization of discrete functions}
>
> -Let use mappers to evaluate a function $f:\Omega\to\mathbb{R}$ for
> +Let us use mappers to evaluate a function $f:\Omega\to\mathbb{R}$ for
> certain entities and store the values in a vector. Then, in order to
> do something useful, we use the vector to produce a graphical
> visualization of the function.
> @@ -1406,7 +1435,7 @@ the second last argument of \lstinline!grape.displayVector!. This argument is th
> degree of the approximation.
> \end{itemize}
>
> -Finally the following listing shows the main programm that drives the
> +Finally the following listing shows the main program that calls the
> two functions just discussed:
>
> \begin{lst}[File dune-grid-howto/visualization.cc] \mbox{}
> @@ -1542,7 +1571,7 @@ can also be parallelized in a straightforward way.
>
> \subsection{Implementation}
>
> -First, we need to specify the problem parameters, i.~e.~initial
> +First, we need to specify the problem parameters, i.e.~initial
> condition, boundary condition and velocity field. This is done by the
> following functions.
>
> @@ -1557,8 +1586,8 @@ condition should also be straightforward now. The function
> \lstinline!initialize! works on a concentration vector \lstinline!c!
> that can be stored in any container type with a vector interface
> (\lstinline!operator[]!, \lstinline!size()! and copy constructor are
> -needed). Moreover the grid and a mapper for element-wise data have to
> -passed as well.
> +needed). Moreover the grid and a mapper for element-wise data have to be passed
> +as well.
>
> \begin{lst}[File dune-grid-howto/initialize.hh] \mbox{}
> \nopagebreak
> @@ -1601,12 +1630,10 @@ with copy constructor and size method.
> The computation of the fluxes is done in lines
> \ref{evh:flux0}-\ref{evh:flux1}. An \lstinline!IntersectionIterator!
> is used to access all intersections $\gamma_{ij}$ of a leaf element
> -$\omega_i$. For a full documentation of the
> -\lstinline!IntersectionIterator! we refer to
> -\begin{center}
> -\href{http://www.dune-project.org/doc/doxygen/dune-grid-html/group\_\_GIIntersectionIterator.html}
> -{{\small\texttt{http://www.dune-project.org/doc/doxygen/dune-grid-html/group\_\_GIIntersectionIterator.html}}}
> -\end{center}
> +$\omega_i$. For a full documentation on the
> +\lstinline!Intersection!\ class, we refer to the doxygen module page on
> +\href{http://www.dune-project.org/doc/doxygen/dune-grid-html/group\_\_GIIntersectionIterator.html}{Intersections}\footnote{\href{http://www.dune-project.org/doc/doxygen/dune-grid-html/group\_\_GIIntersectionIterator.html}
> +{{\small\texttt{http://www.dune-project.org/doc/doxygen/dune-grid-html/group\_\_GIIntersectionIterator.html}}}}
> An Intersection is with another element $\omega_j$ if the
> \lstinline!neighbor()! method of the iterator returns true (line
> \ref{evh:neighbor}) or with the external boundary if
> @@ -1665,6 +1692,14 @@ XML file format.
> numberstyle=\tiny, numbersep=5pt]{../vtkout.hh}
> \end{lst}
>
> +In addition to the snapshots that are produced at each timestep, this function
> +also generates a series file which stores the actual \lstinline!time! of an
> +evolution scheme together with the snapshots' filenames. After executing the
> +shell script \lstinline!writePVD!\ on this series file, we get a
> +\underline{P}ara\underline{v}iew {D}ata (PVD) file with the same name as the
> +snapshots. This file opened with paraview then gives us a neat animation over
> +the time period.
> +
> Finally, the main program:
>
> \begin{lst}[File dune-grid-howto/finitevolume.cc] \mbox{}
> @@ -1770,10 +1805,10 @@ numberstyle=\tiny, numbersep=5pt]{../adaptiveintegration.cc}
>
> The work is done in the function \lstinline!adaptiveintegration!.
> Lines \ref{aic:int0}-\ref{aic:int1} compute the value of the integral
> -on the current mesh. After printing the result the decission whether
> +on the current mesh. After printing the result the decision whether
> to continue or not is done in line \ref{aic:finish}. The extrapolation
> strategy relies on the fact that every element has a father. To ensure
> -this the grid is at least once refined globally in the first step
> +this, the grid is at least once refined globally in the first step
> (line \ref{aic:gr}). Now the refinement threshold $\kappa$ can be
> computed in lines \ref{aic:kappa0}-\ref{aic:kappa1}. Finally the last
> loop in lines \ref{aic:mark0}-\ref{aic:mark1} marks elements for
> @@ -1821,7 +1856,7 @@ heuristic strategy that works as follows:
> \item Compute global maximum and minimum of element concentrations:
> $$\overline{C}=\max_i C_i, \ \ \ \underline{C}=\min_i C_i.$$
> \item As the local indicator in cell $\omega_i$ we define $$\eta_i =
> - \max_{\gamma_{ij}} |C_i-C_j|$$. Here $\gamma_{ij}$ denotes
> + \max_{\gamma_{ij}} |C_i-C_j|.$$ Here $\gamma_{ij}$ denotes
> intersections with other elements in the leaf grid.
> \item If for $\omega_i$ we have
> $\eta_i>\overline{\text{tol}}\cdot (\overline{C}-\underline{C})$
> @@ -1848,8 +1883,8 @@ for an element:
> old mesh as a non-leaf element: Compute the cell value as an average
> of the son elements in the old mesh.
> \item The element is a leaf element in the new mesh and is obtained
> - from through refining some element in the old mesh: Copy the value
> - from this element in the old mesh.
> + through refining some element in the old mesh: Copy the value
> + from the element in the old mesh to the new mesh.
> \end{itemize}
>
> The complete mesh adaptation is done by the function
> @@ -1873,7 +1908,7 @@ grid levels and thereby adding the value in an element to the father
> element. The key into the map is the global id of an element. Thus the
> value is accessible also after mesh modification.
>
> -Now grid can really be modified in line \ref{fah:adapt} by calling the
> +Now the grid can really be modified in line \ref{fah:adapt} by calling the
> \lstinline!adapt()! method on the grid object. The mapper is updated
> to reflect the changes in the grid in line \ref{fah:update} and the
> concentration vector is resized to the new size in line
> @@ -1890,6 +1925,22 @@ Here is the new main program with an adapted \lstinline!timeloop!:
> numberstyle=\tiny, numbersep=5pt]{../adaptivefinitevolume.cc}
> \end{lst}
>
> +The program works analogously to the non adaptive \lstinline!finitevolume!\
> +version from the previous chapter. The only differences are inside the
> +\lstinline!timeloop!\ function. During the initialization of the concentration
> +vector in line \ref{afv:in} and after each time step in line \ref{afv:ad} the
> +function \lstinline!finitevolumeadapt! is called in order to refine the grid.
> +The initial adaptation is repeated $\overline{M}$ times. Note that adaptation
> +after each time steps is deactivated during the compiler phase for unstructured
> +grids with help of the \lstinline!Capabilities!\ class. This is because
> +structured grids do not allow a conforming refinement and are therefore
> +unusable for adaptive schemes. In fact, the \lstinline!adapt! method on a grid
> +of \lstinline!YaspGrid! e.g.~ results in a {\em global} grid refinement.
> +
> +\begin{exc}
> + Compile the program with the gridtype set to \lstinline!ALUGRID_SIMPLEX!
> + and \lstinline!ALUGRID_CONFORM! and compare the results visually.
> +\end{exc}
>
> \chapter{Parallelism}
>
> @@ -1924,8 +1975,8 @@ values:
> Entities of codimension 0 are restricted to the three partition types
> \textit{interior}, \textit{overlap} and \textit{ghost}. Entities of
> codimension greater than 0 may take all partition type values.
> -The codimension 0 entities with partition type \textit{interior} for a
> -non-overlapping decomposition of the entity set, i.~e.~for each
> +The codimension 0 entities with partition type \textit{interior} form a
> +non-overlapping decomposition of the entity set, i.e.~for each
> entity of codimension 0 there is exactly one process where this entity
> has partition type \textit{interior}.
> Moreover, the codimension 0 leaf entities in process number $i$ form a
> @@ -2055,7 +2106,7 @@ the receiving end. Here the data is read from the message buffer and
> stored in the user's data structures. The message buffer is realized
> as an input stream delivering items of type \lstinline!DataType!. In
> the \lstinline!scatter! method it is up to the user how the data is to
> -be processed, e.~g.~one can simply overwrite (as is done here), add or
> +be processed, e.g.~one can simply overwrite (as is done here), add or
> compute a maximum.
>
> \section{Parallel finite volume scheme}
> @@ -2073,7 +2124,7 @@ numberstyle=\tiny, numbersep=5pt]{../parevolve.hh}
>
> The first difference to the sequential version is in line
> \ref{peh:assert} where it is checked that the grid provides an overlap
> -of at least one element. The overlap my be either of partition type
> +of at least one element. The overlap may be either of partition type
> \textit{overlap} or \textit{ghost}. The finite volume scheme itself
> only computes the updates for the elements with partition type
> \textit{interior}.
> @@ -2113,9 +2164,16 @@ Finally, we need a new main program, which is in the following listing:
> numberstyle=\tiny, numbersep=5pt]{../parfinitevolume.cc}
> \end{lst}
>
> -The only essential difference to the sequential program is in line
> +Essential differences to the sequential program are in line
> \ref{pfc:rank0} where the printing of the data of the current time
> -step is restricted to the process with rank 0.
> +step is restricted to the process with rank 0 and in line \ref{pfv:lb} where
> +the method \lstinline!loadBalance!\ is called on the grid. This method
> +re-partitions the grid in a way such that on every partition there is an equal
> +amount of grid elements. Some parallel grids like e.g.~\lstinline!YaspGrid!
> +however do not support load balancing and therefore need to start with a
> +sufficiently fine grid that allows a reasonable partition on all processes.
> +Note that during each global refinement step the overlap region of a
> +\lstinline!YaspGrid!\ grows and therefore the communication overhead increases.
>
>
> % \chapter{Input and Output}
> 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..982879b 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();
> @@ -48,14 +52,14 @@ void elementdata (const G& grid, const F& f)
> Dune::FieldVector<ct,dimworld> global = it->geometry().global(local);
>
> // evaluate functor and store value
> - c[mapper.map(*it)] = f(global); /*@\label{edh:feval}@*/
> + c[mapper.map(*it)] = f(global); /*@\label{edh:feval}@*/
> } /*@\label{edh:loop1}@*/
>
> // 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}@*/
> + vtkwriter.write("elementdata", Dune::VTKOptions::binaryappended); /*@\label{edh:vtk1}@*/
>
> // online visualization with Grape
> #if HAVE_GRAPE /*@\label{edh:grape0}@*/
> @@ -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..a8f81bc 100644
> --- a/finitevolume.cc
> +++ b/finitevolume.cc
> @@ -1,18 +1,18 @@
> -#include"config.h" // know what grids are present
> -#include<iostream> // for input/output to shell
> -#include<fstream> // for input/output to files
> -#include<vector> // STL vector class
> -#include<dune/grid/common/mcmgmapper.hh> // mapper class
> +#include "config.h" // know what grids are present
> +#include <iostream> // for input/output to shell
> +#include <fstream> // for input/output to files
> +#include <vector> // STL vector class
> +#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"transportproblem2.hh"
> -#include"initialize.hh"
> -#include"evolve.hh"
> +#include "vtkout.hh"
> +#include "unitcube.hh"
> +#include "transportproblem2.hh"
> +#include "initialize.hh"
> +#include "evolve.hh"
>
> //===============================================================
> // the time loop function working for all types of grids
> @@ -41,43 +41,44 @@ 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}@*/
> - {
> - // augment time step counter
> - ++k;
> + {
> + // augment time step counter
> + ++k;
>
> - // apply finite volume scheme
> - evolve(grid,mapper,c,t,dt);
> + // apply finite volume scheme
> + evolve(grid,mapper,c,t,dt);
>
> - // augment time
> - t += dt;
> + // augment time
> + t += dt;
>
> - // check if data should be written
> - if (t >= saveStep)
> - {
> - // write data
> - vtkout(grid,c,"concentration",counter);
> + // check if data should be written
> + if (t >= saveStep)
> + {
> + // write data
> + vtkout(grid,c,"concentration",counter,t);
>
> - // increase counter and saveStep for next interval
> - saveStep += saveInterval;
> - ++counter;
> - }
> + // increase counter and saveStep for next interval
> + saveStep += saveInterval;
> + ++counter;
> + }
>
> - // print info about time, timestep size and counter
> - std::cout << "s=" << grid.size(0) << " k=" << k << " t=" << t << " dt=" << dt << std::endl;
> - } /*@\label{fvc:loop1}@*/
> + // print info about time, timestep size and counter
> + 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..11ea169 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,31 +119,35 @@ 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 your map entry
> - IdType idi = idset.id(*it);
> - RestrictedValue& rv = restrictionmap[idi];
> -
> - // put your value in the map
> - if (it->isLeaf())
> - {
> - int indexi = mapper.map(*it);
> - rv.value = c[indexi];
> - rv.count = 1;
> - }
> -
> - // average in father
> - if (it.level()>0)
> - {
> - EntityPointer ep = it->father();
> - IdType idf = idset.id(*ep);
> - RestrictedValue& rvf = restrictionmap[idf];
> - rvf.value += rv.value/rv.count;
> - rvf.count += 1;
> - }
> - } /*@\label{fah:loop5}@*/
> + {
> + // 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);
> + RestrictedValue& rv = restrictionmap[idi];
> +
> + // put your value in the map
> + if (it->isLeaf())
> + {
> + int indexi = mapper.map(*it);
> + rv.value = c[indexi];
> + rv.count = 1;
> + }
> +
> + // average in father
> + if (it.level()>0)
> + {
> + EntityPointer ep = it->father();
> + IdType idf = idset.id(*ep);
> + RestrictedValue& rvf = restrictionmap[idf];
> + rvf.value += rv.value/rv.count;
> + rvf.count += 1;
> + }
> + } /*@\label{fah:loop5}@*/
> + }
> grid.preAdapt();
>
> // adapt mesh and mapper
> @@ -145,46 +157,50 @@ 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)
> - {
> - // get your id
> - IdType idi = idset.id(*it);
> -
> - // check map entry
> - typename std::map<IdType,RestrictedValue>::iterator rit = restrictionmap.find(idi);
> - if (rit!=restrictionmap.end())
> - {
> - // entry is in map, write in leaf
> - if (it->isLeaf())
> - {
> - int indexi = mapper.map(*it);
> - c[indexi] = rit->second.value/rit->second.count;
> - }
> - }
> - else
> - {
> - // value is not in map, interpolate
> - if (it.level()>0)
> - {
> - EntityPointer ep = it->father();
> - IdType idf = idset.id(*ep);
> - RestrictedValue& rvf = restrictionmap[idf];
> - if (it->isLeaf())
> - {
> - int indexi = mapper.map(*it);
> - c[indexi] = rvf.value/rvf.count;
> - }
> - else
> - {
> - // create new entry
> - RestrictedValue& rv = restrictionmap[idi];
> - rv.value = rvf.value/rvf.count;
> - rv.count = 1;
> - }
> - }
> - }
> - } /*@\label{fah:loop7}@*/
> + {
> + 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);
> + if (rit!=restrictionmap.end())
> + {
> + // entry is in map, write in leaf
> + if (it->isLeaf())
> + {
> + int indexi = mapper.map(*it);
> + c[indexi] = rit->second.value/rit->second.count;
> + }
> + }
> + else
> + {
> + // value is not in map, interpolate from father element
> + if (it.level()>0)
> + {
> + EntityPointer ep = it->father();
> + IdType idf = idset.id(*ep);
> + RestrictedValue& rvf = restrictionmap[idf];
> + if (it->isLeaf())
> + {
> + int indexi = mapper.map(*it);
> + c[indexi] = rvf.value/rvf.count;
> + }
> + else
> + {
> + // create new entry
> + RestrictedValue& rv = restrictionmap[idi];
> + rv.value = rvf.value/rvf.count;
> + rv.count = 1;
> + }
> + }
> + }
> + } /*@\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/integrateentity.hh b/integrateentity.hh
> index a155c36..3ab4d73 100644
> --- a/integrateentity.hh
> +++ b/integrateentity.hh
> @@ -1,8 +1,8 @@
> #ifndef DUNE_INTEGRATE_ENTITY_HH
> #define DUNE_INTEGRATE_ENTITY_HH
>
> -#include<dune/common/exceptions.hh>
> -#include<dune/grid/common/quadraturerules.hh>
> +#include <dune/common/exceptions.hh>
> +#include <dune/grid/common/quadraturerules.hh>
>
> //! compute integral of function over entity with given order
> template<class Iterator, class Functor>
> diff --git a/integration.cc b/integration.cc
> index 4bbab9b..00ef25d 100644
> --- a/integration.cc
> +++ b/integration.cc
> @@ -1,15 +1,15 @@
> // $Id$
>
> // Dune includes
> -#include"config.h" // file constructed by ./configure script
> -#include<dune/grid/sgrid.hh> // load sgrid definition
> +#include "config.h" // file constructed by ./configure script
> +#include <dune/grid/sgrid.hh> // load sgrid definition
> #include <dune/common/mpihelper.hh> // include mpi helper class
>
> // checks for defined gridtype and includes appropriate dgfparser implementation
> #include <dune/grid/io/file/dgfparser/dgfgridtype.hh>
>
> -#include"functors.hh"
> -#include"integrateentity.hh"
> +#include "functors.hh"
> +#include "integrateentity.hh"
>
> //! uniform refinement test
> template<class Grid>
> @@ -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..3c5b8fd 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,88 +66,88 @@ 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 Intersection &intersection = *is;
> -
> - // get geometry type of face
> - Dune::GeometryType gtf = intersection.intersectionSelfLocal().type();
> + const IntersectionIterator isend = gridView.iend(*it);
> + for( IntersectionIterator is = gridView.ibegin(*it); is != isend; ++is )
> + {
> + const Intersection &intersection = *is;
>
> - const Dune::ReferenceElement< ct, dim-1 > &refElement
> - = Dune::ReferenceElements< ct, dim-1 >::general( gtf );
> -
> - // center in face's reference element
> - const Dune::FieldVector< ct, dim-1 > &facelocal = refElement.position( 0, 0 );
> + // get geometry type of face
> + Dune::GeometryType gtf = intersection.intersectionSelfLocal().type();
>
> - // get normal vector scaled with volume
> - Dune::FieldVector< ct, dimworld > integrationOuterNormal
> - = intersection.integrationOuterNormal( facelocal );
> - integrationOuterNormal *= refElement.volume();
> + const Dune::ReferenceElement< ct, dim-1 > &refElement
> + = Dune::ReferenceElements< ct, dim-1 >::general( gtf );
>
> - // center of face in global coordinates
> - Dune::FieldVector< ct, dimworld > faceglobal
> - = intersection.intersectionGlobal().global( facelocal );
> + // center in face's reference element
> + const Dune::FieldVector< ct, dim-1 > &facelocal = refElement.position( 0, 0 );
>
> - // evaluate velocity at face center
> - Dune::FieldVector<double,dim> velocity = u(faceglobal,t);
> + // get normal vector scaled with volume
> + Dune::FieldVector< ct, dimworld > integrationOuterNormal
> + = intersection.integrationOuterNormal( facelocal );
> + integrationOuterNormal *= refElement.volume();
>
> - // compute factor occuring in flux formula
> - double factor = velocity*integrationOuterNormal/volume;
> + // center of face in global coordinates
> + Dune::FieldVector< ct, dimworld > faceglobal
> + = intersection.intersectionGlobal().global( facelocal );
>
> - // for time step calculation
> - if (factor>=0) sumfactor += factor;
> + // evaluate velocity at face center
> + Dune::FieldVector<double,dim> velocity = u(faceglobal,t);
>
> - // handle interior face
> - if( intersection.neighbor() )
> - {
> - // access neighbor
> - EntityPointer outside = intersection.outside();
> - int indexj = mapper.map(*outside);
> -
> - const int insideLevel = it->level();
> - const int outsideLevel = outside->level();
> -
> - // handle face from one side
> - if( (insideLevel > outsideLevel)
> - || ((insideLevel == outsideLevel) && (indexi < indexj)) )
> - {
> - // compute factor in neighbor
> - Dune::GeometryType nbgt = outside->type();
> - const Dune::FieldVector<ct,dim>&
> - nblocal = Dune::ReferenceElements<ct,dim>::general(nbgt).position(0,0);
> - double nbvolume = outside->geometry().integrationElement(nblocal)
> - *Dune::ReferenceElements<ct,dim>::general(nbgt).volume();
> - double nbfactor = velocity*integrationOuterNormal/nbvolume;
> -
> - if( factor < 0 ) // inflow
> - {
> - update[indexi] -= c[indexj]*factor;
> - update[indexj] += c[indexj]*nbfactor;
> - }
> - else // outflow
> - {
> - update[indexi] -= c[indexi]*factor;
> - update[indexj] += c[indexi]*nbfactor;
> - }
> - }
> - }
> + // compute factor occuring in flux formula
> + double factor = velocity*integrationOuterNormal/volume;
> +
> + // for time step calculation
> + if (factor>=0) sumfactor += factor;
>
> - // handle boundary face
> - if( intersection.boundary() )
> + // handle interior face
> + if( intersection.neighbor() )
> + {
> + // access neighbor
> + EntityPointer outside = intersection.outside();
> + int indexj = mapper.map(*outside);
> +
> + const int insideLevel = it->level();
> + const int outsideLevel = outside->level();
> +
> + // handle face from one side
> + if( (insideLevel > outsideLevel)
> + || ((insideLevel == outsideLevel) && (indexi < indexj)) )
> {
> - if( factor < 0 ) // inflow, apply boundary condition
> - update[indexi] -= b(faceglobal,t)*factor;
> + // compute factor in neighbor
> + Dune::GeometryType nbgt = outside->type();
> + const Dune::FieldVector<ct,dim>&
> + nblocal = Dune::ReferenceElements<ct,dim>::general(nbgt).position(0,0);
> + double nbvolume = outside->geometry().integrationElement(nblocal)
> + *Dune::ReferenceElements<ct,dim>::general(nbgt).volume();
> + double nbfactor = velocity*integrationOuterNormal/nbvolume;
> +
> + if( factor < 0 ) // inflow
> + {
> + update[indexi] -= c[indexj]*factor;
> + update[indexj] += c[indexj]*nbfactor;
> + }
> else // outflow
> + {
> update[indexi] -= c[indexi]*factor;
> + update[indexj] += c[indexi]*nbfactor;
> + }
> }
> - } // end all intersections
> + }
> +
> + // handle boundary face
> + if( intersection.boundary() )
> + {
> + if( factor < 0 ) // inflow, apply boundary condition
> + update[indexi] -= b(faceglobal,t)*factor;
> + else // outflow
> + update[indexi] -= c[indexi]*factor;
> + }
> + } // end all intersections
>
> - // compute dt restriction
> - if (it->partitionType()==Dune::InteriorEntity) /*@\label{peh:inter}@*/
> - dt = std::min(dt,1.0/sumfactor);
> + // compute dt restriction
> + if (it->partitionType()==Dune::InteriorEntity) /*@\label{peh:inter}@*/
> + dt = std::min(dt,1.0/sumfactor);
>
> - } // end grid traversal
> + } // end grid traversal
>
> // global min over all partitions
> dt = grid.comm().min(dt); /*@\label{peh:min}@*/
> @@ -151,12 +157,12 @@ void parevolve (const G& grid, const M& mapper, V& c, double t, double& dt)
> // exchange update
> VectorExchange<M,V> dh(mapper,update); /*@\label{peh:dist0}@*/
> grid.template
> - communicate<VectorExchange<M,V> >(dh,Dune::InteriorBorder_All_Interface,
> - Dune::ForwardCommunication); /*@\label{peh:dist1}@*/
> + communicate<VectorExchange<M,V> >(dh,Dune::InteriorBorder_All_Interface,
> + Dune::ForwardCommunication); /*@\label{peh:dist1}@*/
>
> // update the concentration vector
> for (unsigned int i=0; i<c.size(); ++i)
> - c[i] += dt*update[i];
> + c[i] += dt*update[i];
>
> return;
> }
> diff --git a/parfinitevolume.cc b/parfinitevolume.cc
> index 9f1f2f2..f32c75b 100644
> --- a/parfinitevolume.cc
> +++ b/parfinitevolume.cc
> @@ -1,16 +1,20 @@
> -#include"config.h" // know what grids are present
> -#include<iostream> // for input/output to shell
> -#include<fstream> // for input/output to files
> -#include<vector> // STL vector class
> -#include<dune/grid/common/mcmgmapper.hh> // mapper class
> +#include "config.h" // know what grids are present
> +#include <iostream> // for input/output to shell
> +#include <fstream> // for input/output to files
> +#include <vector> // STL vector class
> +#include <dune/grid/common/mcmgmapper.hh> // mapper class
> #include <dune/common/mpihelper.hh> // include mpi helper class
>
> -#include"vtkout.hh"
> -#include"unitcube.hh"
> -#include"transportproblem.hh"
> -#include"initialize.hh"
> -#include"parfvdatahandle.hh"
> -#include"parevolve.hh"
> +
> +// checks for defined gridtype and inlcudes appropriate dgfparser implementation
> +#include <dune/grid/io/file/dgfparser/dgfgridtype.hh>
> +
> +#include "vtkout.hh"
> +#include "unitcube.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);
> + 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..d83a5fa 100644
> --- a/traversal.cc
> +++ b/traversal.cc
> @@ -24,23 +24,28 @@ 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)
> - { /*@\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]
> - << std::endl;
> - count++; /*@\label{tc:count}@*/
> - } /*@\label{tc:forel1}@*/
> + 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().corner(0)
> + << std::endl;
> + count++; /*@\label{tc:count}@*/
> + } /*@\label{tc:forel1}@*/
>
> std::cout << "there are/is " << count << " leaf element(s)" << std::endl;
>
> @@ -50,48 +55,57 @@ 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)
> - {
> - Dune::GeometryType gt = it->type();
> - std::cout << "visiting " << gt
> - << " at " << it->geometry()[0]
> - << std::endl;
> - count++;
> - }
> + 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().corner(0)
> + << std::endl;
> + count++;
> + }
> std::cout << "there are/is " << count << " leaf vertices(s)"
> - << std::endl;
> + << std::endl;
>
> // Levelwise traversal of codim 0
> 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++)
> - {
> - count = 0;
> - for (ElementLevelIterator it = grid.template lbegin<0>(level);
> - it!=grid.template lend<0>(level); ++it)
> - {
> - Dune::GeometryType gt = it->type();
> - std::cout << "visiting " << gt
> - << " with first vertex at " << it->geometry()[0]
> - << std::endl;
> - count++;
> - }
> - std::cout << "there are/is " << count << " element(s) on level "
> - << level << std::endl;
> - std::cout << std::endl;
> - } /*@\label{tc:level1}@*/
> -} /*@\label{tc:tra1}@*/
> + {
> + // get the instance of the LeafGridView
> + LevelGridView levelView = grid.levelView(level);
> +
> + count = 0;
> + 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().corner(0)
> + << std::endl;
> + count++;
> + }
> + std::cout << "there are/is " << count << " element(s) on level "
> + << level << std::endl;
> + std::cout << std::endl;
> + } /*@\label{tc:level1}@*/
> +} /*@\label{tc:tra1}@*/
>
>
> int main(int argc, char **argv)
> diff --git a/unitcube.hh b/unitcube.hh
> index 6683610..dd594fa 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,14 +25,17 @@ 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"
> -#include"unitcube_sgrid.hh"
> -#include"unitcube_yaspgrid.hh"
> -#include"unitcube_uggrid.hh"
> -#include"unitcube_albertagrid.hh"
> -#include"unitcube_alugrid.hh"
> +#include "unitcube_onedgrid.hh"
> +#include "unitcube_sgrid.hh"
> +#include "unitcube_yaspgrid.hh"
> +#include "unitcube_uggrid.hh"
> +#include "unitcube_albertagrid.hh"
> +#include "unitcube_alugrid.hh"
>
> #endif
> diff --git a/unitcube_alugrid.hh b/unitcube_alugrid.hh
> index 86f64b6..a1887f3 100644
> --- a/unitcube_alugrid.hh
> +++ b/unitcube_alugrid.hh
> @@ -3,65 +3,89 @@
>
> #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>
> +class UnitCube<Dune::ALUSimplexGrid< 3, 3 >, 1 >
> +: public BasicUnitCube< 3 >
> {
> public:
> - typedef Dune::ALUSimplexGrid<3,3> GridType;
> + 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
> template<>
> -class UnitCube<Dune::ALUSimplexGrid<2,2>,1>
> +class UnitCube<Dune::ALUSimplexGrid<2, 2>,1>
> {
> public:
> - typedef Dune::ALUSimplexGrid<2,2> GridType;
> + typedef Dune::ALUSimplexGrid<2, 2> GridType;
>
> UnitCube () : filename("grids/2dsimplex.alu"), grid_(filename.c_str())
> {}
>
> GridType& grid ()
> {
> - return grid_;
> + return grid_;
> }
>
> -private:
> +private:
> std::string filename;
> GridType grid_;
> };
>
> // ALU3dGrid hexahedra specialization. Note: element type determined by type
> template<>
> -class UnitCube<Dune::ALUCubeGrid<3,3>,1>
> +class UnitCube<Dune::ALUCubeGrid< 3, 3 >, 1 >
> +: public BasicUnitCube< 3 >
> {
> public:
> - typedef Dune::ALUCubeGrid<3,3> GridType;
> + 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..be08005 100644
> --- a/visualization.cc
> +++ b/visualization.cc
> @@ -1,20 +1,20 @@
> // $Id$
>
> -#include"config.h"
> -#include<iostream>
> -#include<iomanip>
> -#include<stdio.h>
> +#include "config.h"
> +#include <iostream>
> +#include <iomanip>
> +#include <stdio.h>
> #include <dune/common/mpihelper.hh> // include mpi helper class
>
>
> -#include"elementdata.hh"
> -#include"vertexdata.hh"
> -#include"functors.hh"
> -#include"unitcube.hh"
> +#include "elementdata.hh"
> +#include "vertexdata.hh"
> +#include "functors.hh"
> +#include "unitcube.hh"
>
>
> #ifdef GRIDDIM
> -const int dimGrid = GRIDDIM;
> +const int dimGrid = 3;//GRIDDIM;
> #endif
>
>
> @@ -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 100755
> 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
> +
More information about the Dune
mailing list