diff --git a/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh b/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh
index 3c9aa18881eb52220f04166d4a2f69926c2e196e..bb8ab19b8acf24527ebd0c8b4efe7dfff7f61e70 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 d45c28026a62d30c2bdb8d96afb66027b31580ad..e881fa6ca8a17fb1ec7e00cea44e179dbf52f57c 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 0fe3d33964737f4e2fb06d4730ccb4d0b2468b2c..86d12253fa33be818f16b5098f6ad8c8f0150512 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 7397f651ff13c11353f5f310bfefe99eccd7d021..4cb5df8efffbc7ffbc0e27ccfbbcf094acf321bf 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 8dc4d7cf921af43cfed45798151404ceb0fd21f9..dbe8be91433e41917ad04b14049904927b5983eb 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 8ad17d525260109d16425f308e7481442e5af3a9..016c08de1f43a4375313a68128fb29243cf7f7b5 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 57cfa8a796192ff3726a59f2c066135b52476269..fb7e34bb5b1b23ad6f80841e41cbd0985c6059dd 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