Merge branch 'feature/fluxcomputation-on-networkgrids' into 'next'

Feature/fluxcomputation on networkgrids

Todo: special treatment of branching points in ...

* [x] ... tpfa ficks law
* [x] ... tpfa fouriers law
* [x] Replace mobilityUpwindWeight / massUpwindWeight by upwindWeight

See merge request !304

dumux/discretization/cellcentered/mpfa/fickslaw.hh

@@ -77,27 +77,27 @@ public:
         const auto& tij = fluxVarsCache.diffusionTij(phaseIdx, compIdx);
 
         // get the scaling factor for the effective diffusive fluxes
-        auto factor = calculateEffectiveDiffusivityFactor(elemVolVars, scvf, phaseIdx);
+        auto effFactor = calculateEffectiveDiffusivityFactor(elemVolVars, scvf, phaseIdx);
 
         // if factor is zero, the flux will end up zero anyway
-        if (factor == 0.0)
+        if (effFactor == 0.0)
             return 0.0;
 
         // lambda functions depending on if we use mole or mass fractions
-        auto xDensity = [useMoles, phaseIdx] (const VolumeVariables& volVars)
-        { return useMoles ? volVars.molarDensity(phaseIdx) : volVars.density(phaseIdx); };
-
-        auto xFraction = [useMoles, phaseIdx, compIdx] (const VolumeVariables& volVars)
+        auto getX = [useMoles, phaseIdx, compIdx] (const VolumeVariables& volVars)
         { return useMoles ? volVars.moleFraction(phaseIdx, compIdx) : volVars.massFraction(phaseIdx, compIdx); };
 
+        auto getRho = [useMoles, phaseIdx] (const VolumeVariables& volVars)
+        { return useMoles ? volVars.molarDensity(phaseIdx) : volVars.density(phaseIdx); };
+
         // calculate Tij*xj
         Scalar flux(0.0);
         unsigned int localIdx = 0;
         for (const auto volVarIdx : volVarsStencil)
-            flux += tij[localIdx++]*xFraction(elemVolVars[volVarIdx]);
+            flux += tij[localIdx++]*getX(elemVolVars[volVarIdx]);
 
         // return effective mass flux
-        return flux*interpolateDensity(elemVolVars, scvf, xDensity)*factor;
+        return flux*interpolateDensity(elemVolVars, scvf, getRho)*effFactor;
     }
 
     static Stencil stencil(const Problem& problem,
@@ -115,21 +115,21 @@ public:
     }
 
 private:
-    template<typename DensityFunction>
+    template<typename GetRhoFunction>
     static Scalar interpolateDensity(const ElementVolumeVariables& elemVolVars,
                                      const SubControlVolumeFace& scvf,
-                                     const DensityFunction& rhoFunction)
+                                     const GetRhoFunction& getRho)
     {
         // use arithmetic mean of the densities around the scvf
         if (!scvf.boundary())
         {
-            Scalar rho = rhoFunction(elemVolVars[scvf.insideScvIdx()]);
+            Scalar rho = getRho(elemVolVars[scvf.insideScvIdx()]);
             for (auto outsideIdx : scvf.outsideScvIndices())
-                rho += rhoFunction(elemVolVars[outsideIdx]);
+                rho += getRho(elemVolVars[outsideIdx]);
             return rho/(scvf.outsideScvIndices().size()+1);
         }
         else
-            return rhoFunction(elemVolVars[scvf.outsideScvIdx()]);
+            return getRho(elemVolVars[scvf.outsideScvIdx()]);
     }
 
     //! Here we want to calculate the factors with which the diffusion coefficient has to be

dumux/discretization/cellcentered/tpfa/darcyslaw.hh

@@ -78,28 +78,28 @@ public:
                        const Element& element,
                        const FVElementGeometry& fvGeometry,
                        const ElementVolumeVariables& elemVolVars,
