From bfde81bdc3c1e31d65ca9148925c6763671a2921 Mon Sep 17 00:00:00 2001
From: Nicolas Schwenck <nicolas.schwenck@iws.uni-stuttgart.de>
Date: Wed, 29 Oct 2014 13:56:22 +0000
Subject: [PATCH] fixed wrong velocity output calculation for non-cube elements
in decoupled models (except MPFA which anyway only works on quadrilateral
grids)
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@13584 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
.../decoupled/1p/diffusion/fv/fvvelocity1p.hh | 40 +++++++++---
.../decoupled/2p/diffusion/fv/fvvelocity2p.hh | 62 ++++++++++++++----
.../2p/diffusion/mimetic/mimeticpressure2p.hh | 63 +++++++++++++++----
.../mimetic/mimeticpressure2padaptive.hh | 63 +++++++++++++++----
.../impes/gridadaptionindicator2plocalflux.hh | 39 +++++++++---
5 files changed, 211 insertions(+), 56 deletions(-)
diff --git a/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh b/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh
index a5679dca9f..3c9aa18881 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 d59ee2803e..d45c28026a 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 58badc2a04..0fe3d33964 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 546ea533a9..7397f651ff 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 ed95557630..8dc4d7cf92 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]);
--
GitLab