[mpfa] change iv index set concept

the dual grid index set now does not contain the data on the global scv indices
that are connected to the faces. The local data is only available in the iv index set.
This should be better for methods other than the mpfa-o method, as in those the local mappings
of the dual grid index set is maybe useless (e.g. for the mpfa-l scheme).

dumux/discretization/cellcentered/mpfa/dualgridindexset.hh

@@ -67,10 +67,7 @@ public:
                 const std::vector<IndexType>& outsideScvIndices)
     {
         // check if it this really hasn't been called yet
-        assert( std::find_if( scvfIndices_.begin(),
-                              scvfIndices_.end(),
-                              scvfIdx
-                            ) == scvfIndices_.end() && "scvf has already been inserted!");
+        assert( std::find(scvfIndices_.begin(), scvfIndices_.end(), scvfIdx ) == scvfIndices_.end() && "scvf has already been inserted!");
 
         // the local index of the scvf data about to be inserted
         const LocalIndexType curScvfLocalIdx = scvfIndices_.size();
@@ -82,12 +79,12 @@ public:
         scvIndices.insert(scvIndices.end(), outsideScvIndices.begin(), outsideScvIndices.end());
 
         // if scvf is on boundary, increase counter
-        if (boundary)
-            numBoundaryScvfs_++;
+        if (boundary) numBoundaryScvfs_++;
 
+        // insert data on the new scv
         scvfIndices_.push_back(scvfIdx);
         scvfIsOnBoundary_.push_back(boundary);
-        scvfNeighborScvIndicesGlobal_.emplace_back( std::move(scvIndices) );
+        scvfNeighborScvIndices_.emplace_back( std::move(scvIndices) );
 
         // if entry for the inside scv exists append scvf local index, create otherwise
         auto it = std::find( scvIndices_.begin(), scvIndices_.end(), insideScvIdx );
@@ -106,28 +103,6 @@ public:
         }
     }
 
-    // This should be called AFTER (!) all scvf data has been inserted
-    void makeLocal()
-    {
-        // make sure the created index set makes sense so far
-        assert(checkGlobalIndexSet_());
-
-        // compute local neighboring scv indices for the scvfs
-        scvfNeighborScvIndices_.resize(numScvfs());
-        for (unsigned int i = 0; i < numScvfs(); ++i)
-        {
-            const auto& neighborsGlobal = scvfNeighborScvIndicesGlobal_[i];
-            const auto numNeighbors = scvfIsOnBoundary_[i] ? 1 : neighborsGlobal.size();
-
-            scvfNeighborScvIndices_[i].resize(numNeighbors);
-            for (unsigned int j = 0; j < numNeighbors; ++j)
-                scvfNeighborScvIndices_[i][j] = findLocalScvIdx_(neighborsGlobal[j]);
-        }
-
-        // delete global neighboring scv data (not used anymore)
-        scvfNeighborScvIndicesGlobal_.clear();
-    }
-
     //! returns the number of scvs
     std::size_t numScvs() const
     { return scvIndices_.size(); }
@@ -148,52 +123,37 @@ public:
     const GlobalIndexContainer& globalScvfIndices() const
     { return scvfIndices_; }
 
-    //! returns a global scv idx for a given local scv index
-    IndexType scvIdxGlobal(LocalIndexType scvIdxLocal) const
-    { return scvIndices_[scvIdxLocal]; }
-
-    //! returns the local indices of the scvfs embedded in a local scv
-    const LocalIndexContainer& localScvfIndicesInScv(LocalIndexType scvIdxLocal) const
-    { return localScvfIndicesInScv_[scvIdxLocal]; }
-
-    //! returns a global scvf idx for a given local scvf index
-    IndexType scvfIdxGlobal(LocalIndexType scvfIdxLocal) const
-    { return scvfIndices_[scvfIdxLocal]; }
+    //! returns the global scv idx of the i-th scv
+    IndexType scvIdxGlobal(unsigned int i) const
+    { return scvIndices_[i]; }
 
