diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dpressure2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dpressure2p.hh
index 8f6f2bad4f7549fceb575a367a08c10e200c1892..4dbeb4dd067fc3374dd11b4032ce3caca0cb6f2e 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dpressure2p.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dpressure2p.hh
@@ -41,12 +41,26 @@ namespace Dumux
  *  At Dirichlet boundaries a two-point flux approximation is used.
  * \f[ \Phi = g \;  \text{on} \; \Gamma_1, \quad \text{and} \quad
  * -\text{div}\, \boldsymbol v_t \cdot \mathbf{n} = J \;  \text{on}  \; \Gamma_2. \f]
+ *
  *  Here, \f$ \Phi_\alpha \f$ denotes the potential of phase \f$ \alpha \f$, \f$ \boldsymbol K \f$ the intrinsic permeability,
  * \f$ \lambda_t \f$ the total mobility, \f$ f_\alpha \f$ the phase fractional flow function.
  *
- * More details on the equations can be found...
+ * More details on the equations can be found in
+ *
+ * Wolff 2013: http://elib.uni-stuttgart.de/opus/volltexte/2013/8661/
+ *
+ * M. Wolff, Y. Cao, B. Flemisch, R. Helmig, and B. Wohlmuth (2013a). Multi-point flux
+ * approximation L-method in 3D: numerical convergence and application to two-phase
+ * flow through porous media. In P. Bastian, J. Kraus, R. Scheichl, and M. Wheeler,
+ * editors, Simulation of Flow in Porous Media - Applications in Energy and Environment. De Gruyter.
+ *
+ * M. Wolff, B. Flemisch, R. Helmig, I. Aavatsmark.
+ * Treatment of tensorial relative permeabilities with multipoint flux approximation.
+ * International Journal of Numerical Analysis and Modeling (9), pp. 725-744, 2012.
  *
- *  Remark: only for 3-d hexahedral grids
+ *  Remark1: only for 2-D quadrilateral grid
+ *
+ *  Remark2: implemented for UGGrid, ALUGrid, or SGrid/YaspGrid
  *
  *\tparam TypeTag The problem Type Tag
  */
@@ -464,12 +478,15 @@ private:
      * this errors due to wrong time-stepping without losing mass conservation. The error term looks as follows:
      * \f[
      *  q_{error} = \begin{cases}
-     *          S < 0 & a_{error} \frac{S}{\Delta t} V \\
-     *          S > 1 & a_{error} \frac{(S - 1)}{\Delta t} V \\
+     *          S < 0 &  \frac{S}{\Delta t} V \\
+     *          S > 1 &  \frac{(S - 1)}{\Delta t} V \\
      *          0 \le S \le 1 & 0
      *      \end{cases}
      *  \f]