-                       const SubControlVolumeFace& scvFace,
+                       const SubControlVolumeFace& scvf,
                        int phaseIdx,
                        const ElementFluxVarsCache& elemFluxVarsCache)
     {
-        const auto& fluxVarsCache = elemFluxVarsCache[scvFace];
+        const auto& fluxVarsCache = elemFluxVarsCache[scvf];
 
         // Get the inside and outside volume variables
-        const auto& insideScv = fvGeometry.scv(scvFace.insideScvIdx());
+        const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
         const auto& insideVolVars = elemVolVars[insideScv];
-        const auto& outsideVolVars = elemVolVars[scvFace.outsideScvIdx()];
+        const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
 
         auto hInside = insideVolVars.pressure(phaseIdx);
-        auto hOutside = scvFace.numOutsideScvs() <= 1 ? outsideVolVars.pressure(phaseIdx)
+        auto hOutside = scvf.numOutsideScvs() <= 1 ? outsideVolVars.pressure(phaseIdx)
                         : branchingFacetPressure_(phaseIdx, problem, element, fvGeometry,
-                                                  elemVolVars, scvFace,
+                                                  elemVolVars, scvf,
                                                   elemFluxVarsCache, hInside);
 
         if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity))
         {
             // do averaging for the density over all neighboring elements
-            const auto rho = scvFace.numOutsideScvs() <= 1 ? (insideVolVars.density(phaseIdx) + outsideVolVars.density(phaseIdx))*0.5
-                             : branchingFacetDensity_(phaseIdx, elemVolVars, scvFace, insideVolVars.density(phaseIdx));
+            const auto rho = scvf.numOutsideScvs() <= 1 ? (insideVolVars.density(phaseIdx) + outsideVolVars.density(phaseIdx))*0.5
+                             : branchingFacetDensity_(phaseIdx, elemVolVars, scvf, insideVolVars.density(phaseIdx));
 
             // ask for the gravitational acceleration in the inside neighbor
             const auto xInside = insideScv.center();
@@ -108,15 +108,15 @@ public:
             hInside -= rho*(gInside*xInside);
 
             // and the outside neighbor
-            if (scvFace.boundary() || scvFace.numOutsideScvs() > 1)
+            if (scvf.boundary() || scvf.numOutsideScvs() > 1)
             {
-                const auto xOutside = scvFace.center();
+                const auto xOutside = scvf.ipGlobal();
                 const auto gOutside = problem.gravityAtPos(xOutside);
                 hOutside -= rho*(gOutside*xOutside);
             }
             else
             {
-                const auto outsideScvIdx = scvFace.outsideScvIdx();
+                const auto outsideScvIdx = scvf.outsideScvIdx();
                 // as we assemble fluxes from the neighbor to our element the outside index
                 // refers to the scv of our element, so we use the scv method
                 const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
@@ -132,19 +132,19 @@ public:
     static Stencil stencil(const Problem& problem,
                            const Element& element,
                            const FVElementGeometry& fvGeometry,
-                           const SubControlVolumeFace& scvFace)
+                           const SubControlVolumeFace& scvf)
     {
-        if (scvFace.boundary())
-            return Stencil({scvFace.insideScvIdx()});
-        else if (scvFace.numOutsideScvs() > 1)
+        if (scvf.boundary())
+            return Stencil({scvf.insideScvIdx()});
+        else if (scvf.numOutsideScvs() > 1)
         {
-            Stencil stencil({scvFace.insideScvIdx()});
-            for (unsigned int i = 0; i < scvFace.numOutsideScvs(); ++i)
-                stencil.push_back(scvFace.outsideScvIdx(i));
+            Stencil stencil({scvf.insideScvIdx()});
+            for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i)
+                stencil.push_back(scvf.outsideScvIdx(i));
             return stencil;
         }
         else
-            return Stencil({scvFace.insideScvIdx(), scvFace.outsideScvIdx()});
+            return Stencil({scvf.insideScvIdx(), scvf.outsideScvIdx()});
     }
 
     // The flux variables cache has to be bound to an element prior to flux calculations
