[3p2cni] Make fluxvariables default constructable

lecture/mm/heavyoil/3p2cni/fluxvariables.hh

@@ -46,6 +46,7 @@ namespace Dumux
 template <class TypeTag>
 class ThreePTwoCNIFluxVariables : public GET_PROP_TYPE(TypeTag, BaseFluxVariables)
 {
+    friend typename GET_PROP_TYPE(TypeTag, BaseFluxVariables); // be friends with base class
     typedef typename GET_PROP_TYPE(TypeTag, BaseFluxVariables) BaseFluxVariables;
     typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
     typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
@@ -81,61 +82,60 @@ class ThreePTwoCNIFluxVariables : public GET_PROP_TYPE(TypeTag, BaseFluxVariable
     };
 
 public:
-    /*
-     * \brief The constructor
+    /*!
+     * \brief Compute / update the flux variables
      *
      * \param problem The problem
      * \param element The finite element
-     * \param fvGeometry The finite-volume geometry in the fully implicit scheme
-     * \param faceIdx The local index of the SCV (sub-control-volume) face
+     * \param fvGeometry The finite-volume geometry
+     * \param fIdx The local index of the SCV (sub-control-volume) face
      * \param elemVolVars The volume variables of the current element
      * \param onBoundary A boolean variable to specify whether the flux variables
      * are calculated for interior SCV faces or boundary faces, default=false
      */
-    ThreePTwoCNIFluxVariables(const Problem &problem,
-                              const Element &element,
-                              const FVElementGeometry &fvGeometry,
-                              const int faceIdx,
-                              const ElementVolumeVariables &elemVolVars,
-                              const bool onBoundary = false)
-    : BaseFluxVariables(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary)
+    void update(const Problem &problem,
+                const Element &element,
+                const FVElementGeometry &fvGeometry,
+                const int fIdx,
+                const ElementVolumeVariables &elemVolVars,
+                const bool onBoundary = false)
     {
-        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
-            density_[phaseIdx] = Scalar(0);
-            molarDensity_[phaseIdx] = Scalar(0);
-            massFractionCompWGrad_[phaseIdx] = Scalar(0);
-            massFractionCompNGrad_[phaseIdx] = Scalar(0);
-            moleFractionCompWGrad_[phaseIdx] = Scalar(0);
-            moleFractionCompNGrad_[phaseIdx] = Scalar(0);
-        }
-
-        calculateGradients_(problem, element, elemVolVars);
-        calculateporousDiffCoeff_(problem, element, elemVolVars);
+        BaseFluxVariables::update(problem, element, fvGeometry, fIdx, elemVolVars, onBoundary);
+        calculatePorousDiffCoeff_(problem, element, elemVolVars);
 
         // The spatial parameters calculates the actual heat flux vector
         DimVector temperatureGrad(0);
         DimVector tmp(0.0);
-        problem.spatialParams().matrixHeatFlux(tmp,
-                                                   *this,
-                                                   elemVolVars,
-                                                   temperatureGrad,
-                                                   element,
-                                                   fvGeometry,
-                                                   faceIdx);
+        problem.spatialParams().matrixHeatFlux(tmp, *this,
+                                               elemVolVars,
+                                               temperatureGrad,
+                                               element,
+                                               fvGeometry,
+                                               fIdx);
         // project the heat flux vector on the face's normal vector
         normalMatrixHeatFlux_ = tmp*this->face().normal;
-    };
+    }
 
 private:
     void calculateGradients_(const Problem &problem,
                              const Element &element,
                              const ElementVolumeVariables &elemVolVars)
     {
+        // initialize to zero
+        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            density_[phaseIdx] = Scalar(0);
+            molarDensity_[phaseIdx] = Scalar(0);
+            massFractionCompWGrad_[phaseIdx] = Scalar(0);
+            massFractionCompNGrad_[phaseIdx] = Scalar(0);
+            moleFractionCompWGrad_[phaseIdx] = Scalar(0);
+            moleFractionCompNGrad_[phaseIdx] = Scalar(0);
+        }
+
+        BaseFluxVariables::calculateGradients_(problem, element, elemVolVars);
+
         // calculate gradients
         DimVector tmp(0.0);
-        for (int idx = 0;
-             idx < this->face().numFap;
-             idx++) // loop over adjacent vertices
+        for (int idx = 0; idx < this->face().numFap; idx++) // loop over adjacent vertices
         {
             // FE gradient at vertex idx
             const DimVector &feGrad = this->face().grad[idx];
@@ -196,39 +196,24 @@ private:
             moleFractionCompNGrad_[gPhaseIdx] += tmp;
 
         }
-            // calculate temperature gradient using finite element
-            // gradients
-            DimVector temperatureGrad(0);
-            for (int idx = 0; idx < this->face().numFap; idx++)
-            {
-                tmp = this->face().grad[idx];
+        // calculate temperature gradient using finite element
+        // gradients
+        DimVector temperatureGrad(0);
+        for (int idx = 0; idx < this->face().numFap; idx++)
+        {
+            tmp = this->face().grad[idx];
 
-                // index for the element volume variables
-                int volVarsIdx = this->face().fapIndices[idx];
+            // index for the element volume variables
+            int volVarsIdx = this->face().fapIndices[idx];
 
-                tmp *= elemVolVars[volVarsIdx].temperature();
-                temperatureGrad += tmp;
-            }
+            tmp *= elemVolVars[volVarsIdx].temperature();
+            temperatureGrad += tmp;
         }
-
-    Scalar rhoFactor_(int phaseIdx, int scvIdx, const ElementVolumeVariables &elemVolVars)
-    {
-        static const Scalar eps = 1e-2;
-        const Scalar sat = elemVolVars[scvIdx].density(phaseIdx);
-        if (sat > eps)
-            return 0.5;
-        if (sat <= 0)
-            return 0;
-
-        static const Dumux::Spline<Scalar> sp(0, eps, // x0, x1
-                                              0, 0.5, // y0, y1
-                                              0, 0); // m0, m1
-        return sp.eval(sat);
     }
 
-    void calculateporousDiffCoeff_(const Problem &problem,
-                               const Element &element,
-                               const ElementVolumeVariables &elemVolVars)
+    void calculatePorousDiffCoeff_(const Problem &problem,
+                                   const Element &element,
+                                   const ElementVolumeVariables &elemVolVars)
     {
 
         const VolumeVariables &volVarsI = elemVolVars[this->face().i];

lecture/mm/heavyoil/3p2cni/localresidual.hh

@@ -176,14 +176,15 @@ public:
      * \param onBoundary A boolean variable to specify whether the flux variables
      *        are calculated for interior SCV faces or boundary faces, default=false
      */
-    void computeFlux(PrimaryVariables &flux, const int faceIdx, const bool onBoundary=false) const
+    void computeFlux(PrimaryVariables &flux, const int fIdx, const bool onBoundary=false) const
     {
-        FluxVariables fluxVars(this->problem_(),
-                           this->element_(),
-                           this->fvGeometry_(),
-                           faceIdx,
-                           this->curVolVars_(),
-                           onBoundary);
+        FluxVariables fluxVars;
+        fluxVars.update(this->problem_(),
+                        this->element_(),
+                        this->fvGeometry_(),
+                        fIdx,
+                        this->curVolVars_(),
+                        onBoundary);
 
         flux = 0;
         asImp_()->computeAdvectiveFlux(flux, fluxVars);