-     *  where \f$a_{error}\f$ is a weighting factor (default: \f$a_{error} = 0.5\f$)
+     *
+     *  \param cellData The IMPES <tt>CellData</tt> object of the current cell.
+     *
+     *  \return The scalar value of the error term.
     */
     Scalar evaluateErrorTerm_(CellData& cellData)
     {
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dpressure2padaptive.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dpressure2padaptive.hh
index 4bf579f4cb8783ee162ef169e1b1d315e93b86a6..c368f6f8598c46005854dbba825deb6106a7f917 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dpressure2padaptive.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dpressure2padaptive.hh
@@ -44,12 +44,23 @@ namespace Dumux
  *  Here, \f$ \Phi_\alpha \f$ denotes the potential of phase \f$ \alpha \f$, \f$ \boldsymbol K \f$ the intrinsic permeability,
  * \f$ \lambda_t \f$ the total mobility, \f$ f_\alpha \f$ the phase fractional flow function.
  *
- * More details on the equations can be found in H. Hoteit, A. Firoozabadi.
- * Numerical modeling of thwo-phase flow in heterogeneous permeable media with different capillarity pressures. Adv Water Res (31), 2008
+ * More details on the equations can be found in
+ *
+ * Wolff 2013: http://elib.uni-stuttgart.de/opus/volltexte/2013/8661/
+ *
+ * M. Wolff, Y. Cao, B. Flemisch, R. Helmig, and B. Wohlmuth (2013a). Multi-point flux
+ * approximation L-method in 3D: numerical convergence and application to two-phase
+ * flow through porous media. In P. Bastian, J. Kraus, R. Scheichl, and M. Wheeler,
+ * editors, Simulation of Flow in Porous Media - Applications in Energy and Environment. De Gruyter.
+ *
+ * M. Wolff, B. Flemisch, R. Helmig, I. Aavatsmark.
+ * Treatment of tensorial relative permeabilities with multipoint flux approximation.
+ * International Journal of Numerical Analysis and Modeling (9), pp. 725-744, 2012.
+ *
  *
  *  Remark1: only for 2-D quadrilateral grid
  *
- *  Remark2: implemented for UGGrid, ALUGrid, or SGrid/YaspGrid
+ *  Remark2: implemented for UGGrid, ALUGrid
  *
  *  Remark3: Allowed difference in grid levels of two neighboring cells: 1
  *
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dpressurevelocity2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dpressurevelocity2p.hh
index 66522762147bc4fea68081147c522ca498a6468f..f53abb72ec6f22d54638c4e2535a99c3027c95a2 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dpressurevelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dpressurevelocity2p.hh
@@ -30,18 +30,10 @@ namespace Dumux
 {
 
 //! \ingroup FVPressure2p
-/*! \brief Determines the velocity from a finite volume solution of the  pressure equation of a sequential model (IMPES).
- * Calculates phase velocities or total velocity from a known pressure field applying a finite volume discretization and a MPFA L-method.
- * At Dirichlet boundaries a two-point flux approximation is used.
- * The pressure has to be given as piecewise constant cell values.
- * The velocities are calculated as
+/*! \brief Class for the calculation of velocities from the  pressure solution of an IMPES scheme using a MPFA L-method.
  *
- *\f[ \boldsymbol v_\alpha = - \lambda_\alpha \boldsymbol K \textbf{grad}\, \Phi_\alpha, \f]
- * and,
- * \f[ \boldsymbol v_t = \boldsymbol v_w + \boldsymbol v_n,\f]
- *
- * where \f$ \Phi_\alpha \f$ denotes the potential of phase \f$ \alpha \f$, \f$ \boldsymbol K \f$ the intrinsic permeability,
- * and \f$ \lambda_\alpha \f$ a phase mobility.
+ * Can be used for calculating the complete velocity field before the solution of the transport equation (more efficient),
+ * or for face-wise velocity calculation directly in the transport solution (less efficient).
  *
  * Remark1: only for 2-D quadrilateral grids!
  *
@@ -133,6 +125,7 @@ public:
     void calculateVelocity(const Intersection&, CellData&);
     void calculateVelocityOnBoundary(const Intersection& intersection, CellData& cellData);
 
+    //! Function for updating the velocity field if iterations are necessary in the transport solution
     void updateVelocity()
     {
         this->updateMaterialLaws();
@@ -220,9 +213,10 @@ private:
 };
 // end of template
 
-/*! \brief Calculates the velocities at a cell-cell interfaces.
+
+/*! \brief Calculates the velocities at a cell-cell interfaces for the entire simulation grid.
  *
- * Calculates the velocities at a cell-cell interfaces from a given pressure field.
+ * Calculates the velocities at a cell-cell interfaces from a given pressure field for the entire simulation grid.
  *
  */
 template<class TypeTag>
@@ -295,7 +289,13 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocity()
     return;
 } // end method calcTotalVelocity
 
-
+/*! \brief Calculates the velocity at a cell-cell interface.
+ *
+ * Calculates the velocity at a cell-cell interface from a given pressure field.
+ *
+ * \param intersection Intersection of two grid cells
+ * \param cellData Object containing all model relevant cell data
+ */
 template<class TypeTag>
 void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocity(const Intersection& intersection, CellData& cellData)
 {
@@ -371,6 +371,13 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocity(const Intersection&
     cellDataJ.fluxData().setVelocityMarker(indexInOutside);
 }
 
+/*! \brief Calculates the velocity at a boundary.
+ *
+ * Calculates the velocity at a boundary from a given pressure field.
+ *
+ * \param intersection Boundary intersection
+ * \param cellData Object containing all model relevant cell data
+ */
 template<class TypeTag>
 void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocityOnBoundary(const Intersection& intersection, CellData& cellData)
 {
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dpressurevelocity2padaptive.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dpressurevelocity2padaptive.hh
index ebb22295398e45f2174687217e5d281959379eb6..9751b5c91bf84c2f7f6fa5421ff3172555fb67ea 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dpressurevelocity2padaptive.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dpressurevelocity2padaptive.hh
@@ -30,22 +30,14 @@ namespace Dumux
 {
 
 //! \ingroup FVPressure2p
-/*! \brief Determines the velocity from a grid adaptive finite volume solution of the  pressure equation of a sequential model (IMPES).
- * Calculates phase velocities or total velocity from a known pressure field applying a grid adaptive finite volume discretization and a MPFA L-method.
- * At Dirichlet boundaries a two-point flux approximation is used.
- * The pressure has to be given as piecewise constant cell values.
- * The velocities are calculated as
+/*! \brief Class for the calculation of velocities from the  pressure solution of an IMPES scheme using a grid adaptive MPFA L-method.
  *
- *\f[ \boldsymbol v_\alpha = - \lambda_\alpha \boldsymbol K \textbf{grad}\, \Phi_\alpha, \f]
- * and,
- * \f[ \boldsymbol v_t = \boldsymbol v_w + \boldsymbol v_n,\f]
- *
- * where \f$ \Phi_\alpha \f$ denotes the potential of phase \f$ \alpha \f$, \f$ \boldsymbol K \f$ the intrinsic permeability,
- * and \f$ \lambda_\alpha \f$ a phase mobility.
+ * Can be used for calculating the complete velocity field before the solution of the transport equation (more efficient),
+ * or for face-wise velocity calculation directly in the transport solution (less efficient).
  *
  * Remark1: only for 2-D quadrilateral grids!
  *
- * Remark2: can use UGGrid, ALUGrid or SGrid/YaspGrid!
+ * Remark2: can use UGGrid, ALUGrid
  *
  * Remark3: Allowed difference in grid levels of two neighboring cells: 1
  *
@@ -125,7 +117,7 @@ template<class TypeTag> class FvMpfaL2dPressureVelocity2pAdaptive: public FvMpfa
     typedef Dune::FieldVector<Scalar, dim> DimVector;
 
 public:
-    //! Constructs a FVMPFAO2pFABoundVelocity2p object
+    //! Constructs a FvMpfaL2dPressureVelocity2pAdaptive object
     /*!
      * \param problem A problem class object
      */
@@ -147,6 +139,7 @@ public:
     void calculateVelocity(const Intersection&, CellData&);
     void calculateVelocityOnBoundary(const Intersection& intersection, CellData& cellData);
 
+    //! Function for updating the velocity field if iterations are necessary in the transport solution
     void updateVelocity()
     {
         this->updateMaterialLaws();
@@ -234,9 +227,9 @@ private:
 };
 // end of template
 
-/*! \brief Calculates the velocities at a cell-cell interfaces.
+/*! \brief Calculates the velocities at a cell-cell interfaces for the entire simulation grid.
  *
- * Calculates the velocities at a cell-cell interfaces from a given pressure field.
+ * Calculates the velocities at a cell-cell interfaces from a given pressure field for the entire simulation grid.
  *
  */
 template<class TypeTag>
@@ -334,6 +327,13 @@ void FvMpfaL2dPressureVelocity2pAdaptive<TypeTag>::calculateVelocity()
     return;
 } // end method calcTotalVelocity
 
+/*! \brief Calculates the velocity at a cell-cell interface.
+ *
+ * Calculates the velocity at a cell-cell interface from a given pressure field.
+ *
+ * \param intersection Intersection of two grid cells
+ * \param cellData Object containing all model relevant cell data
+ */
 template<class TypeTag>
 void FvMpfaL2dPressureVelocity2pAdaptive<TypeTag>::calculateVelocity(const Intersection& intersection, CellData& cellData)
 {
@@ -502,6 +502,13 @@ void FvMpfaL2dPressureVelocity2pAdaptive<TypeTag>::calculateVelocity(const Inter
     }
 }
 
+/*! \brief Calculates the velocity at a boundary.
+ *
+ * Calculates the velocity at a boundary from a given pressure field.
+ *
+ * \param intersection Boundary intersection
+ * \param cellData Object containing all model relevant cell data
+ */
 template<class TypeTag>
 void FvMpfaL2dPressureVelocity2pAdaptive<TypeTag>::calculateVelocityOnBoundary(const Intersection& intersection, CellData& cellData)
 {
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dtransmissibilitycalculator.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dtransmissibilitycalculator.hh
index 70bbffb27f14eae738f34c2b5ecb96a7e9c655c3..fb80cc97c86571fdc14f8c7f33ffcaaad8967776 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dtransmissibilitycalculator.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dtransmissibilitycalculator.hh
@@ -25,13 +25,19 @@
 
 /**
  * @file
- * @brief Provides methods for transmissibility calculation
+ * @brief Provides methods for transmissibility calculation 2-d
  */
 
 namespace Dumux
 {
 //! \ingroup FVPressure2p
-/*! Provides methods for transmissibility calculation.
+/*! \brief Provides methods for transmissibility calculation in 2-d.
+ *
+ *  The transmissibilities are calculated using the MPFA L-method.
+ *
+ *  Aavatsmark et al. A compact multipoint flux calculation method with improved robustness.
+ *  Numerical Methods for Partial Differential Equations 24. 2008
+ *
  */
 template<class TypeTag>
 class FvMpfaL2dTransmissibilityCalculator
@@ -56,13 +62,14 @@ class FvMpfaL2dTransmissibilityCalculator
     typedef typename Dumux::FVMPFALInteractionVolume<TypeTag> InteractionVolume;
 
 public:
-    typedef Dune::FieldMatrix<Scalar, dim, 2 * dim - dim + 1> TransmissibilityType;
+    typedef Dune::FieldMatrix<Scalar, dim, 2 * dim - dim + 1> TransmissibilityType;//!< Type of the transmissibility matrix
 
+    //! return values for the transmissibility functions
     enum
     {
-        leftTriangle = -1,
-        noTransmissibility = 0,
-        rightTriangle = 1
+        leftTriangle = -1,//!< Left L-shape
+        noTransmissibility = 0,//!< No transmissibility calculated
+        rightTriangle = 1 //!< Right L-shape
     };
 
     // Calculates tranmissibility matrix
@@ -84,7 +91,10 @@ public:
     std::vector<DimVector >& lambda,
     int idx1, int idx2, int idx3);
 
-
+    //! Constructs a FvMpfaL2dTransmissibilityCalculator object
+    /**
+     * \param problem A problem class object
+     */
     FvMpfaL2dTransmissibilityCalculator(Problem& problem) :
         problem_(problem), R_(0)
     {
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dvelocity2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dvelocity2p.hh
index d6241a583c0c900f98cdc83b57cf0e41e49fc75e..3874d6737deeb0bdf3ed1ffda493427b554edd50 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dvelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dvelocity2p.hh
@@ -33,6 +33,24 @@
 
 namespace Dumux
 {
+//! \ingroup FVPressure2p
+/*! \brief Class for calculating 2-d velocities from cell-wise constant pressure values.
+ * Calculates phase velocities or total velocity from a known pressure field applying a finite volume discretization and a MPFA L-method.
+ * At Dirichlet boundaries a two-point flux approximation is used.
+ * The pressure has to be given as piecewise constant cell values.
+ * The velocities are calculated as
+ *
+ *\f[ \boldsymbol v_\alpha = - \lambda_\alpha \boldsymbol K \textbf{grad}\, \Phi_\alpha, \f]
+ * and,
+ * \f[ \boldsymbol v_t = \boldsymbol v_w + \boldsymbol v_n,\f]
+ *
+ * where \f$ \Phi_\alpha \f$ denotes the potential of phase \f$ \alpha \f$, \f$ \boldsymbol K \f$ the intrinsic permeability,
+ * and \f$ \lambda_\alpha \f$ a phase mobility.
+ *
+ * Remark1: only for quadrilaterals!
+ *
+ * \tparam TypeTag The problem Type Tag
+ */
 template<class TypeTag> class FvMpfaL2dVelocity2p
 {
     typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
@@ -114,6 +132,10 @@ template<class TypeTag> class FvMpfaL2dVelocity2p
     typedef Dune::FieldVector<Scalar, dim> DimVector;
 
 public:
+    //! Constructs a FvMpfaL2dVelocity2p object
+    /*!
+     * \param problem A problem class object
+     */
     FvMpfaL2dVelocity2p(Problem& problem) :
         problem_(problem), transmissibilityCalculator_(problem), gravity_(problem.gravity())
     {
@@ -125,9 +147,11 @@ public:
         vtkOutputLevel_ = GET_PARAM_FROM_GROUP(TypeTag, int, Vtk, OutputLevel);
     }
 
+    //calculate velocities for flux faces of an interaction volume
     void calculateInnerInteractionVolumeVelocity(InteractionVolume& interactionVolume, CellData& cellData1, CellData& cellData2, CellData& cellData3, CellData& cellData4, InnerBoundaryVolumeFaces& innerBoundaryVolumeFaces);
     void calculateBoundaryInteractionVolumeVelocity(InteractionVolume& interactionVolume, CellData& cellData, int elemIdx);
 
+    //!Initializes the velocity model
     void initialize()
     {
         ElementIterator element = problem_.gridView().template begin<0>();
@@ -145,8 +169,16 @@ public:
         return;
     }
 
-    //! \brief Write data files
-    /*  \param name file name */
+    /*! \brief Adds velocity output to the output file
+     *
+     * Adds the phase velocities or a total velocity (depending on the formulation) to the output.
+     * If the VtkOutputLevel is equal to zero (default) only primary variables are written,
+     * if it is larger than zero also secondary variables are written.
+     *
+     * \tparam MultiWriter Class defining the output writer
+     * \param writer The output writer (usually a <tt>VTKMultiWriter</tt> object)
+     *
+     */
     template<class MultiWriter>
     void addOutputVtkFields(MultiWriter &writer)
     {
@@ -235,6 +267,17 @@ protected:
 };
 // end of template
 
+/*! \brief Calculates the velocities at the flux faces of an interation volume around a vertex which is not a boundary vertex.
+ *
+ *  Calculates the velocities at the flux faces of an interation volume around a vertex which is not a boundary vertex and adds them to the face velocity vectors in the <tt>CellData</tt> objects.
+ *
+ * \param interactionVolume An <tt>InteractionVolume</tt> object including the information for calculating the MPFA transmissibilities
+ * \param cellData1  <tt>CellData</tt> object of an IMPES model for sub-volume 1
+ * \param cellData2  <tt>CellData</tt> object of an IMPES model for sub-volume 2
+ * \param cellData3  <tt>CellData</tt> object of an IMPES model for sub-volume 3
+ * \param cellData4  <tt>CellData</tt> object of an IMPES model for sub-volume 4
+ * \param innerBoundaryVolumeFaces container including information about faces intersecting a boundary
+ */
 template<class TypeTag>
 void FvMpfaL2dVelocity2p<TypeTag>::calculateInnerInteractionVolumeVelocity(InteractionVolume& interactionVolume, CellData& cellData1, CellData& cellData2, CellData& cellData3, CellData& cellData4, InnerBoundaryVolumeFaces& innerBoundaryVolumeFaces)
 {
@@ -677,6 +720,14 @@ void FvMpfaL2dVelocity2p<TypeTag>::calculateInnerInteractionVolumeVelocity(Inter
     cellData4.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(3, 1));
 }
 
+/*! \brief Calculates the velocity at a boundary flux faces.
+ *
+ *  Calculates the velocity at a boundary flux face and adds it to the face velocity vector in the <tt>CellData</tt> object.
+ *
+ * \param interactionVolume An <tt>InteractionVolume</tt> object including the information for calculating the MPFA transmissibilities
+ * \param cellData  <tt>CellData</tt> object of an IMPES model for the sub-volume including the boundary face
+ * \param elemIdx local sub-volume index
+ */
 template<class TypeTag>
 void FvMpfaL2dVelocity2p<TypeTag>::calculateBoundaryInteractionVolumeVelocity(InteractionVolume& interactionVolume, CellData& cellData, int elemIdx)
 {
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dvelocity2padaptive.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dvelocity2padaptive.hh
index 7f3b4fc26b3e055b83b40603b56eeed94721d70a..82090e3346a570a497f7814a3bb8e56b120f453d 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dvelocity2padaptive.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dvelocity2padaptive.hh
@@ -24,11 +24,29 @@
 
 /**
  * @file
- * @brief  Velocity calculation using a 2-d MPFA L-method
+ * @brief  Velocity calculation on adaptive grids using a 2-d MPFA L-method
  */
 
 namespace Dumux
 {
+//! \ingroup FVPressure2p
+/*! \brief Class for calculating 2-d velocities from cell-wise constant pressure values.
+ * Calculates phase velocities or total velocity from a known pressure field applying a finite volume discretization and a MPFA L-method.
+ * At Dirichlet boundaries a two-point flux approximation is used.
+ * The pressure has to be given as piecewise constant cell values.
+ * The velocities are calculated as
+ *
+ *\f[ \boldsymbol v_\alpha = - \lambda_\alpha \boldsymbol K \textbf{grad}\, \Phi_\alpha, \f]
+ * and,
+ * \f[ \boldsymbol v_t = \boldsymbol v_w + \boldsymbol v_n,\f]
+ *
+ * where \f$ \Phi_\alpha \f$ denotes the potential of phase \f$ \alpha \f$, \f$ \boldsymbol K \f$ the intrinsic permeability,
+ * and \f$ \lambda_\alpha \f$ a phase mobility.
+ *
+ * Remark1: only for quadrilaterals!
+ *
+ * \tparam TypeTag The problem Type Tag
+ */
 template<class TypeTag> class FvMpfaL2dVelocity2pAdaptive : public FvMpfaL2dVelocity2p<TypeTag>
 {
     typedef FvMpfaL2dVelocity2p<TypeTag> ParentType;
@@ -111,10 +129,15 @@ template<class TypeTag> class FvMpfaL2dVelocity2pAdaptive : public FvMpfaL2dVelo
     typedef Dune::FieldVector<Scalar, dim> DimVector;
 
 public:
+    //! Constructs a FvMpfaL2dVelocity2pAdaptive object
+    /*!
+     * \param problem A problem class object
+     */
     FvMpfaL2dVelocity2pAdaptive(Problem& problem) :
         ParentType(problem), problem_(problem)
     {}
 
+    //calculate velocities for flux faces of a hanging node interaction volume
     void calculateHangingNodeInteractionVolumeVelocity(InteractionVolume& interactionVolume, CellData& cellData1, CellData& cellData2, CellData& cellData4, InnerBoundaryVolumeFaces& innerBoundaryVolumeFaces);
 
 private:
@@ -122,6 +145,16 @@ private:
 };
 // end of template
 
+/*! \brief Calculates the velocities at the flux faces of an interation volume around a hanging node vertex.
+ *
+ *  Calculates the velocities at the flux faces of an interation volume around a hanging node vertex and adds them to the face velocity vectors in the <tt>CellData</tt> objects.
+ *
+ * \param interactionVolume An <tt>InteractionVolume</tt> object including the information for calculating the MPFA transmissibilities
+ * \param cellData1  <tt>CellData</tt> object of an IMPES model for sub-volume 1
+ * \param cellData2  <tt>CellData</tt> object of an IMPES model for sub-volume 2
+ * \param cellData4  <tt>CellData</tt> object of an IMPES model for sub-volume 4
+ * \param innerBoundaryVolumeFaces container including information about faces intersecting a boundary
+ */
 template<class TypeTag>
 void FvMpfaL2dVelocity2pAdaptive<TypeTag>::calculateHangingNodeInteractionVolumeVelocity(InteractionVolume& interactionVolume, CellData& cellData1, CellData& cellData2, CellData& cellData4, InnerBoundaryVolumeFaces& innerBoundaryVolumeFaces)
 {
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dinteractionvolumecontainer.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dinteractionvolumecontainer.hh
index 98dfdac7b60547b6070e7441d1743d54e85fad72..20fadc1e542c5c144050f78a9d417b561af80e03 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dinteractionvolumecontainer.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dinteractionvolumecontainer.hh
@@ -38,7 +38,14 @@ bool sort_compare(const std::vector<int>& entryI, const std::vector<int>& entryJ
 }
 
 //! \ingroup FVPressure2p
-/*! Interactionvolume container for 3-d MPFA L-method
+/*! \brief Interactionvolume container for 3-d MPFA L-method
+ *
+ * Container class which stores MPFA-interaction-volume information for each vertex of a DUNE grid.
+ * Each <tt>InteractionVolume</tt> object stores the information which is necessary to calculate MPFA transmissibility matrices:
+ *
+ * - relationship and orientation of the elements around a vertex (see doc/docextra/3dmpfa)
+ * - geometric information, such as element/face/edge positions, normals, ...
+ *
  */
 template<class TypeTag>
 class FvMpfaL3dInteractionVolumeContainer
@@ -98,7 +105,7 @@ class FvMpfaL3dInteractionVolumeContainer
         };
 
 public:
-    typedef typename GET_PROP_TYPE(TypeTag, MPFAInteractionVolume) InteractionVolume;
+    typedef typename GET_PROP_TYPE(TypeTag, MPFAInteractionVolume) InteractionVolume; //!< Type for storing an MPFA-interaction-volume. (Usually of type Dumux::FvMpfaL3dInteractionVolume or Dumux::FvMpfaL3dInteractionVolumeAdaptive)
 
 private:
     typedef std::vector<InteractionVolume> GlobalInteractionVolumeVector;
@@ -112,6 +119,10 @@ private:
     void storeInteractionVolumeInfo();
 public:
 
+    /*! \brief Updates the interaction volume container
+     *
+     * Rebuilds and stores the interaction volumes for the entire grid
+     */
     void update()
     {
         interactionVolumes_.clear();
@@ -123,39 +134,58 @@ public:
         asImp_().storeInteractionVolumeInfo();
     }
 
+
+    /*! \brief Initializes the interaction volume container
+     *
+     * Builds and stores the interaction volumes for the entire grid
+     */
     void initialize(bool solveTwice = true)
     {
-        interactionVolumes_.clear();
-        realFluxFaceArea_.clear();
-
-        realFluxFaceArea_.resize(problem_.gridView().size(dim), Dune::FieldVector<Dune::FieldVector<Scalar, 2>, 2 * dim>(Dune::FieldVector<Scalar, 2>(0.0)));
-        interactionVolumes_.resize(problem_.gridView().size(dim));
-
-        asImp_().storeInteractionVolumeInfo();
+        update();
 
         return;
     }
 
+    //! Returns an interaction volume
+    /*!
+     *  \param vertIdx Global index of a vertex in the DUNE grid
+     */
     InteractionVolume& interactionVolume(int vertexIdx)
     {
         return interactionVolumes_[vertexIdx];
     }
 
+    //! Returns an interaction volume
+    /*!
+     *  \param vertIdx Global index of a vertex in the DUNE grid
+     */
     InteractionVolume& interactionVolume(int vertexIdx) const
     {
         return interactionVolumes_[vertexIdx];
     }
 
+    //! Returns the interaction volumes container
     GlobalInteractionVolumeVector& interactionVolumesGlobal()
     {
         return interactionVolumes_;
     }
 
+    //! Returns the interaction volumes container
     GlobalInteractionVolumeVector& interactionVolumesGlobal() const
     {
         return interactionVolumes_;
     }
 
+    //! Returns the area weighting factor for the fluxes
+    /*!
+     *  \param interactionVolume An interaction volume object
+     *  \param elemGlobalIdx Global index of an element in the DUNE grid
+     *  \param elemLocalIdx Local index of an element in the interaction volume
+     *  \param localFaceIdx  Local index of a flux face with respect to an element of the interaction volume
+     *
+     *  \return Ratio of the element face area and the flux face area through which fluxes are calculated by the MPFA method
+     *  (1 if an element does not touches the domain boundary!)
+     */
     Scalar faceAreaFactor(InteractionVolume& interactionVolume, int elemGlobalIdx, int elemLocalIdx, int localFaceIdx)
     {
         Scalar factor = getRealFaceArea(interactionVolume, elemGlobalIdx, elemLocalIdx, localFaceIdx);
@@ -164,6 +194,14 @@ public:
         return factor;
     }
 
+    //! Returns the area weighting factor for the fluxes
+    /*!
+     *  \param elemGlobalIdx Global index of an element in the DUNE grid
+     *  \param indexInInside Local index of the face in the DUNE reference element
+     *
+     *  \return Ratio of the element face area and the flux face area through which fluxes are calculated by the MPFA method
+     *  (1 if an element does not touches the domain boundary!)
+     */
     Scalar faceAreaFactor(int elemGlobalIdx, int indexInInside)
     {
         Scalar factor = getRealFaceArea(elemGlobalIdx, indexInInside);
@@ -172,6 +210,15 @@ public:
         return factor;
     }
 
+    //! Returns the area trough which fluxes are calculated by the MPFA
+    /*!
+     *  \param interactionVolume An interaction volume object
+     *  \param elemGlobalIdx Global index of an element in the DUNE grid
+     *  \param elemLocalIdx Local index of an element in the interaction volume
+     *  \param localFaceIdx  Local index of a flux face with respect to an element of the interaction volume
+     *
+     *  \return flux face area (equal to the element face area if an element does not touches the domain boundary!)
+     */
     Scalar getRealFluxFaceArea(InteractionVolume& interactionVolume, int elemGlobalIdx, int elemLocalIdx, int localFaceIdx)
     {
         Scalar factor = realFluxFaceArea_[elemGlobalIdx][interactionVolume.getIndexOnElement(elemLocalIdx, localFaceIdx)][fluxFaceArea];
@@ -179,6 +226,13 @@ public:
         return factor;
     }
 
+    //! Returns the area trough which fluxes are calculated by the MPFA
+    /*!
+     *  \param elemGlobalIdx Global index of an element in the DUNE grid
+     *  \param indexInInside Local index of the face in the DUNE reference element
+     *
+     *  \return flux face area (equal to the element face area if an element does not touches the domain boundary!)
+     */
     Scalar getRealFluxFaceArea(int elemGlobalIdx, int indexInInside)
     {
         Scalar factor = realFluxFaceArea_[elemGlobalIdx][indexInInside][fluxFaceArea];
@@ -186,6 +240,15 @@ public:
         return factor;
     }
 
+    //! Returns the face area of the element
+    /*!
+     *  \param interactionVolume An interaction volume object
+     *  \param elemGlobalIdx Global index of an element in the DUNE grid
+     *  \param elemLocalIdx Local index of an element in the interaction volume
+     *  \param localFaceIdx  Local index of a flux face with respect to an element of the interaction volume
+     *
+     *  \return the face area of the element
+     */
     Scalar getRealFaceArea(InteractionVolume& interactionVolume, int elemGlobalIdx, int elemLocalIdx, int localFaceIdx)
     {
         Scalar factor = realFluxFaceArea_[elemGlobalIdx][interactionVolume.getIndexOnElement(elemLocalIdx, localFaceIdx)][realFaceArea];
@@ -193,6 +256,13 @@ public:
         return factor;
     }
 
+    //! Returns the face area of the element
+    /*!
+     *  \param elemGlobalIdx Global index of an element in the DUNE grid
+     *  \param indexInInside Local index of the face in the DUNE reference element
+     *
+     *  \return the face area of the element
+     */
     Scalar getRealFaceArea(int elemGlobalIdx, int indexInInside)
     {
         Scalar factor = realFluxFaceArea_[elemGlobalIdx][indexInInside][realFaceArea];
@@ -200,6 +270,10 @@ public:
         return factor;
     }
 
+    //! Constructs a FvMpfaL3dInteractionVolumeContainer object
+    /**
+     * \param problem A problem class object
+     */
     FvMpfaL3dInteractionVolumeContainer(Problem& problem) :
         problem_(problem)
     {
@@ -228,11 +302,20 @@ private:
     Implementation &asImp_()
     { return *static_cast<Implementation *>(this); }
 
-    //! \copydoc Dumux::IMPETProblem::asImp_()
+    //! Returns the implementation of the problem (i.e. static polymorphism)
     const Implementation &asImp_() const
     { return *static_cast<const Implementation *>(this); }
 };
 
+/*! \brief Function for storing the elements of an interaction volume and
+ * constructing a map from a vertex to its surrounding elements
+ *
+ * Stores an element in all interaction volumes it belongs to. Additionally, the global index of an element is stored at the position of its local index according to a DUNE reference element
+ * for every vertex of the element.
+ *
+ * \param element A level 0 Entity of a DUNE grid
+ * \param elemVertMap Vector containing the global vertex-element map
+ */
 template<class TypeTag>
 void FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeSubVolumeElements(const Element& element, std::vector < std::vector<int> >& elemVertMap)
 {
@@ -271,6 +354,14 @@ void FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeSubVolumeElements(const
     elemVertMap[globalVertIdx][0] = globalIdx;
 }
 
+/*! \brief Stores information with respect to DUNE intersections in the interaction volumes
+ *
+ * Stores information with respect to DUNE intersections, such as normals, in the interaction volumes. Assumes a local storage following the DUNE
+ * reference element index, which is performed by the function Dumux::FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeSubVolumeElements(const Element& element, std::vector < std::vector<int> >& elemVertMap).
+ *
+ * \param element A level 0 Entity of a DUNE grid
+ * \param elemVertMap Vector containing the global vertex-element map
+ */
 template<class TypeTag>
 void FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeIntersectionInfo(const Element& element, std::vector < std::vector<int> >& elemVertMap)
 {
@@ -1226,6 +1317,20 @@ void FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeIntersectionInfo(const E
     }
 }
 
+/*! \brief Stores additional information which can be constructed for interaction volumes of non-boundary vertices.
+ *
+ * Stores additional information which can be constructed for interaction volumes of non-boundary vertices:
+ *
+ *  - edge coordinates (coordinates of edge-continuity-points)
+ *  - flux face areas
+ *
+ *  Assumes a local storage following the DUNE reference element index, which is performed by the
+ *  function Dumux::FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeSubVolumeElements(const Element& element, std::vector < std::vector<int> >& elemVertMap).
+ *
+ * \param interactionVolume An interaction volume object
+ * \param vertex The vertex (level dim entity) for which the interaction volume is stored
+ * \param sameLevel Level indicator: true if all elements of an interaction volume are of the same level
+ */
 template<class TypeTag>
 void FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeInnerInteractionVolume(InteractionVolume& interactionVolume, const Vertex& vertex, bool sameLevel)
 {
@@ -1343,6 +1448,23 @@ void FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeInnerInteractionVolume(I
     interactionVolume.setFaceArea(faceArea, 11);
 }
 
+/*! \brief Stores additional information for interaction volumes of boundary vertices.
+ *
+ * Stores additional information for interaction volumes of boundary vertices:
+ *
+ *  - boundary conditions
+ *  - information for flux weighting along boundary faces (see  Wolff 2013: http://elib.uni-stuttgart.de/opus/volltexte/2013/8661/, or
+ *  M. Wolff, Y. Cao, B. Flemisch, R. Helmig, and B. Wohlmuth (2013a). Multi-point flux
+ * approximation L-method in 3D: numerical convergence and application to two-phase
+ * flow through porous media. In P. Bastian, J. Kraus, R. Scheichl, and M. Wheeler,
+ * editors, Simulation of Flow in Porous Media - Applications in Energy and Environment. De Gruyter.)
+ *
+ * Assumes a local storage following the DUNE reference element index, which is performed by the
+ * function Dumux::FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeSubVolumeElements(const Element& element, std::vector < std::vector<int> >& elemVertMap).
+ *
+ * \param interactionVolume An interaction volume object
+ * \param vertex The vertex (level dim entity) for which the interaction volume is stored
+ */
 template<class TypeTag>
 void FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeBoundaryInteractionVolume(InteractionVolume& interactionVolume, const Vertex& vertex)
 {
@@ -1905,12 +2027,13 @@ void FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeBoundaryInteractionVolum
 }
 
 
-// only for 3-D general quadrilateral
+//! \brief Stores interaction volumes for each grid vertex
 template<class TypeTag>
 void FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeInteractionVolumeInfo()
 {
     std::vector < std::vector<int> > elemVertMap(problem_.gridView().size(dim), std::vector<int>(8, -1));
 
+    //Add elements to the interaction volumes and store element-vertex map
     ElementIterator eEndIt = problem_.gridView().template end<0>();
     for (ElementIterator eIt = problem_.gridView().template begin<0>(); eIt != eEndIt; ++eIt)
     {
@@ -1922,17 +2045,15 @@ void FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeInteractionVolumeInfo()
         if (interactionVolumes_[i].getElementNumber() == 0)
             interactionVolumes_[i].printInteractionVolumeInfo();
 
-    // run through all elements
-
+    // Store information related to DUNE intersections for all interaction volumes
     for (ElementIterator eIt = problem_.gridView().template begin<0>(); eIt != eEndIt; ++eIt)
     {
-        // get common geometry information for the following computation
-
         ElementPointer ePtr = *eIt;
         storeIntersectionInfo(*ePtr, elemVertMap);
     }
 
-    // run through all vertices
+    // Complete storage of the interaction volumes using the previously stored information
+    // about the orientation and relationship of the DUNE elements in the interaction volumes (see doc/docextra/3dmpfa)
     VertexIterator vEndIt = problem_.gridView().template end<dim>();
     for (VertexIterator vIt = problem_.gridView().template begin<dim>(); vIt != vEndIt; ++vIt)
     {
@@ -1953,6 +2074,11 @@ void FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeInteractionVolumeInfo()
             DUNE_THROW(Dune::NotImplemented,"Interaction volume is no boundary volume but consists of less than 8 elements");
         }
 
+        // Store information about the MPFA flux face areas for correct flux weighting of
+        // fluxes though faces which intersect the domain boundary (see  M. Wolff, Y. Cao, B. Flemisch, R. Helmig, and B. Wohlmuth (2013a). Multi-point flux
+        // approximation L-method in 3D: numerical convergence and application to two-phase
+        // flow through porous media. In P. Bastian, J. Kraus, R. Scheichl, and M. Wheeler,
+        // editors, Simulation of Flow in Porous Media - Applications in Energy and Environment. De Gruyter.)
         if (!interactionVolume.isBoundaryInteractionVolume())
         {
             ElementPointer& elementPointer1 = interactionVolume.getSubVolumeElement(0);
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dinteractionvolumecontaineradaptive.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dinteractionvolumecontaineradaptive.hh
index 045408ad649ff1533f429b0965dd0f8147c2b490..870e85d871732fc002adbf3bad34f73ab8331f6a 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dinteractionvolumecontaineradaptive.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dinteractionvolumecontaineradaptive.hh
@@ -27,14 +27,21 @@
 
 /**
  * @file
- * @brief  Interactionvolume container for 3-d MPFA L-method
+ * @brief  Interactionvolume container for 3-d MPFA L-method on an h-adaptive grid
  */
 
 namespace Dumux
 {
 
 //! \ingroup FVPressure2p
-/*! Interactionvolume container for 3-d MPFA L-method
+/*! \brief Interactionvolume container for 3-d MPFA L-method on an h-adaptive grid
+ *
+ * Container class which stores MPFA-interaction-volume information for each vertex of a DUNE grid.
+ * Each <tt>InteractionVolume</tt> object stores the information which is necessary to calculate MPFA transmissibility matrices:
+ *
+ * - relationship and orientation of the elements around a vertex (see doc/docextra/3dmpfa)
+ * - geometric information, such as element/face/edge positions, normals, ...
+ *
  */
 template<class TypeTag>
 class FvMpfaL3dInteractionVolumeContainerAdaptive: public FvMpfaL3dInteractionVolumeContainer<TypeTag>
@@ -83,7 +90,7 @@ class FvMpfaL3dInteractionVolumeContainerAdaptive: public FvMpfaL3dInteractionVo
     typedef IndexTranslatorAdaptive IndexTranslator;
 
 public:
-    typedef typename GET_PROP_TYPE(TypeTag, MPFAInteractionVolume) InteractionVolume;
+    typedef typename GET_PROP_TYPE(TypeTag, MPFAInteractionVolume) InteractionVolume;//!< Type for storing an MPFA-interaction-volume. (Usually of type Dumux::FvMpfaL3dInteractionVolume or Dumux::FvMpfaL3dInteractionVolumeAdaptive)
 
 private:
 
@@ -94,11 +101,22 @@ private:
 
 public:
 
+    /*! \brief Returns the set of vertices on an element face
+     *
+     * The DUNE reference elements does not allow to access hanging nodes from a given element face.
+     * However, if a flux through a entire element face has to be calculated, e.g. if single fluxes
+     * have to be updated in an implicit treatment of the transport equation, it is necessary to get
+     * the complete set of vertices on a face: 4 corners + all hanging nodes.
+     */
     std::set<int>& faceVerticeIndices(int globalElemIdx, int faceIdx)
     {
         return faceVertices_[globalElemIdx][faceIdx];
     }
 
+    //! Constructs a FvMpfaL3dInteractionVolumeContainerAdaptive object
+    /**
+     * \param problem A problem class object
+     */
     FvMpfaL3dInteractionVolumeContainerAdaptive(Problem& problem) :
         ParentType(problem), problem_(problem)
     {
@@ -113,11 +131,32 @@ private:
     FaceVerticesVector faceVertices_;
 };
 
+/*! \brief Stores additional information which can be constructed for interaction volumes of non-boundary vertices.
+ *
+ * Stores additional information which can be constructed for interaction volumes of non-boundary vertices:
+ *
+ *  - edge coordinates (coordinates of edge-continuity-points)
+ *  - flux face areas
+ *
+ *  Assumes a local storage following the DUNE reference element index, which is performed by the
+ *  function Dumux::FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeSubVolumeElements(const Element& element, std::vector < std::vector<int> >& elemVertMap).
+ *
+ *  In the case of an adaptive grids with hanging nodes it is important to notice, that the smaller cell of two intersecting cells of different grid level determines the geometric information
+ *    of the flux face (size, position of the center, etc.). This has to be stored correctly for both sides (elements) of
+ *    an intersection.
+ *
+ * \param interactionVolume An interaction volume object
+ * \param vertex The vertex (level dim entity) for which the interaction volume is stored
+ */
 template<class TypeTag>
 void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeInnerInteractionVolume(InteractionVolume& interactionVolume, const Vertex& vertex)
 {
+    // if cells of different grid level appear, the geometric information is constructed going from the
+    // coarsest level to the finest level. This ensures that coarser information ins always overwritten by
+    // finer information.
     if (!interactionVolume.sameLevel())
     {
+        // sort the local element indices according to the grid level
         std::vector < std::vector<int> > levelIdx(8, std::vector<int>(2));
         for (int i = 0; i < 8; i++)
         {
@@ -127,6 +166,10 @@ void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeInnerInteraction
 
         std::sort(levelIdx.begin(), levelIdx.end(), sort_compare);
 
+        // Generate and store the geometric information going from the coarsest to the finest level.
+        // For the calculation we take advantage from the fact that the ordering inside the interaction volume
+        // with respect to the DUNE reference element is known due to the storage process of the elements in
+        // Dumux::FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeSubVolumeElements(const Element& element, std::vector < std::vector<int> >& elemVertMap)
         for (int i = 0; i < 8; i++)
         {
             int idx = levelIdx[i][0];
@@ -236,6 +279,25 @@ void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeInnerInteraction
     }
 }
 
+/*! \brief Stores additional information which can be constructed for interaction volumes around hanging nodes.
+ *
+ *  Stores additional information which can be constructed for interaction volumes around hanging nodes.
+ *
+ *  - missing cells: As hanging nodes cannot be accessed from a cell face using the DUNE reference element,
+ *                   the attached coarser cells do not appear in the interaction volume object after execution
+ *                   of the function Dumux::FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeSubVolumeElements(const Element& element, std::vector < std::vector<int> >& elemVertMap).
+ *                   We take advantage of this fact because it allows us to identify the type of hanging-node interaction volume.
+ *                   If, for example, only two cells are stored, we know that the interaction volume is of type 5 according to
+ *                   Wolff 2013: http://elib.uni-stuttgart.de/opus/volltexte/2013/8661/
+ *                   As first step, the orientation of the stored elements in the local interaction volume indices
+ *                   is rotated such that cells 1 and 2 according to the local index (see doc/docextra/3dmpfa) are always of the smallest existing level. Thus, the relationship between
+ *                   the grid elements in a hanging-node-interaction-volume is always known if the type of interaction volume is known
+ *                   (1-7, see Wolff 2013: http://elib.uni-stuttgart.de/opus/volltexte/2013/8661/). As a second step, the missing cells
+ *                   are added.
+ *  - edge coordinates (coordinates of edge-continuity-points)
+ *  - flux face areas
+ *
+ */
 template<class TypeTag>
 void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeHangingNodeInteractionVolume(InteractionVolume& interactionVolume, const Vertex& vertex)
 {
@@ -243,6 +305,7 @@ void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeHangingNodeInter
 
     interactionVolume.setCenterPosition(centerPos);
 
+    // sort the local element indices according to the grid level
     std::vector < std::vector<int> > levelIdx(8, std::vector<int>(2));
     for (int i = 0; i < 8; i++)
     {
@@ -255,7 +318,10 @@ void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeHangingNodeInter
 
     std::sort(levelIdx.begin(), levelIdx.end(), sort_compare);
 
-
+    // Generate and store the geometric information going from the coarsest to the finest level.
+     // For the calculation we take advantage from the fact that the ordering inside the interaction volume
+     // with respect to the DUNE reference element is known due to the storage process of the elements in
+     // Dumux::FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeSubVolumeElements(const Element& element, std::vector < std::vector<int> >& elemVertMap)
     for (int i = 0; i < 8; i++)
     {
         if (levelIdx[i][1] < 0)
@@ -362,8 +428,12 @@ void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeHangingNodeInter
         }
     }
 
+    // Choose the type of the hanging-node-interaction-volume depending on the number of stored fine
+    // elements (see dissertation M. Wolff, http://elib.uni-stuttgart.de/opus/volltexte/2013/8661/)
+    // and add missing elements and geometric information
     switch (interactionVolume.getElementNumber())
     {
+    // hanging-node interaction volume of type 5 or 7
     case 2:
         {
             InteractionVolume hangingNodeVolume;
@@ -567,6 +637,7 @@ void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeHangingNodeInter
 
             break;
         }
+        //hanging-node interaction volume of type 1, 3 or 4
     case 4:
         {
             InteractionVolume hangingNodeVolume;
@@ -1325,6 +1396,7 @@ void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeHangingNodeInter
 
             break;
         }
+        //hanging-node interaction volume of type 2 or 6
     case 6:
         {
             InteractionVolume hangingNodeVolume;
@@ -1528,12 +1600,13 @@ void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeHangingNodeInter
     }
 }
 
-// only for 3-D general quadrilateral
+//! \brief Stores interaction volumes for each grid vertex
 template<class TypeTag>
 void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeInteractionVolumeInfo()
 {
     std::vector < std::vector<int> > elemVertMap(problem_.gridView().size(dim), std::vector<int>(8, -1));
 
+    //Add elements to the interaction volumes and store element-vertex map
     ElementIterator eEndIt = problem_.gridView().template end<0>();
     for (ElementIterator eIt = problem_.gridView().template begin<0>(); eIt != eEndIt; ++eIt)
     {
@@ -1545,12 +1618,9 @@ void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeInteractionVolum
         if (this->interactionVolumes_[i].getElementNumber() == 0)
             this->interactionVolumes_[i].printInteractionVolumeInfo();
 
-    // run through all elements
-
+    // Store information related to DUNE intersections for all interaction volumes
     for (ElementIterator eIt = problem_.gridView().template begin<0>(); eIt != eEndIt; ++eIt)
     {
-        // get common geometry information for the following computation
-
         ElementPointer ePtr = *eIt;
         this->storeIntersectionInfo(*ePtr, elemVertMap);
     }
@@ -1558,7 +1628,8 @@ void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeInteractionVolum
     faceVertices_.clear();
     faceVertices_.resize(problem_.gridView().size(0));
 
-    // run through all vertices
+    // Complete storage of the interaction volumes using the previously stored information
+    // about the orientation and relationship of the DUNE elements in the interaction volumes (see doc/docextra/3dmpfa)
     VertexIterator vEndIt = problem_.gridView().template end<dim>();
     for (VertexIterator vIt = problem_.gridView().template begin<dim>(); vIt != vEndIt; ++vIt)
     {
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dpressure2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dpressure2p.hh
index 356b3a52b8adb19541a1c48b43678b018ee29e7b..8a6cdcd3aed1083760f6067cf5f80c6b64589b0b 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dpressure2p.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dpressure2p.hh
@@ -27,29 +27,43 @@
 
 /**
  * @file
- * @brief  Finite Volume-MPFAL implementation of a two-phase pressure equation
- * @brief  Remark1: only for 3-D hexahedron grid.
- * @brief  Remark3: number of grid cells in each direction > 1
+ * @brief  3-d finite Volume-MPFAL implementation of a two-phase pressure equation
+ * @brief  Remark1: only for 3-D hexahedrons of quadrilaterals.
  */
 
 namespace Dumux
 {
 //! \ingroup FVPressure2p
-/*! Finite Volume-MPFAL-Implementation of the equation
- * \f$ - \text{div}\, \mathbf{v}_t = - \text{div}\, (\lambda_t \mathbf{K} \text{grad}\, \Phi_w + f_n \lambda_t \mathbf{K} \text{grad}\, \Phi_{cap}   ) = 0, \f$, or
+/*! \brief 3-d finite volume MPFA L-method discretization of a two-phase flow pressure equation of the sequential IMPES model
+ *
+ * Finite Volume-MPFAL-Implementation of the equation
+ * \f$ - \text{div}\, \mathbf{v}_t = - \text{div}\, (\lambda_t \mathbf{K} \text{grad}\, \Phi_w + f_n \lambda_t \mathbf{K} \text{grad}\, \Phi_{cap}   ) = 0, \f$,
+ * or
  * \f$ - \text{div}\, \mathbf{v}_t = - \text{div}\, (\lambda_t \mathbf{K} \text{grad}\, \Phi_n - f_w \lambda_t \mathbf{K} \text{grad}\, \Phi_{cap}   ) = 0, \f$.
- * \f$\Phi = g\f$ on \f$\Gamma_1\f$, and
- * \f$-\text{div}\, \mathbf{v}_t \cdot \mathbf{n} = J\f$
- * on \f$\Gamma_2\f$. Here,
- * \f$Phi_\alpha \f$ denotes the potential of phase \f$\alpha\f$, \f$K\f$ the intrinsic permeability,
- * \f$\lambda_t\f$ the total mobility, \f$this->f_\alpha\f$ the phase fractional flow function.
+ *
+ * \f$ \Phi = g \f$ on \f$ \Gamma_1 \f$, and
+ * \f$ - \text{div} \, \mathbf{v}_t \cdot \mathbf{n} = J \f$
+ * on \f$ \Gamma_2 \f$.
+ *
+ * Here, \f$ \Phi_\alpha \f$ denotes the potential of phase \f$ \alpha \f$, \f$ \mathbf{K} \f$ the intrinsic permeability,
+ * \f$ \lambda_t \f$ the total mobility, \f$ f_\alpha \f$ the phase fractional flow function.
  *
  * More details on the equations can be found in
+ *
+ * Wolff 2013: http://elib.uni-stuttgart.de/opus/volltexte/2013/8661/
+ *
  * M. Wolff, Y. Cao, B. Flemisch, R. Helmig, and B. Wohlmuth (2013a). Multi-point flux
  * approximation L-method in 3D: numerical convergence and application to two-phase
  * flow through porous media. In P. Bastian, J. Kraus, R. Scheichl, and M. Wheeler,
  * editors, Simulation of Flow in Porous Media - Applications in Energy and Environment. De Gruyter.
  *
+ * M. Wolff, B. Flemisch, R. Helmig, I. Aavatsmark.
+ * Treatment of tensorial relative permeabilities with multipoint flux approximation.
+ * International Journal of Numerical Analysis and Modeling (9), pp. 725-744, 2012.
+ *
+ * Remark1: only for 3-D hexahedrons of quadrilaterals.
+ *
+ * \tparam TypeTag The problem Type Tag
  */
 template<class TypeTag>
 class FvMpfaL3dPressure2p: public FVPressure<TypeTag>
@@ -142,8 +156,8 @@ class FvMpfaL3dPressure2p: public FVPressure<TypeTag>
     typedef typename GET_PROP_TYPE(TypeTag, MPFAInteractionVolumeContainer) InteractionVolumeContainer;
     typedef  Dumux::FvMpfaL3dTransmissibilityCalculator<TypeTag> TransmissibilityCalculator;
 public:
-    typedef typename TransmissibilityCalculator::TransmissibilityType TransmissibilityType;
-    typedef typename GET_PROP_TYPE(TypeTag, MPFAInteractionVolume) InteractionVolume;
+    typedef typename TransmissibilityCalculator::TransmissibilityType TransmissibilityType;//!< Type including methods for calculation of MPFA transmissibilities
+    typedef typename GET_PROP_TYPE(TypeTag, MPFAInteractionVolume) InteractionVolume;//!< Type for storing interaction volume information
 protected:
     //initializes the matrix to store the system of equations
     friend class FVPressure<TypeTag>;
@@ -161,6 +175,10 @@ public:
     //constitutive functions are initialized and stored in the variables object
     void updateMaterialLaws();
 
+    /*! \brief Initializes the pressure model
+     *
+     * \copydetails ParentType::initialize()
+     */
     void initialize(bool solveTwice = true)
     {
         ElementIterator element = problem_.gridView().template begin<0>();
@@ -188,6 +206,7 @@ public:
         return;
     }
 
+    //! \brief Globally stores the pressure solution
     void storePressureSolution()
     {
         // iterate through leaf grid an evaluate c0 at cell center
@@ -198,6 +217,10 @@ public:
         }
     }
 
+    /*! \brief Stores the pressure solution of a cell
+     *
+     * \param element Dune grid element
+     */
     void storePressureSolution(const Element& element)
     {
         int globalIdx = problem_.variables().index(element);
@@ -245,7 +268,11 @@ public:
         cellData.fluxData().resetVelocity();
     }
 
-    //! updates the pressure field (analog to update function in Dumux::IMPET)
+    /*! \brief Pressure update
+     *
+     * \copydetails ParentType::update()
+     *
+     */
     void update()
     {
         int size = problem_.gridView().size(0);
@@ -295,6 +322,22 @@ public:
         return;
     }
 
+    /* \brief Volume correction term to correct for unphysical saturation overshoots/undershoots.
+     *
+     * These can occur if the estimated time step for the explicit transport was too large. Correction by an artificial source term allows to correct
+     * this errors due to wrong time-stepping without losing mass conservation. The error term looks as follows:
+     * \f[
+     *  q_{error} = \begin{cases}
+     *          S < 0 & \frac{S}{\Delta t} V \\
+     *          S > 1 & \frac{(S - 1)}{\Delta t} V \\
+     *          0 \le S \le 1 & 0
+     *      \end{cases}
+     *  \f]
+     *
+     *  \param cellData The IMPES <tt>CellData</tt> object of the current cell.
+     *
+     *  \return The scalar value of the error term.
+    */
     Scalar evaluateErrorTerm(CellData& cellData)
     {
         //error term for incompressible models to correct unphysical saturation over/undershoots due to saturation transport
@@ -404,16 +447,25 @@ public:
         }
     }
 
+    //! Returns the global container of the stored interaction volumes
     InteractionVolumeContainer& interactionVolumes()
     {
         return interactionVolumes_;
     }
 
+    //! Returns the transmissibility calculator
+    /*!
+     *  Object including methods for the MPFA transmissibility calculation
+     */
     TransmissibilityCalculator& transmissibilityCalculator()
     {
         return transmissibilityCalculator_;
     }
 
+    //! Constructs a FvMpfaL3dPressure2p object
+    /**
+     * \param problem A problem class object
+     */
     FvMpfaL3dPressure2p(Problem& problem) :
         ParentType(problem), problem_(problem),
         interactionVolumes_(problem), transmissibilityCalculator_(problem),
@@ -453,8 +505,8 @@ private:
     Problem& problem_;
 
 protected:
-    InteractionVolumeContainer interactionVolumes_;
-    TransmissibilityCalculator transmissibilityCalculator_;
+    InteractionVolumeContainer interactionVolumes_;//!<Global container of the stored interaction volumes
+    TransmissibilityCalculator transmissibilityCalculator_;//!<The transmissibility calculator including methods for the MPFA transmissibility calculation
 
 private:
     //! Returns the implementation of the problem (i.e. static polymorphism)
@@ -484,6 +536,7 @@ private:
     static const int velocityType_ = GET_PROP_VALUE(TypeTag, VelocityFormulation);//!< gives kind of velocity used (\f$ 0 = v_w\f$, \f$ 1 = v_n\f$, \f$ 2 = v_t\f$)
 };
 
+//! Initializes the sparse matrix for the pressure solution
 template<class TypeTag>
 void FvMpfaL3dPressure2p<TypeTag>::initializeMatrix()
 {
@@ -493,6 +546,7 @@ void FvMpfaL3dPressure2p<TypeTag>::initializeMatrix()
     this->A_.endindices();
 }
 
+//! Initializes the row size of the sparse matrix for the pressure solution
 template<class TypeTag>
 void FvMpfaL3dPressure2p<TypeTag>::initializeMatrixRowSize()
 {
@@ -532,6 +586,7 @@ void FvMpfaL3dPressure2p<TypeTag>::initializeMatrixRowSize()
     return;
 }
 
+//! Initializes the indices of the sparse matrix for the pressure solution
 template<class TypeTag>
 void FvMpfaL3dPressure2p<TypeTag>::initializeMatrixIndices()
 {
@@ -568,7 +623,7 @@ void FvMpfaL3dPressure2p<TypeTag>::initializeMatrixIndices()
     return;
 }
 
-// only for 3-D general hexahedron
+//! assembles the global matrix and rhs vector for the pressure solution
 template<class TypeTag>
 void FvMpfaL3dPressure2p<TypeTag>::assemble()
 {
@@ -601,6 +656,7 @@ void FvMpfaL3dPressure2p<TypeTag>::assemble()
     return;
 }
 
+//! assembles the matrix entries of one interaction volume into the global matrix
 template<class TypeTag>
 void FvMpfaL3dPressure2p<TypeTag>::assembleInnerInteractionVolume(InteractionVolume& interactionVolume)
 {
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dpressure2padaptive.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dpressure2padaptive.hh
index 45c1cf5516f04e25562023a7650c4ce4c9a8c045..43a39a0be332c7fd3834b5c38457073ce8db959d 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dpressure2padaptive.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dpressure2padaptive.hh
@@ -26,27 +26,47 @@
 
 /**
  * @file
- * @brief  Finite Volume-MPFAL implementation of a two-phase pressure equation
- * @brief  Remark1: only for 3-D hexahedron grid.
- * @brief  Remark3: number of grid cells in each direction > 1
+ * @brief  3-d finite Volume-MPFAL implementation of a two-phase pressure equation on h-adaptive grids
+ * @brief  Remark1: only for 3-D hexahedrons of quadrilaterals.
+ * @brief  Remark2: number of grid cells in each direction > 1
  */
 
 namespace Dumux
 {
 //! \ingroup FVPressure2p
-/*! Finite Volume-MPFAL-Implementation of the equation
- * \f$ - \text{div}\, \mathbf{v}_t = - \text{div}\, (\lambda_t \mathbf{K} \text{grad}\, \Phi_w + f_n \lambda_t \mathbf{K} \text{grad}\, \Phi_{cap}   ) = 0, \f$, or
+/*! \brief 3-d finite volume MPFA L-method discretization of a two-phase flow pressure equation of the sequential IMPES model on h-adaptive grids.
+ *
+ * * Finite Volume-MPFAL-Implementation of the equation
+ *
+ * \f$ - \text{div}\, \mathbf{v}_t = - \text{div}\, (\lambda_t \mathbf{K} \text{grad}\, \Phi_w + f_n \lambda_t \mathbf{K} \text{grad}\, \Phi_{cap}   ) = 0, \f$,
+ * or
  * \f$ - \text{div}\, \mathbf{v}_t = - \text{div}\, (\lambda_t \mathbf{K} \text{grad}\, \Phi_n - f_w \lambda_t \mathbf{K} \text{grad}\, \Phi_{cap}   ) = 0, \f$.
- * \f$\Phi = g\f$ on \f$\Gamma_1\f$, and
- * \f$-\text{div}\, \mathbf{v}_t \cdot \mathbf{n} = J\f$
- * on \f$\Gamma_2\f$. Here,
- * \f$Phi_\alpha \f$ denotes the potential of phase \f$\alpha\f$, \f$K\f$ the intrinsic permeability,
- * \f$\lambda_t\f$ the total mobility, \f$this->f_\alpha\f$ the phase fractional flow function.
+ *
+ * \f$ \Phi = g \f$ on \f$ \Gamma_1 \f$, and
+ * \f$ - \text{div} \, \mathbf{v}_t \cdot \mathbf{n} = J \f$
+ * on \f$ \Gamma_2 \f$.
+ *
+ * Here, \f$ \Phi_\alpha \f$ denotes the potential of phase \f$ \alpha \f$, \f$ \mathbf{K} \f$ the intrinsic permeability,
+ * \f$ \lambda_t \f$ the total mobility, \f$ f_\alpha \f$ the phase fractional flow function.
+ *
+ * More details on the equations can be found in
+ *
+ * Wolff 2013: http://elib.uni-stuttgart.de/opus/volltexte/2013/8661/
  *
  * M. Wolff, Y. Cao, B. Flemisch, R. Helmig, and B. Wohlmuth (2013a). Multi-point flux
  * approximation L-method in 3D: numerical convergence and application to two-phase
  * flow through porous media. In P. Bastian, J. Kraus, R. Scheichl, and M. Wheeler,
  * editors, Simulation of Flow in Porous Media - Applications in Energy and Environment. De Gruyter.
+ *
+ * M. Wolff, B. Flemisch, R. Helmig, I. Aavatsmark.
+ * Treatment of tensorial relative permeabilities with multipoint flux approximation.
+ * International Journal of Numerical Analysis and Modeling (9), pp. 725-744, 2012.
+ *
+ * Remark1: only for 3-D hexahedrons of quadrilaterals.
+ *
+ * Remark2: number of grid cells in each direction > 1
+ *
+ * \tparam TypeTag The problem Type Tag
  */
 template<class TypeTag>
 class FvMpfaL3dPressure2pAdaptive: public FvMpfaL3dPressure2p<TypeTag>
@@ -138,10 +158,11 @@ class FvMpfaL3dPressure2pAdaptive: public FvMpfaL3dPressure2p<TypeTag>
     typedef typename GET_PROP_TYPE(TypeTag, MPFAInteractionVolumeContainer) InteractionVolumeContainer;
     typedef  Dumux::FvMpfaL3dTransmissibilityCalculator<TypeTag> TransmissibilityCalculator;
 public:
-    typedef typename TransmissibilityCalculator::TransmissibilityType TransmissibilityType;
-    typedef typename GET_PROP_TYPE(TypeTag, MPFAInteractionVolume) InteractionVolume;
+    typedef typename TransmissibilityCalculator::TransmissibilityType TransmissibilityType;//!< Type including methods for calculation of MPFA transmissibilities
+    typedef typename GET_PROP_TYPE(TypeTag, MPFAInteractionVolume) InteractionVolume;//!< Type for storing interaction volume information
 
 protected:
+    //initializes the matrix to store the system of equations
     friend class FVPressure<TypeTag>;
     void initializeMatrix();
     void initializeMatrixRowSize();
@@ -152,6 +173,10 @@ protected:
     void assemble();
     void assembleHangingNodeInteractionVolume(InteractionVolume& interactionVolume);
 public:
+    /*! \brief Initializes the pressure model
+     *
+     * \copydetails ParentType::initialize()
+     */
     void initialize(bool solveTwice = true)
     {
         ElementIterator element = problem_.gridView().template begin<0>();
@@ -189,6 +214,10 @@ public:
         ParentType::update();
     }
 
+    //! Constructs a FvMpfaL3dPressure2pAdaptive object
+    /**
+     * \param problem A problem class object
+     */
     FvMpfaL3dPressure2pAdaptive(Problem& problem) :
         ParentType(problem), problem_(problem),
         gravity_(problem.gravity())
@@ -219,7 +248,6 @@ public:
 private:
     Problem& problem_;
 
-private:
     //! Returns the implementation of the problem (i.e. static polymorphism)
     Implementation &asImp_()
     {   return *static_cast<Implementation *>(this);}
@@ -239,6 +267,7 @@ private:
     static const int velocityType_ = GET_PROP_VALUE(TypeTag, VelocityFormulation);//!< gives kind of velocity used (\f$ 0 = v_w\f$, \f$ 1 = v_n\f$, \f$ 2 = v_t\f$)
 };
 
+//! Initializes the sparse matrix for the pressure solution
 template<class TypeTag>
 void FvMpfaL3dPressure2pAdaptive<TypeTag>::initializeMatrix()
 {
@@ -248,6 +277,7 @@ void FvMpfaL3dPressure2pAdaptive<TypeTag>::initializeMatrix()
     this->A_.endindices();
 }
 
+//! Initializes the row size of the sparse matrix for the pressure solution
 template<class TypeTag>
 void FvMpfaL3dPressure2pAdaptive<TypeTag>::initializeMatrixRowSize()
 {
@@ -335,6 +365,7 @@ void FvMpfaL3dPressure2pAdaptive<TypeTag>::initializeMatrixRowSize()
     return;
 }
 
+//! Initializes the indices of the sparse matrix for the pressure solution
 template<class TypeTag>
 void FvMpfaL3dPressure2pAdaptive<TypeTag>::initializeMatrixIndices()
 {
@@ -382,7 +413,7 @@ void FvMpfaL3dPressure2pAdaptive<TypeTag>::initializeMatrixIndices()
     return;
 }
 
-// only for 3-D general hexahedron
+//! assembles the global matrix and rhs vector for the pressure solution
 template<class TypeTag>
 void FvMpfaL3dPressure2pAdaptive<TypeTag>::assemble()
 {
@@ -427,6 +458,7 @@ void FvMpfaL3dPressure2pAdaptive<TypeTag>::assemble()
     return;
 }
 
+//! assembles the matrix entries of one hanging node interaction volume into the global matrix
 template<class TypeTag>
 void FvMpfaL3dPressure2pAdaptive<TypeTag>::assembleHangingNodeInteractionVolume(InteractionVolume& interactionVolume)
 {
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dpressurevelocity2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dpressurevelocity2p.hh
index 49050cb03b4cdaf2074aef035c2d6748bccc6b12..8f59356fa460cfc478ee5d499840ab07c462410f 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dpressurevelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dpressurevelocity2p.hh
@@ -24,11 +24,21 @@
 
 /**
  * @file
- * @brief  Base class for defining an instance of a numerical diffusion model
+ * @brief  3d Velocity Field from a finite volume solution of a pressure equation using a MPFA L-method.
  */
 
 namespace Dumux
 {
+//! \ingroup FVPressure2p
+/*! \brief Class for the calculation of 3d velocities from the  pressure solution of an IMPES scheme using a MPFA L-method.
+ *
+ * Can be used for calculating the complete velocity field before the solution of the transport equation (more efficient),
+ * or for face-wise velocity calculation directly in the transport solution (less efficient).
+ *
+ * Remark1: only for 3-D Hexahedrons of quadrilaterals!
+ *
+ * \tparam TypeTag The problem Type Tag
+ */
 template<class TypeTag> class FvMpfaL3dPressureVelocity2p: public FvMpfaL3dPressure2p<TypeTag>
 {
     typedef FvMpfaL3dPressure2p<TypeTag> ParentType;
@@ -87,6 +97,10 @@ template<class TypeTag> class FvMpfaL3dPressureVelocity2p: public FvMpfaL3dPress
     typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix;
 
 public:
+    //! Constructs a FvMpfaL3dPressureVelocity2p object
+    /*!
+     * \param problem A problem class object
+     */
     FvMpfaL3dPressureVelocity2p(Problem& problem) :
         ParentType(problem), problem_(problem), velocity_(problem)
     {
@@ -106,6 +120,7 @@ public:
     void calculateVelocity(const Intersection&, CellData&);
     void calculateVelocityOnBoundary(const Intersection& intersection, CellData& cellData);
 
+    //! Function for updating the velocity field if iterations are necessary in the transport solution
     void updateVelocity()
     {
         this->updateMaterialLaws();
@@ -116,6 +131,10 @@ public:
         calculateVelocity();
     }
 
+    /*! \brief Initializes pressure and velocity
+     *
+     * \copydetails ParentType::initialize()
+     */
     void initialize(bool solveTwice = true)
     {
         ElementIterator element = problem_.gridView().template begin<0>();
@@ -139,6 +158,11 @@ public:
         return;
     }
 
+    /*! \brief Pressure and velocity update
+     *
+     * \copydetails ParentType::update()
+     *
+     */
     void update()
     {
         ParentType::update();
@@ -156,8 +180,16 @@ public:
         return calcVelocityInTransport_;
     }
 
-    //! \brief Write data files
-    /*  \param name file name */
+    /*! \brief Adds velocity output to the output file
+     *
+     * Adds the phase velocities or a total velocity (depending on the formulation) to the output.
+     * If the VtkOutputLevel is equal to zero (default) only primary variables are written,
+     * if it is larger than zero also secondary variables are written.
+     *
+     * \tparam MultiWriter Class defining the output writer
+     * \param writer The output writer (usually a <tt>VTKMultiWriter</tt> object)
+     *
+     */
     template<class MultiWriter>
     void addOutputVtkFields(MultiWriter &writer)
     {
@@ -178,7 +210,11 @@ private:
 };
 // end of template
 
-// only for 3-D general hexahedron
+/*! \brief Calculates the velocities at a cell-cell interfaces for the entire simulation grid.
+ *
+ * Calculates the velocities at a cell-cell interfaces from a given pressure field for the entire simulation grid.
+ *
+ */
 template<class TypeTag>
 void FvMpfaL3dPressureVelocity2p<TypeTag>::calculateVelocity()
 {
@@ -255,6 +291,13 @@ void FvMpfaL3dPressureVelocity2p<TypeTag>::calculateVelocity()
     return;
 }
 
+/*! \brief Calculates the velocity at a cell-cell interface.
+ *
+ * Calculates the velocity at a cell-cell interface from a given pressure field.
+ *
+ * \param intersection Intersection of two grid cells
+ * \param cellData Object containing all model relevant cell data
+ */
 template<class TypeTag>
 void FvMpfaL3dPressureVelocity2p<TypeTag>::calculateVelocity(const Intersection& intersection, CellData& cellData)
 {
@@ -329,6 +372,13 @@ void FvMpfaL3dPressureVelocity2p<TypeTag>::calculateVelocity(const Intersection&
     cellDataJ.fluxData().setVelocityMarker(indexInOutside);
 }
 
+/*! \brief Calculates the velocity at a boundary.
+ *
+ * Calculates the velocity at a boundary from a given pressure field.
+ *
+ * \param intersection Boundary intersection
+ * \param cellData Object containing all model relevant cell data
+ */
 template<class TypeTag>
 void FvMpfaL3dPressureVelocity2p<TypeTag>::calculateVelocityOnBoundary(const Intersection& intersection, CellData& cellData)
 {
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dpressurevelocity2padaptive.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dpressurevelocity2padaptive.hh
index 71c4a7bfca7e50f8ff99873e0a14fae877be66f7..7f0bf8f134a4ead1c016665db39e939ae07e1d4f 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dpressurevelocity2padaptive.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dpressurevelocity2padaptive.hh
@@ -23,11 +23,23 @@
 
 /**
  * @file
- * @brief  Base class for defining an instance of a numerical diffusion model
+ * @brief  3d Velocity Field from a finite volume solution of a pressure equation using a grid adaptive MPFA L-method.
  */
 
 namespace Dumux
 {
+//! \ingroup FVPressure2p
+/*! \brief Class for the calculation of 3d velocities from the  pressure solution of an IMPES scheme using a grid adaptive MPFA L-method.
+ *
+ * Can be used for calculating the complete velocity field before the solution of the transport equation (more efficient),
+ * or for face-wise velocity calculation directly in the transport solution (less efficient).
+ *
+ * Remark1: only for 3-D Hexahedrons of quadrilaterals!
+ *
+ * Remark2: Allowed difference in grid levels of two neighboring cells: 1
+ *
+ * \tparam TypeTag The problem Type Tag
+ */
 template<class TypeTag> class FvMpfaL3dPressureVelocity2pAdaptive: public FvMpfaL3dPressure2pAdaptive<TypeTag>
 {
     typedef FvMpfaL3dPressure2pAdaptive<TypeTag> ParentType;
@@ -87,6 +99,10 @@ template<class TypeTag> class FvMpfaL3dPressureVelocity2pAdaptive: public FvMpfa
     typedef Dune::FieldVector<Scalar, dim> DimVector;
 
 public:
+    //! Constructs a FvMpfaL3dPressureVelocity2pAdaptive object
+    /*!
+     * \param problem A problem class object
+     */
     FvMpfaL3dPressureVelocity2pAdaptive(Problem& problem) :
         ParentType(problem), problem_(problem), velocity_(problem)
 {
@@ -100,13 +116,11 @@ public:
 
     void calculateVelocity();
 
-public:
-
     // Calculates the velocity at a cell-cell interface.
     void calculateVelocity(const Intersection&, CellData&);
     void calculateVelocityOnBoundary(const Intersection& intersection, CellData& cellData);
 
-
+    //! Function for updating the velocity field if iterations are necessary in the transport solution
     void updateVelocity()
     {
         this->updateMaterialLaws();
@@ -117,6 +131,10 @@ public:
             calculateVelocity();
     }
 
+    /*! \brief Initializes pressure and velocity
+     *
+     * \copydetails ParentType::initialize()
+     */
     void initialize(bool solveTwice = true)
     {
         ElementIterator element = problem_.gridView().template begin<0>();
@@ -140,6 +158,11 @@ public:
         return;
     }
 
+    /*! \brief Pressure and velocity update
+     *
+     * \copydetails ParentType::update()
+     *
+     */
     void update()
     {
         ParentType::update();
@@ -157,8 +180,16 @@ public:
         return calcVelocityInTransport_;
     }
 
-    //! \brief Write data files
-    /*  \param name file name */
+    /*! \brief Adds velocity output to the output file
+     *
+     * Adds the phase velocities or a total velocity (depending on the formulation) to the output.
+     * If the VtkOutputLevel is equal to zero (default) only primary variables are written,
+     * if it is larger than zero also secondary variables are written.
+     *
+     * \tparam MultiWriter Class defining the output writer
+     * \param writer The output writer (usually a <tt>VTKMultiWriter</tt> object)
+     *
+     */
     template<class MultiWriter>
     void addOutputVtkFields(MultiWriter &writer)
     {
@@ -179,7 +210,11 @@ private:
 };
 // end of template
 
-// only for 3-D general hexahedron
+/*! \brief Calculates the velocities at a cell-cell interfaces for the entire simulation grid.
+ *
+ * Calculates the velocities at a cell-cell interfaces from a given pressure field for the entire simulation grid.
+ *
+ */
 template<class TypeTag>
 void FvMpfaL3dPressureVelocity2pAdaptive<TypeTag>::calculateVelocity()
 {
@@ -265,6 +300,13 @@ void FvMpfaL3dPressureVelocity2pAdaptive<TypeTag>::calculateVelocity()
     return;
 }
 
+/*! \brief Calculates the velocity at a cell-cell interface.
+ *
+ * Calculates the velocity at a cell-cell interface from a given pressure field.
+ *
+ * \param intersection Intersection of two grid cells
+ * \param cellData Object containing all model relevant cell data
+ */
 template<class TypeTag>
 void FvMpfaL3dPressureVelocity2pAdaptive<TypeTag>::calculateVelocity(const Intersection& intersection, CellData& cellData)
 {
@@ -431,6 +473,13 @@ void FvMpfaL3dPressureVelocity2pAdaptive<TypeTag>::calculateVelocity(const Inter
     }
 }
 
+/*! \brief Calculates the velocity at a boundary.
+ *
+ * Calculates the velocity at a boundary from a given pressure field.
+ *
+ * \param intersection Boundary intersection
+ * \param cellData Object containing all model relevant cell data
+ */
 template<class TypeTag>
 void FvMpfaL3dPressureVelocity2pAdaptive<TypeTag>::calculateVelocityOnBoundary(const Intersection& intersection, CellData& cellData)
 {
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dtransmissibilitycalculator.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dtransmissibilitycalculator.hh
index 8d57d10604ad3695b9d8f895c4aa82343618a9eb..3e69b9ba83c7a2838d9e76ad81e3a380e37e4cb9 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dtransmissibilitycalculator.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dtransmissibilitycalculator.hh
@@ -24,13 +24,30 @@
 
 /**
  * @file
- * @brief Provides methods for transmissibility calculation
+ * @brief Provides methods for transmissibility calculation in 3-d
  */
 
 namespace Dumux
 {
 //! \ingroup FVPressure2p
-/*! Provides methods for transmissibility calculation.
+/*! \brief Provides methods for transmissibility calculation in 3-d.
+ *
+ *  The transmissibilities are calculated using the MPFA L-method.
+ *
+ * [1] Aavatsmark et al. A new finite-volume approach to efficient discretization on challenging grids. SPE Journal 15, 2010.
+ *
+ * [2] M. Wolff, Y. Cao, B. Flemisch, R. Helmig, and B. Wohlmuth (2013a). Multi-point flux
+ *     approximation L-method in 3D: numerical convergence and application to two-phase
+ *     flow through porous media. In P. Bastian, J. Kraus, R. Scheichl, and M. Wheeler,
+ *     editors, Simulation of Flow in Porous Media - Applications in Energy and Environment. De Gruyter.
+ *
+ *  Various parameters can be defined via an input parameter file or the property system:
+ *
+ * MPFAEnableSimpleLStencil - enables the two centered flux stencils
+ * MPFAEnableComplexLStencil - enables the two non-centered flux stencils
+ * MPFAEnableTPFA - enables the use of TPFA if neighboring cells are of the same grid level
+ * MPFATransmissibilityCriterion - 0: Criterion of [1] (default), 1: Criterion of [2]
+ *
  */
 template<class TypeTag>
 class FvMpfaL3dTransmissibilityCalculator
@@ -61,12 +78,8 @@ class FvMpfaL3dTransmissibilityCalculator
 
 
 public:
-    typedef Dune::FieldMatrix<Scalar, dim, 2 * dim - dim + 1> TransmissibilityType;
+    typedef Dune::FieldMatrix<Scalar, dim, 2 * dim - dim + 1> TransmissibilityType;//!< Type of the transmissibility matrix
 
-    enum
-        {
-            hangingNodeCell = 99
-        };
     int chooseTransmissibility(TransmissibilityType& transmissibilityOne, TransmissibilityType& transmissibilityTwo, int lTypeOne, int lTypeTwo);
 
     int transmissibility(Dune::FieldMatrix<Scalar,dim,2*dim-dim+1>& transmissibility,
@@ -109,7 +122,10 @@ public:
                                  int idx1, int idx2, int idx3, int idx6);
 
 
-
+    //! Constructs a FvMpfaL3dTransmissibilityCalculator object
+    /**
+     * \param problem A problem class object
+     */
     FvMpfaL3dTransmissibilityCalculator(Problem& problem) :
         problem_(problem)
     {
@@ -155,6 +171,18 @@ private:
     int transCriterion_;
 };
 
+/*! \brief Compares two transmissibility matrices according to a L-selection criterion
+ *
+ * Compares two transmissibility matrices according to the L-selection criterion which is chosen via the parameter/property
+ * MPFATransmissibilityCriterion (Criterion of [1], 1: Criterion of [2]) and returns the number of the preferable L-shape (1-4).
+ *
+ * \param transmissibilityOne A first transmissibility matrix
+ * \param transmissibilityTwo A second transmissibility matrix
+ * \param lTypeOne Type of the L-shape of the first matrix (1-4)
+ * \param lTypeTwo Type of the L-shape of the second matrix (1-4)
+ *
+ * \return The type of the preferable L-shape (1-4)
+ */
 template<class TypeTag>
 int FvMpfaL3dTransmissibilityCalculator<TypeTag>::chooseTransmissibility(TransmissibilityType& transmissibilityOne, TransmissibilityType& transmissibilityTwo, int lTypeOne, int lTypeTwo)
 {
@@ -218,6 +246,20 @@ int FvMpfaL3dTransmissibilityCalculator<TypeTag>::chooseTransmissibility(Transmi
 
 }
 
+/*! \brief Calculates the transmissibility matrix for the flux face between the cells with the local index idx1 and idx2
+ *
+ *  \param transmissibility Reference to the local transmissibility matrix
+ *  \param interactionVolume An interaction volume object
+ *  \param lambda A vector containing the mobilities of the interaction volume cells in the order of the local interaction volume index
+ *  \param idx1 Local index of cell 1 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa)
+ *  \param idx2 Local index of cell 2 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa)
+ *  \param idx3 Local index of cell 3 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa)
+ *  \param idx4 Local index of cell 4 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa)
+ *  \param idx5 Local index of cell 5 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa)
+ *  \param idx6 Local index of cell 6 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa)
+ *
+ *  \return The type of the chosen L-shape (1-4)
+ */
 template<class TypeTag>
 int FvMpfaL3dTransmissibilityCalculator<TypeTag>::transmissibility(
                                                                    TransmissibilityType& transmissibility, InteractionVolume& interactionVolume,
@@ -300,6 +342,21 @@ int FvMpfaL3dTransmissibilityCalculator<TypeTag>::transmissibility(
     }
 }
 
+/*! \brief Calculates the transmissibility matrix for the flux face between the cells with the local index idx1 and idx2
+ *
+ *  \param transmissibility Reference to the local transmissibility matrix
+ *  \param interactionVolume An interaction volume object
+ *  \param lambda A vector containing the mobilities of the interaction volume cells in the order of the local interaction volume index
+ *  \param idx1 Local index of cell 1 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa)
+ *  \param idx2 Local index of cell 2 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa)
+ *  \param idx3 Local index of cell 3 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa)
+ *  \param idx4 Local index of cell 4 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa)
+ *  \param idx5 Local index of cell 5 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa)
+ *  \param idx6 Local index of cell 6 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa)
+ *  \param useCases Vector of enabling/disabling single L-shape types (Can be necessary for certain configurations in the case of hanging nodes)
+ *
+ *  \return The type of the chosen L-shape (1-4)
+ */
 template<class TypeTag>
 int FvMpfaL3dTransmissibilityCalculator<TypeTag>::transmissibility(
                                                                    TransmissibilityType& transmissibility, InteractionVolume& interactionVolume,
@@ -454,6 +511,15 @@ int FvMpfaL3dTransmissibilityCalculator<TypeTag>::transmissibility(
     }
 }
 
+/*! \brief Calculates a TPFA transmissibility matrix for the flux face between the cells with the local index idx1 and idx2
+ *
+ *  \param transmissibility Reference to the local transmissibility matrix
+ *  \param interactionVolume An interaction volume object
+ *  \param lambda A vector containing the mobilities of the interaction volume cells in the order of the local interaction volume index
+ *  \param idx1 Local index of cell 1 for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa)
+ *  \param idx2 Local index of cell 2 for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa)
+ *
+ */
 template<class TypeTag>
 int FvMpfaL3dTransmissibilityCalculator<TypeTag>::transmissibilityTPFA(
                                                                        Dune::FieldMatrix<Scalar,dim,2*dim-dim+1>& transmissibility,
@@ -542,6 +608,27 @@ int FvMpfaL3dTransmissibilityCalculator<TypeTag>::transmissibilityTPFA(
     return 1;
 }
 
+/*! \brief Calculates the transmissibility matrix for the flux face between the cells with the local index idx1 and idx2 using the L-shape type 1
+ *
+ * For more details on L-shape type 1 see:
+ *
+ *  Wolff 2013: http://elib.uni-stuttgart.de/opus/volltexte/2013/8661/, or
+ *
+ *  M. Wolff, Y. Cao, B. Flemisch, R. Helmig, and B. Wohlmuth (2013a). Multi-point flux
+ * approximation L-method in 3D: numerical convergence and application to two-phase
+ * flow through porous media. In P. Bastian, J. Kraus, R. Scheichl, and M. Wheeler,
+ * editors, Simulation of Flow in Porous Media - Applications in Energy and Environment. De Gruyter.
+ *
+ *  \param transmissibility Reference to the local transmissibility matrix
+ *  \param interactionVolume An interaction volume object
+ *  \param lambda A vector containing the mobilities of the interaction volume cells in the order of the local interaction volume index
+ *  \param idx1 Local index of cell 1 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa)
+ *  \param idx2 Local index of cell 2 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa)
+ *  \param idx3 Local index of cell 3 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa
+ *  \param idx5 Local index of cell 5 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa)
+ *
+ *  \return 1 (L-shape 1)
+ */
 template<class TypeTag>
 int FvMpfaL3dTransmissibilityCalculator<TypeTag>::transmissibilityCaseOne(
                                                                           Dune::FieldMatrix<Scalar,dim,2*dim-dim+1>& transmissibility,
@@ -906,6 +993,27 @@ int FvMpfaL3dTransmissibilityCalculator<TypeTag>::transmissibilityCaseOne(
     return 1;
 }
 
+/*! \brief Calculates the transmissibility matrix for the flux face between the cells with the local index idx1 and idx2 using the L-shape type 2
+ *
+ * For more details on L-shape type 2 see:
+ *
+ *  Wolff 2013: http://elib.uni-stuttgart.de/opus/volltexte/2013/8661/, or
+ *
+ *  M. Wolff, Y. Cao, B. Flemisch, R. Helmig, and B. Wohlmuth (2013a). Multi-point flux
+ * approximation L-method in 3D: numerical convergence and application to two-phase
+ * flow through porous media. In P. Bastian, J. Kraus, R. Scheichl, and M. Wheeler,
+ * editors, Simulation of Flow in Porous Media - Applications in Energy and Environment. De Gruyter.
+ *
+ *  \param transmissibility Reference to the local transmissibility matrix
+ *  \param interactionVolume An interaction volume object
+ *  \param lambda A vector containing the mobilities of the interaction volume cells in the order of the local interaction volume index
+ *  \param idx1 Local index of cell 1 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa)
+ *  \param idx2 Local index of cell 2 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa)
+ *  \param idx4 Local index of cell 4 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa
+ *  \param idx6 Local index of cell 6 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa)
+ *
+ *  \return 2 (L-shape 2)
+ */
 template<class TypeTag>
 int FvMpfaL3dTransmissibilityCalculator<TypeTag>::transmissibilityCaseTwo(
                                                                           Dune::FieldMatrix<Scalar,dim,2*dim-dim+1>& transmissibility,
@@ -1283,6 +1391,27 @@ int FvMpfaL3dTransmissibilityCalculator<TypeTag>::transmissibilityCaseTwo(
     return 2;
 }
 
+/*! \brief Calculates the transmissibility matrix for the flux face between the cells with the local index idx1 and idx2 using the L-shape type 3
+ *
+ * For more details on L-shape type 3 see:
+ *
+ *  Wolff 2013: http://elib.uni-stuttgart.de/opus/volltexte/2013/8661/, or
+ *
+ *  M. Wolff, Y. Cao, B. Flemisch, R. Helmig, and B. Wohlmuth (2013a). Multi-point flux
+ * approximation L-method in 3D: numerical convergence and application to two-phase
+ * flow through porous media. In P. Bastian, J. Kraus, R. Scheichl, and M. Wheeler,
+ * editors, Simulation of Flow in Porous Media - Applications in Energy and Environment. De Gruyter.
+ *
+ *  \param transmissibility Reference to the local transmissibility matrix
+ *  \param interactionVolume An interaction volume object
+ *  \param lambda A vector containing the mobilities of the interaction volume cells in the order of the local interaction volume index
+ *  \param idx1 Local index of cell 1 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa)
+ *  \param idx2 Local index of cell 2 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa)
+ *  \param idx4 Local index of cell 4 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa
+ *  \param idx5 Local index of cell 5 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa)
+ *
+ *  \return 3 (L-shape 3)
+ */
 template<class TypeTag>
 int FvMpfaL3dTransmissibilityCalculator<TypeTag>::transmissibilityCaseThree(
                                                                             Dune::FieldMatrix<Scalar,dim,2*dim-dim+1>& transmissibility,
@@ -1661,6 +1790,27 @@ int FvMpfaL3dTransmissibilityCalculator<TypeTag>::transmissibilityCaseThree(
     return 3;
 }
 
+/*! \brief Calculates the transmissibility matrix for the flux face between the cells with the local index idx1 and idx2 using the L-shape type 4
+ *
+ * For more details on L-shape type 4 see:
+ *
+ *  Wolff 2013: http://elib.uni-stuttgart.de/opus/volltexte/2013/8661/, or
+ *
+ *  M. Wolff, Y. Cao, B. Flemisch, R. Helmig, and B. Wohlmuth (2013a). Multi-point flux
+ * approximation L-method in 3D: numerical convergence and application to two-phase
+ * flow through porous media. In P. Bastian, J. Kraus, R. Scheichl, and M. Wheeler,
+ * editors, Simulation of Flow in Porous Media - Applications in Energy and Environment. De Gruyter.
+ *
+ *  \param transmissibility Reference to the local transmissibility matrix
+ *  \param interactionVolume An interaction volume object
+ *  \param lambda A vector containing the mobilities of the interaction volume cells in the order of the local interaction volume index
+ *  \param idx1 Local index of cell 1 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa)
+ *  \param idx2 Local index of cell 2 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa)
+ *  \param idx3 Local index of cell 3 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa
+ *  \param idx6 Local index of cell 6 of 6 cells for the calculation of the flux between cell idx1 and idx2  (see doc/docextra/3dmpfa)
+ *
+ *  \return 4 (L-shape 4)
+ */
 template<class TypeTag>
 int FvMpfaL3dTransmissibilityCalculator<TypeTag>::transmissibilityCaseFour(
                                                                            Dune::FieldMatrix<Scalar,dim,2*dim-dim+1>& transmissibility,
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dvelocity2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dvelocity2p.hh
index 65ea09d18bfcf23d766509c5bc25b44fb89b2a93..96d3017f018aa32e512843456f77b4f6538042ed 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dvelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dvelocity2p.hh
@@ -28,11 +28,29 @@
 
 /**
  * @file
- * @brief  Base class for defining an instance of a numerical diffusion model
+ * @brief 2-d velocity calculation using a 3-d MPFA L-method
  */
 
 namespace Dumux
 {
+//! \ingroup FVPressure2p
+/*! \brief Class for calculating 3-d velocities from cell-wise constant pressure values.
+ * Calculates phase velocities or total velocity from a known pressure field applying a finite volume discretization and a MPFA L-method.
+ * At Dirichlet boundaries a two-point flux approximation is used.
+ * The pressure has to be given as piecewise constant cell values.
+ * The velocities are calculated as
+ *
+ *\f[ \boldsymbol v_\alpha = - \lambda_\alpha \boldsymbol K \textbf{grad}\, \Phi_\alpha, \f]
+ * and,
+ * \f[ \boldsymbol v_t = \boldsymbol v_w + \boldsymbol v_n,\f]
+ *
+ * where \f$ \Phi_\alpha \f$ denotes the potential of phase \f$ \alpha \f$, \f$ \boldsymbol K \f$ the intrinsic permeability,
+ * and \f$ \lambda_\alpha \f$ a phase mobility.
+ *
+ * Remark1: only for 3-D Hexahedrons of quadrilaterals!
+ *
+ * \tparam TypeTag The problem Type Tag
+ */
 template<class TypeTag> class FvMpfaL3dVelocity2p
 {
     typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
@@ -129,23 +147,29 @@ template<class TypeTag> class FvMpfaL3dVelocity2p
     typedef Dune::FieldVector<Scalar, dim> DimVector;
 
 public:
+    //! Constructs a FvMpfaL3dVelocity2p object
+    /*!
+     * \param problem A problem class object
+     */
     FvMpfaL3dVelocity2p(Problem& problem) :
         problem_(problem), gravity_(problem.gravity())
-{
+    {
         density_[wPhaseIdx] = 0.;
         density_[nPhaseIdx] = 0.;
         viscosity_[wPhaseIdx] = 0.;
         viscosity_[nPhaseIdx] = 0.;
 
         vtkOutputLevel_ = GET_PARAM_FROM_GROUP(TypeTag, int, Vtk, OutputLevel);
-}
+    }
 
+    //calculate velocities for flux faces of an interaction volume
     void calculateInnerInteractionVolumeVelocity(InteractionVolume& interactionVolume,
             CellData & cellData1,  CellData & cellData2, CellData & cellData3, CellData & cellData4,
             CellData & cellData5, CellData & cellData6, CellData & cellData7, CellData & cellData8,
             InteractionVolumeContainer& interactionVolumes, TransmissibilityCalculator& transmissibilityCalculator, int faceIdx = -1);
     void calculateBoundaryInteractionVolumeVelocity(InteractionVolume& interactionVolume, CellData& cellData, int elemIdx);
 
+    //!Initializes the velocity model
     void initialize(bool solveTwice = true)
     {
         ElementIterator element = problem_.gridView().template begin<0>();
@@ -163,8 +187,16 @@ public:
         return;
     }
 
-    //! \brief Write data files
-    /*  \param name file name */
+    /*! \brief Adds velocity output to the output file
+     *
+     * Adds the phase velocities or a total velocity (depending on the formulation) to the output.
+     * If the VtkOutputLevel is equal to zero (default) only primary variables are written,
+     * if it is larger than zero also secondary variables are written.
+     *
+     * \tparam MultiWriter Class defining the output writer
+     * \param writer The output writer (usually a <tt>VTKMultiWriter</tt> object)
+     *
+     */
     template<class MultiWriter>
     void addOutputVtkFields(MultiWriter &writer)
     {
@@ -252,6 +284,24 @@ private:
 };
 // end of template
 
+/*! \brief Calculates the velocities at the flux faces of an interation volume around a vertex which is not a boundary vertex.
+ *
+ *  Calculates the velocities at the flux faces of an interation volume around a vertex which is not a boundary vertex and adds them to the face velocity vectors in the <tt>CellData</tt> objects.
+ *
+ * \param interactionVolume An <tt>InteractionVolume</tt> object including the information for calculating the MPFA transmissibilities
+ * \param cellData1  <tt>CellData</tt> object of an IMPES model for sub-volume 1
+ * \param cellData2  <tt>CellData</tt> object of an IMPES model for sub-volume 2
+ * \param cellData3  <tt>CellData</tt> object of an IMPES model for sub-volume 3
+ * \param cellData4  <tt>CellData</tt> object of an IMPES model for sub-volume 4
+ * \param cellData5  <tt>CellData</tt> object of an IMPES model for sub-volume 5
+ * \param cellData6  <tt>CellData</tt> object of an IMPES model for sub-volume 6
+ * \param cellData7  <tt>CellData</tt> object of an IMPES model for sub-volume 7
+ * \param cellData8  <tt>CellData</tt> object of an IMPES model for sub-volume 8
+ * \param interactionVolumes Container including the interaction volume information for the complete grid
+ * \param transmissibilityCalculator Object including the methods for calculating the transmissibilities
+ * \param faceIdx Index of the flux face for which the velocity has to be calculated. If no face index is given, <tt>faceIdx</tt> = -1
+ * and velocities for all flux faces in the interaction volume are calculated!
+ */
 template<class TypeTag>
 void FvMpfaL3dVelocity2p<TypeTag>::calculateInnerInteractionVolumeVelocity(InteractionVolume& interactionVolume,
         CellData & cellData1,  CellData & cellData2, CellData & cellData3, CellData & cellData4,
@@ -1868,6 +1918,14 @@ void FvMpfaL3dVelocity2p<TypeTag>::calculateInnerInteractionVolumeVelocity(Inter
     cellData8.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(7, 2));
         }
 
+/*! \brief Calculates the velocity at a boundary flux faces.
+ *
+ *  Calculates the velocity at a boundary flux face and adds it to the face velocity vector in the <tt>CellData</tt> object.
+ *
+ * \param interactionVolume An <tt>InteractionVolume</tt> object including the information for calculating the MPFA transmissibilities
+ * \param cellData  <tt>CellData</tt> object of an IMPES model for the sub-volume including the boundary face
+ * \param elemIdx local sub-volume index
+ */
 template<class TypeTag>
 void FvMpfaL3dVelocity2p<TypeTag>::calculateBoundaryInteractionVolumeVelocity(InteractionVolume& interactionVolume, CellData& cellData, int elemIdx)
 {
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dvelocity2padaptive.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dvelocity2padaptive.hh
index 9c68a75dd4fe22fb313dab449565926c842a33ca..7fc0c520a3ba851b95c6a1768b5ce3a72bf7d10d 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dvelocity2padaptive.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dvelocity2padaptive.hh
@@ -24,12 +24,31 @@
 
 /**
  * @file
- * @brief  Base class for defining an instance of a numerical diffusion model
+ * @brief  3-d velocity calculation on adaptive grids using a 3-d MPFA L-method
  */
 
 namespace Dumux
 {
-
+//! \ingroup FVPressure2p
+/*! \brief Class for calculating 3-d velocities from cell-wise constant pressure values.
+ * Calculates phase velocities or total velocity from a known pressure field applying a finite volume discretization and a MPFA L-method.
+ * At Dirichlet boundaries a two-point flux approximation is used.
+ * The pressure has to be given as piecewise constant cell values.
+ * The velocities are calculated as
+ *
+ *\f[ \boldsymbol v_\alpha = - \lambda_\alpha \boldsymbol K \textbf{grad}\, \Phi_\alpha, \f]
+ * and,
+ * \f[ \boldsymbol v_t = \boldsymbol v_w + \boldsymbol v_n,\f]
+ *
+ * where \f$ \Phi_\alpha \f$ denotes the potential of phase \f$ \alpha \f$, \f$ \boldsymbol K \f$ the intrinsic permeability,
+ * and \f$ \lambda_\alpha \f$ a phase mobility.
+ *
+ * Remark1: only for 3-D Hexahedrons of quadrilaterals!
+ *
+ * Remark2: Allowed difference in grid levels of two neighboring cells: 1
+ *
+ * \tparam TypeTag The problem Type Tag
+ */
 template<class TypeTag> class FvMpfaL3dVelocity2pAdaptive: public FvMpfaL3dVelocity2p<TypeTag>
 {
     typedef FvMpfaL3dVelocity2p<TypeTag> ParentType;
@@ -134,6 +153,10 @@ template<class TypeTag> class FvMpfaL3dVelocity2pAdaptive: public FvMpfaL3dVeloc
     }
 
 public:
+    //! Constructs a FvMpfaL3dVelocity2pAdaptive object
+    /*!
+     * \param problem A problem class object
+     */
     FvMpfaL3dVelocity2pAdaptive(Problem& problem) :
         ParentType(problem), problem_(problem), gravity_(problem.gravity())
 {
@@ -143,11 +166,12 @@ public:
         viscosity_[nPhaseIdx] = 0.;
 }
 
+    //calculate velocities for flux faces of a hanging node interaction volume
     void calculateHangingNodeInteractionVolumeVelocity(InteractionVolume& interactionVolume,
             CellData & cellData1,  CellData & cellData2, CellData & cellData3, CellData & cellData4,
             CellData & cellData5, CellData & cellData6, CellData & cellData7, CellData & cellData8, int faceIdx = -1);
 
-
+    //!Initializes the velocity model
     void initialize()
     {
         ParentType::initialize();
@@ -181,7 +205,22 @@ private:
 };
 // end of template
 
-// only for 3-D general hexahedron
+/*! \brief Calculates the velocities at the flux faces of an interation volume around a hanging node vertex.
+ *
+ *  Calculates the velocities at the flux faces of an interation volume around  a hanging node vertex and adds them to the face velocity vectors in the <tt>CellData</tt> objects.
+ *
+ * \param interactionVolume An <tt>InteractionVolume</tt> object including the information for calculating the MPFA transmissibilities
+ * \param cellData1  <tt>CellData</tt> object of an IMPES model for sub-volume 1
+ * \param cellData2  <tt>CellData</tt> object of an IMPES model for sub-volume 2
+ * \param cellData3  <tt>CellData</tt> object of an IMPES model for sub-volume 3
+ * \param cellData4  <tt>CellData</tt> object of an IMPES model for sub-volume 4
+ * \param cellData5  <tt>CellData</tt> object of an IMPES model for sub-volume 5
+ * \param cellData6  <tt>CellData</tt> object of an IMPES model for sub-volume 6
+ * \param cellData7  <tt>CellData</tt> object of an IMPES model for sub-volume 7
+ * \param cellData8  <tt>CellData</tt> object of an IMPES model for sub-volume 8
+ * \param faceIdx Index of the flux face for which the velocity has to be calculated. If no face index is given, <tt>faceIdx</tt> = -1
+ * and velocities for all flux faces in the interaction volume are calculated!
+ */
 template<class TypeTag>
 void FvMpfaL3dVelocity2pAdaptive<TypeTag>::calculateHangingNodeInteractionVolumeVelocity(InteractionVolume& interactionVolume,
         CellData & cellData1,  CellData & cellData2, CellData & cellData3, CellData & cellData4,
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2dpressure2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2dpressure2p.hh
index bfbf29d6246405dcb9522ef0e3511f7a35135733..331ffcf911ce4ad34b8bef5a407207ce51216dc0 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2dpressure2p.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2dpressure2p.hh
@@ -43,8 +43,16 @@ namespace Dumux
  *  Here, \f$ \Phi_\alpha \f$ denotes the potential of phase \f$ \alpha \f$, \f$ \boldsymbol K \f$ the intrinsic permeability,
  * \f$ \lambda_t \f$ the total mobility, \f$ f_\alpha \f$ the phase fractional flow function.
  *
- * More details on the equations can be found in H. Hoteit, A. Firoozabadi.
- * Numerical modeling of thwo-phase flow in heterogeneous permeable media with different capillarity pressures. Adv Water Res (31), 2008
+ * More details on the equations can be found in
+ *
+ * M. Wolff, Y. Cao, B. Flemisch, R. Helmig, and B. Wohlmuth (2013a). Multi-point flux
+ * approximation L-method in 3D: numerical convergence and application to two-phase
+ * flow through porous media. In P. Bastian, J. Kraus, R. Scheichl, and M. Wheeler,
+ * editors, Simulation of Flow in Porous Media - Applications in Energy and Environment. De Gruyter.
+ *
+ * M. Wolff, B. Flemisch, R. Helmig, I. Aavatsmark.
+ * Treatment of tensorial relative permeabilities with multipoint flux approximation.
+ * International Journal of Numerical Analysis and Modeling (9), pp. 725-744, 2012.
  *
  *  Remark1: only for 2-D quadrilateral grid
  *
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2dpressurevelocity2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2dpressurevelocity2p.hh
index 27156a2b8b3ca5e97426e45061fa77f96254039e..0854e4ffdcdff3ac8d5df48dcee248bdb57b10c4 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2dpressurevelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2dpressurevelocity2p.hh
@@ -30,18 +30,10 @@ namespace Dumux
 {
 
 //! \ingroup FVPressure2p
-/*! \brief Determines the velocity from a finite volume solution of the  pressure equation of a sequential model (IMPES).
- * Calculates phase velocities or total velocity from a known pressure field applying a finite volume discretization and a MPFA O-method.
- * At Dirichlet boundaries a two-point flux approximation is used.
- * The pressure has to be given as piecewise constant cell values.
- * The velocities are calculated as
+/*! \brief Class for the calculation of velocities from the  pressure solution of an IMPES scheme using a MPFA O-method.
  *
- *\f[ \boldsymbol v_\alpha = - \lambda_\alpha \boldsymbol K \textbf{grad}\, \Phi_\alpha, \f]
- * and,
- * \f[ \boldsymbol v_t = \boldsymbol v_w + \boldsymbol v_n,\f]
- *
- * where \f$ \Phi_\alpha \f$ denotes the potential of phase \f$ \alpha \f$, \f$ \boldsymbol K \f$ the intrinsic permeability,
- * and \f$ \lambda_\alpha \f$ a phase mobility.
+ * Can be used for calculating the complete velocity field before the solution of the transport equation (more efficient),
+ * or for face-wise velocity calculation directly in the transport solution (less efficient).
  *
  * Remark1: only for 2-D quadrilateral grids!
  *
@@ -134,6 +126,7 @@ public:
     void calculateVelocity(const Intersection&, CellData&);
     void calculateVelocityOnBoundary(const Intersection& intersection, CellData& cellData);
 
+    //! Function for updating the velocity field if iterations are necessary in the transport solution
     void updateVelocity()
     {
         this->updateMaterialLaws();
@@ -208,6 +201,11 @@ public:
             velocity_.addOutputVtkFields(writer);
     }
 
+    FvMpfaO2dVelocity2P<TypeTag> velocity()
+    {
+        return velocity_;
+    }
+
 private:
     Problem& problem_;
     FvMpfaO2dVelocity2P<TypeTag> velocity_;
@@ -219,9 +217,9 @@ private:
 };
 // end of template
 
-/*! \brief Calculates the velocities at a cell-cell interfaces.
+/*! \brief Calculates the velocities at a cell-cell interfaces for the entire simulation grid.
  *
- * Calculates the velocities at a cell-cell interfaces from a given pressure field.
+ * Calculates the velocities at a cell-cell interfaces from a given pressure field for the entire simulation grid.
  *
  */
 template<class TypeTag>
@@ -291,6 +289,13 @@ void FvMpfaO2dPressureVelocity2p<TypeTag>::calculateVelocity()
     return;
 } // end method calcTotalVelocity
 
+/*! \brief Calculates the velocity at a cell-cell interface.
+ *
+ * Calculates the velocity at a cell-cell interface from a given pressure field.
+ *
+ * \param intersection Intersection of two grid cells
+ * \param cellData Object containing all model relevant cell data
+ */
 template<class TypeTag>
 void FvMpfaO2dPressureVelocity2p<TypeTag>::calculateVelocity(const Intersection& intersection, CellData& cellData)
 {
@@ -362,10 +367,19 @@ void FvMpfaO2dPressureVelocity2p<TypeTag>::calculateVelocity(const Intersection&
         }
     }
     }
+
+
     cellData.fluxData().setVelocityMarker(indexInInside);
     cellDataJ.fluxData().setVelocityMarker(indexInOutside);
 }
 
+/*! \brief Calculates the velocity at a boundary.
+ *
+ * Calculates the velocity at a boundary from a given pressure field.
+ *
+ * \param intersection Boundary intersection
+ * \param cellData Object containing all model relevant cell data
+ */
 template<class TypeTag>
 void FvMpfaO2dPressureVelocity2p<TypeTag>::calculateVelocityOnBoundary(const Intersection& intersection, CellData& cellData)
 {
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2dvelocity2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2dvelocity2p.hh
index 914fff2d425eed268c56cdc9fe8ad9c63268e9bd..6378c3f402b500d928fdde5e76e742323ebc94b9 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2dvelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2dvelocity2p.hh
@@ -32,6 +32,24 @@
 
 namespace Dumux
 {
+//! \ingroup FVPressure2p
+/*! \brief Class for calculating velocities from cell-wise constant pressure values.
+ * Calculates phase velocities or total velocity from a known pressure field applying a finite volume discretization and a MPFA O-method.
+ * At Dirichlet boundaries a two-point flux approximation is used.
+ * The pressure has to be given as piecewise constant cell values.
+ * The velocities are calculated as
+ *
+ *\f[ \boldsymbol v_\alpha = - \lambda_\alpha \boldsymbol K \textbf{grad}\, \Phi_\alpha, \f]
+ * and,
+ * \f[ \boldsymbol v_t = \boldsymbol v_w + \boldsymbol v_n,\f]
+ *
+ * where \f$ \Phi_\alpha \f$ denotes the potential of phase \f$ \alpha \f$, \f$ \boldsymbol K \f$ the intrinsic permeability,
+ * and \f$ \lambda_\alpha \f$ a phase mobility.
+ *
+ * Remark1: only for 2-D quadrilateral grids!
+ *
+ * \tparam TypeTag The problem Type Tag
+ */
 template<class TypeTag> class FvMpfaO2dVelocity2P
 {
     typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
@@ -112,6 +130,10 @@ template<class TypeTag> class FvMpfaO2dVelocity2P
     typedef Dune::FieldVector<Scalar, dim> DimVector;
 
 public:
+    //! Constructs a FvMpfaO2dVelocity2P object
+    /*!
+     * \param problem A problem class object
+     */
     FvMpfaO2dVelocity2P(Problem& problem) :
         problem_(problem), gravity_(problem.gravity())
     {
@@ -123,9 +145,11 @@ public:
         vtkOutputLevel_ = GET_PARAM_FROM_GROUP(TypeTag, int, Vtk, OutputLevel);
     }
 
+    //calculate velocities for all flux faces of an interaction volume
     void calculateInnerInteractionVolumeVelocity(InteractionVolume& interactionVolume, CellData& cellData1, CellData& cellData2, CellData& cellData3, CellData& cellData4, InnerBoundaryVolumeFaces& innerBoundaryVolumeFaces);
     void calculateBoundaryInteractionVolumeVelocity(InteractionVolume& interactionVolume, CellData& cellData, int elemIdx);
 
+    //!Initializes the velocity model
     void initialize()
     {
         ElementIterator element = problem_.gridView().template begin<0>();
@@ -143,8 +167,16 @@ public:
         return;
     }
 
-    //! \brief Write data files
-    /*  \param name file name */
+    /*! \brief Adds velocity output to the output file
+     *
+     * Adds the phase velocities or a total velocity (depending on the formulation) to the output.
+     * If the VtkOutputLevel is equal to zero (default) only primary variables are written,
+     * if it is larger than zero also secondary variables are written.
+     *
+     * \tparam MultiWriter Class defining the output writer
+     * \param writer The output writer (usually a <tt>VTKMultiWriter</tt> object)
+     *
+     */
     template<class MultiWriter>
     void addOutputVtkFields(MultiWriter &writer)
     {
@@ -230,6 +262,17 @@ private:
 };
 // end of template
 
+/*! \brief Calculates the velocities at the flux faces of an interation volume around a vertex which is not a boundary vertex.
+ *
+ *  Calculates the velocities at the flux faces of an interation volume around a vertex which is not a boundary vertex and adds them to the face velocity vectors in the <tt>CellData</tt> objects.
+ *
+ * \param interactionVolume An <tt>InteractionVolume</tt> object including the information for calculating the MPFA transmissibilities
+ * \param cellData1  <tt>CellData</tt> object of an IMPES model for sub-volume 1
+ * \param cellData2  <tt>CellData</tt> object of an IMPES model for sub-volume 2
+ * \param cellData3  <tt>CellData</tt> object of an IMPES model for sub-volume 3
+ * \param cellData4  <tt>CellData</tt> object of an IMPES model for sub-volume 4
+ * \param innerBoundaryVolumeFaces container including information about faces intersecting a boundary
+ */
 template<class TypeTag>
 void FvMpfaO2dVelocity2P<TypeTag>::calculateInnerInteractionVolumeVelocity(InteractionVolume& interactionVolume, CellData& cellData1, CellData& cellData2, CellData& cellData3, CellData& cellData4, InnerBoundaryVolumeFaces& innerBoundaryVolumeFaces)
 {
@@ -513,6 +556,14 @@ void FvMpfaO2dVelocity2P<TypeTag>::calculateInnerInteractionVolumeVelocity(Inter
     cellData4.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(3, 1));
 }
 
+/*! \brief Calculates the velocity at a boundary flux faces.
+ *
+ *  Calculates the velocity at a boundary flux face and adds it to the face velocity vector in the <tt>CellData</tt> object.
+ *
+ * \param interactionVolume An <tt>InteractionVolume</tt> object including the information for calculating the MPFA transmissibilities
+ * \param cellData  <tt>CellData</tt> object of an IMPES model for the sub-volume including the boundary face
+ * \param elemIdx local sub-volume index
+ */
 template<class TypeTag>
 void FvMpfaO2dVelocity2P<TypeTag>::calculateBoundaryInteractionVolumeVelocity(InteractionVolume& interactionVolume, CellData& cellData, int elemIdx)
 {