Commit be2d6649 authored by Martin Nolte's avatar Martin Nolte

adapt to new intersection syntax

[[Imported from SVN: r215]]
parent f263c5bc
......@@ -24,6 +24,7 @@ othergrids
parfinitevolume
traversal
visualization
*.swp
*.vtu
dependencies.m4
dune-grid-howto-*.*.tar.gz
......
......@@ -58,7 +58,7 @@ void evolve (const G& grid, const M& mapper, V& c, double t, double& dt)
for (IntersectionIterator is = it->ileafbegin(); is!=isend; ++is)
{
// get geometry type of face
Dune::GeometryType gtf = is.intersectionSelfLocal().type();
Dune::GeometryType gtf = is->intersectionSelfLocal().type();
// center in face's reference element
const Dune::FieldVector<ct,dim-1>&
......@@ -66,13 +66,13 @@ void evolve (const G& grid, const M& mapper, V& c, double t, double& dt)
// get normal vector scaled with volume
Dune::FieldVector<ct,dimworld> integrationOuterNormal
= is.integrationOuterNormal(facelocal);
= is->integrationOuterNormal(facelocal);
integrationOuterNormal
*= Dune::ReferenceElements<ct,dim-1>::general(gtf).volume();
// center of face in global coordinates
Dune::FieldVector<ct,dimworld>
faceglobal = is.intersectionGlobal().global(facelocal);
faceglobal = is->intersectionGlobal().global(facelocal);
// evaluate velocity at face center
Dune::FieldVector<double,dim> velocity = u(faceglobal,t);
......@@ -84,10 +84,10 @@ 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 /*@\label{evh:neighbor}@*/
if (is->neighbor()) // "correct" version /*@\label{evh:neighbor}@*/
{
// access neighbor
EntityPointer outside = is.outside();
EntityPointer outside = is->outside();
int indexj = mapper.map(*outside);
// compute flux from one side only
......@@ -117,11 +117,13 @@ void evolve (const G& grid, const M& mapper, V& c, double t, double& dt)
}
// handle boundary face
if (is.boundary()) /*@\label{evh:bndry}@*/
if (is->boundary()) /*@\label{evh:bndry}@*/
{
if (factor<0) // inflow, apply boundary condition
update[indexi] -= b(faceglobal,t)*factor;
else // outflow
update[indexi] -= c[indexi]*factor;
}
} // end all intersections /*@\label{evh:flux1}@*/
// compute dt restriction
......
......@@ -27,7 +27,8 @@ bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax, int k)
typedef typename G::template Codim<0>::LeafIterator LeafIterator;
typedef typename G::template Codim<0>::LevelIterator LevelIterator;
// entity pointer
// entity and entity pointer
typedef typename G::template Codim<0>::Entity Entity;
typedef typename G::template Codim<0>::EntityPointer EntityPointer;
// intersection iterator type
......@@ -53,22 +54,26 @@ bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax, int k)
IntersectionIterator isend = it->ileafend();
for (IntersectionIterator is = it->ileafbegin(); is!=isend; ++is)
if (is.neighbor())
{
const typename IntersectionIterator::Intersection &intersection = *is;
if( !intersection.neighbor() )
continue;
// access neighbor
const EntityPointer pOutside = intersection.outside();
const Entity &outside = *pOutside;
int indexj = mapper.map( outside );
// handle face from one side only
if ( it.level() > outside.level() ||
(it.level() == outside.level() && indexi<indexj) )
{
// access neighbor
EntityPointer outside = is.outside();
int indexj = mapper.map(*outside);
// handle face from one side only
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);
}
} /*@\label{fah:loop1}@*/
}
} /*@\label{fah:loop1}@*/
// mark cells for refinement/coarsening
double globaldelta = globalmax-globalmin;
......@@ -79,21 +84,30 @@ bool finitevolumeadapt (G& grid, M& mapper, V& c, int lmin, int lmax, int k)
if (indicator[mapper.map(*it)]>refinetol*globaldelta
&& (it.level()<lmax || !it->isRegular()))
{
grid.mark(1,it);
marked++;
IntersectionIterator isend = it->ileafend();
for (IntersectionIterator is = it->ileafbegin(); is!=isend; ++is)
if (is.neighbor())
if (is.outside().level()<lmax || !is.outside()->isRegular())
grid.mark(1,is.outside());
const Entity &entity = *it;
grid.mark( 1, entity );
++marked;
IntersectionIterator isend = entity.ileafend();
for( IntersectionIterator is = entity.ileafbegin(); is != isend; ++is )
{
const typename IntersectionIterator::Intersection intersection = *is;
if( !intersection.neighbor() )
continue;
const EntityPointer pOutside = intersection.outside();
const Entity &outside = *pOutside;
if( (outside.level() < lmax) || !outside.isRegular() )
grid.mark( 1, outside );
}
}
if (indicator[mapper.map(*it)]<coarsentol*globaldelta && it.level()>lmin)
{
grid.mark(-1,it);
marked++;
grid.mark( -1, *it );
++marked;
}
} /*@\label{fah:loop3}@*/
if (marked==0) return false;
} /*@\label{fah:loop3}@*/
if( marked==0 )
return false;
// restrict to coarse elements
std::map<IdType,RestrictedValue> restrictionmap; // restricted concentration /*@\label{fah:loop4}@*/
......
......@@ -22,6 +22,9 @@ void parevolve (const G& grid, const M& mapper, V& c, double t, double& dt)
// intersection iterator type
typedef typename G::template Codim<0>::LeafIntersectionIterator IntersectionIterator;
// type of intersection
typedef typename IntersectionIterator::Intersection Intersection;
// entity pointer type
typedef typename G::template Codim<0>::EntityPointer EntityPointer;
......@@ -59,25 +62,28 @@ 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
IntersectionIterator isend = it->ileafend();
for (IntersectionIterator is = it->ileafbegin(); is!=isend; ++is)
const IntersectionIterator isend = it->ileafend();
for( IntersectionIterator is = it->ileafbegin(); is != isend; ++is )
{
const Intersection &intersection = *is;
// get geometry type of face
Dune::GeometryType gtf = is.intersectionSelfLocal().type();
Dune::GeometryType gtf = intersection.intersectionSelfLocal().type();
const Dune::ReferenceElement< ct, dim-1 > &refElement
= Dune::ReferenceElements< ct, dim-1 >::general( gtf );
// center in face's reference element
const Dune::FieldVector<ct,dim-1>&
facelocal = Dune::ReferenceElements<ct,dim-1>::general(gtf).position(0,0);
const Dune::FieldVector< ct, dim-1 > &facelocal = refElement.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
= intersection.integrationOuterNormal( facelocal );
integrationOuterNormal *= refElement.volume();
// center of face in global coordinates
Dune::FieldVector<ct,dimworld>
faceglobal = is.intersectionGlobal().global(facelocal);
Dune::FieldVector< ct, dimworld > faceglobal
= intersection.intersectionGlobal().global( facelocal );
// evaluate velocity at face center
Dune::FieldVector<double,dim> velocity = u(faceglobal,t);
......@@ -89,15 +95,18 @@ 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())
if( intersection.neighbor() )
{
// access neighbor
EntityPointer outside = is.outside();
EntityPointer outside = intersection.outside();
int indexj = mapper.map(*outside);
const int insideLevel = it->level();
const int outsideLevel = outside->level();
// handle face from one side
if ( it->level()>outside->level() ||
(it->level()==outside->level() && indexi<indexj) )
if( (insideLevel > outsideLevel)
|| ((insideLevel == outsideLevel) && (indexi < indexj)) )
{
// compute factor in neighbor
Dune::GeometryType nbgt = outside->type();
......@@ -107,12 +116,12 @@ void parevolve (const G& grid, const M& mapper, V& c, double t, double& dt)
*Dune::ReferenceElements<ct,dim>::general(nbgt).volume();
double nbfactor = velocity*integrationOuterNormal/nbvolume;
if (factor<0) // inflow
if( factor < 0 ) // inflow
{
update[indexi] -= c[indexj]*factor;
update[indexj] += c[indexj]*nbfactor;
}
else // outflow
else // outflow
{
update[indexi] -= c[indexi]*factor;
update[indexj] += c[indexi]*nbfactor;
......@@ -121,12 +130,14 @@ void parevolve (const G& grid, const M& mapper, V& c, double t, double& dt)
}
// handle boundary face
if (is.boundary())
if (factor<0) // inflow, apply boundary condition
if( intersection.boundary() )
{
if( factor < 0 ) // inflow, apply boundary condition
update[indexi] -= b(faceglobal,t)*factor;
else // outflow
else // outflow
update[indexi] -= c[indexi]*factor;
} // end all intersections
}
} // end all intersections
// compute dt restriction
if (it->partitionType()==Dune::InteriorEntity) /*@\label{peh:inter}@*/
......
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