-    //! returns whether or not an scvf touches the boundary
-    bool scvfIsOnBoundary(LocalIndexType scvfIdxLocal) const
-    { return scvfIsOnBoundary_[scvfIdxLocal]; }
-
-    //! returns the local indices of the neighboring scvs of an scvf
-    const LocalIndexContainer& localNeighboringScvIndices(LocalIndexType scvfIdxLocal) const
-    { return scvfNeighborScvIndices_[scvfIdxLocal]; }
-
-private:
-    //! returns the node-local scv index to a given global scv index
-    unsigned int findLocalScvIdx_(IndexType globalScvIdx) const
+    //! returns the global index of the j-th scvf embedded in the i-th scv
+    IndexType scvfIdxGlobal(unsigned int i, unsigned int j) const
     {
-        auto it = std::find( scvIndices_.begin(), scvIndices_.end(), globalScvIdx );
-        assert(it != scvIndices_.end() && "Global scv index not found in local container!");
-        return std::distance(scvIndices_.begin(), it);
+        assert(j < dim); // only dim faces can be embedded in an scv
+        return scvfIndices_[ localScvfIndicesInScv_[i][j] ];
     }
 
-    //! checks whether or not the inserted index set makes sense
-    bool checkGlobalIndexSet_() const
+    //! returns the nodel-local index of the j-th scvf embedded in the i-th scv
+    IndexType scvfIdxLocal(unsigned int i, unsigned int j) const
     {
-        //! do all scvs have dim embedded scvfs?
-        for (const auto& scvfIndices : localScvfIndicesInScv_)
-            if (scvfIndices.size() != dim)
-                DUNE_THROW(Dune::InvalidStateException, "Number of scv faces found for an scv != dimension");
-
-        //! are all scvfs unique?
-        for (const auto& scvfIdx : scvfIndices_)
-            for (const auto& scvfIdx2 : scvfIndices_)
-                if (scvfIdx == scvfIdx2)
-                    DUNE_THROW(Dune::InvalidStateException, "Inserted scvfs seem to not be unique");
-
-        return true;
+        assert(j < dim); // only dim faces can be embedded in an scv
+        return localScvfIndicesInScv_[i][j];
     }
 
+    //! returns the index of the i-th scvf
+    IndexType scvfIdxGlobal(unsigned int i) const
+    { return scvfIndices_[i]; }
+
+    //! returns whether or not the i-th scvf touches the boundary
+    bool scvfIsOnBoundary(unsigned int i) const
+    { return scvfIsOnBoundary_[i]; }
+
+    //! returns the indices of the neighboring scvs of the i-th scvf
+    const GlobalIndexContainer& neighboringScvIndices(unsigned int i) const
+    { return scvfNeighborScvIndices_[i]; }
+
+private:
     //! the indices of the scvs around a dual grid node
     GlobalIndexContainer scvIndices_;
     //! maps to each scv a list of scvf indices embedded in it
@@ -203,12 +163,9 @@ private:
     GlobalIndexContainer scvfIndices_;
     //! maps to each scvf a boolean to indicate if it is on the boundary
     std::vector<bool> scvfIsOnBoundary_;
-    //! maps to each scvf a list of neighbouring local scv indices
+    //! maps to each scvf a list of neighbouring scv indices
     //! ordering: 0 - inside scv idx; 1..n - outside scv indices
-    std::vector<LocalIndexContainer> scvfNeighborScvIndices_;
-    //! maps to each scvf a list of neighbouring global scv indices
-    //! This container is destroyed when makeLocal() is called
-    std::vector<GlobalIndexContainer> scvfNeighborScvIndicesGlobal_;
+    std::vector<GlobalIndexContainer> scvfNeighborScvIndices_;
     //! stores how many boundary scvfs are embedded in this dual grid node
     std::size_t numBoundaryScvfs_;
 };
@@ -236,14 +193,6 @@ public:
         (*this).resize(gridView.size(dim));
     }
 
