From 1ba4ef46c9423a739af45f920477300ea7fb53d0 Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Fri, 17 Aug 2012 09:56:11 +0000
Subject: [PATCH] decoupled element velocity calculations: generalize to
 arbitrary dimension

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@8891 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
 dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh    |  4 ++--
 dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh    | 14 ++++----------
 .../fvmpfa/lmethod/fvmpfal2pfaboundvelocity2p.hh   |  8 ++++----
 .../lmethod/fvmpfal2pfaboundvelocity2padaptive.hh  |  8 ++++----
 .../fvmpfa/omethod/fvmpfao2pfaboundvelocity2p.hh   |  8 ++++----
 5 files changed, 18 insertions(+), 24 deletions(-)

diff --git a/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh b/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh
index 928f04308f..089b699afd 100644
--- a/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh
+++ b/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh
@@ -142,8 +142,8 @@ public:
             }
 
             Dune::FieldVector<Scalar, dim> refVelocity(0);
-            refVelocity[0] = 0.5 * (flux[1] - flux[0]);
-            refVelocity[1] = 0.5 * (flux[3] - flux[2]);
+            for (int i = 0; i < dim; i++)
+                refVelocity[i] = 0.5 * (flux[2*i + 1] - flux[2*i]);
 
             typedef Dune::GenericReferenceElements<Scalar, dim> ReferenceElements;
             const Dune::FieldVector<Scalar, dim>& localPos = ReferenceElements::general(eIt->geometry().type()).position(0,
diff --git a/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh b/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
index 5b061207e9..02b36c210e 100644
--- a/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
@@ -207,11 +207,8 @@ public:
             }
 
             Dune::FieldVector < Scalar, dim > refVelocity(0);
-            refVelocity[0] = 0.5 * (fluxW[1] - fluxW[0]);
-            if (dim > 1)
-            refVelocity[1] = 0.5 * (fluxW[3] - fluxW[2]);
-            if (dim == 3)
-            refVelocity[2] = 0.5 * (fluxW[5] - fluxW[4]);
+            for (int i = 0; i < dim; i++)
+                refVelocity[i] = 0.5 * (fluxW[2*i + 1] - fluxW[2*i]);
 
             const Dune::FieldVector<Scalar, dim>& localPos =
             ReferenceElementContainer::general(eIt->geometry().type()).position(0, 0);
@@ -229,11 +226,8 @@ public:
             velocity[globalIdx] = elementVelocity;
 
             refVelocity = 0;
-            refVelocity[0] = 0.5 * (fluxNW[1] - fluxNW[0]);
-            if (dim > 1)
-            refVelocity[1] = 0.5 * (fluxNW[3] - fluxNW[2]);
-            if (dim == 3)
-            refVelocity[2] = 0.5 * (fluxNW[5] - fluxNW[4]);
+            for (int i = 0; i < dim; i++)
+                refVelocity[i] = 0.5 * (fluxNW[2*i + 1] - fluxNW[2*i]);
 
             // calculate the element velocity by the Piola transformation
             elementVelocity = 0;
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundvelocity2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundvelocity2p.hh
index 242111aaee..0d31e740b0 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundvelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundvelocity2p.hh
@@ -185,8 +185,8 @@ public:
             }
 
             DimVector refVelocity(0);
-            refVelocity[0] = 0.5 * (fluxW[1] - fluxW[0]);
-            refVelocity[1] = 0.5 * (fluxW[3] - fluxW[2]);
+            for (int i = 0; i < dim; i++)
+                refVelocity[i] = 0.5 * (fluxW[2*i + 1] - fluxW[2*i]);
 
             const DimVector& localPos =
                     ReferenceElementContainer::general(eIt->geometry().type()).position(0, 0);
@@ -202,8 +202,8 @@ public:
             velocityWetting[globalIdx] = elementVelocity;
 
             refVelocity = 0;
-            refVelocity[0] = 0.5 * (fluxNW[1] - fluxNW[0]);
-            refVelocity[1] = 0.5 * (fluxNW[3] - fluxNW[2]);
+            for (int i = 0; i < dim; i++)
+                refVelocity[i] = 0.5 * (fluxNW[2*i + 1] - fluxNW[2*i]);
 
             // calculate the element velocity by the Piola transformation
             elementVelocity = 0;
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundvelocity2padaptive.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundvelocity2padaptive.hh
index 8c82617881..387ed571ec 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundvelocity2padaptive.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundvelocity2padaptive.hh
@@ -191,8 +191,8 @@ public:
              }
 
              DimVector refVelocity(0);
-             refVelocity[0] = 0.5 * (fluxW[1] - fluxW[0]);
-             refVelocity[1] = 0.5 * (fluxW[3] - fluxW[2]);
+            for (int i = 0; i < dim; i++)
+                refVelocity[i] = 0.5 * (fluxW[2*i + 1] - fluxW[2*i]);
 
              const DimVector& localPos =
                      ReferenceElementContainer::general(eIt->geometry().type()).position(0, 0);
@@ -208,8 +208,8 @@ public:
              velocityWetting[globalIdx] = elementVelocity;
 
              refVelocity = 0;
-             refVelocity[0] = 0.5 * (fluxNW[1] - fluxNW[0]);
-             refVelocity[1] = 0.5 * (fluxNW[3] - fluxNW[2]);
+            for (int i = 0; i < dim; i++)
+                refVelocity[i] = 0.5 * (fluxNW[2*i + 1] - fluxNW[2*i]);
 
              // calculate the element velocity by the Piola transformation
              elementVelocity = 0;
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2pfaboundvelocity2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2pfaboundvelocity2p.hh
index fd14abb154..787fe0eb46 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2pfaboundvelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2pfaboundvelocity2p.hh
@@ -185,8 +185,8 @@ public:
             }
 
             DimVector refVelocity(0);
-            refVelocity[0] = 0.5 * (fluxW[1] - fluxW[0]);
-            refVelocity[1] = 0.5 * (fluxW[3] - fluxW[2]);
+            for (int i = 0; i < dim; i++)
+                refVelocity[i] = 0.5 * (fluxW[2*i + 1] - fluxW[2*i]);
 
             const DimVector& localPos =
                     ReferenceElementContainer::general(eIt->geometry().type()).position(0, 0);
@@ -202,8 +202,8 @@ public:
             velocityWetting[globalIdx] = elementVelocity;
 
             refVelocity = 0;
-            refVelocity[0] = 0.5 * (fluxNW[1] - fluxNW[0]);
-            refVelocity[1] = 0.5 * (fluxNW[3] - fluxNW[2]);
+            for (int i = 0; i < dim; i++)
+                refVelocity[i] = 0.5 * (fluxNW[2*i + 1] - fluxNW[2*i]);
 
             // calculate the element velocity by the Piola transformation
             elementVelocity = 0;
-- 
GitLab