From 6d2f5eaf2555be2f612779de4e2d9d61ec5251ef Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Thu, 8 Dec 2016 18:44:06 +0100
Subject: [PATCH] [staggeredGrid][fvGeometry] Add map (eIdx, localFaceIdx) ->
 scvf

*Add localFaceIdx in geometry helper
*Comment out code which is potentially superfluous (to be deleted later)
---
 .../staggered/fvelementgeometry.hh            |  5 +++
 .../staggered/globalfvgeometry.hh             | 21 ++++++++--
 .../staggered/staggeredgeometryhelper.hh      | 41 +++++++++++++++++++
 3 files changed, 64 insertions(+), 3 deletions(-)

diff --git a/dumux/discretization/staggered/fvelementgeometry.hh b/dumux/discretization/staggered/fvelementgeometry.hh
index 7e12e9fc43..5f7e1db8c5 100644
--- a/dumux/discretization/staggered/fvelementgeometry.hh
+++ b/dumux/discretization/staggered/fvelementgeometry.hh
@@ -83,6 +83,11 @@ public:
         return globalFvGeometry().scvf(scvfIdx);
     }
 
+    const SubControlVolumeFace& scvf(IndexType eIdx ,IndexType localScvfIdx) const
+    {
+        return globalFvGeometry().scvf(eIdx, localScvfIdx);
+    }
+
     //! iterator range for sub control volumes. Iterates over
     //! all scvs of the bound element (not including neighbor scvs)
     //! This is a free function found by means of ADL
diff --git a/dumux/discretization/staggered/globalfvgeometry.hh b/dumux/discretization/staggered/globalfvgeometry.hh
index d7e16024e5..f9218b4707 100644
--- a/dumux/discretization/staggered/globalfvgeometry.hh
+++ b/dumux/discretization/staggered/globalfvgeometry.hh
@@ -123,6 +123,7 @@ public:
         scvs_.resize(numScvs);
         scvfs_.reserve(numScvf);
         scvfIndicesOfScv_.resize(numScvs);
+        localToGlobalScvfIndices_.resize(numScvs);
 
         // Build the scvs and scv faces
         IndexType scvfIdx = 0;
@@ -131,10 +132,10 @@ public:
         {
             auto eIdx = problem.elementMapper().index(element);
 
-            // get the element geometry
-//             auto elementGeometry = element.geometry();
-//             GeometryHelper geometryHelper(elementGeometry);
 
+            // reserve memory for the localToGlobalScvfIdx map
+            auto numLocalFaces = element.subEntities(1);
+            localToGlobalScvfIndices_[eIdx].resize(numLocalFaces);
 
             scvs_[eIdx] = SubControlVolume(element.geometry(), eIdx);
 
@@ -158,6 +159,7 @@ public:
                                         std::vector<IndexType>({eIdx, nIdx}),
                                         geometryHelper
                                         );
+                    localToGlobalScvfIndices_[eIdx][intersection.indexInInside()] = scvfIdx;
                     scvfsIndexSet.push_back(scvfIdx++);
                 }
                 // boundary sub control volume faces
@@ -169,6 +171,7 @@ public:
                                         std::vector<IndexType>({eIdx, gridView_.size(0) + numBoundaryScvf_++}),
                                         geometryHelper
                                         );
+                    localToGlobalScvfIndices_[eIdx][intersection.indexInInside()] = scvfIdx;
                     scvfsIndexSet.push_back(scvfIdx++);
                 }
             }
@@ -206,6 +209,17 @@ public:
         return scvfIndicesOfScv_[scvIdx];
     }
 
+    const auto localToGlobalScvfIndex(IndexType eIdx, IndexType localScvfIdx) const
+    {
+        return localToGlobalScvfIndices_[eIdx][localScvfIdx];
+    }
+
+    const SubControlVolumeFace& scvf(IndexType eIdx ,IndexType localScvfIdx) const
+    {
+        return scvf(localToGlobalScvfIndex(eIdx, localScvfIdx));
+    }
+
+
 private:
     const Problem& problem_() const
     { return *problemPtr_; }