-    //! after the global data has been inserted, this
-    //! method should be called to set up the local index sets
-    void makeLocalIndexSets()
-    {
-        for (auto& node : (*this))
-            node.makeLocal();
-    }
-
     //! access with an scvf
     const NodalIndexSet& operator[] (const SubControlVolumeFace& scvf) const
     { return (*this)[scvf.vertexIndex()]; }

dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh

@@ -361,9 +361,6 @@ public:
             }
         }
 
-        // make the local index sets of the dual grid
-        dualIdSet.makeLocalIndexSets();
-
         // calculate the total number of vertices using the primary/secondary interaction volumes
         numVertUsingPrimaryIV_ = 0;
         numVertUsingSecondaryIV_ = 0;
@@ -714,9 +711,6 @@ public:
             neighborVolVarIndices_[eIdx] = neighborVolVarIndexSet;
         }
 
-        // make the local index sets of the dual grid
-        dualIdSet.makeLocalIndexSets();
-
         // calculate the total number of vertices using the primary/secondary interaction volumes
         numVertUsingPrimaryIV_ = 0;
         numVertUsingSecondaryIV_ = 0;

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

@@ -175,7 +175,7 @@ public:
         }
 
         // // set vol vars stencil & positions pointer in handle
-        dataHandle.setVolVarsStencilPointer(indexSet().nodalIndexSet().globalScvIndices());
+        dataHandle.setVolVarsStencilPointer(indexSet().globalScvIndices());
         dataHandle.setDirichletDataPointer(dirichletData_);
 
         // on surface grids, additionally prepare the outside transmissibilities
@@ -293,7 +293,7 @@ private:
         numFaces_ = indexSet().numFaces();
 
         //! number of interaction-volume-local (and node-local) scvs
-        const auto& scvIndices = indexSet().nodalIndexSet().globalScvIndices();
+        const auto& scvIndices = indexSet().globalScvIndices();
         const auto numLocalScvs = scvIndices.size();
 
         //! number of global scvfs appearing in this interaction volume
