diff --git a/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh b/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh
index a5679dca9f3713179547d2aaae062cecccabe465..3c9aa18881eb52220f04166d4a2f69926c2e196e 100644
--- a/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh
+++ b/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh
@@ -123,7 +123,15 @@ public:
CellData& cellData = problem_.variables().cellData(globalIdx);
- Dune::FieldVector<Scalar, 2*dim> flux(0);
+ const typename Element::Geometry& geometry = eIt->geometry();
+ // get corresponding reference element
+ typedef Dune::ReferenceElements<Scalar, dim> ReferenceElements;
+ const Dune::ReferenceElement< Scalar , dim > & refElement =
+ ReferenceElements::general( geometry.type() );
+ const int numberOfFaces=refElement.size(1);
+
+ std::vector<Scalar> flux(numberOfFaces,0);
+
// run through all intersections with neighbors and boundary
IntersectionIterator
isEndIt = problem_.gridView().iend(*eIt);
@@ -136,19 +144,31 @@ public:
flux[isIndex] = isIt->geometry().volume()
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(isIndex));
}
-
Dune::FieldVector<Scalar, dim> refVelocity(0);
- for (int i = 0; i < dim; i++)
- refVelocity[i] = 0.5 * (flux[2*i + 1] - flux[2*i]);
-
- const typename Element::Geometry& geometry = eIt->geometry();
+ //2D simplices
+ if ( dim==2 && geometry.corners() == 3 ) {
+ refVelocity[0] = 1.0/3.0*(2.0*flux[1] - flux[2]);
+ refVelocity[1] = 1.0/3.0*(2.0*flux[0] - flux[2]);
+ }
+ //3D tetrahedra
+ else if ( dim==3 && geometry.corners() == 4 ){
+ DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
+ }
+ //2D and 3D cubes
+ else if ( refElement.type().isCube() ){
+ for (int i = 0; i < dim; i++)
+ refVelocity[i] = 0.5 * (flux[2*i + 1] - flux[2*i]);
+ }
+ //3D prism and pyramids
+ else {
+ DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
+ }
- typedef Dune::ReferenceElements<Scalar, dim> ReferenceElements;
- const Dune::FieldVector<Scalar, dim>& localPos
- = ReferenceElements::general(geometry.type()).position(0, 0);
+ const Dune::FieldVector<Scalar, dim>& localPos
+ = refElement.position(0, 0);
// get the transposed Jacobian of the element mapping
- const typename Element::Geometry::JacobianTransposed& jacobianT =
+ const typename Element::Geometry::JacobianTransposed& jacobianT =
geometry.jacobianTransposed(localPos);
// calculate the element velocity by the Piola transformation
diff --git a/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh b/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
index d59ee2803e68e448c3b482fcd6d65d2529ccb1ea..d45c28026a62d30c2bdb8d96afb66027b31580ad 100644
--- a/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
@@ -198,8 +198,16 @@ public:
CellData& cellData = problem_.variables().cellData(globalIdx);
- Dune::FieldVector < Scalar, 2 * dim > fluxW(0);
- Dune::FieldVector < Scalar, 2 * dim > fluxNw(0);
+ const typename Element::Geometry& geometry = eIt->geometry();
+ // get corresponding reference element
+ typedef Dune::ReferenceElements<Scalar, dim> ReferenceElements;
+ const Dune::ReferenceElement< Scalar , dim > & refElement =
+ ReferenceElements::general( geometry.type() );
+ const int numberOfFaces=refElement.size(1);
+
+ std::vector<Scalar> fluxW(numberOfFaces,0);
+ 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)
@@ -213,31 +221,61 @@ public:
}
Dune::FieldVector < Scalar, dim > refVelocity(0);
- for (int i = 0; i < dim; i++)
- refVelocity[i] = 0.5 * (fluxW[2*i + 1] - fluxW[2*i]);
-
+ //2D simplices
+ if ( dim==2 && geometry.corners() == 3 ) {
+ refVelocity[0] = 1.0/3.0*(2.0*fluxW[1] - fluxW[2]);
+ refVelocity[1] = 1.0/3.0*(2.0*fluxW[0] - fluxW[2]);
+ }
+ //3D tetrahedra
+ else if ( dim==3 && geometry.corners() == 4 ){
+ DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
+ }
+ //2D and 3D cubes
+ else if ( refElement.type().isCube() ){
+ for (int i = 0; i < dim; i++)
+ refVelocity[i] = 0.5 * (fluxW[2*i + 1] - fluxW[2*i]);
+ }
+ //3D prism and pyramids
+ else {
+ DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
+ }
const Dune::FieldVector<Scalar, dim> localPos =
- ReferenceElements::general(eIt->geometry().type()).position(0, 0);
+ refElement.position(0, 0);
// get the transposed Jacobian of the element mapping
const JacobianTransposed jacobianT =
- eIt->geometry().jacobianTransposed(localPos);
+ geometry.jacobianTransposed(localPos);
// calculate the element velocity by the Piola transformation
Dune::FieldVector < Scalar, dim > elementVelocity(0);
jacobianT.umtv(refVelocity, elementVelocity);
- elementVelocity /= eIt->geometry().integrationElement(localPos);
+ elementVelocity /= geometry.integrationElement(localPos);
velocity[globalIdx] = elementVelocity;
refVelocity = 0;
- for (int i = 0; i < dim; i++)
- refVelocity[i] = 0.5 * (fluxNw[2*i + 1] - fluxNw[2*i]);
-
+ //2D simplices
+ if ( dim==2 && geometry.corners() == 3 ) {
+ refVelocity[0] = 1.0/3.0*(2.0*fluxNw[1] - fluxNw[2]);
+ refVelocity[1] = 1.0/3.0*(2.0*fluxNw[0] - fluxNw[2]);
+ }
+ //3D tetrahedra
+ else if ( dim==3 && geometry.corners() == 4 ){
+ DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
+ }
+ //2D and 3D cubes
+ else if ( refElement.type().isCube() ){
+ for (int i = 0; i < dim; i++)
+ refVelocity[i] = 0.5 * (fluxNw[2*i + 1] - fluxNw[2*i]);
+ }
+ //3D prism and pyramids
+ else {
+ DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
+ }
// calculate the element velocity by the Piola transformation
elementVelocity = 0;
jacobianT.umtv(refVelocity, elementVelocity);
- elementVelocity /= eIt->geometry().integrationElement(localPos);
+ elementVelocity /= geometry.integrationElement(localPos);
velocitySecondPhase[globalIdx] = elementVelocity;
}
diff --git a/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2p.hh b/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2p.hh
index 58badc2a04c417d9533d953c10528669bccbe4e8..0fe3d33964737f4e2fb06d4730ccb4d0b2468b2c 100644
--- a/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2p.hh
+++ b/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2p.hh
@@ -280,8 +280,15 @@ public:
(*pressureSecond)[globalIdx] = cellData.pressure(wPhaseIdx);
}
- Dune::FieldVector < Scalar, 2 * dim > fluxW(0);
- Dune::FieldVector < Scalar, 2 * dim > fluxNw(0);
+ const typename Element::Geometry& geometry = eIt->geometry();
+ // get corresponding reference element
+ typedef Dune::ReferenceElements<Scalar, dim> ReferenceElements;
+ const Dune::ReferenceElement< Scalar , dim > & refElement =
+ ReferenceElements::general( geometry.type() );
+ const int numberOfFaces=refElement.size(1);
+
+ std::vector<Scalar> fluxW(numberOfFaces,0);
+ std::vector<Scalar> fluxNw(numberOfFaces,0);
// run through all intersections with neighbors and boundary
IntersectionIterator isEndIt = problem_.gridView().iend(*eIt);
@@ -295,29 +302,59 @@ public:
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
}
- DimVector refVelocity;
- for (int k = 0; k < dim; k++)
- refVelocity[k] = 0.5 * (fluxW[2*k+1] - fluxW[2*k]);
-
- const DimVector& localPos = ReferenceElements::general(eIt->geometry().type()).position(0, 0);
+ DimVector refVelocity(0.0);
+ //2D simplices
+ if ( dim==2 && geometry.corners() == 3 ) {
+ refVelocity[0] = 1.0/3.0*(2.0*fluxW[1] - fluxW[2]);
+ refVelocity[1] = 1.0/3.0*(2.0*fluxW[0] - fluxW[2]);
+ }
+ //3D tetrahedra
+ else if ( dim==3 && geometry.corners() == 4 ){
+ DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
+ }
+ //2D and 3D cubes
+ else if ( refElement.type().isCube() ){
+ for (int k = 0; k < dim; k++)
+ refVelocity[k] = 0.5 * (fluxW[2*k+1] - fluxW[2*k]);
+ }
+ //3D prism and pyramids
+ else {
+ DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
+ }
+ const DimVector& localPos = refElement.position(0, 0);
// get the transposed Jacobian of the element mapping
- const JacobianTransposed jacobianT = eIt->geometry().jacobianTransposed(localPos);
+ const JacobianTransposed jacobianT = geometry.jacobianTransposed(localPos);
// calculate the element velocity by the Piola transformation
DimVector elementVelocity(0);
jacobianT.umtv(refVelocity, elementVelocity);
- elementVelocity /= eIt->geometry().integrationElement(localPos);
+ elementVelocity /= geometry.integrationElement(localPos);
(*velocityWetting)[globalIdx] = elementVelocity;
- for (int k = 0; k < dim; k++)
- refVelocity[k] = 0.5 * (fluxNw[2*k+1] - fluxNw[2*k]);
-
+ //2D simplices
+ if ( dim==2 && geometry.corners() == 3 ) {
+ refVelocity[0] = 1.0/3.0*(2.0*fluxNw[1] - fluxNw[2]);
+ refVelocity[1] = 1.0/3.0*(2.0*fluxNw[0] - fluxNw[2]);
+ }
+ //3D tetrahedra
+ else if ( dim==3 && geometry.corners() == 4 ){
+ DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
+ }
+ //2D and 3D cubes
+ else if ( refElement.type().isCube() ){
+ for (int k = 0; k < dim; k++)
+ refVelocity[k] = 0.5 * (fluxNw[2*k+1] - fluxNw[2*k]);
+ }
+ //3D prism and pyramids
+ else {
+ DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
+ }
// calculate the element velocity by the Piola transformation
elementVelocity = 0;
jacobianT.umtv(refVelocity, elementVelocity);
- elementVelocity /= eIt->geometry().integrationElement(localPos);
+ elementVelocity /= geometry.integrationElement(localPos);
(*velocityNonwetting)[globalIdx] = elementVelocity;
}
diff --git a/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2padaptive.hh b/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2padaptive.hh
index 546ea533a9b7a890341915313fff69f4fb53277b..7397f651ff13c11353f5f310bfefe99eccd7d021 100644
--- a/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2padaptive.hh
+++ b/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2padaptive.hh
@@ -292,8 +292,15 @@ public:
(*pressureSecond)[globalIdx] = cellData.pressure(wPhaseIdx);
}
- Dune::FieldVector < Scalar, 2 * dim > fluxW(0);
- Dune::FieldVector < Scalar, 2 * dim > fluxNw(0);
+ const typename Element::Geometry& geometry = eIt->geometry();
+ // get corresponding reference element
+ typedef Dune::ReferenceElements<Scalar, dim> ReferenceElements;
+ const Dune::ReferenceElement< Scalar , dim > & refElement =
+ ReferenceElements::general( geometry.type() );
+ const int numberOfFaces=refElement.size(1);
+
+ std::vector<Scalar> fluxW(numberOfFaces,0);
+ std::vector<Scalar> fluxNw(numberOfFaces,0);
// run through all intersections with neighbors and boundary
IntersectionIterator isEndIt = problem_.gridView().iend(*eIt);
@@ -307,29 +314,59 @@ public:
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
}
- DimVector refVelocity;
- for (int k = 0; k < dim; k++)
- refVelocity[k] = 0.5 * (fluxW[2*k+1] - fluxW[2*k]);
-
- const DimVector& localPos = ReferenceElements::general(eIt->geometry().type()).position(0, 0);
+ Dune::FieldVector < Scalar, dim > refVelocity(0);
+ //2D simplices
+ if ( dim==2 && geometry.corners() == 3 ) {
+ refVelocity[0] = 1.0/3.0*(2.0*fluxW[1] - fluxW[2]);
+ refVelocity[1] = 1.0/3.0*(2.0*fluxW[0] - fluxW[2]);
+ }
+ //3D tetrahedra
+ else if ( dim==3 && geometry.corners() == 4 ){
+ DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
+ }
+ //2D and 3D cubes
+ else if ( refElement.type().isCube() ){
+ for (int k = 0; k < dim; k++)
+ refVelocity[k] = 0.5 * (fluxW[2*k+1] - fluxW[2*k]);
+ }
+ //3D prism and pyramids
+ else {
+ DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
+ }
+ const DimVector& localPos = refElement.position(0, 0);
// get the transposed Jacobian of the element mapping
- const JacobianTransposed jacobianT = eIt->geometry().jacobianTransposed(localPos);
+ const JacobianTransposed jacobianT = geometry.jacobianTransposed(localPos);
// calculate the element velocity by the Piola transformation
DimVector elementVelocity(0);
jacobianT.umtv(refVelocity, elementVelocity);
- elementVelocity /= eIt->geometry().integrationElement(localPos);
+ elementVelocity /= geometry.integrationElement(localPos);
(*velocityWetting)[globalIdx] = elementVelocity;
- for (int k = 0; k < dim; k++)
- refVelocity[k] = 0.5 * (fluxNw[2*k+1] - fluxNw[2*k]);
-
+ //2D simplices
+ if ( dim==2 && geometry.corners() == 3 ) {
+ refVelocity[0] = 1.0/3.0*(2.0*fluxNw[1] - fluxNw[2]);
+ refVelocity[1] = 1.0/3.0*(2.0*fluxNw[0] - fluxNw[2]);
+ }
+ //3D tetrahedra
+ else if ( dim==3 && geometry.corners() == 4 ){
+ DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
+ }
+ //2D and 3D cubes
+ else if ( refElement.type().isCube() ){
+ for (int k = 0; k < dim; k++)
+ refVelocity[k] = 0.5 * (fluxNw[2*k+1] - fluxNw[2*k]);
+ }
+ //3D prism and pyramids
+ else {
+ DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
+ }
// calculate the element velocity by the Piola transformation
elementVelocity = 0;
jacobianT.umtv(refVelocity, elementVelocity);
- elementVelocity /= eIt->geometry().integrationElement(localPos);
+ elementVelocity /= geometry.integrationElement(localPos);
(*velocityNonwetting)[globalIdx] = elementVelocity;
}
diff --git a/dumux/decoupled/2p/impes/gridadaptionindicator2plocalflux.hh b/dumux/decoupled/2p/impes/gridadaptionindicator2plocalflux.hh
index ed9555763093d2e20a272b1fbf0d93e04220a2d4..8dc4d7cf921af43cfed45798151404ceb0fd21f9 100644
--- a/dumux/decoupled/2p/impes/gridadaptionindicator2plocalflux.hh
+++ b/dumux/decoupled/2p/impes/gridadaptionindicator2plocalflux.hh
@@ -211,8 +211,16 @@ public:
break;
}
- Dune::FieldVector < Scalar, 2 * dim > flux(0);
- // Berechne Verfeinerungsindikator an allen Zellen
+ const typename Element::Geometry& geometry = eIt->geometry();
+ // get corresponding reference element
+ typedef Dune::ReferenceElements<Scalar, dim> ReferenceElements;
+ const Dune::ReferenceElement< Scalar , dim > & refElement =
+ ReferenceElements::general( geometry.type() );
+ const int numberOfFaces=refElement.size(1);
+
+ std::vector<Scalar> flux(numberOfFaces,0);
+
+ // Calculate refinement indicator on all cells
IntersectionIterator isItend = problem_.gridView().iend(*eIt);
for (IntersectionIterator isIt = problem_.gridView().ibegin(*eIt); isIt != isItend; ++isIt)
{
@@ -331,19 +339,34 @@ public:
if (useFluxInd_ || usePercentileFlux_)
{
DimVector refVelocity(0);
- for (int i = 0; i < dim; i++)
- refVelocity[i] = 0.5 * (flux[2*i + 1] - flux[2*i]);
+ //2D simplices
+ if ( dim==2 && geometry.corners() == 3 ) {
+ refVelocity[0] = 1.0/3.0*(2.0*flux[1] - flux[2]);
+ refVelocity[1] = 1.0/3.0*(2.0*flux[0] - flux[2]);
+ }
+ //3D tetrahedra
+ else if ( dim==3 && geometry.corners() == 4 ){
+ DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
+ }
+ //2D and 3D cubes
+ else if ( refElement.typ().isCube() ){
+ for (int i = 0; i < dim; i++)
+ refVelocity[i] = 0.5 * (flux[2*i + 1] - flux[2*i]);
+ }
+ //3D prism and pyramids
+ else {
+ DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
+ }
- const DimVector& localPos =
- ReferenceElements::general(eIt->geometry().type()).position(0, 0);
+ const DimVector& localPos = refElement.position(0, 0);
// get the transposed Jacobian of the element mapping
- const DimMatrix jacobianT = eIt->geometry().jacobianTransposed(localPos);
+ const DimMatrix jacobianT = geometry.jacobianTransposed(localPos);
// calculate the element velocity by the Piola transformation
DimVector elementVelocity(0);
jacobianT.umtv(refVelocity, elementVelocity);
- elementVelocity /= eIt->geometry().integrationElement(localPos);
+ elementVelocity /= geometry.integrationElement(localPos);
Scalar velNorm = elementVelocity.two_norm();
indicatorVectorFlux_[globalIdxI] = std::max(velNorm, indicatorVectorFlux_[globalIdxI]);