@@ -153,25 +153,25 @@ public:
                                               const Element& element,
                                               const FVElementGeometry& fvGeometry,
                                               const ElementVolumeVariables& elemVolVars,
-                                              const SubControlVolumeFace& scvFace)
+                                              const SubControlVolumeFace& scvf)
     {
         Scalar tij;
 
-        const auto insideScvIdx = scvFace.insideScvIdx();
+        const auto insideScvIdx = scvf.insideScvIdx();
         const auto& insideScv = fvGeometry.scv(insideScvIdx);
         const auto& insideVolVars = elemVolVars[insideScvIdx];
         const auto insideK = problem.spatialParams().intrinsicPermeability(insideScv, insideVolVars);
-        Scalar ti = calculateOmega_(problem, scvFace, insideK, element, insideScv);
+        Scalar ti = calculateOmega_(problem, scvf, insideK, element, insideScv);
 
         // for the boundary (dirichlet) or at branching points we only need ti
-        if (scvFace.boundary() || scvFace.numOutsideScvs() > 1)
+        if (scvf.boundary() || scvf.numOutsideScvs() > 1)
         {
-            tij = scvFace.area()*ti;
+            tij = scvf.area()*ti;
         }
         // otherwise we compute a tpfa harmonic mean
         else
         {
-            const auto outsideScvIdx = scvFace.outsideScvIdx();
+            const auto outsideScvIdx = scvf.outsideScvIdx();
             // as we assemble fluxes from the neighbor to our element the outside index
             // refers to the scv of our element, so we use the scv method
             const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
@@ -181,15 +181,15 @@ public:
             Scalar tj;
             if (dim == dimWorld)
                 // assume the normal vector from outside is anti parallel so we save flipping a vector
-                tj = -1.0*calculateOmega_(problem, scvFace, outsideK, outsideElement, outsideScv);
+                tj = -1.0*calculateOmega_(problem, scvf, outsideK, outsideElement, outsideScv);
             else
-                tj = calculateOmega_(problem, fvGeometry.flipScvf(scvFace.index()), outsideK, outsideElement, outsideScv);
+                tj = calculateOmega_(problem, fvGeometry.flipScvf(scvf.index()), outsideK, outsideElement, outsideScv);
 
             // check for division by zero!
             if (ti*tj <= 0.0)
                 tij = 0;
             else
-                tij = scvFace.area()*(ti * tj)/(ti + tj);
+                tij = scvf.area()*(ti * tj)/(ti + tj);
         }
 
         return tij;
@@ -198,15 +198,15 @@ public:
 private:
     //! compute the transmissibility ti, overload for tensor permeabilites
     static Scalar calculateOmega_(const Problem& problem,
-                                  const SubControlVolumeFace& scvFace,
+                                  const SubControlVolumeFace& scvf,
                                   const DimWorldMatrix &K,
                                   const Element& element,
                                   const SubControlVolume &scv)
     {
         GlobalPosition Knormal;
-        K.mv(scvFace.unitOuterNormal(), Knormal);
+        K.mv(scvf.unitOuterNormal(), Knormal);
 
-        auto distanceVector = scvFace.center();
+        auto distanceVector = scvf.ipGlobal();
         distanceVector -= scv.center();
         distanceVector /= distanceVector.two_norm2();
 
@@ -218,16 +218,16 @@ private:
 
     //! compute the transmissibility ti, overload for scalar permeabilites
     static Scalar calculateOmega_(const Problem& problem,
-                                  const SubControlVolumeFace& scvFace,
+                                  const SubControlVolumeFace& scvf,
                                   const Scalar K,
                                   const Element& element,
                                   const SubControlVolume &scv)
     {
-        auto distanceVector = scvFace.center();
+        auto distanceVector = scvf.ipGlobal();
         distanceVector -= scv.center();
         distanceVector /= distanceVector.two_norm2();
 
-        Scalar omega = K * (distanceVector * scvFace.unitOuterNormal());
+        Scalar omega = K * (distanceVector * scvf.unitOuterNormal());
         omega *= problem.boxExtrusionFactor(element, scv);
 
         return omega;
@@ -239,18 +239,18 @@ private:
                                           const Element& element,
                                           const FVElementGeometry& fvGeometry,
                                           const ElementVolumeVariables& elemVolVars,
-                                          const SubControlVolumeFace& scvFace,
+                                          const SubControlVolumeFace& scvf,
                                           const ElementFluxVarsCache& elemFluxVarsCache,
                                           Scalar insideP)
     {
-        const auto& insideFluxVarsCache = elemFluxVarsCache[scvFace];
+        const auto& insideFluxVarsCache = elemFluxVarsCache[scvf];
         Scalar sumTi(insideFluxVarsCache.tij());
         Scalar sumPTi(insideFluxVarsCache.tij()*insideP);
 
-        for (unsigned int i = 0; i < scvFace.numOutsideScvs(); ++i)
+        for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i)
         {
-            const auto outsideScvIdx = scvFace.outsideScvIdx(i);
-            const auto& flippedScvf = fvGeometry.flipScvf(scvFace.index(), i);
+            const auto outsideScvIdx = scvf.outsideScvIdx(i);
+            const auto& flippedScvf = fvGeometry.flipScvf(scvf.index(), i);
             const auto& outsideVolVars = elemVolVars[outsideScvIdx];
             const auto& outsideFluxVarsCache = elemFluxVarsCache[flippedScvf];
             sumTi += outsideFluxVarsCache.tij();
@@ -262,17 +262,17 @@ private:
     //! compute the density at branching facets for network grids as arithmetic mean
     static Scalar branchingFacetDensity_(int phaseIdx,
                                          const ElementVolumeVariables& elemVolVars,
-                                         const SubControlVolumeFace& scvFace,
+                                         const SubControlVolumeFace& scvf,
                                          Scalar insideRho)
     {
         Scalar rho(insideRho);
-        for (unsigned int i = 0; i < scvFace.numOutsideScvs(); ++i)
+        for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i)
         {
-            const auto outsideScvIdx = scvFace.outsideScvIdx(i);
+            const auto outsideScvIdx = scvf.outsideScvIdx(i);
             const auto& outsideVolVars = elemVolVars[outsideScvIdx];
             rho += outsideVolVars.density(phaseIdx);
         }
-        return rho/(scvFace.numOutsideScvs()+1);
+        return rho/(scvf.numOutsideScvs()+1);
     }
 };
 

dumux/discretization/cellcentered/tpfa/fickslaw.hh

@@ -39,7 +39,6 @@ namespace Properties
 {
 // forward declaration of properties
 NEW_PROP_TAG(NumPhases);
-NEW_PROP_TAG(FluidSystem);
 NEW_PROP_TAG(EffectiveDiffusivityModel);
 }
 
@@ -52,7 +51,6 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCTpfa >
 {
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
-    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     using EffDiffModel = typename GET_PROP_TYPE(TypeTag, EffectiveDiffusivityModel);
     using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
     using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
@@ -60,13 +58,14 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCTpfa >
     using IndexType = typename GridView::IndexSet::IndexType;
     using Stencil = typename std::vector<IndexType>;
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
     using Element = typename GridView::template Codim<0>::Entity;
     using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
 
-    enum { dim = GridView::dimension} ;
-    enum { dimWorld = GridView::dimensionworld} ;
-    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)} ;
+    static const int dim = GridView::dimension;
+    static const int dimWorld = GridView::dimensionworld;
+    static const int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
 
     using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
     using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
