Commit 1c3a5de5 authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

[range-for] adapt intersection loops in dumux/decoupled

For the adaptive decoupled models, several loops could not be
adapted. In particular, incrementation of iterators is used not
only for loop steering. Moreover, iterators are explicitly stored.
parent ca320997
......@@ -47,21 +47,20 @@ template<class TypeTag>
class FVVelocity1P
{
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
typedef typename GET_PROP_TYPE(TypeTag, Fluid) Fluid;
typedef typename GET_PROP_TYPE(TypeTag, Fluid) Fluid;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
typedef typename SolutionTypes::PrimaryVariables PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::Traits::template Codim<0>::Entity Element;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename GridView::Traits::template Codim<0>::Entity Element;
typedef typename GridView::Intersection Intersection;
enum
......@@ -132,16 +131,12 @@ public:
std::vector<Scalar> flux(numberOfFaces,0);
// run through all intersections with neighbors and boundary
IntersectionIterator
isEndIt = problem_.gridView().iend(*eIt);
for (IntersectionIterator
isIt = problem_.gridView().ibegin(*eIt); isIt
!=isEndIt; ++isIt)
for (const auto& intersection : Dune::intersections(problem_.gridView(), *eIt))
{
int isIndex = isIt->indexInInside();
int isIndex = intersection.indexInInside();
flux[isIndex] = isIt->geometry().volume()
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(isIndex));
flux[isIndex] = intersection.geometry().volume()
* (intersection.centerUnitOuterNormal() * cellData.fluxData().velocity(isIndex));
}
// calculate velocity on reference element as the Raviart-Thomas-0
......
......@@ -82,7 +82,6 @@ template<class TypeTag> class FVPressure2PAdaptive: public FVPressure2P<TypeTag>
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::Intersection Intersection;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldMatrix<Scalar, dim, dim> FieldMatrix;
......@@ -306,16 +305,15 @@ void FVPressure2PAdaptive<TypeTag>::getFlux(EntryType& entry, const Intersection
// 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)
for (const auto& intersectionI : Dune::intersections(problem_.gridView(), elementI))
{
if (isItI->neighbor())
if (intersectionI.neighbor())
{
auto neighbor2 = isItI->outside();
auto neighbor2 = intersectionI.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 (neighbor2 != elementJ && isItI->indexInInside() == isIndexI)
if (neighbor2 != elementJ && intersectionI.indexInInside() == isIndexI)
{
globalIdxK = problem_.variables().index(neighbor2);
elementK = neighbor2;
......
......@@ -77,7 +77,6 @@ class FVVelocity2P
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::Traits::template Codim<0>::Entity Element;
typedef typename GridView::Intersection Intersection;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename Element::Geometry Geometry;
typedef typename Geometry::JacobianTransposed JacobianTransposed;
......@@ -209,15 +208,14 @@ public:
std::vector<Scalar> fluxNw(numberOfFaces,0);
// run through all intersections with neighbors and boundary
IntersectionIterator isEndIt = problem_.gridView().iend(*eIt);
for (IntersectionIterator isIt = problem_.gridView().ibegin(*eIt); isIt != isEndIt; ++isIt)
for (const auto& intersection : Dune::intersections(problem_.gridView(), *eIt))
{
int isIndex = isIt->indexInInside();
int isIndex = intersection.indexInInside();
fluxW[isIndex] += isIt->geometry().volume()
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(wPhaseIdx, isIndex));
fluxNw[isIndex] += isIt->geometry().volume()
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
fluxW[isIndex] += intersection.geometry().volume()
* (intersection.centerUnitOuterNormal() * cellData.fluxData().velocity(wPhaseIdx, isIndex));
fluxNw[isIndex] += intersection.geometry().volume()
* (intersection.centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
}
// calculate velocity on reference element as the Raviart-Thomas-0
......
......@@ -52,7 +52,6 @@ class FVVelocity2PAdaptive: public FVVelocity2P<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename GridView::Intersection Intersection;
enum
......@@ -209,20 +208,19 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters
// 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)
for (const auto& intersectionI : Dune::intersections(problem_.gridView(), elementI))
{
if (isItI->neighbor())
if (intersectionI.neighbor())
{
auto neighbor2 = isItI->outside();
auto neighbor2 = intersectionI.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 (neighbor2 != elementJ && isItI->indexInInside() == isIndexI)
if (neighbor2 != elementJ && intersectionI.indexInInside() == isIndexI)
{
globalIdxK = problem_.variables().index(neighbor2);
elementK = neighbor2;
faceAreaSum += isItI->geometry().volume();
faceAreaSum += intersectionI.geometry().volume();
break;
}
......
......@@ -85,7 +85,6 @@ template<class TypeTag> class FvMpfaL2dVelocity2p
typedef typename GridView::IndexSet IndexSet;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename Element::Geometry Geometry;
typedef typename Geometry::JacobianTransposed JacobianTransposed;
......@@ -197,15 +196,14 @@ public:
Dune::FieldVector < Scalar, 2 * dim > fluxNw(0);
// run through all intersections with neighbors and boundary
IntersectionIterator isEndIt = problem_.gridView().iend(*eIt);
for (IntersectionIterator isIt = problem_.gridView().ibegin(*eIt); isIt != isEndIt; ++isIt)
for (const auto& intersection : Dune::intersections(problem_.gridView(), *eIt))
{
int isIndex = isIt->indexInInside();
int isIndex = intersection.indexInInside();
fluxW[isIndex] += isIt->geometry().volume()
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(wPhaseIdx, isIndex));
fluxNw[isIndex] += isIt->geometry().volume()
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
fluxW[isIndex] += intersection.geometry().volume()
* (intersection.centerUnitOuterNormal() * cellData.fluxData().velocity(wPhaseIdx, isIndex));
fluxNw[isIndex] += intersection.geometry().volume()
* (intersection.centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
}
DimVector refVelocity(0);
......
......@@ -84,7 +84,6 @@ template<class TypeTag> class FvMpfaL3dVelocity2p
typedef typename GridView::IndexSet IndexSet;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename Element::Geometry Geometry;
typedef typename Geometry::JacobianTransposed JacobianTransposed;
......@@ -212,15 +211,14 @@ public:
Dune::FieldVector < Scalar, 2 * dim > fluxNw(0);
// run through all intersections with neighbors and boundary
IntersectionIterator isEndIt = problem_.gridView().iend(*eIt);
for (IntersectionIterator isIt = problem_.gridView().ibegin(*eIt); isIt != isEndIt; ++isIt)
for (const auto& intersection : Dune::intersections(problem_.gridView(), *eIt))
{
int isIndex = isIt->indexInInside();
int isIndex = intersection.indexInInside();
fluxW[isIndex] += isIt->geometry().volume()
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(wPhaseIdx, isIndex));
fluxNw[isIndex] += isIt->geometry().volume()
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
fluxW[isIndex] += intersection.geometry().volume()
* (intersection.centerUnitOuterNormal() * cellData.fluxData().velocity(wPhaseIdx, isIndex));
fluxNw[isIndex] += intersection.geometry().volume()
* (intersection.centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
}
DimVector refVelocity(0);
......
......@@ -84,7 +84,6 @@ template<class TypeTag> class FvMpfaO2dVelocity2P
typedef typename GridView::IndexSet IndexSet;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename Element::Geometry Geometry;
typedef typename Geometry::JacobianTransposed JacobianTransposed;
......@@ -195,15 +194,14 @@ public:
Dune::FieldVector < Scalar, 2 * dim > fluxNw(0);
// run through all intersections with neighbors and boundary
IntersectionIterator isEndIt = problem_.gridView().iend(*eIt);
for (IntersectionIterator isIt = problem_.gridView().ibegin(*eIt); isIt != isEndIt; ++isIt)
for (const auto& intersection : Dune::intersections(problem_.gridView(), *eIt))
{
int isIndex = isIt->indexInInside();
int isIndex = intersection.indexInInside();
fluxW[isIndex] += isIt->geometry().volume()
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(wPhaseIdx, isIndex));
fluxNw[isIndex] += isIt->geometry().volume()
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
fluxW[isIndex] += intersection.geometry().volume()
* (intersection.centerUnitOuterNormal() * cellData.fluxData().velocity(wPhaseIdx, isIndex));
fluxNw[isIndex] += intersection.geometry().volume()
* (intersection.centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
}
DimVector refVelocity(0);
......
......@@ -100,7 +100,6 @@ class MimeticTwoPLocalStiffness: public LocalStiffness<TypeTag, 1>
typedef typename GridView::Grid Grid;
typedef typename GridView::Traits::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
......@@ -270,10 +269,9 @@ public:
int eIdxGlobal = problem_.variables().index(element);
Dune::FieldVector<Scalar, 2 * dim> faceVol(0);
IntersectionIterator isEndIt = gridView_.iend(element);
for (IntersectionIterator isIt = gridView_.ibegin(element); isIt != isEndIt; ++isIt)
for (const auto& intersection : Dune::intersections(gridView_, element))
{
faceVol[isIt->indexInInside()] = isIt->geometry().volume();
faceVol[intersection.indexInInside()] = intersection.geometry().volume();
}
vel = 0;
......@@ -288,10 +286,9 @@ public:
int eIdxGlobal = problem_.variables().index(element);
Dune::FieldVector<Scalar, 2 * dim> faceVol(0);
IntersectionIterator isEndIt = gridView_.iend(element);
for (IntersectionIterator isIt = gridView_.ibegin(element); isIt != isEndIt; ++isIt)
for (const auto& intersection : Dune::intersections(gridView_, element))
{
faceVol[isIt->indexInInside()] = isIt->geometry().volume();
faceVol[intersection.indexInInside()] = intersection.geometry().volume();
}
vel = 0;
......@@ -306,13 +303,12 @@ public:
Dune::FieldVector<Scalar, 2 * dim> c(0);
Dune::FieldMatrix<Scalar, 2 * dim, 2 * dim> Pi(0);
IntersectionIterator isEndIt = gridView_.iend(element);
for (IntersectionIterator isIt = gridView_.ibegin(element); isIt != isEndIt; ++isIt)
for (const auto& intersection : Dune::intersections(gridView_, element))
{
// local number of facet
int idx = isIt->indexInInside();
int idx = intersection.indexInInside();
Scalar faceVol = isIt->geometry().volume();
Scalar faceVol = intersection.geometry().volume();
// Corresponding to the element under consideration,
// calculate the part of the matrix C coupling velocities and element pressures.
......@@ -477,17 +473,16 @@ void MimeticTwoPLocalStiffness<TypeTag>::assembleElementMatrices(const Element&
Scalar gravPot = (problem_.bBoxMax() - centerGlobal) * problem_.gravity() * (density_[nPhaseIdx] - density_[wPhaseIdx]);
int i = -1;
IntersectionIterator isEndIt = gridView_.iend(element);
for (IntersectionIterator isIt = gridView_.ibegin(element); isIt != isEndIt; ++isIt)
for (const auto& intersection : Dune::intersections(gridView_, element))
{
// local number of facet
i = isIt->indexInInside();
i = intersection.indexInInside();
Dune::FieldVector<Scalar, dim> faceGlobal = isIt->geometry().center();
faceVol[i] = isIt->geometry().volume();
Dune::FieldVector<Scalar, dim> faceGlobal = intersection.geometry().center();
faceVol[i] = intersection.geometry().volume();
// get normal vector
const Dune::FieldVector<Scalar, dim>& unitOuterNormal = isIt->centerUnitOuterNormal();
const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
N[i] = unitOuterNormal;
......@@ -613,9 +608,9 @@ void MimeticTwoPLocalStiffness<TypeTag>::assembleElementMatrices(const Element&
// std::cout << "Pi = \dim" << Pi << "c = " << c << ", F = " << F << std::endl;
//accumulate fluxes due to capillary potential (pc + gravity!)
for (IntersectionIterator isIt = gridView_.ibegin(element); isIt != isEndIt; ++isIt)
for (const auto& intersection : Dune::intersections(gridView_, element))
{
int idx = isIt->indexInInside();
int idx = intersection.indexInInside();
Scalar fracFlow = 0;
......@@ -624,9 +619,9 @@ void MimeticTwoPLocalStiffness<TypeTag>::assembleElementMatrices(const Element&
flux += W_[eIdxGlobal][idx][j] * faceVol[j] * (gravPot - gravPotFace[j]);
//it is enough to evaluate the capillary/gravity flux for neighbors -> not needed anyway at the boundaries!
if (isIt->neighbor())
if (intersection.neighbor())
{
int eIdxGlobalNeighbor = problem_.variables().index(isIt->outside());
int eIdxGlobalNeighbor = problem_.variables().index(intersection.outside());
if (flux >= 0.)
{
switch (pressureType)
......@@ -648,10 +643,10 @@ void MimeticTwoPLocalStiffness<TypeTag>::assembleElementMatrices(const Element&
rhs_[eIdxGlobalNeighbor] += (faceVol[idx] * fracFlow * flux);
}
}
else if (isIt->boundary())
else if (intersection.boundary())
{
BoundaryTypes bctype;
problem_.boundaryTypes(bctype, *isIt);
problem_.boundaryTypes(bctype, intersection);
if (bctype.isDirichlet(pressureEqIdx))
{
......@@ -674,7 +669,7 @@ void MimeticTwoPLocalStiffness<TypeTag>::assembleElementMatrices(const Element&
else if (flux < 0. && bctype.isDirichlet(satEqIdx))
{
PrimaryVariables boundValues(0.0);
problem_.dirichlet(boundValues, *isIt);
problem_.dirichlet(boundValues, intersection);
Scalar krw = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element),
boundValues[saturationIdx]);
......@@ -714,32 +709,29 @@ template<class TypeTag>
void MimeticTwoPLocalStiffness<TypeTag>::assembleBC(const Element& element, int k)
{
// evaluate boundary conditions via intersection iterator
typedef typename GridView::IntersectionIterator IntersectionIterator;
IntersectionIterator endIsIt = gridView_.iend(element);
for (IntersectionIterator isIt = gridView_.ibegin(element); isIt != endIsIt; ++isIt)
for (const auto& intersection : Dune::intersections(gridView_, element))
{
int faceIndex = isIt->indexInInside();
if (isIt->boundary())
int faceIndex = intersection.indexInInside();
if (intersection.boundary())
{
problem_.boundaryTypes(this->bctype[faceIndex], *isIt);
problem_.boundaryTypes(this->bctype[faceIndex], intersection);
PrimaryVariables boundValues(0.0);
if (this->bctype[faceIndex].isNeumann(pressureEqIdx))
{
problem_.neumann(boundValues, *isIt);
problem_.neumann(boundValues, intersection);
Scalar J = (boundValues[wPhaseIdx]/density_[wPhaseIdx] + boundValues[nPhaseIdx]/density_[nPhaseIdx]);
this->b[faceIndex] -= J * isIt->geometry().volume();
this->b[faceIndex] -= J * intersection.geometry().volume();
}
else if (this->bctype[faceIndex].isDirichlet(pressureEqIdx))
{
problem_.dirichlet(boundValues, *isIt);
problem_.dirichlet(boundValues, intersection);
if (pressureType == pw)
this->b[faceIndex] = boundValues[pressureIdx] +
(problem_.bBoxMax() - isIt->geometry().center()) * problem_.gravity() * density_[wPhaseIdx];
(problem_.bBoxMax() - intersection.geometry().center()) * problem_.gravity() * density_[wPhaseIdx];
else if (pressureType == pn)
this->b[faceIndex] = boundValues[pressureIdx] +
(problem_.bBoxMax() - isIt->geometry().center()) * problem_.gravity() * density_[nPhaseIdx];
(problem_.bBoxMax() - intersection.geometry().center()) * problem_.gravity() * density_[nPhaseIdx];
else
this->b[faceIndex] = boundValues[pressureIdx];
}
......
......@@ -100,7 +100,6 @@ class MimeticTwoPLocalStiffnessAdaptive: public LocalStiffness<TypeTag, 1>
typedef typename GridView::Grid Grid;
typedef typename GridView::Traits::template Codim<0>::Entity Element;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
......@@ -296,11 +295,10 @@ public:
Dune::FieldVector<Scalar, 2*dim> faceVolumeReal(0.0);
int idx = 0;
IntersectionIterator isEndIt = gridView_.iend(element);
for (IntersectionIterator isIt = gridView_.ibegin(element); isIt != isEndIt; ++isIt)
for (const auto& intersection : Dune::intersections(gridView_, element))
{
faceVol[idx] = isIt->geometry().volume();
int indexInInside = isIt->indexInInside();
faceVol[idx] = intersection.geometry().volume();
int indexInInside = intersection.indexInInside();
faceVolumeReal[indexInInside] += faceVol[idx];
idx++;
......@@ -329,10 +327,9 @@ public:
Dune::DynamicMatrix<Scalar> Pi(numFaces, numFaces);
int idx = 0;
IntersectionIterator isEndIt = gridView_.iend(element);
for (IntersectionIterator isIt = gridView_.ibegin(element); isIt != isEndIt; ++isIt)
for (const auto& intersection : Dune::intersections(gridView_, element))
{
Scalar faceVol = isIt->geometry().volume();
Scalar faceVol = intersection.geometry().volume();
// Corresponding to the element under consideration,
// calculate the part of the matrix C coupling velocities and element pressures.
......@@ -498,16 +495,15 @@ void MimeticTwoPLocalStiffnessAdaptive<TypeTag>::assembleElementMatrices(const E
* (density_[nPhaseIdx] - density_[wPhaseIdx]);
int idx = 0;
IntersectionIterator isEndIt = gridView_.iend(element);
for (IntersectionIterator isIt = gridView_.ibegin(element); isIt != isEndIt; ++isIt)
for (const auto& intersection : Dune::intersections(gridView_, element))
{
// local number of facet
Dune::FieldVector<Scalar, dim> faceGlobal(isIt->geometry().center());
faceVol[idx] = isIt->geometry().volume();
Dune::FieldVector<Scalar, dim> faceGlobal(intersection.geometry().center());
faceVol[idx] = intersection.geometry().volume();
// get normal vector
const Dune::FieldVector<Scalar, dim>& unitOuterNormal = isIt->centerUnitOuterNormal();
const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
for (int k = 0; k < dim; k++)
......@@ -640,7 +636,7 @@ void MimeticTwoPLocalStiffnessAdaptive<TypeTag>::assembleElementMatrices(const E
//accumulate fluxes due to capillary potential (pc + gravity!)
idx = 0;
for (IntersectionIterator isIt = gridView_.ibegin(element); isIt != isEndIt; ++isIt)
for (const auto& intersection : Dune::intersections(gridView_, element))
{
Scalar fracFlow = 0;
......@@ -648,10 +644,10 @@ void MimeticTwoPLocalStiffnessAdaptive<TypeTag>::assembleElementMatrices(const E
for (int j = 0; j < numFaces; j++)
flux += W_[eIdxGlobal][idx][j] * faceVol[j] * (pcPot - pcPotFace[j]);
if (isIt->neighbor())
if (intersection.neighbor())
{
int eIdxGlobalNeighbor = problem_.variables().index(isIt->outside());
int eIdxGlobalNeighbor = problem_.variables().index(intersection.outside());
if (flux > 0.)
{
switch (pressureType)
......@@ -672,10 +668,10 @@ void MimeticTwoPLocalStiffnessAdaptive<TypeTag>::assembleElementMatrices(const E
rhs_[eIdxGlobalNeighbor] += (faceVol[idx] * fracFlow * flux);
}
}
else if (isIt->boundary())
else if (intersection.boundary())
{
BoundaryTypes bctype;
problem_.boundaryTypes(bctype, *isIt);
problem_.boundaryTypes(bctype, intersection);
if (bctype.isDirichlet(pressureEqIdx))
{
......@@ -698,7 +694,7 @@ void MimeticTwoPLocalStiffnessAdaptive<TypeTag>::assembleElementMatrices(const E
else if (flux < 0. && bctype.isDirichlet(satEqIdx))
{
PrimaryVariables boundValues(0.0);
problem_.dirichlet(boundValues, *isIt);
problem_.dirichlet(boundValues, intersection);
Scalar krw = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element),
boundValues[saturationIdx]);
......@@ -742,33 +738,32 @@ void MimeticTwoPLocalStiffnessAdaptive<TypeTag>::assembleBC(const Element& eleme
typedef typename GridView::IntersectionIterator IntersectionIterator;
unsigned int faceIndex = 0;
IntersectionIterator endIsIt = gridView_.iend(element);
for (IntersectionIterator isIt = gridView_.ibegin(element); isIt != endIsIt; ++isIt)
for (const auto& intersection : Dune::intersections(gridView_, element))
{
if (isIt->neighbor())
if (intersection.neighbor())
{
}
else if (isIt->boundary())
else if (intersection.boundary())
{
problem_.boundaryTypes(this->bctype[faceIndex], *isIt);
problem_.boundaryTypes(this->bctype[faceIndex], intersection);
PrimaryVariables boundValues(0.0);
if (this->bctype[faceIndex].isNeumann(pressureEqIdx))
{
problem_.neumann(boundValues, *isIt);
problem_.neumann(boundValues, intersection);
Scalar J = (boundValues[wPhaseIdx]/density_[wPhaseIdx] + boundValues[nPhaseIdx]/density_[nPhaseIdx]);
this->b[faceIndex] -= J * isIt->geometry().volume();
this->b[faceIndex] -= J * intersection.geometry().volume();
}
else if (this->bctype[faceIndex].isDirichlet(pressureEqIdx))
{
problem_.dirichlet(boundValues, *isIt);
problem_.dirichlet(boundValues, intersection);
if (pressureType == pw)
this->b[faceIndex] = boundValues[pressureIdx] + (problem_.bBoxMax()
- isIt->geometry().center()) * problem_.gravity() * density_[wPhaseIdx];
- intersection.geometry().center()) * problem_.gravity() * density_[wPhaseIdx];
else if (pressureType == pn)
this->b[faceIndex] = boundValues[pressureIdx] + (problem_.bBoxMax()
- isIt->geometry().center()) * problem_.gravity() * density_[nPhaseIdx];
- intersection.geometry().center()) * problem_.gravity() * density_[nPhaseIdx];
else
this->b[faceIndex] = boundValues[pressureIdx];
}
......
......@@ -60,7 +60,6 @@ class MimeticOperatorAssemblerTwoP: public CROperatorAssemblerTwoP<TypeTag>
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
......@@ -129,17 +128,15 @@ public:
FieldVector globalPos = eIt->geometry().center();