Fix an error in the DG skeleton flux term

The upwind numeric flux already contains the upwind conductivity factor

dune/dorie/solver/operator_DG.hh

@@ -380,7 +380,7 @@ public:
         }
 
         // compute upwind flux
-        RF numFlux_upwind = numFlux;
+        RF numFlux_upwind;
         RF satCond_upwind;
         if (upwinding == RichardsDGUpwinding::fullUpwind) {
           if (numFlux > 0) {
@@ -397,10 +397,27 @@ public:
         // diffusion term
         // consistency term
         // + penalty term
-        for (unsigned int i = 0; i<lfsv_s.size(); i++)
-          r_s.accumulate(lfsv_s,i,   relCond_upwind * numFlux_upwind * phi_s[i] * factor);
-        for (unsigned int i = 0; i<lfsv_n.size(); i++)
-          r_n.accumulate(lfsv_n,i, - relCond_upwind * numFlux_upwind * phi_n[i] * factor);
+        if (upwinding == RichardsDGUpwinding::none)
+        {
+          for (unsigned int i = 0; i < lfsv_s.size(); i++)
+            r_s.accumulate(lfsv_s, i,   numFlux * phi_s[i] * factor);
+          for (unsigned int i = 0; i < lfsv_n.size(); i++)
+            r_n.accumulate(lfsv_n, i, - numFlux * phi_n[i] * factor);
+        }
+        else if (upwinding == RichardsDGUpwinding::semiUpwind)
+        {
+          for (unsigned int i = 0; i<lfsv_s.size(); i++)
+            r_s.accumulate(lfsv_s,i,   relCond_upwind * numFlux * phi_s[i] * factor);
+          for (unsigned int i = 0; i<lfsv_n.size(); i++)
+            r_n.accumulate(lfsv_n,i, - relCond_upwind * numFlux * phi_n[i] * factor);
+        }
+        else if (upwinding == RichardsDGUpwinding::fullUpwind)
+        {
+          for (unsigned int i = 0; i < lfsv_s.size(); i++)
+            r_s.accumulate(lfsv_s, i,   numFlux_upwind * phi_s[i] * factor);
+          for (unsigned int i = 0; i < lfsv_n.size(); i++)
+            r_n.accumulate(lfsv_n, i, - numFlux_upwind * phi_n[i] * factor);
+        }
 
         // update residual (symmetry term)
         // (non-)symmetric IP term
@@ -426,7 +443,6 @@ public:
           for (unsigned int i = 0; i < lfsv_n.size(); i++)
             r_n.accumulate(lfsv_n, i, theta * (tgradphi_n[i] * normal) * relCond_upwind * satCond_upwind * jump * factor);
         }
-
       }
     }