Commit a4d822ae authored by Peter Bastian's avatar Peter Bastian

works with UG now :-)

still missing: Alberta intersection iterator conforming to interface

[[Imported from SVN: r13]]
parent 32bb34f3
......@@ -11,7 +11,7 @@
#include "initialize.hh"
#include "evolve.hh"
#include "finitevolumeadapt.hh"
#include "vtkconcentration.hh"
#include "vtkout.hh"
//===============================================================
// the time loop function working for all types of grids
......@@ -45,9 +45,10 @@ void timeloop (G& grid, double tend, int lmin, int lmax)
initialize(grid,mapper,c);
for (int i=grid.maxLevel(); i<lmax; i++)
{
finitevolumeadapt(grid,mapper,c,lmin,lmax);
finitevolumeadapt(grid,mapper,c,lmin,lmax,0);
initialize(grid,mapper,c);
}
vtkout(grid,c,"concentration",0);
double dt, t=0;
int k=0;
......@@ -56,13 +57,12 @@ void timeloop (G& grid, double tend, int lmin, int lmax)
k++;
evolve(grid,mapper,c,t,dt);
t += dt;
if (k%20==0) vtkout(grid,c,"concentration",k/20);
std::cout << "s=" << grid.size(0) << " k=" << k
<< " t=" << t << " dt=" << dt << std::endl;
finitevolumeadapt(grid,mapper,c,lmin,lmax);
finitevolumeadapt(grid,mapper,c,lmin,lmax,k);
}
// output results
vtkconcentration(grid,c,0);
vtkout(grid,c,"concentration",k/20);
}
//===============================================================
......@@ -77,15 +77,15 @@ int main (int argc , char ** argv)
UnitCube<Dune::YaspGrid<2,2>,1> uc;
#if HAVE_UG
UnitCube<Dune::UGGrid<2,2>,2> uc2;
UnitCube<Dune::UGGrid<2,2>,1> uc2;
#endif
#if HAVE_ALBERTA
UnitCube<Dune::AlbertaGrid<2,2>,1> uc3;
#endif
uc3.grid().globalRefine(8);
timeloop(uc3.grid(),0.5,8,17);
// uc2.grid().globalRefine(5);
// timeloop(uc2.grid(),0.5,5,9);
timeloop(uc3.grid(),0.5,8,18);
// uc2.grid().globalRefine(4);
// timeloop(uc2.grid(),0.5,4,9);
#if HAVE_MPI
MPI_Finalize();
......
......@@ -65,9 +65,12 @@ Oliver Sander$^\ddagger$}
\date{\today}
\publishers{%
\bigskip
\includegraphics[width=0.43\textwidth]{EPS/conc.eps}\\
\bigskip
\vspace{10mm}
\includegraphics[width=0.32\textwidth]{EPS/alberta2d-view2.eps}\hfill
\includegraphics[width=0.32\textwidth]{EPS/ug2dtri-view2.eps}\hfill
\includegraphics[width=0.32\textwidth]{EPS/ug2dquad-view2.eps}
\vspace{10mm}
{\normalsize $^\ast$Interdisziplinäres Zentrum für Wissenschaftliches Rechnen,
Universität Heidelberg,\\
Im Neuenheimer Feld 368, D-69120 Heidelberg, Germany}\\
......@@ -1052,10 +1055,10 @@ numberstyle=\tiny, numbersep=5pt]{../initialize.hh}
numberstyle=\tiny, numbersep=5pt]{../evolve.hh}
\end{lst}
\begin{lst}[File dune-grid-howto/vtkconcentration.hh] \mbox{}
\begin{lst}[File dune-grid-howto/vtkout.hh] \mbox{}
\nopagebreak
\lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
numberstyle=\tiny, numbersep=5pt]{../vtkconcentration.hh}
numberstyle=\tiny, numbersep=5pt]{../vtkout.hh}
\end{lst}
\begin{lst}[File dune-grid-howto/finitevolume.cc] \mbox{}
......
......@@ -79,57 +79,64 @@ void evolve (const G& grid, const M& mapper, V& c, double t, double& dt)
if (factor>=0) sumfactor += factor;
// handle interior face
// if (is.neighbor())
// {
// // access neighbor
// EntityPointer outside = is.outside();
// // handle face from correct side
// if (outside->isLeaf())
// {
// int indexj = mapper.map(*outside);
// // std::cout << "leaf neighbor i: " << indexi << "." << it->level()
// // << " j=" << indexj << "." << outside->level()
// // << " => " << (it->level()-outside->level()) << std::endl;
// if ( it->level()>outside->level() ||
// (it->level()==outside->level() && indexi<indexj) )
// {
// // compute factor in neighbor
// Dune::GeometryType nbgt = outside->geometry().type();
// const Dune::FieldVector<ct,dim>&
// nblocal = Dune::ReferenceElements<ct,dim>::general(nbgt).position(0,0);
// double nbvolume = outside->geometry().integrationElement(nblocal)
// *Dune::ReferenceElements<ct,dim>::general(nbgt).volume();
// double nbfactor = velocity*integrationOuterNormal/nbvolume;
// if (factor<0) // inflow
// {
// update[indexi] -= c[indexj]*factor;
// update[indexj] += c[indexj]*nbfactor;
// }
// else // outflow
// {
// update[indexi] -= c[indexi]*factor;
// update[indexj] += c[indexi]*nbfactor;
// }
// }
// }
// }
if (is.neighbor()) // Alberta
if (0)
{
// access neighbor
EntityPointer outside = is.outside();
if (is.neighbor()) // "correct" version
{
// access neighbor
EntityPointer outside = is.outside();
// handle face from correct side
if (outside->isLeaf())
{
int indexj = mapper.map(*outside);
// std::cout << "leaf neighbor i: " << indexi << "." << it->level()
// << " j=" << indexj << "." << outside->level()
// << " => " << (it->level()-outside->level()) << std::endl;
if ( it->level()>outside->level() ||
(it->level()==outside->level() && indexi<indexj) )
{
// compute factor in neighbor
Dune::GeometryType nbgt = outside->geometry().type();
const Dune::FieldVector<ct,dim>&
nblocal = Dune::ReferenceElements<ct,dim>::general(nbgt).position(0,0);
double nbvolume = outside->geometry().integrationElement(nblocal)
*Dune::ReferenceElements<ct,dim>::general(nbgt).volume();
double nbfactor = velocity*integrationOuterNormal/nbvolume;
if (factor<0) // inflow
{
update[indexi] -= c[indexj]*factor;
update[indexj] += c[indexj]*nbfactor;
}
else // outflow
{
update[indexi] -= c[indexi]*factor;
update[indexj] += c[indexi]*nbfactor;
}
}
}
}
}
// neighbors index
int indexj = mapper.map(*outside);
if (factor<0) // inflow
update[indexi] -= c[indexj]*factor;
else // outflow
update[indexi] -= c[indexi]*factor;
if (1)
{
if (is.neighbor()) // Alberta
{
// access neighbor
EntityPointer outside = is.outside();
// neighbors index
int indexj = mapper.map(*outside);
if (factor<0) // inflow
update[indexi] -= c[indexj]*factor;
else // outflow
update[indexi] -= c[indexi]*factor;
}
}
......
......@@ -10,7 +10,7 @@
#include "transportproblem.hh"
#include "initialize.hh"
#include "evolve.hh"
#include "vtkconcentration.hh"
#include "vtkout.hh"
//===============================================================
// the time loop function working for all types of grids
......@@ -55,7 +55,7 @@ void timeloop (const G& grid, double tend)
}
// output results
vtkconcentration(grid,c,k);
vtkout(grid,c,"concentration",k);
}
//===============================================================
......
......@@ -14,11 +14,11 @@ struct RestrictedValue
};
template<class G, class M, class V>
bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax)
bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax, int k)
{
// tol value for refinement strategy
const double refinetol = 0.02;
const double coarsentol = 0.0002;
const double refinetol = 0.05;
const double coarsentol = 0.001;
// type used for coordinates in the grid
typedef typename G::ctype ct;
......@@ -38,32 +38,68 @@ bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax)
typedef typename IdSet::IdType IdType;
// compute cell indicators
V indicator(c.size());
double globaldelta = -1E100;
V indicator(c.size(),-1E100);
double globalmax = -1E100;
double globalmin = 1E100;
for (LeafIterator it = grid.template leafbegin<0>();
it!=grid.template leafend<0>(); ++it)
{
// my index
int indexi = mapper.map(*it);
double localdelta = -1E100;
// global min/max
globalmax = std::max(globalmax,c[indexi]);
globalmin = std::min(globalmin,c[indexi]);
IntersectionIterator isend = it->iend();
for (IntersectionIterator is = it->ibegin(); is!=isend; ++is)
if (is.neighbor())
{
int indexj = mapper.map(*(is.outside()));
localdelta = std::max(localdelta,std::abs(c[indexj]-c[indexi]));
// access neighbor
EntityPointer outside = is.outside();
if (0) // correct
{
// handle face from correct side
if (outside->isLeaf())
{
int indexj = mapper.map(*outside);
if ( it.level()>outside->level() ||
(it.level()==outside->level() && indexi<indexj) )
{
double localdelta = std::abs(c[indexj]-c[indexi]);
indicator[indexi] = std::max(indicator[indexi],localdelta);
indicator[indexj] = std::max(indicator[indexj],localdelta);
}
}
}
if (1) // alberta
{
int indexj = mapper.map(*(is.outside()));
double localdelta = std::abs(c[indexj]-c[indexi]);
indicator[indexi] = std::max(indicator[indexi],localdelta);
}
}
indicator[indexi] = localdelta;
globaldelta = std::max(globaldelta,localdelta);
}
// mark cells for refinement/coarsening
double globaldelta = globalmax-globalmin;
V mark(c.size(),0.0);
for (LeafIterator it = grid.template leafbegin<0>();
it!=grid.template leafend<0>(); ++it)
{
if (indicator[mapper.map(*it)]>refinetol*globaldelta && it.level()<lmax)
if (indicator[mapper.map(*it)]>refinetol*globaldelta && (it.level()<lmax || it->isRegular()==false))
{
grid.mark(1,it);
mark[mapper.map(*it)] = 1;
}
if (indicator[mapper.map(*it)]<coarsentol*globaldelta && it.level()>lmin)
{
grid.mark(-1,it);
mark[mapper.map(*it)] = -1;
}
}
// restrict to coarse elements
......
......@@ -2,12 +2,15 @@
// vi: set et ts=4 sw=2 sts=2:
#include "dune/io/file/vtk/vtkwriter.hh"
#include "dune/disc/functions/p0function.hh"
#include <stdio.h>
template<class G, class V>
void vtkconcentration (const G& grid, const V& c, int k)
void vtkout (const G& grid, const V& c, char* name, int k)
{
Dune::VTKWriter<G> vtkwriter(grid);
Dune::LeafP0FunctionWrapper<G,std::vector<double> > cc(grid,c);
vtkwriter.addCellData(cc,"concentration");
vtkwriter.write("concentration",Dune::VTKOptions::binaryappended);
char fname[128];
sprintf(fname,"%s-%05d",name,k);
vtkwriter.addCellData(cc,"celldata");
vtkwriter.write(fname,Dune::VTKOptions::binaryappended);
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment