[Dune] dune-grid-howto broken
Christian Engwer
christi at uni-hd.de
Fri Mar 6 17:53:29 CET 2009
Hi Martin,
On Fri, Mar 06, 2009 at 05:42:53PM +0100, Martin Drohmann wrote:
> Hello Christian,
>
> sorry for the whitespace changes. I built a new patch and excluded all lines
> with "whitespace-only" changes.
>
> Additionally I resorted the big patch, made small ones out of it and
> documented them individually. These patches are in the attached
> small_patches.tgz. In the header of each patch you will find a description
> of what I changed. I hope this helps.
On mall question is left concerning
0005-updated-the-parfinitevolume-example-now-using-dgf-f.patch
Is there a reason for having both version, UnitCube and DGFParser in
the example? Or can I just remove the UnitCube part?
Christian
>
> Bye,
> Martin
>
> 2009/3/6 Christian Engwer <christi at uni-hd.de>
>
> > Dear Martin,
> >
> > thank you for your work on the dune-grid-howto. The changes in the tex
> > file include some references to a basicunitcube.hh. The code you sent
> > include a lot of white space changes and it is not obvious what all
> > the changes are.
> >
> > Please update your patches to avoid whitespace changes and explain a
> > little bit more about the source-code changes.
> >
> > In the meantime I'll remove the paragraph quoting basicunitcube, in
> > order to make the howto work again.
> >
> > Thank you
> > Christian
> >
> >
>
> diff --git a/adaptivefinitevolume.cc b/adaptivefinitevolume.cc
> index a9f2241..e1425d4 100644
> --- a/adaptivefinitevolume.cc
> +++ b/adaptivefinitevolume.cc
> @@ -32,43 +32,6 @@ struct P0Layout
> };
>
> template<class G>
> -void gnuplot (G& grid, std::vector<double>& c)
> -{
> - // first we extract the dimensions of the grid
> - const int dim = G::dimension;
> - const int dimworld = G::dimensionworld;
> -
> - // type used for coordinates in the grid
> - // such a type is exported by every grid implementation
> - typedef typename G::ctype ct;
> -
> - // the grid has an iterator providing the access to
> - // all elements (better codim 0 entities) which are leafs
> - // of the refinement tree.
> - // Note the use of the typename keyword and the traits class
> - typedef typename G::template Codim<0>::LeafIterator ElementLeafIterator;
> -
> - // make a mapper for codim 0 entities in the leaf grid
> - Dune::LeafMultipleCodimMultipleGeomTypeMapper<G,P0Layout>
> - mapper(grid);
> -
> - // iterate through all entities of codim 0 at the leafs
> - int count = 0;
> - for (ElementLeafIterator it = grid.template leafbegin<0>();
> - it!=grid.template leafend<0>(); ++it)
> - {
> - Dune::GeometryType gt = it->type();
> - const Dune::FieldVector<ct,dim>&
> - local = Dune::ReferenceElements<ct,dim>::general(gt).position(0,0);
> - Dune::FieldVector<ct,dimworld>
> - global = it->geometry().global(local);
> - std::cout << global[0] << " " << c[mapper.map(*it)] << std::endl;
> - count++;
> - }
> -}
> -
> -
> -template<class G>
> void timeloop (G& grid, double tend, int lmin, int lmax)
> {
> // make a mapper for codim 0 entities in the leaf grid
> @@ -83,12 +46,12 @@ void timeloop (G& grid, double tend, int lmin, int lmax)
> for (int i=grid.maxLevel(); i<lmax; i++)
> {
> if (grid.maxLevel()>=lmax) break;
> - finitevolumeadapt(grid,mapper,c,lmin,lmax,0);
> + finitevolumeadapt(grid,mapper,c,lmin,lmax,0); /*@\label{afv:in}@*/
> initialize(grid,mapper,c);
> }
>
> // write initial data
> - vtkout(grid,c,"concentration",0);
> + vtkout(grid,c,"concentration",0,0);
>
> // variables for time, timestep etc.
> double dt, t=0;
> @@ -113,7 +76,7 @@ void timeloop (G& grid, double tend, int lmin, int lmax)
> if (t >= saveStep)
> {
> // write data
> - vtkout(grid,c,"concentration",counter);
> + vtkout(grid,c,"concentration",counter,t);
>
> // increase counter and saveStep for next interval
> saveStep += saveInterval;
> @@ -121,17 +84,20 @@ void timeloop (G& grid, double tend, int lmin, int lmax)
> }
>
> // print info about time, timestep size and counter
> - std::cout << "s=" << grid.size(0) << " k=" << k << " t=" << t << " dt=" << dt << std::endl;
> + std::cout << "s=" << grid.size(0)
> + << " k=" << k << " t=" << t << " dt=" << dt << std::endl;
>
> // for unstructured grids call adaptation algorithm
> if( Dune :: Capabilities :: IsUnstructured<G> :: v )
> {
> - finitevolumeadapt(grid,mapper,c,lmin,lmax,k);
> + finitevolumeadapt(grid,mapper,c,lmin,lmax,k); /*@\label{afv:ad}@*/
> }
> }
>
> // write last time step
> - vtkout(grid,c,"concentration",counter);
> + vtkout(grid,c,"concentration",counter,tend);
> +
> + // write
> }
>
> //===============================================================
> diff --git a/adaptiveintegration.cc b/adaptiveintegration.cc
> index ded8d9c..5396c44 100644
> --- a/adaptiveintegration.cc
> +++ b/adaptiveintegration.cc
> @@ -19,8 +19,13 @@
> template<class Grid, class Functor>
> void adaptiveintegration (Grid& grid, const Functor& f)
> {
> + // get grid view type for leaf grid part
> + typedef typename Grid::LeafGridView GridView;
> // get iterator type
> - typedef typename Grid::template Codim<0>::LeafIterator ElementLeafIterator;
> + typedef typename GridView::template Codim<0>::Iterator ElementLeafIterator;
> +
> + // get grid view on leaf part
> + GridView gridView = grid.leafView();
>
> // algorithm parameters
> const double tol=1E-8;
> @@ -33,8 +38,8 @@ void adaptiveintegration (Grid& grid, const Functor& f)
> {
> // compute integral on current mesh
> double value=0; /*@\label{aic:int0}@*/
> - for (ElementLeafIterator it = grid.template leafbegin<0>();
> - it!=grid.template leafend<0>(); ++it)
> + for (ElementLeafIterator it = gridView.template begin<0>();
> + it!=gridView.template end<0>(); ++it)
> value += integrateentity(it,f,highorder); /*@\label{aic:int1}@*/
>
> // print result
> @@ -87,8 +92,8 @@ void adaptiveintegration (Grid& grid, const Functor& f)
> double kappa = std::min(maxextrapolatederror,0.5*maxerror); /*@\label{aic:kappa1}@*/
>
> // mark elements for refinement
> - for (ElementLeafIterator it = grid.template leafbegin<0>(); /*@\label{aic:mark0}@*/
> - it!=grid.template leafend<0>(); ++it)
> + for (ElementLeafIterator it = gridView.template begin<0>(); /*@\label{aic:mark0}@*/
> + it!=gridView.template end<0>(); ++it)
> {
> double lowresult=integrateentity(it,f,loworder);
> double highresult=integrateentity(it,f,highorder);
> @@ -103,7 +108,7 @@ void adaptiveintegration (Grid& grid, const Functor& f)
> }
>
> // write grid in VTK format
> - Dune::VTKWriter<typename Grid::LeafGridView> vtkwriter(grid.leafView());
> + Dune::VTKWriter<typename Grid::LeafGridView> vtkwriter(gridView);
> vtkwriter.write("adaptivegrid",Dune::VTKOptions::binaryappended);
> }
>
> diff --git a/basicunitcube.hh b/basicunitcube.hh
> new file mode 100644
> index 0000000..65a30dd
> --- /dev/null
> +++ b/basicunitcube.hh
> @@ -0,0 +1,111 @@
> +#ifndef BASICUNITCUBE_HH
> +#define BASICUNITCUBE_HH
> +
> +#include <dune/grid/common/gridfactory.hh>
> +
> +// declaration of a basic unit cube that uses the GridFactory
> +template< int dim >
> +class BasicUnitCube;
> +
> +// unit cube in two dimensions with 2 variants: triangle and rectangle elements
> +template<>
> +class BasicUnitCube< 2 >
> +{
> +protected:
> + template< class Grid >
> + static void insertVertices ( Dune::GridFactory< Grid > &factory )
> + {
> + Dune::FieldVector<double,2> pos;
> +
> + pos[0] = 0; pos[1] = 0;
> + factory.insertVertex(pos); /*@\label{uc:iv}@*/
> +
> + pos[0] = 1; pos[1] = 0;
> + factory.insertVertex(pos);
> +
> + pos[0] = 0; pos[1] = 1;
> + factory.insertVertex(pos);
> +
> + pos[0] = 1; pos[1] = 1;
> + factory.insertVertex(pos);
> + }
> +
> + template< class Grid >
> + static void insertSimplices ( Dune::GridFactory< Grid > &factory )
> + {
> + const Dune::GeometryType type( Dune::GeometryType::simplex, 2 );
> + std::vector< unsigned int > cornerIDs( 3 );
> +
> + cornerIDs[0] = 0; cornerIDs[1] = 1; cornerIDs[2] = 2;
> + factory.insertElement( type, cornerIDs ); /*@\label{uc:ie}@*/
> +
> + cornerIDs[0] = 2; cornerIDs[1] = 1; cornerIDs[2] = 3;
> + factory.insertElement( type, cornerIDs );
> + }
> +
> + template< class Grid >
> + static void insertCubes ( Dune::GridFactory< Grid > &factory )
> + {
> + const Dune::GeometryType type( Dune::GeometryType::cube, 2 );
> + std::vector< unsigned int > cornerIDs( 4 );
> + for( int i = 0; i < 4; ++i )
> + cornerIDs[ i ] = i;
> + factory.insertElement( type, cornerIDs );
> + }
> +};
> +
> +// unit cube in 3 dimensions with two variants: tetraheda and hexahedra
> +template<>
> +class BasicUnitCube< 3 >
> +{
> +protected:
> + template< class Grid >
> + static void insertVertices ( Dune::GridFactory< Grid > &factory )
> + {
> + Dune::FieldVector< double, 3 > pos;
> +
> + pos[0] = 0; pos[1] = 0; pos[2] = 0; factory.insertVertex(pos);
> + pos[0] = 1; pos[1] = 0; pos[2] = 0; factory.insertVertex(pos);
> + pos[0] = 0; pos[1] = 1; pos[2] = 0; factory.insertVertex(pos);
> + pos[0] = 1; pos[1] = 1; pos[2] = 0; factory.insertVertex(pos);
> + pos[0] = 0; pos[1] = 0; pos[2] = 1; factory.insertVertex(pos);
> + pos[0] = 1; pos[1] = 0; pos[2] = 1; factory.insertVertex(pos);
> + pos[0] = 0; pos[1] = 1; pos[2] = 1; factory.insertVertex(pos);
> + pos[0] = 1; pos[1] = 1; pos[2] = 1; factory.insertVertex(pos);
> + }
> +
> + template< class Grid >
> + static void insertSimplices ( Dune::GridFactory< Grid > &factory )
> + {
> + const Dune::GeometryType type( Dune::GeometryType::simplex, 3 );
> + std::vector< unsigned int > cornerIDs( 4 );
> +
> + cornerIDs[0] = 0; cornerIDs[1] = 1; cornerIDs[2] = 2; cornerIDs[3] = 4;
> + factory.insertElement( type, cornerIDs );
> +
> + cornerIDs[0] = 1; cornerIDs[1] = 3; cornerIDs[2] = 2; cornerIDs[3] = 7;
> + factory.insertElement( type, cornerIDs );
> +
> + cornerIDs[0] = 1; cornerIDs[1] = 7; cornerIDs[2] = 2; cornerIDs[3] = 4;
> + factory.insertElement( type, cornerIDs );
> +
> + cornerIDs[0] = 1; cornerIDs[1] = 7; cornerIDs[2] = 4; cornerIDs[3] = 5;
> + factory.insertElement( type, cornerIDs );
> +
> + cornerIDs[0] = 4; cornerIDs[1] = 7; cornerIDs[2] = 2; cornerIDs[3] = 6;
> + factory.insertElement( type, cornerIDs );
> + }
> +
> + template< class Grid >
> + static void insertCubes ( Dune::GridFactory< Grid > &factory )
> + {
> + const Dune::GeometryType type( Dune::GeometryType::cube, 3 );
> + std::vector< unsigned int > cornerIDs( 8 );
> + for( int i = 0; i < 8; ++i )
> + cornerIDs[ i ] = i;
> + factory.insertElement( type, cornerIDs );
> + }
> +};
> +
> +#endif /*BASICUNITCUBE_HH*/
> +
> diff --git a/dune.module b/dune.module
> index ef05236..c5e2876 100644
> --- a/dune.module
> +++ b/dune.module
> @@ -1,5 +1,5 @@
> # parameters for dune control
> Module: dune-grid-howto
> Maintainer: dune at dune-project.org
> -Version: 1.2svn
> +Version: 1.2
> Depends: dune-common dune-grid
> diff --git a/elementdata.hh b/elementdata.hh
> index 1696aa9..a4c8b88 100644
> --- a/elementdata.hh
> +++ b/elementdata.hh
> @@ -24,7 +24,11 @@ void elementdata (const G& grid, const F& f)
> const int dim = G::dimension;
> const int dimworld = G::dimensionworld;
> typedef typename G::ctype ct;
> - typedef typename G::template Codim<0>::LeafIterator ElementLeafIterator;
> + typedef typename G::LeafGridView GridView;
> + typedef typename GridView::template Codim<0>::Iterator ElementLeafIterator;
> +
> + // get grid view on leaf part
> + GridView gridView = grid.leafView();
>
> // make a mapper for codim 0 entities in the leaf grid
> Dune::LeafMultipleCodimMultipleGeomTypeMapper<G,P0Layout>
> @@ -34,8 +38,8 @@ void elementdata (const G& grid, const F& f)
> std::vector<double> c(mapper.size()); /*@\label{edh:c}@*/
>
> // iterate through all entities of codim 0 at the leafs
> - for (ElementLeafIterator it = grid.template leafbegin<0>(); /*@\label{edh:loop0}@*/
> - it!=grid.template leafend<0>(); ++it)
> + for (ElementLeafIterator it = gridView.template begin<0>(); /*@\label{edh:loop0}@*/
> + it!=gridView.template end<0>(); ++it)
> {
> // cell geometry type
> Dune::GeometryType gt = it->type();
> @@ -53,7 +57,7 @@ void elementdata (const G& grid, const F& f)
>
> // generate a VTK file
> // Dune::LeafP0Function<G,double> cc(grid,c);
> - Dune::VTKWriter<typename G::LeafGridView> vtkwriter(grid.leafView()); /*@\label{edh:vtk0}@*/
> + Dune::VTKWriter<typename G::LeafGridView> vtkwriter(gridView); /*@\label{edh:vtk0}@*/
> vtkwriter.addCellData(c,"data");
> vtkwriter.write("elementdata",Dune::VTKOptions::binaryappended); /*@\label{edh:vtk1}@*/
>
> @@ -67,7 +71,7 @@ void elementdata (const G& grid, const F& f)
> // display data
> grape.displayVector("concentration", // name of data that appears in grape
> c, // data vector
> - grid.leafIndexSet(), // used index set
> + gridView.indexSet(), // used index set
> polynomialOrder, // polynomial order of data
> dimRange); // dimRange of data
> }
> diff --git a/evolve.hh b/evolve.hh
> index 1b4262a..af0aac4 100644
> --- a/evolve.hh
> +++ b/evolve.hh
> @@ -10,15 +10,21 @@ void evolve (const G& grid, const M& mapper, V& c, double t, double& dt)
> // type used for coordinates in the grid
> typedef typename G::ctype ct;
>
> - // iterator type
> - typedef typename G::template Codim<0>::LeafIterator LeafIterator;
> + // type of grid view on leaf part
> + typedef typename G::LeafGridView GridView;
> +
> + // element iterator type
> + typedef typename GridView::template Codim<0>::Iterator LeafIterator;
>
> // intersection iterator type
> - typedef typename G::template Codim<0>::LeafIntersectionIterator IntersectionIterator;
> + typedef typename GridView::IntersectionIterator IntersectionIterator;
>
> // entity pointer type
> typedef typename G::template Codim<0>::EntityPointer EntityPointer;
>
> + // get grid view on leaf part
> + GridView gridView = grid.leafView();
> +
> // allocate a temporary vector for the update
> V update(c.size()); /*@\label{evh:update}@*/
> for (typename V::size_type i=0; i<c.size(); i++) update[i] = 0;
> @@ -27,8 +33,8 @@ void evolve (const G& grid, const M& mapper, V& c, double t, double& dt)
> dt = 1E100;
>
> // compute update vector and optimum dt in one grid traversal
> - LeafIterator endit = grid.template leafend<0>(); /*@\label{evh:loop0}@*/
> - for (LeafIterator it = grid.template leafbegin<0>(); it!=endit; ++it)
> + LeafIterator endit = gridView.template end<0>(); /*@\label{evh:loop0}@*/
> + for (LeafIterator it = gridView.template begin<0>(); it!=endit; ++it)
> {
> // cell geometry type
> Dune::GeometryType gt = it->type();
> @@ -52,8 +58,8 @@ void evolve (const G& grid, const M& mapper, V& c, double t, double& dt)
> double sumfactor = 0.0;
>
> // run through all intersections with neighbors and boundary
> - IntersectionIterator isend = it->ileafend(); /*@\label{evh:flux0}@*/
> - for (IntersectionIterator is = it->ileafbegin(); is!=isend; ++is)
> + IntersectionIterator isend = gridView.iend(*it); /*@\label{evh:flux0}@*/
> + for (IntersectionIterator is = gridView.ibegin(*it); is!=isend; ++is)
> {
> // get geometry type of face
> Dune::GeometryType gtf = is->intersectionSelfLocal().type();
> diff --git a/finitevolume.cc b/finitevolume.cc
> index 315c72e..be9ab43 100644
> --- a/finitevolume.cc
> +++ b/finitevolume.cc
> @@ -41,14 +41,14 @@ void timeloop (const G& grid, double tend)
>
> // initialize concentration with initial values
> initialize(grid,mapper,c); /*@\label{fvc:init}@*/
> - vtkout(grid,c,"concentration",0);
> + vtkout(grid,c,"concentration",0,0.0);
>
> // now do the time steps
> double t=0,dt;
> int k=0;
> const double saveInterval = 0.1;
> double saveStep = 0.1;
> - int counter = 0;
> + int counter = 1;
>
> while (t<tend) /*@\label{fvc:loop0}@*/
> {
> @@ -65,7 +65,7 @@ void timeloop (const G& grid, double tend)
> if (t >= saveStep)
> {
> // write data
> - vtkout(grid,c,"concentration",counter);
> + vtkout(grid,c,"concentration",counter,t);
>
> // increase counter and saveStep for next interval
> saveStep += saveInterval;
> @@ -73,11 +73,12 @@ void timeloop (const G& grid, double tend)
> }
>
> // print info about time, timestep size and counter
> - std::cout << "s=" << grid.size(0) << " k=" << k << " t=" << t << " dt=" << dt << std::endl;
> + std::cout << "s=" << grid.size(0)
> + << " k=" << k << " t=" << t << " dt=" << dt << std::endl;
> } /*@\label{fvc:loop1}@*/
>
> // output results
> - vtkout(grid,c,"concentration",counter); /*@\label{fvc:file}@*/
> + vtkout(grid,c,"concentration",counter,tend); /*@\label{fvc:file}@*/
> }
>
> //===============================================================
> @@ -93,7 +94,7 @@ int main (int argc , char ** argv)
> try {
> using namespace Dune;
>
> - // use unitcube from grids
> + // use unitcube from dgf grids
> std::stringstream dgfFileName;
> dgfFileName << "grids/unitcube" << GridType :: dimension << ".dgf";
>
> diff --git a/finitevolumeadapt.hh b/finitevolumeadapt.hh
> index 42b0f35..dd79b49 100644
> --- a/finitevolumeadapt.hh
> +++ b/finitevolumeadapt.hh
> @@ -21,27 +21,35 @@ bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax, int k)
> // type used for coordinates in the grid
> typedef typename G::ctype ct;
>
> + // grid view types
> + typedef typename G::LeafGridView LeafGridView;
> + typedef typename G::LevelGridView LevelGridView;
> +
> // iterator types
> - typedef typename G::template Codim<0>::LeafIterator LeafIterator;
> - typedef typename G::template Codim<0>::LevelIterator LevelIterator;
> + typedef typename LeafGridView::template Codim<0>::Iterator LeafIterator;
> + typedef typename LevelGridView::template Codim<0>::Iterator LevelIterator;
>
> // entity and entity pointer
> typedef typename G::template Codim<0>::Entity Entity;
> typedef typename G::template Codim<0>::EntityPointer EntityPointer;
>
> // intersection iterator type
> - typedef typename G::template Codim<0>::LeafIntersectionIterator IntersectionIterator;
> + typedef typename LeafGridView::IntersectionIterator LeafIntersectionIterator;
>
> - // global id set types
> + // global id set types, local means that the numbering is unique in a single process only.
> typedef typename G::template Codim<0>::LocalIdSet IdSet;
> + // type for the index set, note that this is _not_ an integer
> typedef typename IdSet::IdType IdType;
>
> + // get grid view on leaf grid
> + LeafGridView leafView = grid.leafView();
> +
> // compute cell indicators
> V indicator(c.size(),-1E100);
> double globalmax = -1E100;
> double globalmin = 1E100;
> - for (LeafIterator it = grid.template leafbegin<0>(); /*@\label{fah:loop0}@*/
> - it!=grid.template leafend<0>(); ++it)
> + for (LeafIterator it = leafView.template begin<0>(); /*@\label{fah:loop0}@*/
> + it!=leafView.template end<0>(); ++it)
> {
> // my index
> int indexi = mapper.map(*it);
> @@ -50,10 +58,10 @@ bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax, int k)
> globalmax = std::max(globalmax,c[indexi]);
> globalmin = std::min(globalmin,c[indexi]);
>
> - IntersectionIterator isend = it->ileafend();
> - for (IntersectionIterator is = it->ileafbegin(); is!=isend; ++is)
> + LeafIntersectionIterator isend = leafView.iend(*it);
> + for (LeafIntersectionIterator is = leafView.ibegin(*it); is!=isend; ++is)
> {
> - const typename IntersectionIterator::Intersection &intersection = *is;
> + const typename LeafIntersectionIterator::Intersection &intersection = *is;
> if( !intersection.neighbor() )
> continue;
>
> @@ -76,8 +84,8 @@ bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax, int k)
> // mark cells for refinement/coarsening
> double globaldelta = globalmax-globalmin;
> int marked=0;
> - for (LeafIterator it = grid.template leafbegin<0>(); /*@\label{fah:loop2}@*/
> - it!=grid.template leafend<0>(); ++it)
> + for (LeafIterator it = leafView.template begin<0>(); /*@\label{fah:loop2}@*/
> + it!=leafView.template end<0>(); ++it)
> {
> if (indicator[mapper.map(*it)]>refinetol*globaldelta
> && (it.level()<lmax || !it->isRegular()))
> @@ -85,10 +93,10 @@ bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax, int k)
> const Entity &entity = *it;
> grid.mark( 1, entity );
> ++marked;
> - IntersectionIterator isend = entity.ileafend();
> - for( IntersectionIterator is = entity.ileafbegin(); is != isend; ++is )
> + LeafIntersectionIterator isend = leafView.iend(entity);
> + for( LeafIntersectionIterator is = leafView.ibegin(entity); is != isend; ++is )
> {
> - const typename IntersectionIterator::Intersection intersection = *is;
> + const typename LeafIntersectionIterator::Intersection intersection = *is;
> if( !intersection.neighbor() )
> continue;
>
> @@ -111,8 +119,11 @@ bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax, int k)
> std::map<IdType,RestrictedValue> restrictionmap; // restricted concentration /*@\label{fah:loop4}@*/
> const IdSet& idset = grid.localIdSet();
> for (int level=grid.maxLevel(); level>=0; level--)
> - for (LevelIterator it = grid.template lbegin<0>(level);
> - it!=grid.template lend<0>(level); ++it)
> + {
> + // get grid view on level grid
> + LevelGridView levelView = grid.levelView(level);
> + for (LevelIterator it = levelView.template begin<0>();
> + it!=levelView.template end<0>(); ++it)
> {
> // get your map entry
> IdType idi = idset.id(*it);
> @@ -136,6 +147,7 @@ bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax, int k)
> rvf.count += 1;
> }
> } /*@\label{fah:loop5}@*/
> + }
> grid.preAdapt();
>
> // adapt mesh and mapper
> @@ -145,14 +157,17 @@ bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax, int k)
>
> // interpolate new cells, restrict coarsened cells
> for (int level=0; level<=grid.maxLevel(); level++) /*@\label{fah:loop6}@*/
> - for (LevelIterator it = grid.template lbegin<0>(level);
> - it!=grid.template lend<0>(level); ++it)
> + {
> + LevelGridView levelView = grid.levelView(level);
> + for (LevelIterator it = levelView.template begin<0>();
> + it!=levelView.template end<0>(); ++it)
> {
> // get your id
> IdType idi = idset.id(*it);
>
> // check map entry
> - typename std::map<IdType,RestrictedValue>::iterator rit = restrictionmap.find(idi);
> + typename std::map<IdType,RestrictedValue>::iterator rit
> + = restrictionmap.find(idi);
> if (rit!=restrictionmap.end())
> {
> // entry is in map, write in leaf
> @@ -164,7 +179,7 @@ bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax, int k)
> }
> else
> {
> - // value is not in map, interpolate
> + // value is not in map, interpolate from father element
> if (it.level()>0)
> {
> EntityPointer ep = it->father();
> @@ -185,6 +200,7 @@ bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax, int k)
> }
> }
> } /*@\label{fah:loop7}@*/
> + }
> grid.postAdapt();
>
> return rv;
> diff --git a/grids/3dgrid.al b/grids/3dgrid.al
> deleted file mode 100644
> index 2bde72f..0000000
> --- a/grids/3dgrid.al
> +++ /dev/null
> @@ -1,39 +0,0 @@
> -DIM: 3
> -DIM_OF_WORLD: 3
> -
> -number of vertices: 8
> -number of elements: 6
> -
> -vertex coordinates:
> - 0.0 0.0 0.0
> - 1.0 0.0 0.0
> - 0.0 0.0 1.0
> - 1.0 0.0 1.0
> - 1.0 1.0 0.0
> - 1.0 1.0 1.0
> - 0.0 1.0 0.0
> - 0.0 1.0 1.0
> -
> -element vertices:
> - 0 5 4 1
> - 0 5 3 1
> - 0 5 3 2
> - 0 5 4 6
> - 0 5 7 6
> - 0 5 7 2
> -
> -element boundaries:
> - 1 1 0 0
> - 1 1 0 0
> - 1 1 0 0
> - 1 1 0 0
> - 1 1 0 0
> - 1 1 0 0
> -
> -element neighbours:
> - -1 -1 1 3
> - -1 -1 0 2
> - -1 -1 5 1
> - -1 -1 4 0
> - -1 -1 3 5
> - -1 -1 2 4
> diff --git a/grids/Makefile.am b/grids/Makefile.am
> index 6527044..8e7c6a0 100644
> --- a/grids/Makefile.am
> +++ b/grids/Makefile.am
> @@ -1,4 +1,4 @@
> -EXTRA_DIST = 2dgrid.al 2dsimplex.alu 3dgrid.al cube.hexa cube.tetra \
> +EXTRA_DIST = 2dgrid.al 2dsimplex.alu \
> unitcube1.dgf unitcube2.dgf unitcube3.dgf
>
> griddir = $(datadir)/doc/dune-grid-howto/examples/grids
> diff --git a/grids/cube.hexa b/grids/cube.hexa
> deleted file mode 100644
> index 2085520..0000000
> --- a/grids/cube.hexa
> +++ /dev/null
> @@ -1,31 +0,0 @@
> -!Hexahedra
> -
> -8
> -0.000000 0.000000 0.000000
> -1.000000 0.000000 0.000000
> -1.000000 1.000000 0.000000
> -0.000000 1.000000 0.000000
> -0.000000 0.000000 1.000000
> -1.000000 0.000000 1.000000
> -1.000000 1.000000 1.000000
> -0.000000 1.000000 1.000000
> -
> -1
> -0 1 2 3 4 5 6 7
> -
> -6
> --2 4 0 3 7 4
> --3 4 1 5 6 2
> --1 4 0 4 5 1
> --1 4 3 2 6 7
> --1 4 0 1 2 3
> --1 4 5 4 7 6
> -
> -0 -1
> -1 -1
> -2 -1
> -3 -1
> -4 -1
> -5 -1
> -6 -1
> -7 -1
> diff --git a/grids/cube.tetra b/grids/cube.tetra
> deleted file mode 100644
> index c9eeb98..0000000
> --- a/grids/cube.tetra
> +++ /dev/null
> @@ -1,41 +0,0 @@
> -!Tetrahedra
> -8
> -0.000000 0.000000 0.000000
> -1.000000 0.000000 0.000000
> -0.000000 1.000000 1.000000
> -0.000000 0.000000 1.000000
> -0.000000 1.000000 0.000000
> -1.000000 0.000000 1.000000
> -1.000000 1.000000 0.000000
> -1.000000 1.000000 1.000000
> -
> -6
> -0 1 2 3
> -0 2 1 4
> -1 5 2 3
> -1 6 4 2
> -1 6 2 7
> -1 7 2 5
> -
> -12
> --2 3 3 2 0
> --2 3 1 3 0
> --1 3 4 1 0
> --2 3 2 4 0
> --1 3 2 3 5
> --2 3 5 3 1
> --2 3 4 2 6
> --1 3 4 6 1
> --2 3 2 7 6
> --2 3 6 7 1
> --1 3 2 5 7
> --2 3 7 5 1
> -
> -0 -1
> -1 -1
> -2 -1
> -3 -1
> -4 -1
> -5 -1
> -6 -1
> -7 -1
> diff --git a/initialize.hh b/initialize.hh
> index cc4872a..96bbbd2 100644
> --- a/initialize.hh
> +++ b/initialize.hh
> @@ -11,12 +11,18 @@ void initialize (const G& grid, const M& mapper, V& c)
> // type used for coordinates in the grid
> typedef typename G::ctype ct;
>
> + // type of grid view on leaf part
> + typedef typename G::LeafGridView GridView;
> +
> // leaf iterator type
> - typedef typename G::template Codim<0>::LeafIterator LeafIterator;
> + typedef typename GridView::template Codim<0>::Iterator LeafIterator;
> +
> + // get grid view on leaf part
> + GridView gridView = grid.leafView();
>
> // iterate through leaf grid an evaluate c0 at cell center
> - LeafIterator endit = grid.template leafend<0>();
> - for (LeafIterator it = grid.template leafbegin<0>(); it!=endit; ++it)
> + LeafIterator endit = gridView.template end<0>();
> + for (LeafIterator it = gridView.template begin<0>(); it!=endit; ++it)
> {
> // get geometry type
> Dune::GeometryType gt = it->type();
> diff --git a/integration.cc b/integration.cc
> index 4bbab9b..8de5237 100644
> --- a/integration.cc
> +++ b/integration.cc
> @@ -18,8 +18,14 @@ void uniformintegration (Grid& grid)
> // function to integrate
> Exp<typename Grid::ctype,Grid::dimension> f;
>
> + // get GridView on leaf grid - type
> + typedef typename Grid :: LeafGridView GridView;
> +
> + // get GridView instance
> + GridView gridView = grid.leafView();
> +
> // get iterator type
> - typedef typename Grid::template Codim<0>::LeafIterator LeafIterator;
> + typedef typename GridView :: template Codim<0> :: Iterator LeafIterator;
>
> // loop over grid sequence
> double oldvalue=1E100;
> @@ -27,8 +33,8 @@ void uniformintegration (Grid& grid)
> {
> // compute integral with some order
> double value = 0.0;
> - LeafIterator eendit = grid.template leafend<0>();
> - for (LeafIterator it = grid.template leafbegin<0>(); it!=eendit; ++it)
> + LeafIterator eendit = gridView.template end<0>();
> + for (LeafIterator it = gridView.template begin<0>(); it!=eendit; ++it)
> value += integrateentity(it,f,1); /*@\label{ic:call}@*/
>
> // print result and error estimate
> diff --git a/parevolve.hh b/parevolve.hh
> index c2c0f96..6d724c9 100644
> --- a/parevolve.hh
> +++ b/parevolve.hh
> @@ -13,12 +13,15 @@ void parevolve (const G& grid, const M& mapper, V& c, double t, double& dt)
> // type used for coordinates in the grid
> typedef typename G::ctype ct;
>
> + // type for grid view on leaf part
> + typedef typename G::LeafGridView GridView;
> +
> // iterator type
> - typedef typename G::template Codim<0>::
> - template Partition<Dune::All_Partition>::LeafIterator LeafIterator; /*@\label{peh:pit}@*/
> + typedef typename GridView::template Codim<0>::
> + template Partition<Dune::All_Partition>::Iterator LeafIterator; /*@\label{peh:pit}@*/
>
> // intersection iterator type
> - typedef typename G::template Codim<0>::LeafIntersectionIterator IntersectionIterator;
> + typedef typename GridView::IntersectionIterator IntersectionIterator;
>
> // type of intersection
> typedef typename IntersectionIterator::Intersection Intersection;
> @@ -33,10 +36,13 @@ void parevolve (const G& grid, const M& mapper, V& c, double t, double& dt)
> // initialize dt very large
> dt = 1E100;
>
> + // get grid view instance on leaf grid
> + GridView gridView = grid.leafView();
> +
> // compute update vector and optimum dt in one grid traversal
> // iterate over all entities, but update is only used on interior entities
> - LeafIterator endit = grid.template leafend<0,Dune::All_Partition>(); /*@\label{peh:end}@*/
> - for (LeafIterator it = grid.template leafbegin<0,Dune::All_Partition>(); it!=endit; ++it) /*@\label{peh:begin}@*/
> + LeafIterator endit = gridView.template end<0,Dune::All_Partition>(); /*@\label{peh:end}@*/
> + for (LeafIterator it = gridView.template begin<0,Dune::All_Partition>(); it!=endit; ++it) /*@\label{peh:begin}@*/
> {
> // cell geometry type
> Dune::GeometryType gt = it->type();
> @@ -60,8 +66,8 @@ void parevolve (const G& grid, const M& mapper, V& c, double t, double& dt)
> double sumfactor = 0.0;
>
> // run through all intersections with neighbors and boundary
> - const IntersectionIterator isend = it->ileafend();
> - for( IntersectionIterator is = it->ileafbegin(); is != isend; ++is )
> + const IntersectionIterator isend = gridView.iend(*it);
> + for( IntersectionIterator is = gridView.ibegin(*it); is != isend; ++is )
> {
> const Intersection &intersection = *is;
>
> diff --git a/parfinitevolume.cc b/parfinitevolume.cc
> index 9f1f2f2..cc801de 100644
> --- a/parfinitevolume.cc
> +++ b/parfinitevolume.cc
> @@ -5,9 +5,13 @@
> #include<dune/grid/common/mcmgmapper.hh> // mapper class
> #include <dune/common/mpihelper.hh> // include mpi helper class
>
> +
> +// checks for defined gridtype and inlcudes appropriate dgfparser implementation
> +#include <dune/grid/io/file/dgfparser/dgfgridtype.hh>
> +
> #include"vtkout.hh"
> #include"unitcube.hh"
> -#include"transportproblem.hh"
> +#include "transportproblem2.hh"
> #include"initialize.hh"
> #include"parfvdatahandle.hh"
> #include"parevolve.hh"
> @@ -40,21 +44,41 @@ void partimeloop (const G& grid, double tend)
>
> // initialize concentration with initial values
> initialize(grid,mapper,c);
> - vtkout(grid,c,"pconc",0);
> + vtkout(grid,c,"pconc",0,0.0,grid.comm().rank());
>
> // now do the time steps
> double t=0,dt;
> int k=0;
> + const double saveInterval = 0.1;
> + double saveStep = 0.1;
> + int counter = 1;
> while (t<tend)
> {
> + // augment time step counter
> k++;
> +
> + // apply finite volume scheme
> parevolve(grid,mapper,c,t,dt);
> +
> + // augment time
> t += dt;
> +
> + // check if data should be written
> + if (t >= saveStep)
> + {
> + // write data
> + vtkout(grid,c,"pconc",counter,t,grid.comm().rank());
> +
> + //increase counter and saveStep for next interval
> + saveStep += saveInterval;
> + ++counter;
> + }
> +
> + // print info about time, timestep size and counter
> if (grid.comm().rank()==0) /*@\label{pfc:rank0}@*/
> std::cout << "k=" << k << " t=" << t << " dt=" << dt << std::endl;
> - if (k%20==0) vtkout(grid,c,"pconc",k/20);
> }
> - vtkout(grid,c,"pconc",k/20);
> + vtkout(grid,c,"pconc",counter,tend,grid.comm().rank());
> }
>
> //===============================================================
> @@ -68,9 +92,36 @@ int main (int argc , char ** argv)
>
> // start try/catch block to get error messages from dune
> try {
> - UnitCube<Dune::YaspGrid<2>,64> uc;
> + using namespace Dune;
> +
> + /*
> + UnitCube<YaspGrid<2>,64> uc;
> uc.grid().globalRefine(2);
> partimeloop(uc.grid(),0.5);
> + */
> +
> + // use unitcube from dgf grids
> + std::stringstream dgfFileName;
> + dgfFileName << "grids/unitcube" << GridType :: dimension << ".dgf";
> +
> + // create grid pointer, GridType is defined by gridtype.hh
> + GridPtr<GridType> gridPtr( dgfFileName.str() );
> +
> + // grid reference
> + GridType& grid = *gridPtr;
> +
> + // half grid width 4 times
> + int level = 6 * DGFGridInfo<GridType>::refineStepsForHalf();
> +
> + // refine grid until upper limit of level
> + grid.globalRefine(level);
> +
> + // re-partition grid for better load balancing
> + grid.loadBalance(); /*@\label{pfv:lb}@*/
> +
> + // do time loop until end time 0.5
> + partimeloop(grid, 0.5);
> +
> }
> catch (std::exception & e) {
> std::cout << "STL ERROR: " << e.what() << std::endl;
> diff --git a/traversal.cc b/traversal.cc
> index de444bb..f061529 100644
> --- a/traversal.cc
> +++ b/traversal.cc
> @@ -24,20 +24,25 @@ void traversal (G& grid)
> // Leaf Traversal
> std::cout << "*** Traverse codim 0 leaves" << std::endl;
>
> - // the grid has an iterator providing the access to
> - // all elements (better codim 0 entities) which are leafs
> - // of the refinement tree.
> - // Note the use of the typename keyword and the traits class
> - typedef typename G::template Codim<0>::LeafIterator ElementLeafIterator; /*@\label{tc:ittype}@*/
> + // type of the GridView used for traversal
> + // every grid exports a LeafGridView and a LevelGridView
> + typedef typename G :: LeafGridView LeafGridView; /*@\label{tc:lfgv}@*/
> +
> + // get the instance of the LeafGridView
> + LeafGridView leafView = grid.leafView(); /*@\label{tc:lfv}@*/
> +
> + // Get the iterator type
> + // Note the use of the typename and template keywords
> + typedef typename LeafGridView::template Codim<0>::Iterator ElementLeafIterator; /*@\label{tc:ittype}@*/
>
> // iterate through all entities of codim 0 at the leafs
> int count = 0;
> - for (ElementLeafIterator it = grid.template leafbegin<0>(); /*@\label{tc:forel}@*/
> - it!=grid.template leafend<0>(); ++it)
> + for (ElementLeafIterator it = leafView.template begin<0>(); /*@\label{tc:forel}@*/
> + it!=leafView.template end<0>(); ++it)
> { /*@\label{tc:forel0}@*/
> Dune::GeometryType gt = it->type(); /*@\label{tc:reftype}@*/
> std::cout << "visiting leaf " << gt /*@\label{tc:print}@*/
> - << " with first vertex at " << it->geometry()[0]
> + << " with first vertex at " << it->geometry().corner(0)
> << std::endl;
> count++; /*@\label{tc:count}@*/
> } /*@\label{tc:forel1}@*/
> @@ -50,16 +55,17 @@ void traversal (G& grid)
>
> // Get the iterator type
> // Note the use of the typename and template keywords
> - typedef typename G::template Codim<dim>::LeafIterator VertexLeafIterator; /*@\label{tc:vertit}@*/
> + typedef typename LeafGridView :: template Codim<dim>
> + :: Iterator VertexLeafIterator; /*@\label{tc:vertit}@*/
>
> // iterate through all entities of codim 0 on the given level
> count = 0;
> - for (VertexLeafIterator it = grid.template leafbegin<dim>(); /*@\label{tc:forve}@*/
> - it!=grid.template leafend<dim>(); ++it)
> + for (VertexLeafIterator it = leafView.template begin<dim>(); /*@\label{tc:forve}@*/
> + it!=leafView.template end<dim>(); ++it)
> {
> Dune::GeometryType gt = it->type();
> std::cout << "visiting " << gt
> - << " at " << it->geometry()[0]
> + << " at " << it->geometry().corner(0)
> << std::endl;
> count++;
> }
> @@ -70,20 +76,28 @@ void traversal (G& grid)
> std::cout << std::endl;
> std::cout << "*** Traverse codim 0 level-wise" << std::endl;
>
> + // type of the GridView used for traversal
> + // every grid exports a LeafGridView and a LevelGridView
> + typedef typename G :: LevelGridView LevelGridView; /*@\label{tc:level0}@*/
> +
> // Get the iterator type
> // Note the use of the typename and template keywords
> - typedef typename G::template Codim<0>::LevelIterator ElementLevelIterator; /*@\label{tc:level0}@*/
> + typedef typename LevelGridView :: template Codim<0>
> + :: Iterator ElementLevelIterator;
>
> // iterate through all entities of codim 0 on the given level
> for (int level=0; level<=grid.maxLevel(); level++)
> {
> + // get the instance of the LeafGridView
> + LevelGridView levelView = grid.levelView(level);
> +
> count = 0;
> - for (ElementLevelIterator it = grid.template lbegin<0>(level);
> - it!=grid.template lend<0>(level); ++it)
> + for (ElementLevelIterator it = levelView.template begin<0>();
> + it!=levelView.template end<0>(); ++it)
> {
> Dune::GeometryType gt = it->type();
> std::cout << "visiting " << gt
> - << " with first vertex at " << it->geometry()[0]
> + << " with first vertex at " << it->geometry().corner(0)
> << std::endl;
> count++;
> }
> diff --git a/unitcube.hh b/unitcube.hh
> index 6683610..c1ef749 100644
> --- a/unitcube.hh
> +++ b/unitcube.hh
> @@ -4,116 +4,8 @@
> #include <dune/common/exceptions.hh>
> #include <dune/common/fvector.hh>
>
> -#include <dune/grid/common/gridfactory.hh>
> -
> -// UGGrid 3d, variant 2 (tetrahedra) specialization
> -template< int dim >
> -class BasicUnitCube;
> -
> -
> -template<>
> -class BasicUnitCube< 2 >
> -{
> -protected:
> - template< class Grid >
> - static void insertVertices ( Dune::GridFactory< Grid > &factory )
> - {
> - Dune::FieldVector<double,2> pos;
> -
> - pos[0] = 0; pos[1] = 0;
> - factory.insertVertex(pos);
> -
> - pos[0] = 1; pos[1] = 0;
> - factory.insertVertex(pos);
> -
> - pos[0] = 0; pos[1] = 1;
> - factory.insertVertex(pos);
> -
> - pos[0] = 1; pos[1] = 1;
> - factory.insertVertex(pos);
> - }
> -
> - template< class Grid >
> - static void insertSimplices ( Dune::GridFactory< Grid > &factory )
> - {
> - const Dune::GeometryType type( Dune::GeometryType::simplex, 2 );
> - std::vector< unsigned int > cornerIDs( 3 );
> -
> - cornerIDs[0] = 0; cornerIDs[1] = 1; cornerIDs[2] = 2;
> - factory.insertElement( type, cornerIDs );
> -
> - cornerIDs[0] = 2; cornerIDs[1] = 1; cornerIDs[2] = 3;
> - factory.insertElement( type, cornerIDs );
> - }
> -
> - template< class Grid >
> - static void insertCubes ( Dune::GridFactory< Grid > &factory )
> - {
> - const Dune::GeometryType type( Dune::GeometryType::cube, 2 );
> - std::vector< unsigned int > cornerIDs( 4 );
> - for( int i = 0; i < 4; ++i )
> - cornerIDs[ i ] = i;
> - factory.insertElement( type, cornerIDs );
> - }
> -};
> -
> -
> -template<>
> -class BasicUnitCube< 3 >
> -{
> -protected:
> - template< class Grid >
> - static void insertVertices ( Dune::GridFactory< Grid > &factory )
> - {
> - Dune::FieldVector< double, 3 > pos;
> -
> - pos[0] = 0; pos[1] = 0; pos[2] = 0; factory.insertVertex(pos);
> - pos[0] = 1; pos[1] = 0; pos[2] = 0; factory.insertVertex(pos);
> - pos[0] = 0; pos[1] = 1; pos[2] = 0; factory.insertVertex(pos);
> - pos[0] = 1; pos[1] = 1; pos[2] = 0; factory.insertVertex(pos);
> - pos[0] = 0; pos[1] = 0; pos[2] = 1; factory.insertVertex(pos);
> - pos[0] = 1; pos[1] = 0; pos[2] = 1; factory.insertVertex(pos);
> - pos[0] = 0; pos[1] = 1; pos[2] = 1; factory.insertVertex(pos);
> - pos[0] = 1; pos[1] = 1; pos[2] = 1; factory.insertVertex(pos);
> - }
> -
> - template< class Grid >
> - static void insertSimplices ( Dune::GridFactory< Grid > &factory )
> - {
> - const Dune::GeometryType type( Dune::GeometryType::simplex, 3 );
> - std::vector< unsigned int > cornerIDs( 4 );
> -
> - cornerIDs[0] = 0; cornerIDs[1] = 1; cornerIDs[2] = 2; cornerIDs[3] = 4;
> - factory.insertElement( type, cornerIDs );
> -
> - cornerIDs[0] = 1; cornerIDs[1] = 3; cornerIDs[2] = 2; cornerIDs[3] = 7;
> - factory.insertElement( type, cornerIDs );
> -
> - cornerIDs[0] = 1; cornerIDs[1] = 7; cornerIDs[2] = 2; cornerIDs[3] = 4;
> - factory.insertElement( type, cornerIDs );
> -
> - cornerIDs[0] = 1; cornerIDs[1] = 7; cornerIDs[2] = 4; cornerIDs[3] = 5;
> - factory.insertElement( type, cornerIDs );
> -
> - cornerIDs[0] = 4; cornerIDs[1] = 7; cornerIDs[2] = 2; cornerIDs[3] = 6;
> - factory.insertElement( type, cornerIDs );
> - }
> -
> - template< class Grid >
> - static void insertCubes ( Dune::GridFactory< Grid > &factory )
> - {
> - const Dune::GeometryType type( Dune::GeometryType::cube, 3 );
> - std::vector< unsigned int > cornerIDs( 8 );
> - for( int i = 0; i < 8; ++i )
> - cornerIDs[ i ] = i;
> - factory.insertElement( type, cornerIDs );
> - }
> -};
> -
> -
> -
> // default implementation for any template parameter
> -template<typename T, int variant>
> +template<typename T, int variant> /*@\label{uc:uc0}@*/
> class UnitCube
> {
> public:
> @@ -133,7 +25,10 @@ public:
> private:
> // the constructed grid object
> T grid_;
> -};
> +}; /*@\label{uc:uc1}@*/
> +
> +// include basic unitcube using GridFactory concept
> +#include "basicunitcube.hh"
>
> // include specializations
> #include"unitcube_onedgrid.hh"
> diff --git a/unitcube_alugrid.hh b/unitcube_alugrid.hh
> index 86f64b6..9e2efb1 100644
> --- a/unitcube_alugrid.hh
> +++ b/unitcube_alugrid.hh
> @@ -3,25 +3,38 @@
>
> #if HAVE_ALUGRID
> #include <dune/grid/alugrid.hh>
> +#include <dune/grid/alugrid/3d/alu3dgridfactory.hh>
>
> -// ALU3dGrid tetrahedra specialization. Note: element type determined by type
> +// ALU3dGrid and ALU2dGrid simplex specialization.
> +// Note: element type determined by type
> template<>
> class UnitCube<Dune::ALUSimplexGrid<3,3>,1>
> +: public BasicUnitCube< 3 >
> {
> public:
> typedef Dune::ALUSimplexGrid<3,3> GridType;
>
> - UnitCube () : filename("grids/cube.tetra"), grid_(filename.c_str())
> - {}
> +private:
> + GridType * grid_;
>
> - GridType& grid ()
> +public:
> + UnitCube ()
> {
> - return grid_;
> + Dune::GridFactory< GridType > factory;
> + BasicUnitCube< 3 >::insertVertices( factory );
> + BasicUnitCube< 3 >::insertSimplices( factory );
> + grid_ = factory.createGrid( );
> }
>
> -private:
> - std::string filename;
> - GridType grid_;
> + ~UnitCube()
> + {
> + delete grid_;
> + }
> +
> + GridType &grid ()
> + {
> + return *grid_;
> + }
> };
>
> // ALU2SimplexGrid 2d specialization. Note: element type determined by type
> @@ -47,21 +60,32 @@ private:
> // ALU3dGrid hexahedra specialization. Note: element type determined by type
> template<>
> class UnitCube<Dune::ALUCubeGrid<3,3>,1>
> +: public BasicUnitCube< 3 >
> {
> public:
> typedef Dune::ALUCubeGrid<3,3> GridType;
>
> - UnitCube () : filename("grids/cube.hexa"), grid_(filename.c_str())
> - {}
> +private:
> + GridType * grid_;
>
> - GridType& grid ()
> +public:
> + UnitCube ()
> {
> - return grid_;
> + Dune::GridFactory< GridType > factory;
> + BasicUnitCube< 3 >::insertVertices( factory );
> + BasicUnitCube< 3 >::insertCubes( factory );
> + grid_ = factory.createGrid( );
> }
>
> -private:
> - std::string filename;
> - GridType grid_;
> + ~UnitCube()
> + {
> + delete grid_;
> + }
> +
> + GridType &grid ()
> + {
> + return *grid_;
> + }
> };
> #endif
>
> diff --git a/unitcube_uggrid.hh b/unitcube_uggrid.hh
> index 79c6757..bc7c02b 100644
> --- a/unitcube_uggrid.hh
> +++ b/unitcube_uggrid.hh
> @@ -24,7 +24,8 @@ public:
> else if( variant == 2 )
> BasicUnitCube< dim >::insertSimplices( factory );
> else
> - DUNE_THROW( Dune::NotImplemented, "Variant " << variant << " of UG unit cube not implemented." );
> + DUNE_THROW( Dune::NotImplemented, "Variant "
> + << variant << " of UG unit cube not implemented." );
> grid_ = factory.createGrid();
> }
>
> diff --git a/vertexdata.hh b/vertexdata.hh
> index 8fecd47..c4bb8bb 100644
> --- a/vertexdata.hh
> +++ b/vertexdata.hh
> @@ -23,8 +23,12 @@ void vertexdata (const G& grid, const F& f)
> // get dimension and coordinate type from Grid
> const int dim = G::dimension;
> typedef typename G::ctype ct;
> + typedef typename G::LeafGridView GridView;
> // dertermine type of LeafIterator for codimension = dimension
> - typedef typename G::template Codim<dim>::LeafIterator VertexLeafIterator;
> + typedef typename GridView::template Codim<dim>::Iterator VertexLeafIterator;
> +
> + // get grid view on the leaf part
> + GridView gridView = grid.leafView();
>
> // make a mapper for codim 0 entities in the leaf grid
> Dune::LeafMultipleCodimMultipleGeomTypeMapper<G,P1Layout>
> @@ -34,8 +38,8 @@ void vertexdata (const G& grid, const F& f)
> std::vector<double> c(mapper.size());
>
> // iterate through all entities of codim 0 at the leafs
> - for (VertexLeafIterator it = grid.template leafbegin<dim>();
> - it!=grid.template leafend<dim>(); ++it)
> + for (VertexLeafIterator it = gridView.template begin<dim>();
> + it!=gridView.template end<dim>(); ++it)
> {
> // evaluate functor and store value
> c[mapper.map(*it)] = f(it->geometry().corner(0));
> @@ -57,7 +61,7 @@ void vertexdata (const G& grid, const F& f)
> // display data
> grape.displayVector("concentration", // name of data that appears in grape
> c, // data vector
> - grid.leafIndexSet(), // used index set
> + gridView.indexSet(), // used index set
> polynomialOrder, // polynomial order of data
> dimRange); // dimRange of data
> }
> diff --git a/visualization.cc b/visualization.cc
> index 3ccf52b..c3af5d6 100644
> --- a/visualization.cc
> +++ b/visualization.cc
> @@ -43,28 +43,33 @@ int main(int argc, char **argv)
> {
> /*
> UnitCube<Dune::OneDGrid,1> uc0;
> - UnitCube<Dune::YaspGrid<3>,1> uc1;
> - UnitCube<Dune::YaspGrid<2>,1> uc2;
> - */
> + UnitCube<Dune::YaspGrid<dimGrid>,1> uc1;
> +
> #if HAVE_UG
> - UnitCube< Dune::UGGrid< dimGrid >, 2 > uc6;
> - dowork( uc6.grid(), 3 );
> + UnitCube< Dune::UGGrid< dimGrid >, 2 > uc2;
> + dowork( uc2.grid(), 3 );
> #endif
>
> #ifdef GRIDDIM
> #if HAVE_ALBERTA
> - UnitCube< Dune::AlbertaGrid< dimGrid, dimGrid >, 1 > uc7;
> + UnitCube< Dune::AlbertaGrid< dimGrid, dimGrid >, 1 > uc3;
> // note: The 3d cube cannot be bisected recursively
> - dowork( uc7.grid(), (dimGrid < 3 ? 6 : 0 ) );
> + dowork( uc3.grid(), (dimGrid < 3 ? 6 : 0 ) );
> #endif
> #endif
> + */
>
> UnitCube< Dune::SGrid< dimGrid, dimGrid >, 1 > uc4;
> dowork( uc4.grid(), 3 );
>
> #if HAVE_ALUGRID
> - UnitCube< Dune::ALUSimplexGrid< dimGrid, dimGrid > , 1 > uc8;
> - dowork( uc8.grid(), 3 );
> + UnitCube< Dune::ALUSimplexGrid< dimGrid, dimGrid > , 1 > uc5;
> + dowork( uc5.grid(), 3 );
> +
> +#if GRIDDIM == 3
> + UnitCube< Dune::ALUCubeGrid< dimGrid, dimGrid > , 1 > uc6;
> + dowork( uc6.grid(), 3 );
> +#endif
> #endif
> }
> catch (std::exception & e) {
> diff --git a/vtkout.hh b/vtkout.hh
> index 562cb71..a810ec9 100644
> --- a/vtkout.hh
> +++ b/vtkout.hh
> @@ -2,11 +2,20 @@
> #include <stdio.h>
>
> template<class G, class V>
> -void vtkout (const G& grid, const V& c, const char* name, int k)
> +void vtkout (const G& grid, const V& c, const char* name, int k, double time=0.0, int rank=0)
> {
> Dune::VTKWriter<typename G::LeafGridView> vtkwriter(grid.leafView());
> char fname[128];
> + char sername[128];
> sprintf(fname,"%s-%05d",name,k);
> + sprintf(sername,"%s.series",name);
> vtkwriter.addCellData(c,"celldata");
> vtkwriter.write(fname,Dune::VTKOptions::ascii);
> +
> + if ( rank == 0)
> + {
> + std::ofstream serstream(sername, (k==0 ? std::ios_base::out : std::ios_base::app));
> + serstream << k << " " << fname << ".vtu " << time << std::endl;
> + serstream.close();
> + }
> }
> diff --git a/writePVD b/writePVD
> new file mode 100644
> index 0000000..14525da
> --- /dev/null
> +++ b/writePVD
> @@ -0,0 +1,64 @@
> +#! /bin/bash
> +
> +# read parameters
> +if [ -z "$1" ]; then
> + echo "Usage: writePVD [file_prefix]"
> + exit
> +fi
> +
> +param=$1
> +SERIESFILE="${param##*/}"
> +DIR="${param%$SERIESFILE}"
> +DIR="${DIR:-.}"
> +FILEPREFIX="${SERIESFILE%.*}"
> +SERIESFILE=$1
> +
> +echo $SERIESFILE $DIR $FILEPREFIX
> +
> +THISDIR=$PWD
> +# cd $DIR
> +
> +# output file name
> +OUTPUTFILE=$FILEPREFIX.pvd
> +
> +#write header
> +echo "Writing file: $OUTPUTFILE"
> +echo "<?xml version=\"1.0\"?>" > $OUTPUTFILE
> +echo "<VTKFile type=\"Collection\" version=\"0.1\" byte_order=\"LittleEndian\">" >> $OUTPUTFILE
> +echo " <Collection>" >> $OUTPUTFILE
> +echo "" >> $OUTPUTFILE
> +
> +# Set loop separator to end of line
> +BAKIFS=$IFS
> +IFS=$(echo -en "\n\b")
> +# the following applied to the series file
> +exec 3<&0
> +exec 0<$SERIESFILE
> +while read line
> +do
> + # file name and time stamp
> + NAME=$(echo $line | cut -d " " -f 2)
> + TIME=$(echo $line | cut -d " " -f 3)
> + # strip unnecessary path from file name
> + while ! test -e $NAME
> + do
> + NAME="${NAME#*/}"
> + if [ $NAME == "" ];
> + then # file does not exists
> + NAME=$(echo $line | cut -d " " -f 2)
> + echo $NAME " not found!"
> + exit
> + fi
> + done
> + # write line into output file
> + echo "<DataSet timestep=\"$TIME\" group=\"\" part=\"0\" file=\"$NAME\"/>" >> $OUTPUTFILE
> +done
> +exec 0<&3
> +# restore $IFS which was used to determine what the field separators are
> +IFS=$BAKIFS
> +
> +# write the end of file
> +echo "" >> $OUTPUTFILE
> +echo " </Collection>" >> $OUTPUTFILE
> +echo "</VTKFile>" >> $OUTPUTFILE
> +
More information about the Dune
mailing list