diff --git a/dumux/common/fvproblem.hh b/dumux/common/fvproblem.hh
index a9857b41d0d03bea6798584faa01cd715af67d9e..dd8b8956a86594b09e94b5d4655f05172ac4836b 100644
--- a/dumux/common/fvproblem.hh
+++ b/dumux/common/fvproblem.hh
@@ -25,6 +25,7 @@
 #define DUMUX_FV_PROBLEM_HH
 
 #include <memory>
+#include <map>
 
 #include <dune/common/version.hh>
 #include <dune/common/fvector.hh>
@@ -75,6 +76,9 @@ class FVProblem
 
     static constexpr bool isBox = GET_PROP_VALUE(TypeTag, DiscretizationMethod) == DiscretizationMethods::Box;
 
+    using PointSourceMap = std::map<std::pair<std::size_t, std::size_t>,
+                                    std::vector<PointSource> >;
+
 public:
     /*!
      * \brief Constructor
@@ -418,7 +422,7 @@ public:
                 // call convienience problem interface function
                 asImp_().pointSource(pointSource, element, fvGeometry, elemVolVars, scv);
                 // at last take care about multiplying with the correct volume
-                pointSource /= volume;
+                pointSource /= volume*pointSource.embeddings();
                 // add the point source values to the local residual
                 source += pointSource.values();
             }
@@ -431,7 +435,7 @@ public:
      * \brief Compute the point source map, i.e. which scvs have point source contributions
      * \note Call this on the problem before assembly if you want to enable point sources set
      *       via the addPointSources member function.
-    */
+     */
     void computePointSourceMap()
     {
         // clear the given point source maps in case it's not empty
@@ -451,10 +455,16 @@ public:
         }
     }
 
+    /*!
+     * \brief Get the point source map. It stores the point sources per scv
+     */
+    const PointSourceMap& pointSourceMap() const
+    { return pointSourceMap_; }
+
     /*!
      * \brief Applies the initial solution for all degrees of freedom of the grid.
      * \param sol the initial solution vector
-    */
+     */
     void applyInitialSolution(SolutionVector& sol) const
     {
         // set the initial values by forwarding to a specialized method
@@ -517,8 +527,8 @@ public:
      */
     Scalar extrusionFactorAtPos(const GlobalPosition &globalPos) const
     {
-        //! As a default, i.e. if the user's problem does not overload any extrusion factor method
-        //! return 1.0
+        // As a default, i.e. if the user's problem does not overload
+        // any extrusion factor method, return 1.0
         return 1.0;
     }
 
@@ -569,8 +579,7 @@ private:
     std::string problemName_;
 
     //! A map from an scv to a vector of point sources
-    std::map<std::pair<unsigned int, unsigned int>,
-             std::vector<PointSource> > pointSourceMap_;
+    PointSourceMap pointSourceMap_;
 };
 
 } // end namespace Dumux
diff --git a/dumux/common/pointsource.hh b/dumux/common/pointsource.hh
index ac4efb84f9b836800c49ba0fd715fffa4d894384..271228fe8e1c6a34e92ca007e659d641025e9216 100644
--- a/dumux/common/pointsource.hh
+++ b/dumux/common/pointsource.hh
@@ -115,14 +115,8 @@ public:
     }
 
     //! return the source values
-    // don't forget to call this when it's overloaded
-    // in the derived class
     PrimaryVariables values() const
-    {
-        auto values = PrimaryVariables(values_);
-        values /= embeddings_;
-        return values;
-    }
+    { return values_; }
 
     //! return the source position
     const GlobalPosition& position() const
@@ -144,7 +138,17 @@ public:
         embeddings_ = embeddings;
     }
 
-    //! get the number of embeddings for this point source
+    /*!
+     * \brief get the number of embeddings for this point source
+     * \note A point source might be located on the intersection between several scvs.
+     *       If so, there are point sources for every neighboring scv with the same position.
+     *       `embeddings` returns the number of neighboring scvs.
+     * Example: If I want to inject 1kg/s at a location that is on the inner face of an scv
+     *          the point source exists in both scvs. Both have a value of 1kg/s.
+     *          We then divide the value by the number of embeddings to not inject 2kg/s but 1kg/s.
+     * \note This division is done in the problem.scvPointSources() if this behaviour is not explicitly
+     *       changed by e.g. overloading this function in the problem implementation.
+     */
     std::size_t embeddings() const
     {
         return embeddings_;
@@ -286,9 +290,10 @@ class BoundingBoxTreePointSourceHelper
 
 public:
     //! calculate a DOF index to point source map from given vector of point sources
+    template<class PointSourceMap>
     static void computePointSourceMap(const FVGridGeometry& fvGridGeometry,
                                       std::vector<PointSource>& sources,
-                                      std::map<std::pair<unsigned int, unsigned int>, std::vector<PointSource> >& pointSourceMap)
+                                      PointSourceMap& pointSourceMap)
     {
         const auto& boundingBoxTree = fvGridGeometry.boundingBoxTree();