@@ -79,91 +78,141 @@ public:
                        const Element& element,
                        const FVElementGeometry& fvGeometry,
                        const ElementVolumeVariables& elemVolVars,
-                       const SubControlVolumeFace& scvFace,
+                       const SubControlVolumeFace& scvf,
                        int phaseIdx, int compIdx,
                        const ElementFluxVariablesCache& elemFluxVarsCache,
                        bool useMoles = true)
     {
         // diffusion tensors are always solution dependent
-        Scalar tij = calculateTransmissibility_(problem, element, fvGeometry, elemVolVars, scvFace, phaseIdx, compIdx);
+        Scalar tij = calculateTransmissibility_(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, compIdx);
 
-        // Get the inside volume variables
-        const auto& insideScv = fvGeometry.scv(scvFace.insideScvIdx());
-        const auto& insideVolVars = elemVolVars[insideScv];
+        // get inside/outside volume variables
+        const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
+        const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
 
-        // and the outside volume variables
-        const auto& outsideVolVars = elemVolVars[scvFace.outsideScvIdx()];
+        // lambdas to get mole/mass fractions & densities
+        auto getX = [useMoles, phaseIdx, compIdx] (const VolumeVariables& volVars)
+        { return useMoles ? volVars.moleFraction(phaseIdx, compIdx) : volVars.massFraction(phaseIdx, compIdx); };
 
-        // compute the diffusive flux using mole fractions
-        if (useMoles)
-        {
-            const auto xInside = insideVolVars.moleFraction(phaseIdx, compIdx);
-            const auto xOutside = outsideVolVars.moleFraction(phaseIdx, compIdx);
-            const auto rho = 0.5*(insideVolVars.molarDensity(phaseIdx) + outsideVolVars.molarDensity(phaseIdx));
+        auto getRho = [useMoles, phaseIdx](const VolumeVariables& volVars)
+        { return useMoles? volVars.molarDensity(phaseIdx) : volVars.density(phaseIdx); };
 
-            return rho*tij*(xInside - xOutside);
-        }
-        // compute the diffusive flux using mass fractions
-        else
-        {
-            const auto xInside = insideVolVars.massFraction(phaseIdx, compIdx);
-            const auto xOutside = outsideVolVars.massFraction(phaseIdx, compIdx);
-            const auto rho = 0.5*(insideVolVars.density(phaseIdx) + outsideVolVars.density(phaseIdx));
+        // interpolate density
+        const auto rho = scvf.numOutsideScvs() == 1 ? 0.5*(getRho(insideVolVars)+ getRho(outsideVolVars))
+                         : branchingFacetDensity_(elemVolVars, scvf, getRho, getRho(insideVolVars));
 
-            return rho*tij*(xInside - xOutside);
-        }
+        // the inside and outside mole/mass fractions
+        auto xInside = getX(insideVolVars);
+        auto xOutside = scvf.numOutsideScvs() == 1 ? getX(outsideVolVars)
+                        : branchingFacetX_(problem, element, fvGeometry, elemVolVars, scvf, getX, xInside, tij, phaseIdx, compIdx);
+
+        return rho*tij*(xInside - xOutside);
     }
 
     static Stencil stencil(const Problem& problem,
                            const Element& element,
                            const FVElementGeometry& fvGeometry,
-                           const SubControlVolumeFace& scvFace)
+                           const SubControlVolumeFace& scvf)
     {
-        if (!scvFace.boundary())
-            return Stencil({scvFace.insideScvIdx(), scvFace.outsideScvIdx()});
+        if (!scvf.boundary())
+            return Stencil({scvf.insideScvIdx(), scvf.outsideScvIdx()});
         else
-            return Stencil({scvFace.insideScvIdx()});
+            return Stencil({scvf.insideScvIdx()});
     }
 
 private:
 
