From a762c728a96b0ed2077f59e68909476c0a5e798e Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Wed, 5 Nov 2014 08:20:24 +0000
Subject: [PATCH] [decoupled] generalize velocity output to handle also
simplices
So far, velocity output for the decoupled models has been restricted to
cubes. It now also works for simplices such as triangles and
tetrahedrons. For each element, it evaluates the Raviart-Thomas-0
interpolant of the fluxes in the barycenter.
Brought to attention by, discussed with and reviewed by Nicolas.
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@13610 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
.../decoupled/1p/diffusion/fv/fvvelocity1p.hh | 29 +++++----
.../decoupled/2p/diffusion/fv/fvvelocity2p.hh | 57 ++++++++++-------
.../2p/diffusion/mimetic/mimeticpressure2p.hh | 64 +++++++++++--------
.../mimetic/mimeticpressure2padaptive.hh | 64 +++++++++++--------
.../impes/gridadaptionindicator2plocalflux.hh | 28 ++++----
test/decoupled/1p/resultevaluation.hh | 27 ++++++--
test/decoupled/1p/resultevaluation3d.hh | 27 ++++++--
7 files changed, 180 insertions(+), 116 deletions(-)
diff --git a/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh b/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh
index 3c9aa18881..bb8ab19b8a 100644
--- a/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh
+++ b/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh
@@ -144,22 +144,27 @@ public:
flux[isIndex] = isIt->geometry().volume()
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(isIndex));
}
- Dune::FieldVector<Scalar, dim> refVelocity(0);
- //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");
+
+ // calculate velocity on reference element as the Raviart-Thomas-0
+ // interpolant of the fluxes
+ Dune::FieldVector<Scalar, dim> refVelocity;
+ // simplices
+ if (refElement.type().isSimplex()) {
+ for (int dimIdx = 0; dimIdx < dim; dimIdx++)
+ {
+ refVelocity[dimIdx] = -flux[dim - 1 - dimIdx];
+ for (int faceIdx = 0; faceIdx < dim + 1; faceIdx++)
+ {
+ refVelocity[dimIdx] += flux[faceIdx]/(dim + 1);
+ }
+ }
}
- //2D and 3D cubes
- else if ( refElement.type().isCube() ){
+ // 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
+ // 3D prism and pyramids
else {
DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
}
diff --git a/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh b/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
index d45c28026a..e881fa6ca8 100644
--- a/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
@@ -220,25 +220,30 @@ public:
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
}
- 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]);
+ // calculate velocity on reference element as the Raviart-Thomas-0
+ // interpolant of the fluxes
+ Dune::FieldVector<Scalar, dim> refVelocity;
+ // simplices
+ if (refElement.type().isSimplex()) {
+ for (int dimIdx = 0; dimIdx < dim; dimIdx++)
+ {
+ refVelocity[dimIdx] = -fluxW[dim - 1 - dimIdx];
+ for (int faceIdx = 0; faceIdx < dim + 1; faceIdx++)
+ {
+ refVelocity[dimIdx] += fluxW[faceIdx]/(dim + 1);
+ }
+ }
}
- //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() ){
+ // 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
+ // 3D prism and pyramids
else {
DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
}
+
const Dune::FieldVector<Scalar, dim> localPos =
refElement.position(0, 0);
@@ -253,25 +258,29 @@ public:
velocity[globalIdx] = elementVelocity;
- refVelocity = 0;
- //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]);
+ // calculate velocity on reference element as the Raviart-Thomas-0
+ // interpolant of the fluxes
+ // simplices
+ if (refElement.type().isSimplex()) {
+ for (int dimIdx = 0; dimIdx < dim; dimIdx++)
+ {
+ refVelocity[dimIdx] = -fluxNw[dim - 1 - dimIdx];
+ for (int faceIdx = 0; faceIdx < dim + 1; faceIdx++)
+ {
+ refVelocity[dimIdx] += fluxNw[faceIdx]/(dim + 1);
+ }
+ }
}
- //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() ){
+ // 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
+ // 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);
diff --git a/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2p.hh b/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2p.hh
index 0fe3d33964..86d12253fa 100644
--- a/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2p.hh
+++ b/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2p.hh
@@ -302,25 +302,30 @@ public:
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
}
- 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]);
+ // calculate velocity on reference element as the Raviart-Thomas-0
+ // interpolant of the fluxes
+ Dune::FieldVector<Scalar, dim> refVelocity;
+ // simplices
+ if (refElement.type().isSimplex()) {
+ for (int dimIdx = 0; dimIdx < dim; dimIdx++)
+ {
+ refVelocity[dimIdx] = -fluxW[dim - 1 - dimIdx];
+ for (int faceIdx = 0; faceIdx < dim + 1; faceIdx++)
+ {
+ refVelocity[dimIdx] += fluxW[faceIdx]/(dim + 1);
+ }
+ }
}
- //3D tetrahedra
- else if ( dim==3 && geometry.corners() == 4 ){
- DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
+ // cubes
+ else if (refElement.type().isCube()){
+ for (int i = 0; i < dim; i++)
+ refVelocity[i] = 0.5 * (fluxW[2*i + 1] - fluxW[2*i]);
}
- //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
+ // 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
@@ -333,24 +338,29 @@ public:
(*velocityWetting)[globalIdx] = elementVelocity;
- //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]);
+ // calculate velocity on reference element as the Raviart-Thomas-0
+ // interpolant of the fluxes
+ // simplices
+ if (refElement.type().isSimplex()) {
+ for (int dimIdx = 0; dimIdx < dim; dimIdx++)
+ {
+ refVelocity[dimIdx] = -fluxNw[dim - 1 - dimIdx];
+ for (int faceIdx = 0; faceIdx < dim + 1; faceIdx++)
+ {
+ refVelocity[dimIdx] += fluxNw[faceIdx]/(dim + 1);
+ }
+ }
}
- //3D tetrahedra
- else if ( dim==3 && geometry.corners() == 4 ){
- DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
+ // cubes
+ else if (refElement.type().isCube()){
+ for (int i = 0; i < dim; i++)
+ refVelocity[i] = 0.5 * (fluxNw[2*i + 1] - fluxNw[2*i]);
}
- //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
+ // 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);
diff --git a/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2padaptive.hh b/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2padaptive.hh
index 7397f651ff..4cb5df8eff 100644
--- a/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2padaptive.hh
+++ b/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2padaptive.hh
@@ -314,25 +314,30 @@ public:
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
}
- 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]);
+ // calculate velocity on reference element as the Raviart-Thomas-0
+ // interpolant of the fluxes
+ Dune::FieldVector<Scalar, dim> refVelocity;
+ // simplices
+ if (refElement.type().isSimplex()) {
+ for (int dimIdx = 0; dimIdx < dim; dimIdx++)
+ {
+ refVelocity[dimIdx] = -fluxW[dim - 1 - dimIdx];
+ for (int faceIdx = 0; faceIdx < dim + 1; faceIdx++)
+ {
+ refVelocity[dimIdx] += fluxW[faceIdx]/(dim + 1);
+ }
+ }
}
- //3D tetrahedra
- else if ( dim==3 && geometry.corners() == 4 ){
- DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
+ // cubes
+ else if (refElement.type().isCube()){
+ for (int i = 0; i < dim; i++)
+ refVelocity[i] = 0.5 * (fluxW[2*i + 1] - fluxW[2*i]);
}
- //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
+ // 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
@@ -345,24 +350,29 @@ public:
(*velocityWetting)[globalIdx] = elementVelocity;
- //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]);
+ // calculate velocity on reference element as the Raviart-Thomas-0
+ // interpolant of the fluxes
+ // simplices
+ if (refElement.type().isSimplex()) {
+ for (int dimIdx = 0; dimIdx < dim; dimIdx++)
+ {
+ refVelocity[dimIdx] = -fluxNw[dim - 1 - dimIdx];
+ for (int faceIdx = 0; faceIdx < dim + 1; faceIdx++)
+ {
+ refVelocity[dimIdx] += fluxNw[faceIdx]/(dim + 1);
+ }
+ }
}
- //3D tetrahedra
- else if ( dim==3 && geometry.corners() == 4 ){
- DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
+ // cubes
+ else if (refElement.type().isCube()){
+ for (int i = 0; i < dim; i++)
+ refVelocity[i] = 0.5 * (fluxNw[2*i + 1] - fluxNw[2*i]);
}
- //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
+ // 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);
diff --git a/dumux/decoupled/2p/impes/gridadaptionindicator2plocalflux.hh b/dumux/decoupled/2p/impes/gridadaptionindicator2plocalflux.hh
index 8dc4d7cf92..dbe8be9143 100644
--- a/dumux/decoupled/2p/impes/gridadaptionindicator2plocalflux.hh
+++ b/dumux/decoupled/2p/impes/gridadaptionindicator2plocalflux.hh
@@ -338,22 +338,26 @@ public:
if (useFluxInd_ || usePercentileFlux_)
{
- DimVector refVelocity(0);
- //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");
+ // calculate velocity on reference element as the Raviart-Thomas-0
+ // interpolant of the fluxes
+ Dune::FieldVector<Scalar, dim> refVelocity;
+ // simplices
+ if (refElement.type().isSimplex()) {
+ for (int dimIdx = 0; dimIdx < dim; dimIdx++)
+ {
+ refVelocity[dimIdx] = -flux[dim - 1 - dimIdx];
+ for (int faceIdx = 0; faceIdx < dim + 1; faceIdx++)
+ {
+ refVelocity[dimIdx] += flux[faceIdx]/(dim + 1);
+ }
+ }
}
- //2D and 3D cubes
- else if ( refElement.typ().isCube() ){
+ // 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
+ // 3D prism and pyramids
else {
DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
}
diff --git a/test/decoupled/1p/resultevaluation.hh b/test/decoupled/1p/resultevaluation.hh
index 8ad17d5252..016c08de1f 100644
--- a/test/decoupled/1p/resultevaluation.hh
+++ b/test/decoupled/1p/resultevaluation.hh
@@ -200,15 +200,28 @@ public:
}
}
- // calculate velocity on reference element
- Dune::FieldVector<Scalar,dim> refVelocity;
- if (geometry.corners() == 3) {
- refVelocity[0] = 1.0/3.0*(fluxVector[0] + fluxVector[2] - 2.0*fluxVector[1]);
- refVelocity[1] = 1.0/3.0*(fluxVector[0] + fluxVector[1] - 2.0*fluxVector[2]);
+ // calculate velocity on reference element as the Raviart-Thomas-0
+ // interpolant of the fluxes
+ Dune::FieldVector<Scalar, dim> refVelocity;
+ // simplices
+ if (geometry.type().isSimplex()) {
+ for (int dimIdx = 0; dimIdx < dim; dimIdx++)
+ {
+ refVelocity[dimIdx] = -fluxVector[dim - 1 - dimIdx];
+ for (int faceIdx = 0; faceIdx < dim + 1; faceIdx++)
+ {
+ refVelocity[dimIdx] += fluxVector[faceIdx]/(dim + 1);
+ }
+ }
+ }
+ // cubes
+ else if (geometry.type().isCube()){
+ for (int i = 0; i < dim; i++)
+ refVelocity[i] = 0.5 * (fluxVector[2*i + 1] - fluxVector[2*i]);
}
+ // 3D prism and pyramids
else {
- refVelocity[0] = 0.5*(fluxVector[1] - fluxVector[0]);
- refVelocity[1] = 0.5*(fluxVector[3] - fluxVector[2]);
+ DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
}
// get the transposed Jacobian of the element mapping
diff --git a/test/decoupled/1p/resultevaluation3d.hh b/test/decoupled/1p/resultevaluation3d.hh
index 57cfa8a796..fb7e34bb5b 100644
--- a/test/decoupled/1p/resultevaluation3d.hh
+++ b/test/decoupled/1p/resultevaluation3d.hh
@@ -303,15 +303,28 @@ public:
}
}
- // calculate velocity on reference element
- Dune::FieldVector<ct,dim> refVelocity;
- if (geometry.corners() == 8) {
- refVelocity[0] = 0.5*(fluxVector[1] - fluxVector[0]);
- refVelocity[1] = 0.5*(fluxVector[3] - fluxVector[2]);
- refVelocity[2] = 0.5*(fluxVector[5] - fluxVector[4]);
+ // calculate velocity on reference element as the Raviart-Thomas-0
+ // interpolant of the fluxes
+ Dune::FieldVector<double, dim> refVelocity;
+ // simplices
+ if (geometry.type().isSimplex()) {
+ for (int dimIdx = 0; dimIdx < dim; dimIdx++)
+ {
+ refVelocity[dimIdx] = -fluxVector[dim - 1 - dimIdx];
+ for (int faceIdx = 0; faceIdx < dim + 1; faceIdx++)
+ {
+ refVelocity[dimIdx] += fluxVector[faceIdx]/(dim + 1);
+ }
+ }
+ }
+ // cubes
+ else if (geometry.type().isCube()){
+ for (int i = 0; i < dim; i++)
+ refVelocity[i] = 0.5 * (fluxVector[2*i + 1] - fluxVector[2*i]);
}
+ // 3D prism and pyramids
else {
- DUNE_THROW(Dune::NotImplemented, "grid shape is not supported!");
+ DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
}
// get the transposed Jacobian of the element mapping
--
GitLab