[mpfa] restructure global interaction volume seeds

dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseeds.hh

@@ -33,14 +33,7 @@ namespace Dumux
 //! By default we simply inherit from the base class
 //! Actual implementations for other methods have to be provided below
 template<class TypeTag, MpfaMethods method>
-class CCMpfaGlobalInteractionVolumeSeedsImplementation : public CCMpfaGlobalInteractionVolumeSeedsBase<TypeTag>
-{
-    using ParentType = CCMpfaGlobalInteractionVolumeSeedsBase<TypeTag>;
-    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
-
-public:
-    CCMpfaGlobalInteractionVolumeSeedsImplementation(const GridView gridView) : ParentType(gridView)  {}
-};
+class CCMpfaGlobalInteractionVolumeSeedsImplementation;
 
 /*!
  * \ingroup Mpfa
@@ -52,6 +45,8 @@ using CCMpfaGlobalInteractionVolumeSeeds = CCMpfaGlobalInteractionVolumeSeedsImp
 } // end namespace
 
 // the specializations of this class differing from the default have to be included here
+#include <dumux/discretization/cellcentered/mpfa/omethod/globalinteractionvolumeseeds.hh>
 #include <dumux/discretization/cellcentered/mpfa/lmethod/globalinteractionvolumeseeds.hh>
+#include <dumux/discretization/cellcentered/mpfa/omethodfps/globalinteractionvolumeseeds.hh>
 
 #endif

dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseedsbase.hh

@@ -36,8 +36,6 @@ template<class TypeTag>
 class CCMpfaGlobalInteractionVolumeSeedsBase
 {
     using Implementation = typename GET_PROP_TYPE(TypeTag, GlobalInteractionVolumeSeeds);
-    // the actual implementation needs to be friend
-    friend Implementation;
 
     using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
@@ -51,21 +49,18 @@ class CCMpfaGlobalInteractionVolumeSeedsBase
     using IndexType = typename GridView::IndexSet::IndexType;
 
 public:
-    CCMpfaGlobalInteractionVolumeSeedsBase(const GridView gridView) : gridView_(gridView) {}
+    CCMpfaGlobalInteractionVolumeSeedsBase(const GridView& gridView) : gridView_(gridView) {}
 
     // initializes the interaction volumes or the seeds
-    void update(const Problem& problem, std::vector<bool>& boundaryVertices)
+    void update(const Problem& p, const std::vector<bool>& boundaryVertices)
     {
-        problemPtr_ = &problem;
-        seeds_.clear();
-        boundarySeeds_.clear();
-
-        auto numScvf = problem_().model().globalFvGeometry().numScvf();
-        scvfIndexMap_.resize(numScvf);
-        std::vector<bool> isFaceHandled(numScvf, false);
+        problemPtr_ = &p;
 
         // initialize the seeds according to the mpfa method
-        asImp_().initializeSeeds_(boundaryVertices, isFaceHandled);
+        asImp_().initializeSeeds(boundaryVertices,
+                                 scvfIndexMap_,
+                                 seeds_,
+                                 boundarySeeds_);
     }
 
     const InteractionVolumeSeed& seed(const SubControlVolumeFace& scvf) const
@@ -74,70 +69,13 @@ public:
     const BoundaryInteractionVolumeSeed& boundarySeed(const SubControlVolumeFace& scvf) const
     { return boundarySeeds_[scvfIndexMap_[scvf.index()]]; }
 
-private:
+    const Problem& problem() const
+    { return *problemPtr_; }
 
-    void initializeSeeds_(std::vector<bool>& boundaryVertices,
-                          std::vector<bool>& isFaceHandled)
-    {
-        const auto numScvf = problem_().model().globalFvGeometry().numScvf();
-        const auto numBoundaryScvf = problem_().model().globalFvGeometry().numBoundaryScvf();
-
-        // reserve memory
-        seeds_.reserve(numScvf - numBoundaryScvf);
-        boundarySeeds_.reserve(numBoundaryScvf);
-
-        IndexType boundarySeedIndex = 0;
-        IndexType seedIndex = 0;
-        for (const auto& element : elements(gridView_))
-        {
-            auto fvGeometry = localView(problem_().model().globalFvGeometry());
-            fvGeometry.bindElement(element);
-            for (const auto& scvf : scvfs(fvGeometry))
-            {
-                // skip the rest if we already handled this face
-                if (isFaceHandled[scvf.index()])
-                    continue;
-
-                if (boundaryVertices[scvf.vertexIndex()])
-                {
-                    // make the boundary interaction volume seed
-                    boundarySeeds_.emplace_back(Helper::makeBoundaryInteractionVolumeSeed(problem_(), element, fvGeometry, scvf));
-
-                    // update the index map entries for the global scv faces in the interaction volume
-                    for (auto scvfIdxGlobal : boundarySeeds_.back().globalScvfIndices())
-                    {
-                        scvfIndexMap_[scvfIdxGlobal] = boundarySeedIndex;
-                        isFaceHandled[scvfIdxGlobal] = true;
-                    }
-
-                    // increment counter
-                    boundarySeedIndex++;
-                }
-                else
-                {
-                    // make the inner interaction volume seed
-                    seeds_.emplace_back(Helper::makeInnerInteractionVolumeSeed(problem_(), element, fvGeometry, scvf));
-
-                    // update the index map entries for the global scv faces in the interaction volume
-                    for (auto scvfIdxGlobal : seeds_.back().globalScvfIndices())
-                    {
-                        scvfIndexMap_[scvfIdxGlobal] = seedIndex;
-                        isFaceHandled[scvfIdxGlobal] = true;
-                    }
-
-                    // increment counter
-                    seedIndex++;
-                }
-            }
-        }
-
-        // shrink vectors to actual size
-        seeds_.shrink_to_fit();
-        boundarySeeds_.shrink_to_fit();
-    }
+    const GridView& gridView() const
+    { return gridView_; }
 
-    const Problem& problem_() const
-    { return *problemPtr_; }
+private:
 
     const Implementation& asImp_() const
     { return *static_cast<Implementation*>(this); }

dumux/discretization/cellcentered/mpfa/lmethod/globalinteractionvolumeseeds.hh

@@ -33,11 +33,10 @@ namespace Dumux
  * \brief Specialization of the class for the mpfa-l method.
  */
 template<class TypeTag>