+    //! compute the mole/mass fraction at branching facets for network grids
+    template<typename GetXFunction>
+    static Scalar branchingFacetX_(const Problem& problem,
+                                   const Element& element,
+                                   const FVElementGeometry& fvGeometry,
+                                   const ElementVolumeVariables& elemVolVars,
+                                   const SubControlVolumeFace& scvf,
+                                   const GetXFunction& getX,
+                                   Scalar insideX, Scalar insideTi,
+                                   int phaseIdx, int compIdx)
+    {
+        Scalar sumTi(insideTi);
+        Scalar sumXTi(insideTi*insideX);
+
+        for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i)
+        {
+            const auto outsideScvIdx = scvf.outsideScvIdx(i);
+            const auto& outsideVolVars = elemVolVars[outsideScvIdx];
+            const auto outsideElement = fvGeometry.globalFvGeometry().element(outsideScvIdx);
+            const auto& flippedScvf = fvGeometry.flipScvf(scvf.index(), i);
+
+            auto outsideTi = calculateTransmissibility_(problem, outsideElement, fvGeometry, elemVolVars, flippedScvf, phaseIdx, compIdx);
+            sumTi += outsideTi;
+            sumXTi += outsideTi*getX(outsideVolVars);
+        }
+        return sumXTi/sumTi;
+    }
+
+    //! compute the density at branching facets for network grids as arithmetic mean
+    template<typename GetRhoFunction>
+    static Scalar branchingFacetDensity_(const ElementVolumeVariables& elemVolVars,
+                                         const SubControlVolumeFace& scvf,
+                                         const GetRhoFunction& getRho,
+                                         Scalar insideRho)
+    {
+        Scalar rho(insideRho);
+        for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i)
+        {
+            const auto outsideScvIdx = scvf.outsideScvIdx(i);
+            const auto& outsideVolVars = elemVolVars[outsideScvIdx];
+            rho += getRho(outsideVolVars);
+        }
+        return rho/(scvf.numOutsideScvs()+1);
+    }
+
 
     static Scalar calculateTransmissibility_(const Problem& problem,
                                              const Element& element,
                                              const FVElementGeometry& fvGeometry,
                                              const ElementVolumeVariables& elemVolVars,
-                                             const SubControlVolumeFace& scvFace,
-                                             const int phaseIdx, const int compIdx)
+                                             const SubControlVolumeFace& scvf,
+                                             int phaseIdx, int compIdx)
     {
         Scalar tij;
 
-        const auto insideScvIdx = scvFace.insideScvIdx();
+        const auto insideScvIdx = scvf.insideScvIdx();
         const auto& insideScv = fvGeometry.scv(insideScvIdx);
         const auto& insideVolVars = elemVolVars[insideScvIdx];
 
         auto insideD = insideVolVars.diffusionCoefficient(phaseIdx, compIdx);
         insideD = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), insideVolVars.saturation(phaseIdx), insideD);
