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