@@ -327,7 +327,7 @@ private:
         for (LocalIndexType faceIdxLocal = 0; faceIdxLocal < numFaces_; ++faceIdxLocal)
         {
             const auto scvfIdxGlobal = indexSet().scvfIdxGlobal(faceIdxLocal);
-            const auto& neighborScvIndicesLocal = indexSet().localNeighboringScvIndices(faceIdxLocal);
+            const auto& neighborScvIndicesLocal = indexSet().neighboringLocalScvIndices(faceIdxLocal);
             const auto insideLocalScvIdx = neighborScvIndicesLocal[0];
 
             // we have to use the "inside" scv face here
@@ -363,10 +363,9 @@ private:
                     // loop over scvfs in outside scv until we find the one coinciding with current scvf
                     for (int coord = 0; coord < dim; ++coord)
                     {
-                        if (indexSet().localScvfIndexInScv(outsideLocalScvIdx, coord) == faceIdxLocal)
+                        if (indexSet().scvfIdxLocal(outsideLocalScvIdx, coord) == faceIdxLocal)
                         {
-                            const auto nodeLocalScvfIdx = indexSet().nodalIndexSet().localScvfIndicesInScv(outsideLocalScvIdx)[coord];
-                            const auto globalScvfIdx = indexSet().nodalIndexSet().scvfIdxGlobal(nodeLocalScvfIdx);
+                            const auto globalScvfIdx = indexSet().nodalIndexSet().scvfIdxGlobal(outsideLocalScvIdx, coord);
                             const auto& flipScvf = fvGeometry().scvf(globalScvfIdx);
                             globalLocalScvfPairedData_.emplace_back(&flipScvf, LocalFaceData(faceIdxLocal, outsideLocalScvIdx, true));
                             globalScvfIndices_.push_back(flipScvf.index());

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

@@ -52,12 +52,14 @@ public:
         ivToNodeScvf_.reserve(numFaces_);
         nodeToIvScvf_.resize(numNodalScvfs);
 
+        //! the local neighboring scv indices of the faces
+        scvfNeighborScvLocalIndices_.reserve(numFaces_);
+
         //! keeps track of which nodal scvfs have been handled already
         std::vector<bool> isHandled(numNodalScvfs, false);
 
-        //! go over faces in nodal index set, check if iv-local face
-        //! has been inserted already for this scvf and if not, insert
-        //! the index transformation
+        //! go over faces in nodal index set, check if iv-local face has been
+        //! inserted already for this scvf and if not, insert index mapping
         unsigned int findCounter = 0;
         for (LocalIndexType i = 0; i < numNodalScvfs; ++i)
         {
@@ -78,7 +80,7 @@ public:
             //! to this face as well by comparing the set of neighboring scv indices.
             const auto scvIndices = [&nodalIndexSet, i] ()
                                     {
-                                        auto tmp = nodalIndexSet.localNeighboringScvIndices(i);
+                                        auto tmp = nodalIndexSet.neighboringScvIndices(i);
                                         std::sort(tmp.begin(), tmp.end());
                                         return tmp;
                                     } ();
@@ -92,7 +94,7 @@ public:
                 {
                     const auto scvIndices2 = [&nodalIndexSet, j] ()
                                              {
-                                                 auto tmp = nodalIndexSet.localNeighboringScvIndices(j);
+                                                 auto tmp = nodalIndexSet.neighboringScvIndices(j);
                                                  std::sort(tmp.begin(), tmp.end());
                                                  return tmp;
                                              } ();
@@ -122,6 +124,18 @@ public:
 
         //! ensure we found as many faces as anticipated
         assert(findCounter == numFaces_ && "Couldn't find as many faces as anticipated");
+
+        // compute local neighboring scv indices for the iv-local scvfs
+        scvfNeighborScvLocalIndices_.resize(numFaces());
+        for (unsigned int i = 0; i < numFaces(); ++i)
+        {
+            const auto& neighborsGlobal = nodalIndexSet_.neighboringScvIndices(ivToNodeScvf_[i]);
+            const auto numNeighbors = nodalIndexSet_.scvfIsOnBoundary(ivToNodeScvf_[i]) ? 1 : neighborsGlobal.size();
+
+            scvfNeighborScvLocalIndices_[i].resize(numNeighbors);
+            for (unsigned int j = 0; j < numNeighbors; ++j)
+                scvfNeighborScvLocalIndices_[i][j] = findLocalScvIdx_(neighborsGlobal[j]);
+        }
     }
 
     //! returns the corresponding nodal index set
@@ -132,17 +146,17 @@ public:
     GlobalIndexType scvfIdxGlobal(LocalIndexType ivLocalScvfIdx) const
     { return nodalIndexSet_.scvfIdxGlobal(ivToNodeScvf_[ivLocalScvfIdx]); }
 
-    //! returns the iv-local scvf idx of the i-th (coordDir) scvfs embedded in a local scv
-    LocalIndexType localScvfIndexInScv(LocalIndexType scvIdxLocal, unsigned int coordDir) const
-    { return nodeToIvScvf_[nodalIndexSet_.localScvfIndicesInScv(scvIdxLocal)[coordDir]]; }
+    //! returns the iv-local scvf idx of the i-th scvfs embedded in a local scv
+    LocalIndexType scvfIdxLocal(LocalIndexType scvIdxLocal, unsigned int i) const
+    { return nodeToIvScvf_[nodalIndexSet_.scvfIdxLocal(scvIdxLocal, i)]; }
 
     //! returns whether or not an scvf touches the boundary
     bool scvfTouchesBoundary(LocalIndexType ivLocalScvfIdx) const
     { return nodalIndexSet_.scvfTouchesBoundary(ivToNodeScvf_[ivLocalScvfIdx]); }
 
     //! returns the local indices of the neighboring scvs of an scvf
-    const LocalIndexContainer& localNeighboringScvIndices(LocalIndexType ivLocalScvfIdx) const
-    { return nodalIndexSet_.localNeighboringScvIndices(ivToNodeScvf_[ivLocalScvfIdx]); }
+    const LocalIndexContainer& neighboringLocalScvIndices(LocalIndexType ivLocalScvfIdx) const
+    { return scvfNeighborScvLocalIndices_[ivLocalScvfIdx]; }
 
     //! returns the number of faces in the interaction volume
     std::size_t numFaces() const
@@ -152,12 +166,31 @@ public:
     std::size_t numScvs() const
     { return nodalIndexSet_.numScvs(); }
 
+    //! returns the global scv indices connected to this dual grid node
+    const GlobalIndexContainer& globalScvIndices() const
+    { return nodalIndexSet_.globalScvIndices(); }
+
+    //! returns the global scvf indices connected to this dual grid node
+    const GlobalIndexContainer& globalScvfIndices() const
+    { return nodalIndexSet_.globalScvfIndices(); }
+
 private:
+    //! returns the local scv index to a given global scv index
+    unsigned int findLocalScvIdx_(GlobalIndexType globalScvIdx) const
+    {
+        auto it = std::find( nodalIndexSet_.globalScvIndices().begin(), nodalIndexSet_.globalScvIndices().end(), globalScvIdx );
+        assert(it != nodalIndexSet_.globalScvIndices().end() && "Global scv index not found in local container!");
+        return std::distance(nodalIndexSet_.globalScvIndices().begin(), it);
+    }
+
     const DualGridNodalIndexSet& nodalIndexSet_;
 
     std::size_t numFaces_;
     LocalIndexContainer ivToNodeScvf_;
     LocalIndexContainer nodeToIvScvf_;
+    //! maps to each scvf a list of neighbouring scv indices
+    //! ordering: 0 - inside scv idx; 1..n - outside scv indices
+    std::vector<LocalIndexContainer> scvfNeighborScvLocalIndices_;
 };
 } // end namespace
 

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

@@ -55,14 +55,14 @@ public:
     , globalScvIndex_(scv.dofIndex())
     , localDofIndex_(localIndex)
     {
-        const auto& nodeLocalScvfIndices = indexSet.nodalIndexSet().localScvfIndicesInScv(localIndex);
+        // center of the global scv
         const auto center = scv.center();
 
         // set up local basis
         LocalBasis localBasis;
         for (unsigned int coordIdx = 0; coordIdx < dim; ++coordIdx)
         {
-            const auto scvfIdx = indexSet.nodalIndexSet().scvfIdxGlobal(nodeLocalScvfIndices[coordIdx]);
+            const auto scvfIdx = indexSet.nodalIndexSet().scvfIdxGlobal(localDofIndex_, coordIdx);
             const auto& scvf = fvGeometry.scvf(scvfIdx);
             localBasis[coordIdx] = scvf.ipGlobal();
             localBasis[coordIdx] -= center;
@@ -78,7 +78,7 @@ public:
     LocalIndexType scvfIdxLocal(const unsigned int coordDir) const
     {
         assert(coordDir < dim);
-        return indexSet_.localScvfIndexInScv(localDofIndex_, coordDir);
+        return indexSet_.scvfIdxLocal(localDofIndex_, coordDir);
     }
 
     //! returns the index in the set of cell unknowns of the iv

dumux/discretization/fluxstencil.hh

@@ -119,9 +119,9 @@ public:
 
         // return the scv (element) indices in the interaction region
         if (fvGridGeometry.vertexUsesSecondaryInteractionVolume(scvf.vertexIndex()))
-            return fvGridGeometry.gridInteractionVolumeIndexSets().secondaryIndexSet(scvf).nodalIndexSet().globalScvIndices();
+            return fvGridGeometry.gridInteractionVolumeIndexSets().secondaryIndexSet(scvf).globalScvIndices();
         else
-            return fvGridGeometry.gridInteractionVolumeIndexSets().primaryIndexSet(scvf).nodalIndexSet().globalScvIndices();
+            return fvGridGeometry.gridInteractionVolumeIndexSets().primaryIndexSet(scvf).globalScvIndices();
     }
 };