[wmpfa][neumann] Update handling of Neumann faces

dumux/discretization/cellcentered/wmpfa/facedatahandle.hh

@@ -90,8 +90,8 @@ public:
         clear();
     }
 
-    template<class TF, class EG, class SCVF>
-    void decompose(const Interpolator& intOp, const TF& tensor, const EG& fvGeometry, const SCVF& scvf)
+    template<class Problem, class TF, class EG, class SCVF>
+    void decompose(const Problem& problem, const Interpolator& intOp, const TF& tensor, const EG& fvGeometry, const SCVF& scvf)
     {
         boundaryFace_ = scvf.boundary();
         const auto coNormal = mv(tensor(scvf.insideScvIdx()), scvf.unitOuterNormal());
@@ -102,26 +102,41 @@ public:
         else
         {
             WMpfaHelper::eraseZeros(coeff, indices);
-            updateEntries(intOp, indices, coeff, fvGeometry, scvf);
+            updateEntries(problem, intOp, indices, coeff, fvGeometry, scvf);
         }
     }
 
-    template<class Indices, class Coeff, class EG, class SCVF>
-    void updateEntries(const Interpolator& intOp, const Indices& indices, const Coeff& coeff, const EG& fvGeometry, const SCVF& scvf)
+    template<class Problem, class Indices, class Coeff, class EG, class SCVF>
+    void updateEntries(const Problem& problem, const Interpolator& intOp, const Indices& indices, const Coeff& coeff, const EG& fvGeometry, const SCVF& scvf)
     {
         const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
         entries_.push_back({0.0, scv.dofIndex(), scv.center()});
         for(std::size_t i=0; i<indices.size(); ++i)
         {
             const auto& intData = intOp.getInterpolationData(indices[i]);
-            for(const auto& e : intData.entries())
+            bool skipEntries = false;
+            if (fvGeometry.hasBoundaryScvf())
             {
-                if(e.dofIndex() != scvf.insideScvIdx())
+                const auto& scvfJ = fvGeometry.scvf(intData.scvfIdx());
+                if(scvfJ.boundary())
                 {
-                    auto c = coeff[i]*e.weight();
-                    entries_[0].coefficient += c;
-                    //ToDo pos is not needed for each entry!!!
-                    entries_.push_back({-1.0*c, e.dofIndex(), intData.position()});
+                    const auto& bcTypes = problem.boundaryTypes(fvGeometry.gridGeometry().element(scvfJ.insideScvIdx()), scvfJ);
+                    if(bcTypes.hasNeumann())
+                        skipEntries = true;
+                }
+            }
+
+            if(!skipEntries)
+            {
+                for(const auto& e : intData.entries())
+                {
+                    if(e.dofIndex() != scvf.insideScvIdx())
+                    {
+                        auto c = coeff[i]*e.weight();
+                        entries_[0].coefficient += c;
+                        //ToDo pos is not needed for each entry!!!
+                        entries_.push_back({-1.0*c, e.dofIndex(), intData.position()});
+                    }
                 }
             }
         }

dumux/discretization/cellcentered/wmpfa/interpolationoperator.hh

@@ -73,8 +73,9 @@ class HapInterpolatorBase
     {
         LocalInterpolationData() = default;
 
-        LocalInterpolationData(Position p, Entry e1, Entry e2)
+        LocalInterpolationData(GridIndexType idx, Position p, Entry e1, Entry e2)
         {
+            scvfIdx_ = idx;
             position_ = p;
             entries_ = {e1,e2};
         }
@@ -85,7 +86,11 @@ class HapInterpolatorBase
         Position position() const
         {return position_;}
 
+        GridIndexType scvfIdx() const
+        {return scvfIdx_;}
+
         private:
+            GridIndexType scvfIdx_ = {};
             Position position_ = {};
             std::array<Entry, 2> entries_ = {Entry(), Entry()};
     };
@@ -133,11 +138,11 @@ public:
                                 + (distL*distK / (distL*tauK + distK*tauL))
                                    * mv(tensor(insideVolVars) - tensor(outsideVolVars), scvf.unitOuterNormal());
 
-                interpolationData_[scvf.localIndex()] = LocalInterpolationData(position, Entry(insideScvIdx, omegaK), Entry(outsideScvIdx, omegaL) );
+                interpolationData_[scvf.localIndex()] = LocalInterpolationData(scvf.index(), position, Entry(insideScvIdx, omegaK), Entry(outsideScvIdx, omegaL) );
             }
             else
             {
-                interpolationData_[scvf.localIndex()] = LocalInterpolationData(scvf.center(), Entry(insideScvIdx, 0.0), Entry(outsideScvIdx, 1.0) );
+                interpolationData_[scvf.localIndex()] = LocalInterpolationData(scvf.index(), scvf.center(), Entry(insideScvIdx, 0.0), Entry(outsideScvIdx, 1.0) );
             }
         }
         isUpdated_ = true;

dumux/porousmediumflow/fluxvariablescachefiller.hh

@@ -814,7 +814,7 @@ private:
     void prepareHandle_(Handle& handle, const IntOP& intOp, const TensorFunction& tensor, const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf)
     {
         handle.prepare();
-        handle.decompose(intOp, tensor, fvGeometry, scvf);
+        handle.decompose(problem(), intOp, tensor, fvGeometry, scvf);
     }
 
     //! method to fill the advective quantities