From da908103dd4ce00bec89d3b15747a5e08719331a Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Thu, 12 Nov 2015 11:47:23 +0100
Subject: [PATCH] [pointsource] Map from element and scv to pointsource

The mapping from dofIdx to pointsource was inconvenient since the
source functions in dumux are called for each element for each
scv. Now the map maps from a elementIdx scvIdx pair to the source
belonging to the scv. DofSource is thus not the right name anymore
and is now called scvPointSources.
---
 dumux/common/pointsource.hh              | 29 +++++++++++++-----------
 dumux/implicit/1p/1plocalresidual.hh     | 10 ++++----
 dumux/implicit/common/implicitproblem.hh | 20 ++++++++--------
 3 files changed, 31 insertions(+), 28 deletions(-)

diff --git a/dumux/common/pointsource.hh b/dumux/common/pointsource.hh
index e176831ee4..203c2620e9 100644
--- a/dumux/common/pointsource.hh
+++ b/dumux/common/pointsource.hh
@@ -101,7 +101,7 @@ public:
     static void computePointSourceMap(const Problem& problem,
                                       const std::shared_ptr<BoundingBoxTree>& boundingBoxTree,
                                       std::vector<PointSource>& sources,
-                                      std::map<unsigned int, std::vector<PointSource> >& pointSourceMap)
+                                      std::map<std::pair<unsigned int, unsigned int>, std::vector<PointSource> >& pointSourceMap)
     {
         for (auto&& source : sources)
         {
@@ -120,31 +120,34 @@ public:
                     fvGeometry.update(problem.gridView(), element);
                     const auto globalPos = source.position();
                     // loop over all sub control volumes and check if the point source is inside
-                    std::vector<unsigned int> vertices;
+                    std::vector<unsigned int> scvs;
                     for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
                     {
                         auto geometry = fvGeometry.subContVolGeometries[scvIdx];
                         if (BoundingBoxTreeHelper<dimworld>::pointInGeometry(geometry, globalPos))
-                            vertices.push_back(problem.model().dofMapper().subIndex(element, scvIdx, dofCodim));
+                            scvs.push_back(scvIdx);
                     }
-                    for (unsigned int vIdx : vertices)
+                    // for all scvs that where tested positiv add the point sources
+                    // to the element/scv to point source map
+                    for (unsigned int scvIdx : scvs)
                     {
-                        // add the pointsource to the DOF map
-                        if (pointSourceMap.count(vIdx))
-                            pointSourceMap.at(vIdx).push_back(source);
+                        const auto key = std::make_pair(eIdx, scvIdx);
+                        if (pointSourceMap.count(key))
+                            pointSourceMap.at(key).push_back(source);
                         else
-                            pointSourceMap.insert({vIdx, {source}});
-                        // split equally on the number of vertices
-                        pointSourceMap.at(vIdx).back() /= vertices.size();
+                            pointSourceMap.insert({key, {source}});
+                        // split equally on the number of matched scvs
+                        pointSourceMap.at(key).back() /= scvs.size();
                     }
                 }
                 else
                 {
                     // add the pointsource to the DOF map
-                    if (pointSourceMap.count(eIdx))
-                        pointSourceMap.at(eIdx).push_back(source);
+                    const auto key = std::make_pair(eIdx, /*scvIdx=*/ 0);
+                    if (pointSourceMap.count(key))
+                        pointSourceMap.at(key).push_back(source);
                     else
-                        pointSourceMap.insert({eIdx, {source}});
+                        pointSourceMap.insert({key, {source}});
                 }
             }
         }
diff --git a/dumux/implicit/1p/1plocalresidual.hh b/dumux/implicit/1p/1plocalresidual.hh
index b9bf5c915d..4e82d27fe7 100644
--- a/dumux/implicit/1p/1plocalresidual.hh
+++ b/dumux/implicit/1p/1plocalresidual.hh
@@ -161,11 +161,11 @@ public:
                                      this->curVolVars_());
 
         // add contribution from possible point sources
-        this->problem_().dofSources(source,
-                                    this->element_(),
-                                    this->fvGeometry_(),
-                                    scvIdx,
-                                    this->curVolVars_());
+        this->problem_().scvPointSources(source,
+                                         this->element_(),
+                                         this->fvGeometry_(),
+                                         scvIdx,
+                                         this->curVolVars_());
     }
 
     /*!
diff --git a/dumux/implicit/common/implicitproblem.hh b/dumux/implicit/common/implicitproblem.hh
index b53797d3d7..5dbc6033b0 100644
--- a/dumux/implicit/common/implicitproblem.hh
+++ b/dumux/implicit/common/implicitproblem.hh
@@ -951,20 +951,20 @@ public:
     }
 
     /*!
-     * \brief Adds contribution of point sources for a specific DOF
+     * \brief Adds contribution of point sources for a specific sub control volume
      *        to the values.
      */
-    void dofSources(PrimaryVariables &values,
-                    const Element &element,
-                    const FVElementGeometry &fvGeometry,
-                    const int scvIdx,
-                    const ElementVolumeVariables &elemVolVars) const
+    void scvPointSources(PrimaryVariables &values,
+                         const Element &element,
+                         const FVElementGeometry &fvGeometry,
+                         const int scvIdx,
+                         const ElementVolumeVariables &elemVolVars) const
     {
-        unsigned int dofGlobalIdx = model_.dofMapper().subIndex(element, scvIdx, dofCodim);
-        if (pointSourceMap_.count(dofGlobalIdx))
+        auto key = std::make_pair(this->elementMapper().index(element), scvIdx);
+        if (pointSourceMap_.count(key))
         {
             // call the solDependent function. Herein the user might fill/add values to the point sources
-            auto pointSources = pointSourceMap_.at(dofGlobalIdx);
+            auto pointSources = pointSourceMap_.at(key);
             asImp_().solDependentPointSources(pointSources, element, fvGeometry, scvIdx, elemVolVars);
 
             // Add the contributions to the dof source values
@@ -1059,7 +1059,7 @@ private:
     std::shared_ptr<GridAdaptModel> gridAdapt_;
 
     std::shared_ptr<BoundingBoxTree> boundingBoxTree_;
-    std::map<unsigned int, std::vector<PointSource> > pointSourceMap_;
+    std::map<std::pair<unsigned int, unsigned int>, std::vector<PointSource> > pointSourceMap_;
 };
 } // namespace Dumux
 
-- 
GitLab