-        Scalar ti = calculateOmega_(problem, element, scvFace, insideD, insideScv);
+        Scalar ti = calculateOmega_(problem, element, scvf, insideD, insideScv);
 
-        if (!scvFace.boundary())
+        // for the boundary (dirichlet) or at branching points we only need ti
+        if (scvf.boundary() || scvf.numOutsideScvs() > 1)
+        {
+            tij = scvf.area()*ti;
+        }
+        // otherwise we compute a tpfa harmonic mean
+        else
         {
-            const auto outsideScvIdx = scvFace.outsideScvIdx();
+            const auto outsideScvIdx = scvf.outsideScvIdx();
             const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
             const auto& outsideVolVars = elemVolVars[outsideScvIdx];
+            const auto outsideElement = fvGeometry.globalFvGeometry().element(outsideScvIdx);
 
             auto outsideD = outsideVolVars.diffusionCoefficient(phaseIdx, compIdx);
             outsideD = EffDiffModel::effectiveDiffusivity(outsideVolVars.porosity(), outsideVolVars.saturation(phaseIdx), outsideD);
-            Scalar tj = -1.0*calculateOmega_(problem, element, scvFace, outsideD, outsideScv);
+
+            Scalar tj;
+            if (dim == dimWorld)
+                // assume the normal vector from outside is anti parallel so we save flipping a vector
+                tj = -1.0*calculateOmega_(problem, outsideElement, scvf, outsideD, outsideScv);
+            else
+                tj = calculateOmega_(problem, outsideElement, fvGeometry.flipScvf(scvf.index()), outsideD, outsideScv);
 
             // check if we are dividing by zero!
             if (ti*tj <= 0.0)
                 tij = 0;
             else
-                tij = scvFace.area()*(ti * tj)/(ti + tj);
-        }
-        else
-        {
-            tij = scvFace.area()*ti;
+                tij = scvf.area()*(ti * tj)/(ti + tj);
         }
 
         return tij;