@@ -217,6 +231,7 @@ private:
     std::vector<SubControlVolume> scvs_;
     std::vector<SubControlVolumeFace> scvfs_;
     std::vector<std::vector<IndexType>> scvfIndicesOfScv_;
+    std::vector<std::vector<IndexType>> localToGlobalScvfIndices_;
     IndexType numBoundaryScvf_;
 };
 
diff --git a/dumux/discretization/staggered/staggeredgeometryhelper.hh b/dumux/discretization/staggered/staggeredgeometryhelper.hh
index 70e52a9730..de715aaf13 100644
--- a/dumux/discretization/staggered/staggeredgeometryhelper.hh
+++ b/dumux/discretization/staggered/staggeredgeometryhelper.hh
@@ -40,9 +40,12 @@ struct PairData
     int outerParallelFaceDofIdx;
     int outerParallelElementDofIdx;
     std::pair<int,int> normalPair;
+    int localNormalFaceIdx;
     int globalCommonEntIdx;
     Scalar parallelDistance;
     Scalar normalDistance;
+//     Scalar normalFaceArea;
+//     GlobalPosition normalFaceDirection;
     GlobalPosition virtualOuterParallelFaceDofPos;
 };
 
@@ -273,6 +276,7 @@ private:
         }
 
         asImp_().treatVirtualOuterParallelFaceDofs_();
+//         asImp_().setNormalFaceInformation_();
     }
 
 protected:
@@ -409,6 +413,9 @@ private:
         this->pairData_[1].normalPair.first = this->gridView_.indexSet().subIndex(this->intersection_.inside(), indices.normalLocalDofIdx2, dim-1) + this->offset_;
         this->pairData_[0].globalCommonEntIdx = this->gridView_.indexSet().subIndex(this->intersection_.inside(), indices.localCommonEntIdx1, codimCommonEntity);
         this->pairData_[1].globalCommonEntIdx = this->gridView_.indexSet().subIndex(this->intersection_.inside(), indices.localCommonEntIdx2, codimCommonEntity);
+
+        this->pairData_[0].localNormalFaceIdx = indices.normalLocalDofIdx1;
+        this->pairData_[1].localNormalFaceIdx = indices.normalLocalDofIdx2;
     }
 
      /*!
@@ -461,6 +468,40 @@ private:
         }
     }
 
+//     void setNormalFaceInformation_()
+//     {
+//         // get the local index of the facet we are dealing with in this class
+//         const int localIntersectionIdx = this->intersection_.indexInInside();
+//
+//         // iterate over all intersections of the element, this facet is part of
+//         for(const auto& is : intersections(this->gridView_, this->element_))
+//         {
+//             if(this->facetIsNormal_(localIntersectionIdx, is.indexInInside()))
+//             {
+//                 int index;
+//                 switch(localIntersectionIdx)
+//                 {
+//                     case 0:
+//                         index = (is.indexInInside() == 3) ? 0 : 1;
+//                         break;
+//                     case 1:
+//                         index = (is.indexInInside() == 2) ? 0 : 1;
+//                         break;
+//                     case 2:
+//                         index = (is.indexInInside() == 0) ? 0 : 1;
+//                         break;
+//                     case 3:
+//                         index = (is.indexInInside() == 1) ? 0 : 1;
+//                         break;
+//                     default:
+//                         DUNE_THROW(Dune::InvalidStateException, "Something went terribly wrong");
+//                 }
+//                 this->pairData_[index].normalFaceDirection = std::move(is.centerUnitOuterNormal());
+//                 this->pairData_[index].normalFaceArea = std::move(0.5 * is.geometry().volume());
+//             }
+//         }
+//     }
+
 };
 
 template<class GridView>
-- 
GitLab