Skip to content

GitLab

Commit a379dc18 authored by Timo Koch's avatar Timo Koch
Browse files

[decoupled] Remove all occurences of ElementPointer 

Interactionvolumes now store EntitySeeds and return Entities.
parent a74b35d5
Expand all Hide whitespace changes
Inline Side-by-side
dumux/decoupled/1p/diffusion/fv/fvpressure1p.hh
View file @ a379dc18
......@@ -86,7 +86,6 @@ template<class TypeTag> class FVPressure1P: public FVPressure<TypeTag>
typedef typename GridView::Traits::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
typedef typename GridView::Intersection Intersection;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
......@@ -244,12 +243,12 @@ template<class TypeTag>
void FVPressure1P<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entry, const Intersection& intersection
, const CellData& cellData, const bool first)
{
ElementPointer elementI = intersection.inside();
ElementPointer elementJ = intersection.outside();
auto elementI = intersection.inside();
auto elementJ = intersection.outside();
// get global coordinates of cell centers
const GlobalPosition& globalPosI = elementI->geometry().center();
const GlobalPosition& globalPosJ = elementJ->geometry().center();
const GlobalPosition& globalPosI = elementI.geometry().center();
const GlobalPosition& globalPosJ = elementJ.geometry().center();
//get face normal
const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
......@@ -266,8 +265,8 @@ void FVPressure1P<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entry, const I
// compute vectorized permeabilities
DimMatrix meanPermeability(0);
problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(*elementI),
problem_.spatialParams().intrinsicPermeability(*elementJ));
problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(elementI),
problem_.spatialParams().intrinsicPermeability(elementJ));
Dune::FieldVector<Scalar, dim> permeability(0);
meanPermeability.mv(unitOuterNormal, permeability);
......@@ -294,10 +293,10 @@ template<class TypeTag>
void FVPressure1P<TypeTag>::getFluxOnBoundary(Dune::FieldVector<Scalar, 2>& entry,
const Intersection& intersection, const CellData& cellData, const bool first)
{
ElementPointer element = intersection.inside();
auto element = intersection.inside();
// get global coordinates of cell centers
const GlobalPosition& globalPosI = element->geometry().center();
const GlobalPosition& globalPosI = element.geometry().center();
// center of face in global coordinates
const GlobalPosition& globalPosJ = intersection.geometry().center();
......@@ -327,7 +326,7 @@ const Intersection& intersection, const CellData& cellData, const bool first)
DimMatrix meanPermeability(0);
problem_.spatialParams().meanK(meanPermeability,
problem_.spatialParams().intrinsicPermeability(*element));
problem_.spatialParams().intrinsicPermeability(element));
Dune::FieldVector<Scalar, dim> permeability(0);
meanPermeability.mv(unitOuterNormal, permeability);
......
dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh
View file @ a379dc18
......@@ -63,7 +63,6 @@ class FVVelocity1P
typedef typename GridView::Traits::template Codim<0>::Entity Element;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename GridView::Intersection Intersection;
typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
enum
{
......@@ -205,16 +204,16 @@ private:
template<class TypeTag>
void FVVelocity1P<TypeTag>::calculateVelocity(const Intersection& intersection, CellData& cellData)
{
ElementPointer elementI = intersection.inside();
ElementPointer elementJ = intersection.outside();
auto elementI = intersection.inside();
auto elementJ = intersection.outside();
int eIdxGlobalJ = problem_.variables().index(*elementJ);
int eIdxGlobalJ = problem_.variables().index(elementJ);
CellData& cellDataJ = problem_.variables().cellData(eIdxGlobalJ);
// get global coordinates of cell centers
const GlobalPosition& globalPosI = elementI->geometry().center();
const GlobalPosition& globalPosJ = elementJ->geometry().center();
const GlobalPosition& globalPosI = elementI.geometry().center();
const GlobalPosition& globalPosJ = elementJ.geometry().center();
//get face index
int isIndexI = intersection.indexInInside();
......@@ -232,8 +231,8 @@ void FVVelocity1P<TypeTag>::calculateVelocity(const Intersection& intersection,
// compute vectorized permeabilities
DimMatrix meanPermeability(0);
problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(*elementI),
problem_.spatialParams().intrinsicPermeability(*elementJ));
problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(elementI),
problem_.spatialParams().intrinsicPermeability(elementJ));
Dune::FieldVector < Scalar, dim > permeability(0);
meanPermeability.mv(unitOuterNormal, permeability);
......@@ -277,7 +276,7 @@ void FVVelocity1P<TypeTag>::calculateVelocity(const Intersection& intersection,
template<class TypeTag>
void FVVelocity1P<TypeTag>::calculateVelocityOnBoundary(const Intersection& intersection, CellData& cellData)
{
ElementPointer element = intersection.inside();
auto element = intersection.inside();
//get face index
int isIndex = intersection.indexInInside();
......@@ -295,7 +294,7 @@ void FVVelocity1P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte
problem_.dirichlet(boundValues, intersection);
// get global coordinates of cell centers
const GlobalPosition& globalPosI = element->geometry().center();
const GlobalPosition& globalPosI = element.geometry().center();
// center of face in global coordinates
const GlobalPosition& globalPosJ = intersection.geometry().center();
......@@ -310,7 +309,7 @@ void FVVelocity1P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte
// compute vectorized permeabilities
DimMatrix meanPermeability(0);
problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(*element));
problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(element));
//multiply with normal vector at the boundary
Dune::FieldVector < Scalar, dim > permeability(0);
......
dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh
View file @ a379dc18
......@@ -143,7 +143,6 @@ template<class TypeTag> class FVPressure2P: public FVPressure<TypeTag>
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::Traits::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
typedef typename GridView::Intersection Intersection;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
......@@ -684,14 +683,14 @@ template<class TypeTag>
void FVPressure2P<TypeTag>::getFlux(EntryType& entry, const Intersection& intersection
, const CellData& cellData, const bool first)
{
ElementPointer elementI = intersection.inside();
ElementPointer elementJ = intersection.outside();
auto elementI = intersection.inside();
auto elementJ = intersection.outside();
const CellData& cellDataJ = problem_.variables().cellData(problem_.variables().index(*elementJ));
const CellData& cellDataJ = problem_.variables().cellData(problem_.variables().index(elementJ));
// get global coordinates of cell centers
const GlobalPosition& globalPosI = elementI->geometry().center();
const GlobalPosition& globalPosJ = elementJ->geometry().center();
const GlobalPosition& globalPosI = elementI.geometry().center();
const GlobalPosition& globalPosJ = elementJ.geometry().center();
// get mobilities and fractional flow factors
Scalar lambdaWI = cellData.mobility(wPhaseIdx);
......@@ -718,8 +717,8 @@ void FVPressure2P<TypeTag>::getFlux(EntryType& entry, const Intersection& inters
// compute vectorized permeabilities
DimMatrix meanPermeability(0);
problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(*elementI),
problem_.spatialParams().intrinsicPermeability(*elementJ));
problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(elementI),
problem_.spatialParams().intrinsicPermeability(elementJ));
Dune::FieldVector<Scalar, dim> permeability(0);
meanPermeability.mv(unitOuterNormal, permeability);
......@@ -809,10 +808,10 @@ template<class TypeTag>
void FVPressure2P<TypeTag>::getFluxOnBoundary(EntryType& entry,
const Intersection& intersection, const CellData& cellData, const bool first)
{
ElementPointer element = intersection.inside();
auto element = intersection.inside();
// get global coordinates of cell centers
const GlobalPosition& globalPosI = element->geometry().center();
const GlobalPosition& globalPosI = element.geometry().center();
// center of face in global coordinates
const GlobalPosition& globalPosJ = intersection.geometry().center();
......@@ -854,7 +853,7 @@ const Intersection& intersection, const CellData& cellData, const bool first)
DimMatrix meanPermeability(0);
problem_.spatialParams().meanK(meanPermeability,
problem_.spatialParams().intrinsicPermeability(*element));
problem_.spatialParams().intrinsicPermeability(element));
Dune::FieldVector<Scalar, dim> permeability(0);
meanPermeability.mv(unitOuterNormal, permeability);
......@@ -885,13 +884,13 @@ const Intersection& intersection, const CellData& cellData, const bool first)
satW = cellData.saturation(wPhaseIdx);
satNw = cellData.saturation(nPhaseIdx);
}
Scalar temperature = problem_.temperature(*element);
Scalar temperature = problem_.temperature(element);
//get dirichlet pressure boundary condition
Scalar pressBound = boundValues[pressureIdx];
//calculate consitutive relations depending on the kind of saturation used
Scalar pcBound = MaterialLaw::pc(problem_.spatialParams().materialLawParams(*element), satW);
Scalar pcBound = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
//determine phase pressures from primary pressure variable
Scalar pressW = 0;
......@@ -932,9 +931,9 @@ const Intersection& intersection, const CellData& cellData, const bool first)
rhoMeanNw = 0.5 * (cellData.density(nPhaseIdx) + densityNwBound);
}
Scalar lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*element), satW)
Scalar lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
/ viscosityWBound;
Scalar lambdaNwBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*element), satW)
Scalar lambdaNwBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
/ viscosityNwBound;
Scalar fractionalWBound = lambdaWBound / (lambdaWBound + lambdaNwBound);
......
dumux/decoupled/2p/diffusion/fv/fvpressure2padaptive.hh
View file @ a379dc18
......@@ -80,7 +80,6 @@ template<class TypeTag> class FVPressure2PAdaptive: public FVPressure2P<TypeTag>
rhs = ParentType::rhs, matrix = ParentType::matrix
};
typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::Intersection Intersection;
typedef typename GridView::IntersectionIterator IntersectionIterator;
......@@ -251,30 +250,30 @@ template<class TypeTag>
void FVPressure2PAdaptive<TypeTag>::getFlux(EntryType& entry, const Intersection& intersection
, const CellData& cellData, const bool first)
{
ElementPointer elementI = intersection.inside();
ElementPointer elementJ = intersection.outside();
auto elementI = intersection.inside();
auto elementJ = intersection.outside();
if (elementI->level() == elementJ->level() || dim == 3)
if (elementI.level() == elementJ.level() || dim == 3)
{
ParentType::getFlux(entry, intersection, cellData, first);
//add the entry only once in case the VisitFacesOnlyOnce option is enabled!!!
if (GET_PROP_VALUE(TypeTag, VisitFacesOnlyOnce) && elementI->level() < elementJ->level())
if (GET_PROP_VALUE(TypeTag, VisitFacesOnlyOnce) && elementI.level() < elementJ.level())
{
entry = 0.;
}
}
// hanging node situation: neighbor has higher level
else if (elementJ->level() == elementI->level() + 1)
else if (elementJ.level() == elementI.level() + 1)
{
const CellData& cellDataJ = problem_.variables().cellData(problem_.variables().index(*elementJ));
const CellData& cellDataJ = problem_.variables().cellData(problem_.variables().index(elementJ));
// get global coordinates of cell centers
const GlobalPosition& globalPosI = elementI->geometry().center();
const GlobalPosition& globalPosJ = elementJ->geometry().center();
const GlobalPosition& globalPosI = elementI.geometry().center();
const GlobalPosition& globalPosJ = elementJ.geometry().center();
int globalIdxI = problem_.variables().index(*elementI);
int globalIdxJ = problem_.variables().index(*elementJ);
int globalIdxI = problem_.variables().index(elementI);
int globalIdxJ = problem_.variables().index(elementJ);
// get mobilities and fractional flow factors
Scalar lambdaWI = cellData.mobility(wPhaseIdx);
......@@ -302,24 +301,24 @@ void FVPressure2PAdaptive<TypeTag>::getFlux(EntryType& entry, const Intersection
// Count number of hanging nodes
// not really necessary
int globalIdxK = 0;
ElementPointer elementK = intersection.outside();
auto elementK = intersection.outside();
// We are looking for two things:
// IsIndexJ, the index of the interface from the neighbor-cell point of view
// GlobalIdxK, the index of the third cell
// Intersectioniterator around cell I
IntersectionIterator isItEndI = problem_.gridView().iend(*elementI);
for (IntersectionIterator isItI = problem_.gridView().ibegin(*elementI); isItI != isItEndI; ++isItI)
IntersectionIterator isItEndI = problem_.gridView().iend(elementI);
for (IntersectionIterator isItI = problem_.gridView().ibegin(elementI); isItI != isItEndI; ++isItI)
{
if (isItI->neighbor())
{
ElementPointer neighborPointer2 = isItI->outside();
auto neighbor2 = isItI->outside();
// make sure we do not choose elemntI as third element
// -> faces with hanging node have more than one intersection but only one face index!
if (neighborPointer2 != elementJ && isItI->indexInInside() == isIndexI)
if (neighbor2 != elementJ && isItI->indexInInside() == isIndexI)
{
globalIdxK = problem_.variables().index(*neighborPointer2);
elementK = neighborPointer2;
globalIdxK = problem_.variables().index(neighbor2);
elementK = neighbor2;
break;
}
}
......@@ -341,11 +340,11 @@ void FVPressure2PAdaptive<TypeTag>::getFlux(EntryType& entry, const Intersection
FieldMatrix permeabilityK(0);
problem_.spatialParams().meanK(permeabilityI,
problem_.spatialParams().intrinsicPermeability(*elementI));
problem_.spatialParams().intrinsicPermeability(elementI));
problem_.spatialParams().meanK(permeabilityJ,
problem_.spatialParams().intrinsicPermeability(*elementJ));
problem_.spatialParams().intrinsicPermeability(elementJ));
problem_.spatialParams().meanK(permeabilityK,
problem_.spatialParams().intrinsicPermeability(*elementK));
problem_.spatialParams().intrinsicPermeability(elementK));
// Calculate permeablity component normal to interface
Scalar kI, kJ, kK, kMean, ng;
......
dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
View file @ a379dc18
......@@ -76,7 +76,6 @@ class FVVelocity2P
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::Traits::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
typedef typename GridView::Intersection Intersection;
typedef typename GridView::IntersectionIterator IntersectionIterator;
......@@ -332,16 +331,16 @@ private:
template<class TypeTag>
void FVVelocity2P<TypeTag>::calculateVelocity(const Intersection& intersection, CellData& cellData)
{
ElementPointer elementI = intersection.inside();
ElementPointer elementJ = intersection.outside();
auto elementI = intersection.inside();
auto elementJ = intersection.outside();
int eIdxGlobalJ = problem_.variables().index(*elementJ);
int eIdxGlobalJ = problem_.variables().index(elementJ);
CellData& cellDataJ = problem_.variables().cellData(eIdxGlobalJ);
// get global coordinates of cell centers
const GlobalPosition& globalPosI = (*elementI).geometry().center();
const GlobalPosition& globalPosJ = (*elementJ).geometry().center();
const GlobalPosition& globalPosI = (elementI).geometry().center();
const GlobalPosition& globalPosJ = (elementJ).geometry().center();
// get mobilities and fractional flow factors
Scalar lambdaWI = cellData.mobility(wPhaseIdx);
......@@ -369,8 +368,8 @@ void FVVelocity2P<TypeTag>::calculateVelocity(const Intersection& intersection,
// compute vectorized permeabilities
DimMatrix meanPermeability(0);
problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(*elementI),
problem_.spatialParams().intrinsicPermeability(*elementJ));
problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(elementI),
problem_.spatialParams().intrinsicPermeability(elementJ));
Dune::FieldVector<Scalar, dim> permeability(0);
meanPermeability.mv(unitOuterNormal, permeability);
......@@ -491,7 +490,7 @@ void FVVelocity2P<TypeTag>::calculateVelocity(const Intersection& intersection,
template<class TypeTag>
void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& intersection, CellData& cellData)
{
ElementPointer element = intersection.inside();
auto element = intersection.inside();
//get face index
int isIndex = intersection.indexInInside();
......@@ -509,7 +508,7 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte
problem_.dirichlet(boundValues, intersection);
// get global coordinates of cell centers
const GlobalPosition& globalPosI = (*element).geometry().center();
const GlobalPosition& globalPosI = element.geometry().center();
// center of face in global coordinates
const GlobalPosition& globalPosJ = intersection.geometry().center();
......@@ -533,7 +532,7 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte
// compute vectorized permeabilities
DimMatrix meanPermeability(0);
problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(*element));
problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(element));
Dune::FieldVector<Scalar, dim> permeability(0);
meanPermeability.mv(unitOuterNormal, permeability);
......@@ -566,7 +565,7 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte
}
Scalar pressBound = boundValues[pressureIdx];
Scalar pcBound = MaterialLaw::pc(problem_.spatialParams().materialLawParams(*element), satW);
Scalar pcBound = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
//determine phase pressures from primary pressure variable
Scalar pressWBound = 0;
......@@ -583,7 +582,7 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte
}
//get temperature at current position
Scalar temperature = problem_.temperature(*element);
Scalar temperature = problem_.temperature(element);
Scalar densityWBound = density_[wPhaseIdx];
Scalar densityNwBound = density_[nPhaseIdx];
......@@ -605,9 +604,9 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte
viscosityNwBound = FluidSystem::viscosity(fluidState, nPhaseIdx) / densityNwBound;
}
Scalar lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*element), satW)
Scalar lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
/ viscosityWBound;
Scalar lambdaNwBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*element), satW)
Scalar lambdaNwBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
/ viscosityNwBound;
Scalar potentialDiffW = cellData.fluxData().upwindPotential(wPhaseIdx, isIndex);
......
dumux/decoupled/2p/diffusion/fv/fvvelocity2padaptive.hh
View file @ a379dc18
......@@ -54,7 +54,6 @@ class FVVelocity2PAdaptive: public FVVelocity2P<TypeTag>
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename GridView::Intersection Intersection;
typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
enum
{
......@@ -160,20 +159,20 @@ private:
template<class TypeTag>
void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& intersection, CellData& cellData)
{
ElementPointer elementI = intersection.inside();
ElementPointer elementJ = intersection.outside();
auto elementI = intersection.inside();
auto elementJ = intersection.outside();
if (elementI->level() == elementJ->level())
if (elementI.level() == elementJ.level())
{
ParentType::calculateVelocity(intersection, cellData);
}
else if (elementJ->level() == elementI->level() + 1 && dim == 2)
else if (elementJ.level() == elementI.level() + 1 && dim == 2)
{
CellData& cellDataJ = problem_.variables().cellData(problem_.variables().index(*elementJ));
CellData& cellDataJ = problem_.variables().cellData(problem_.variables().index(elementJ));
// get global coordinates of cell centers
const GlobalPosition& globalPosI = elementI->geometry().center();
const GlobalPosition& globalPosJ = elementJ->geometry().center();
const GlobalPosition& globalPosI = elementI.geometry().center();
const GlobalPosition& globalPosJ = elementJ.geometry().center();
// get mobilities and fractional flow factors
Scalar lambdaWI = cellData.mobility(wPhaseIdx);
......@@ -203,26 +202,26 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters
int isIndexJ = intersection.indexInOutside();
int globalIdxK = 0;
ElementPointer elementK = intersection.outside();
auto elementK = intersection.outside();
// We are looking for two things:
// IsIndexJ, the index of the interface from the neighbor-cell point of view
// GlobalIdxK, the index of the third cell
// for efficienty this is done in one IntersectionIterator-Loop
// Intersectioniterator around cell I
IntersectionIterator isItEndI = problem_.gridView().iend(*elementI);
for (IntersectionIterator isItI = problem_.gridView().ibegin(*elementI); isItI != isItEndI; ++isItI)
IntersectionIterator isItEndI = problem_.gridView().iend(elementI);
for (IntersectionIterator isItI = problem_.gridView().ibegin(elementI); isItI != isItEndI; ++isItI)
{
if (isItI->neighbor())
{
ElementPointer neighborPointer2 = isItI->outside();
auto neighbor2 = isItI->outside();
// make sure we do not choose elemntI as third element
// -> faces with hanging node have more than one intersection but only one face index!
if (neighborPointer2 != elementJ && isItI->indexInInside() == isIndexI)
if (neighbor2 != elementJ && isItI->indexInInside() == isIndexI)
{
globalIdxK = problem_.variables().index(*neighborPointer2);
elementK = neighborPointer2;
globalIdxK = problem_.variables().index(neighbor2);
elementK = neighbor2;
faceAreaSum += isItI->geometry().volume();
break;
......@@ -246,11 +245,11 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters
DimMatrix permeabilityK(0);
problem_.spatialParams().meanK(permeabilityI,
problem_.spatialParams().intrinsicPermeability(*elementI));
problem_.spatialParams().intrinsicPermeability(elementI));
problem_.spatialParams().meanK(permeabilityJ,
problem_.spatialParams().intrinsicPermeability(*elementJ));
problem_.spatialParams().intrinsicPermeability(elementJ));
problem_.spatialParams().meanK(permeabilityK,
problem_.spatialParams().intrinsicPermeability(*elementK));
problem_.spatialParams().intrinsicPermeability(elementK));
// Calculate permeablity component normal to interface
Scalar kI, kJ, kK, ng, kMean; //, gI, gJ, gK;
......@@ -449,18 +448,18 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters
cellData.fluxData().addVelocity(wPhaseIdx, isIndexI, velocityW);
cellData.fluxData().addVelocity(nPhaseIdx, isIndexI, velocityNw);
}
else if (elementI->level() > elementJ->level() && dim == 3)
else if (elementI.level() > elementJ.level() && dim == 3)
{
ElementPointer elementI = intersection.inside();
ElementPointer elementJ = intersection.outside();
auto elementI = intersection.inside();
auto elementJ = intersection.outside();
int globalIdxJ = problem_.variables().index(*elementJ);
int globalIdxJ = problem_.variables().index(elementJ);
CellData& cellDataJ = problem_.variables().cellData(globalIdxJ);
// get global coordinates of cell centers
const GlobalPosition& globalPosI = elementI->geometry().center();
const GlobalPosition& globalPosJ = elementJ->geometry().center();
const GlobalPosition& globalPosI = elementI.geometry().center();
const GlobalPosition& globalPosJ = elementJ.geometry().center();
// get mobilities and fractional flow factors
Scalar lambdaWI = cellData.mobility(wPhaseIdx);
......@@ -488,8 +487,8 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters
// compute vectorized permeabilities
DimMatrix meanPermeability(0);
problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(*elementI),
problem_.spatialParams().intrinsicPermeability(*elementJ));
problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(elementI),
problem_.spatialParams().intrinsicPermeability(elementJ));
Dune::FieldVector<Scalar, dim> permeability(0);
meanPermeability.mv(unitOuterNormal, permeability);
......@@ -593,7 +592,7 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters
cellData.fluxData().setVelocity(nPhaseIdx, isIndexI, velocityNw);
cellData.fluxData().setVelocityMarker(isIndexI);
Scalar weightingFactor = std::pow(0.5, (dim - 1)*(elementI->level() - elementJ->level()));
Scalar weightingFactor = std::pow(0.5, (dim - 1)*(elementI.level() - elementJ.level()));
velocityW *= weightingFactor;
velocityNw *= weightingFactor;
......
dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dpressurevelocity2p.hh
View file @ a379dc18
......@@ -76,7 +76,6 @@ template<class TypeTag> class FvMpfaL2dPressureVelocity2p: public FvMpfaL2dPress
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
typedef typename GridView::Intersection Intersection;
typedef typename ParentType::InteractionVolume InteractionVolume;
......@@ -231,16 +230,16 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocity()
if (interactionVolume.isInnerVolume())
{
ElementPointer & elementPointer1 = interactionVolume.getSubVolumeElement(0);
ElementPointer & elementPointer2 = interactionVolume.getSubVolumeElement(1);
ElementPointer & elementPointer3 = interactionVolume.getSubVolumeElement(2);
ElementPointer & elementPointer4 = interactionVolume.getSubVolumeElement(3);
auto element1 = interactionVolume.getSubVolumeElement(0);
auto element2 = interactionVolume.getSubVolumeElement(1);
auto element3 = interactionVolume.getSubVolumeElement(2);
auto element4 = interactionVolume.getSubVolumeElement(3);
// cell index
int eIdxGlobal1 = problem_.variables().index(*elementPointer1);
int eIdxGlobal2 = problem_.variables().index(*elementPointer2);
int eIdxGlobal3 = problem_.variables().index(*elementPointer3);
int eIdxGlobal4 = problem_.variables().index(*elementPointer4);
int eIdxGlobal1 = problem_.variables().index(element1);
int eIdxGlobal2 = problem_.variables().index(element2);
int eIdxGlobal3 = problem_.variables().index(element3);
int eIdxGlobal4 = problem_.variables().index(element4);
//get the cell Data
CellData& cellData1 = problem_.variables().cellData(eIdxGlobal1);
......@@ -271,10 +270,8 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocity()
{
continue;
}
ElementPointer & elementPointer = interactionVolume.getSubVolumeElement(elemIdx);