From 741bbd775b3f9669593356420d78f3a5bc1d1849 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Tue, 7 Mar 2017 21:14:45 +0100
Subject: [PATCH] [box] Simplify Darcys law using axpy operator

---
 dumux/discretization/box/darcyslaw.hh | 25 +++++++------------------
 1 file changed, 7 insertions(+), 18 deletions(-)

diff --git a/dumux/discretization/box/darcyslaw.hh b/dumux/discretization/box/darcyslaw.hh
index f4b1bc04b9..9b1790b131 100644
--- a/dumux/discretization/box/darcyslaw.hh
+++ b/dumux/discretization/box/darcyslaw.hh
@@ -113,23 +113,13 @@ public:
                 rho += volVars.density(phaseIdx)*shapeValues[scv.index()][0];
 
             // the global shape function gradient
-            DimVector gradI;
-            jacInvT.mv(shapeJacobian[scv.index()][0], gradI);
-
-            gradI *= volVars.pressure(phaseIdx);
-            gradP += gradI;
+            DimVector gradN;
+            jacInvT.mv(shapeJacobian[scv.index()][0], gradN);
+            gradP.axpy(volVars.pressure(phaseIdx), gradN);
         }
-        if (enableGravity)
-        {
-            // gravitational acceleration
-            DimVector g(problem.gravityAtPos(scvf.center()));
-
-            // turn gravity into a force
-            g *= rho;
 
-            // subtract from pressure gradient
-            gradP -= g;
-        }
+        if (enableGravity)
+            gradP.axpy(-rho, problem.gravityAtPos(scvf.center()));
 
         // apply the permeability and return the flux
         auto KGradP = applyPermeability_(K, gradP);
@@ -137,15 +127,14 @@ public:
     }
 
 private:
-    static DimVector applyPermeability_(const DimWorldMatrix& K, const DimVector& gradI)
+    inline static DimVector applyPermeability_(const DimWorldMatrix& K, const DimVector& gradI)
     {
         DimVector result(0.0);
         K.mv(gradI, result);
-
         return result;
     }
 
-    static DimVector applyPermeability_(const Scalar k, const DimVector& gradI)
+    inline static DimVector applyPermeability_(const Scalar k, const DimVector& gradI)
     {
         DimVector result(gradI);
         result *= k;
-- 
GitLab