-class CCMpfaGlobalInteractionVolumeSeedsImplementation<TypeTag, MpfaMethods::lMethod> : public CCMpfaGlobalInteractionVolumeSeedsBase<TypeTag>
+class CCMpfaGlobalInteractionVolumeSeedsImplementation<TypeTag, MpfaMethods::lMethod>
+       : public CCMpfaGlobalInteractionVolumeSeedsBase<TypeTag>
 {
     using ParentType = CCMpfaGlobalInteractionVolumeSeedsBase<TypeTag>;
-    // the parent needs to be friend to access the private initialization method
-    friend ParentType;
 
     using InteractionVolume = typename GET_PROP_TYPE(TypeTag, InteractionVolume);
     using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
@@ -50,25 +49,37 @@ class CCMpfaGlobalInteractionVolumeSeedsImplementation<TypeTag, MpfaMethods::lMe
     using Element = typename GridView::template Codim<0>::Entity;
 
 public:
-    CCMpfaGlobalInteractionVolumeSeedsImplementation(const GridView gridView) : ParentType(gridView)  {}
+    CCMpfaGlobalInteractionVolumeSeedsImplementation(const GridView& gridView) : ParentType(gridView) {}
 
-private:
-    void initializeSeeds_(std::vector<bool>& boundaryVertices, std::vector<bool>& isFaceHandled)
+    template<typename SeedVector, typename BoundarySeedVector>
+    void initializeSeeds(const std::vector<bool>& boundaryVertices,
+                         std::vector<IndexType>& scvfIndexMap,
+                         SeedVector& seeds,
+                         BoundarySeedVector& boundarySeeds)
     {
-        const auto numScvf = this->problem_().model().globalFvGeometry().numScvf();
-        const auto numBoundaryScvf = this->problem_().model().globalFvGeometry().numBoundaryScvf();
+        seeds.clear();
+        boundarySeeds.clear();
+        scvfIndexMap.clear();
 
         // reserve memory
-        this->seeds_.reserve(numScvf - numBoundaryScvf);
-        this->boundarySeeds_.reserve(numBoundaryScvf);
+        const auto numScvf = this->problem().model().globalFvGeometry().numScvf();
+        const auto numBoundaryScvf = this->problem().model().globalFvGeometry().numBoundaryScvf();
+
+        seeds.reserve( std::size_t((numScvf-numBoundaryScvf)/2) );
+        boundarySeeds.reserve(numBoundaryScvf);
+        scvfIndexMap.resize(numScvf);
+
+        // Keep track of which faces have been handled already
+        std::vector<bool> isFaceHandled(numScvf, false);
 
         IndexType boundarySeedIndex = 0;
         IndexType seedIndex = 0;
-        for (const auto& element : elements(this->gridView_))
+        for (const auto& element : elements(this->gridView()))
         {
-            auto fvGeometry = localView(this->problem_().model().globalFvGeometry());
+            auto fvGeometry = localView(this->problem().model().globalFvGeometry());
             fvGeometry.bindElement(element);
-            for (const auto& scvf : scvfs(fvGeometry))
+
+            for (auto&& scvf : scvfs(fvGeometry))
             {
                 // skip the rest if we already handled this face
                 if (isFaceHandled[scvf.index()])
@@ -77,12 +88,15 @@ private:
                 if (boundaryVertices[scvf.vertexIndex()])
                 {
                     // make the boundary interaction volume seed
-                    this->boundarySeeds_.emplace_back(Helper::makeBoundaryInteractionVolumeSeed(this->problem_(), element, fvGeometry, scvf));
+                    boundarySeeds.emplace_back(Helper::makeBoundaryInteractionVolumeSeed(this->problem(),
+                                                                                         element,
+                                                                                         fvGeometry,
+                                                                                         scvf));
 
                     // update the index map entries for the global scv faces in the interaction volume
-                    for (auto scvfIdxGlobal : this->boundarySeeds_.back().globalScvfIndices())
+                    for (auto scvfIdxGlobal : boundarySeeds.back().globalScvfIndices())
                     {
-                        this->scvfIndexMap_[scvfIdxGlobal] = boundarySeedIndex;
+                        scvfIndexMap[scvfIdxGlobal] = boundarySeedIndex;
                         isFaceHandled[scvfIdxGlobal] = true;
                     }
 
@@ -92,14 +106,17 @@ private:
                 else
                 {
                     // make the inner interaction volume seed only if we are on highest level of all connected elements
-                    if (isLocalMaxLevel(element, scvf))
+                    if (isLocalMaxLevel_(element, scvf))
                     {
-                        this->seeds_.emplace_back(Helper::makeInnerInteractionVolumeSeed(this->problem_(), element, fvGeometry, scvf));
+                        seeds.emplace_back(Helper::makeInnerInteractionVolumeSeed(this->problem(),
+                                                                                  element,
+                                                                                  fvGeometry,
+                                                                                  scvf));
 
                         // update the index map entries for the global scv faces in the interaction volume
-                        for (auto scvfIdxGlobal : this->seeds_.back().globalScvfIndices())
+                        for (auto scvfIdxGlobal : seeds.back().globalScvfIndices())
                         {
-                            this->scvfIndexMap_[scvfIdxGlobal] = seedIndex;
+                            scvfIndexMap[scvfIdxGlobal] = seedIndex;
                             isFaceHandled[scvfIdxGlobal] = true;
                         }
 
@@ -109,18 +126,13 @@ private:
                 }
             }
         }
-
-        // shrink vectors to actual size
-        this->seeds_.shrink_to_fit();
-        this->boundarySeeds_.shrink_to_fit();
     }
 
-    bool isLocalMaxLevel(const Element& element, const SubControlVolumeFace& scvf) const
+    bool isLocalMaxLevel_(const Element& element, const SubControlVolumeFace& scvf) const
     {
         auto inLevel = element.level();
-        for (auto outsideIdx : scvf.outsideScvIndices())
-            if (this->problem_().model().globalFvGeometry().element(outsideIdx).level() > inLevel)
-                return false;
+        if (this->problem().model().globalFvGeometry().element(scvf.outsideScvIdx()).level() > inLevel)
+            return false;
         return true;
     }
 };

dumux/discretization/cellcentered/mpfa/omethod/globalinteractionvolumeseeds.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Base class for the global interaction volumes of the mpfa-o method.
+ */
+#ifndef DUMUX_DISCRETIZATION_MPFA_O_GLOBALINTERACTIONVOLUMESEEDS_HH
+#define DUMUX_DISCRETIZATION_MPFA_O_GLOBALINTERACTIONVOLUMESEEDS_HH
+
+#include <dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseedsbase.hh>
+#include <dumux/discretization/cellcentered/mpfa/methods.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup Mpfa
+ * \brief Specialization of the class for the mpfa-o method.
+ */
+template<class TypeTag>
+class CCMpfaGlobalInteractionVolumeSeedsImplementation<TypeTag, MpfaMethods::oMethod>
+       : public CCMpfaGlobalInteractionVolumeSeedsBase<TypeTag>
+{
+    using ParentType = CCMpfaGlobalInteractionVolumeSeedsBase<TypeTag>;
+
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using Helper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
+    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+    using InteractionVolume = typename GET_PROP_TYPE(TypeTag, InteractionVolume);
+    using InteractionVolumeSeed = typename InteractionVolume::Seed;
+    using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
+    using BoundaryInteractionVolumeSeed = typename BoundaryInteractionVolume::Seed;
+
+    using IndexType = typename GridView::IndexSet::IndexType;
+
+    static const int dim = GridView::dimension;
+
+public:
+    CCMpfaGlobalInteractionVolumeSeedsImplementation(const GridView& gridView) : ParentType(gridView) {}
+
+    template<typename SeedVector, typename BoundarySeedVector>
+    void initializeSeeds(const std::vector<bool>& boundaryVertices,
+                         std::vector<IndexType>& scvfIndexMap,
+                         SeedVector& seeds,
+                         BoundarySeedVector& boundarySeeds)
+    {
+        seeds.clear();
+        boundarySeeds.clear();
+        scvfIndexMap.clear();
+
+        // reserve memory
+        const auto numScvf = this->problem().model().globalFvGeometry().numScvf();
+        const auto numBoundaryVertices = this->problem().model().globalFvGeometry().numBoundaryVertices();
+        const auto numInteriorVertices = this->gridView().size(dim) - numBoundaryVertices;
+
+        if (numInteriorVertices > 0)
+            seeds.reserve(numInteriorVertices);
+        boundarySeeds.reserve(numBoundaryVertices);
+        scvfIndexMap.resize(numScvf);
+
+        // Keep track of which faces have been handled already
+        std::vector<bool> isFaceHandled(numScvf, false);
+
+        IndexType boundarySeedIndex = 0;
+        IndexType seedIndex = 0;
+        for (const auto& element : elements(this->gridView()))
+        {
+            auto fvGeometry = localView(this->problem().model().globalFvGeometry());
+            fvGeometry.bindElement(element);
+
+            for (auto&& scvf : scvfs(fvGeometry))
+            {
+                // skip the rest if we already handled this face
+                if (isFaceHandled[scvf.index()])
+                    continue;
+
+                if (boundaryVertices[scvf.vertexIndex()])
+                {
+                    // make the boundary interaction volume seed
+                    boundarySeeds.emplace_back(Helper::makeBoundaryInteractionVolumeSeed(this->problem(),
+                                                                                         element,
+                                                                                         fvGeometry,
+                                                                                         scvf));
+
+                    // update the index map entries for the global scv faces in the interaction volume
+                    for (auto scvfIdxGlobal : boundarySeeds.back().globalScvfIndices())
+                    {
+                        scvfIndexMap[scvfIdxGlobal] = boundarySeedIndex;
+                        isFaceHandled[scvfIdxGlobal] = true;
+                    }
+
+                    // increment counter
+                    boundarySeedIndex++;
+                }
+                else
+                {
+                    // make the inner interaction volume seed
+                    seeds.emplace_back(Helper::makeInnerInteractionVolumeSeed(this->problem(),
+                                                                              element,
+                                                                              fvGeometry,
+                                                                              scvf));
+
+                    // update the index map entries for the global scv faces in the interaction volume
+                    for (auto scvfIdxGlobal : seeds.back().globalScvfIndices())
+                    {
+                        scvfIndexMap[scvfIdxGlobal] = seedIndex;
+                        isFaceHandled[scvfIdxGlobal] = true;
+                    }
+
+                    // increment counter
+                    seedIndex++;
+                }
+            }
+        }
+    }
+};
+} // end namespace
+
+
+#endif

dumux/discretization/cellcentered/mpfa/omethodfps/globalinteractionvolumeseeds.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Base class for the global interaction volumes of the mpfa-o-fps method.
+ */
+#ifndef DUMUX_DISCRETIZATION_MPFA_O_FPS_GLOBALINTERACTIONVOLUMESEEDS_HH
+#define DUMUX_DISCRETIZATION_MPFA_O_FPS_GLOBALINTERACTIONVOLUMESEEDS_HH
+
+#include <dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseedsbase.hh>
+#include <dumux/discretization/cellcentered/mpfa/methods.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup Mpfa
+ * \brief Specialization of the class for the mpfa-o-fps method.
+ */
+template<class TypeTag>
+class CCMpfaGlobalInteractionVolumeSeedsImplementation<TypeTag, MpfaMethods::oMethodFps>
+      : public CCMpfaGlobalInteractionVolumeSeedsImplementation<TypeTag, MpfaMethods::oMethod>
+{
+    using ParentType = CCMpfaGlobalInteractionVolumeSeedsImplementation<TypeTag, MpfaMethods::oMethod>;
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+public:
+    CCMpfaGlobalInteractionVolumeSeedsImplementation(const GridView& gridView) : ParentType(gridView) {}
+};
+} // end namespace
+
+
+#endif