@@ -171,14 +220,14 @@ private:
 
     static Scalar calculateOmega_(const Problem& problem,
                                   const Element& element,
-                                  const SubControlVolumeFace& scvFace,
+                                  const SubControlVolumeFace& scvf,
                                   const DimWorldMatrix &D,
                                   const SubControlVolume &scv)
     {
         GlobalPosition Dnormal;
-        D.mv(scvFace.unitOuterNormal(), Dnormal);
+        D.mv(scvf.unitOuterNormal(), Dnormal);
 
-        auto distanceVector = scvFace.center();
+        auto distanceVector = scvf.ipGlobal();
         distanceVector -= scv.center();
         distanceVector /= distanceVector.two_norm2();
 
@@ -190,15 +239,15 @@ private:
 
     static Scalar calculateOmega_(const Problem& problem,
                                   const Element& element,
-                                  const SubControlVolumeFace& scvFace,
+                                  const SubControlVolumeFace& scvf,
                                   Scalar D,
                                   const SubControlVolume &scv)
     {
-        auto distanceVector = scvFace.center();
+        auto distanceVector = scvf.ipGlobal();
         distanceVector -= scv.center();
         distanceVector /= distanceVector.two_norm2();
 
-        Scalar omega = D * (distanceVector * scvFace.unitOuterNormal());
+        Scalar omega = D * (distanceVector * scvf.unitOuterNormal());
         omega *= problem.boxExtrusionFactor(element, scv);
 
         return omega;

dumux/discretization/cellcentered/tpfa/fourierslaw.hh

@@ -37,8 +37,6 @@ namespace Dumux
 namespace Properties
 {
 // forward declaration of properties
-NEW_PROP_TAG(NumPhases);
-NEW_PROP_TAG(FluidSystem);
 NEW_PROP_TAG(ThermalConductivityModel);
 }
 
@@ -51,7 +49,6 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCTpfa>
 {
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
-    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
     using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
     using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
@@ -62,9 +59,8 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCTpfa>
     using Element = typename GridView::template Codim<0>::Entity;
     using ElementFluxVarsCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
 
-    enum { dim = GridView::dimension} ;
-    enum { dimWorld = GridView::dimensionworld} ;
-    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)} ;
+    static const int dim = GridView::dimension;
+    static const int dimWorld = GridView::dimensionworld;
 
     using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
     using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
@@ -79,22 +75,16 @@ public:
                        const Element& element,
                        const FVElementGeometry& fvGeometry,
                        const ElementVolumeVariables& elemVolVars,
-                       const SubControlVolumeFace& scvFace,
+                       const SubControlVolumeFace& scvf,
                        const ElementFluxVarsCache& elemFluxVarsCache)
     {
         // heat conductivities are always solution dependent (?)
-        Scalar tij = calculateTransmissibility_(problem, element, fvGeometry, elemVolVars, scvFace);
+        Scalar tij = calculateTransmissibility_(problem, element, fvGeometry, elemVolVars, scvf);