Commit 7071c3f9 authored by Peter Bastian's avatar Peter Bastian

changes due to neighbor() -> leafNeighbor()

[[Imported from SVN: r38]]
parent cc99ba2e
......@@ -76,19 +76,19 @@ int main (int argc , char ** argv)
UnitCube<Dune::SGrid<1,1>,1> uc1;
UnitCube<Dune::YaspGrid<2,2>,1> uc;
#if HAVE_UG
UnitCube<Dune::UGGrid<3,3>,2> uc2;
UnitCube<Dune::UGGrid<2,2>,2> uc2;
#endif
#if HAVE_ALBERTA
#if DUNE_PROBLEM_DIM==2
UnitCube<Dune::AlbertaGrid<2,2>,1> uc3;
#endif
#endif
// uc3.grid().globalRefine(8);
// timeloop(uc3.grid(),0.5,8,18);
// uc2.grid().globalRefine(3);
// timeloop(uc2.grid(),0.5,3,6);
uc0.grid().globalRefine(4);
timeloop(uc0.grid(),0.5,4,9);
uc3.grid().globalRefine(8);
timeloop(uc3.grid(),0.5,8,18);
// uc2.grid().globalRefine(4);
// timeloop(uc2.grid(),0.5,4,8);
// uc0.grid().globalRefine(4);
// timeloop(uc0.grid(),0.5,4,9);
#if HAVE_MPI
MPI_Finalize();
......
......@@ -2081,13 +2081,15 @@ where the printing of the data of the current time step is restricted
to the process with rank 0.
% \chapter{Input and Output}
\chapter{Input and Output}
% \section{Visualization with Grape}
\section{Visualization with Grape}
% \section{Visualization with VTK}
\section{Visualization with VTK}
% \section{Dune portable grid format}
\section{Dune portable grid format}
\section{Amiramesh output}
% \chapter{Iterative Solver Template Library}
......@@ -2111,7 +2113,36 @@ to the process with rank 0.
% \end{lst}
\chapter{Outlook}
\begin{figure}
\begin{center}
\includegraphics[width=0.32\textwidth]{./EPS/alberta2d.eps}\
\includegraphics[width=0.32\textwidth]{./EPS/ugsimplex2d.eps}\
\includegraphics[width=0.32\textwidth]{./EPS/ugcube2d.eps}
\includegraphics[width=0.32\textwidth]{./EPS/alberta3d.eps}\
\includegraphics[width=0.32\textwidth]{./EPS/alusimplex3d.eps}\
\includegraphics[width=0.32\textwidth]{./EPS/alucube3d.eps}
\includegraphics[width=0.32\textwidth]{./EPS/ugsimplex3d.eps}\
\includegraphics[width=0.32\textwidth]{./EPS/ugcube3d.eps}\
\includegraphics[width=0.32\textwidth]{./EPS/iso.eps}
\end{center}
\caption{Adaptive solution of an elliptic model problem with $P_1$
conforming finite elements and residual based error
estimator. Illustrates that adaptive finite element algorithm can be formulated
independent of dimension, element type and refinement scheme. From
top to bottom, left to right: Alberta (bisection, 2d), UG (red/green
on triangles), UG (red/green on quadrilaterals), Alberta (bisection,
3d), ALU (hanging nodes on tetrahedra), ALU (hanging nodes on
hexahedra), UG (red/green on tetrahedra), UG (red/green on hexahedra,
pyramids and tetrahedra), isosurfaces of solution.}
\label{Fig:OutlookP1}
\end{figure}
......
......@@ -65,8 +65,10 @@ void evolve (const G& grid, const M& mapper, V& c, double t, double& dt)
facelocal = Dune::ReferenceElements<ct,dim-1>::general(gtf).position(0,0);
// get normal vector scaled with volume
Dune::FieldVector<ct,dimworld> integrationOuterNormal = is.integrationOuterNormal(facelocal);
integrationOuterNormal *= Dune::ReferenceElements<ct,dim-1>::general(gtf).volume();
Dune::FieldVector<ct,dimworld> integrationOuterNormal
= is.integrationOuterNormal(facelocal);
integrationOuterNormal
*= Dune::ReferenceElements<ct,dim-1>::general(gtf).volume();
// center of face in global coordinates
Dune::FieldVector<ct,dimworld>
......@@ -82,36 +84,34 @@ 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()) // "correct" version
if (is.leafNeighbor()) // "correct" version
{
// access neighbor
EntityPointer outside = is.outside();
int indexj = mapper.map(*outside);
// handle face from correct side
if (outside->isLeaf())
// compute flux from one side only
// this should become easier with the new IntersectionIterator functionality!
if ( it->level()>outside->level() ||
(it->level()==outside->level() && indexi<indexj) )
{
int indexj = mapper.map(*outside);
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
{
// 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;
}
update[indexi] -= c[indexj]*factor;
update[indexj] += c[indexj]*nbfactor;
}
else // outflow
{
update[indexi] -= c[indexi]*factor;
update[indexj] += c[indexi]*nbfactor;
}
}
}
......
......@@ -53,22 +53,19 @@ bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax, int k)
IntersectionIterator isend = it->iend();
for (IntersectionIterator is = it->ibegin(); is!=isend; ++is)
if (is.neighbor())
if (is.leafNeighbor())
{
// access neighbor
EntityPointer outside = is.outside();
int indexj = mapper.map(*outside);
// handle face from correct side
if (outside->isLeaf())
// handle face from one side only
if ( it.level()>outside->level() ||
(it.level()==outside->level() && indexi<indexj) )
{
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);
}
double localdelta = std::abs(c[indexj]-c[indexi]);
indicator[indexi] = std::max(indicator[indexi],localdelta);
indicator[indexj] = std::max(indicator[indexj],localdelta);
}
}
}
......@@ -84,9 +81,8 @@ bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax, int k)
grid.mark(1,it);
IntersectionIterator isend = it->iend();
for (IntersectionIterator is = it->ibegin(); is!=isend; ++is)
if (is.neighbor())
if (is.outside()->isLeaf() &&
(is.outside().level()<lmax || !is.outside()->isRegular()))
if (is.leafNeighbor())
if (is.outside().level()<lmax || !is.outside()->isRegular())
grid.mark(1,is.outside());
}
if (indicator[mapper.map(*it)]<coarsentol*globaldelta && it.level()>lmin)
......
......@@ -70,8 +70,10 @@ void parevolve (const G& grid, const M& mapper, V& c, double t, double& dt)
facelocal = Dune::ReferenceElements<ct,dim-1>::general(gtf).position(0,0);
// get normal vector scaled with volume
Dune::FieldVector<ct,dimworld> integrationOuterNormal = is.integrationOuterNormal(facelocal);
integrationOuterNormal *= Dune::ReferenceElements<ct,dim-1>::general(gtf).volume();
Dune::FieldVector<ct,dimworld> integrationOuterNormal
= is.integrationOuterNormal(facelocal);
integrationOuterNormal
*= Dune::ReferenceElements<ct,dim-1>::general(gtf).volume();
// center of face in global coordinates
Dune::FieldVector<ct,dimworld>
......@@ -87,36 +89,33 @@ void parevolve (const G& grid, const M& mapper, V& c, double t, double& dt)
if (factor>=0) sumfactor += factor;
// handle interior face
if (is.neighbor()) // "correct" version
if (is.leafNeighbor())
{
// access neighbor
EntityPointer outside = is.outside();
int indexj = mapper.map(*outside);
// handle face from correct side
if (outside->isLeaf())
// handle face from one side
if ( it->level()>outside->level() ||
(it->level()==outside->level() && indexi<indexj) )
{
int indexj = mapper.map(*outside);
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
{
// 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;
}
update[indexi] -= c[indexj]*factor;
update[indexj] += c[indexj]*nbfactor;
}
else // outflow
{
update[indexi] -= c[indexi]*factor;
update[indexj] += c[indexi]*nbfactor;
}
}
}
......
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