diff --git a/dune/fufem/constraints/affineconstraints.hh b/dune/fufem/constraints/affineconstraints.hh
index 4310053efee00d3933809a71ea899770d2a169fa..4825268b732f064417b384627f57ba07d5cdcd4e 100644
--- a/dune/fufem/constraints/affineconstraints.hh
+++ b/dune/fufem/constraints/affineconstraints.hh
@@ -27,40 +27,55 @@ namespace Dune::Fufem {
  * \brief Class to handle affine constraints on a subset of DOFs
  * \ingroup Constraints
  *
- * This class handles affine constraints on a subset of DOFs of the form
+ * This class handles affine constraints of the form
  *
- * \f[ x_i = \sum_{j} L_{ij} x_j + c_i\f]
+ * \f[ x_i = \sum_{j\in \mathcal{I}} L_{ij} x_j + c_i \qquad i \in \mathcal{I}\f]
  *
- * or in short \f$ x = L x + c\f$ for a square matrix
- * \f$ L  \in \mathbb{R}^{n \times n} \f$ and a vector \f$c \in \mathbb{R}^n\f$.
+ * or in short \f$ x = L x + c\f$ for a square matrix \f$ L \f$
+ * where all indices are from the index set \f$\mathcal{I}\f$.
+ * For exposition purposes we assume that \f$ \mathcal{I}=\{1,...,n\} \f$
+ * while the implementation supports multi-indices as provided by the
+ * global bases in dune-functions.
+ * Under this assumption we have \f$L \in \mathbb{R}^{n \times n} \f$ and \f$c \in \mathbb{R}^n\f$.
  * The main purpose is to simplify the handling
  * of problems in the constrained affine subspace
  *
  * \f[ V_{L,c} = \{ x\,|\, x = L x+c\}.\f]
  *
- * We assume that the indices \f$ i \f$ can be classified
- * into non-constrained ones, where (using the Kronecker delta)
- * \f$ L_{ij} = \delta_{ij} \forall j \f$
- * and \f$c_i=0\f$ and constrained ones, where \f$L_{ii}=L_{ij}=0 \f$
- * if \f$j\f$ is another constrained index. Then any constrained \f$x_i\f$
- * can be interpolated from \f$c_i\f$ and \f$ x_j \f$
- * for the non-constraints indices \f$j\neq i\f$.
- * Under these assumptions we have \f$L^2 = L\f$ and \f$L c=0\f$
- * such that the affine map \f$ \Psi:\mathbb{R}^n \to V_{L,c}\f$
- * given by \f$v \mapsto \Psi(v) = L v +c\f$ is surjective projection
- * into \f$V_{L,c}\f$.
+ * We assume that the index set can be decomposed
+ * according to \f$\mathcal{I}= \mathcal{I}_C \cup \mathcal{I}_N\f$ into disjoint subsets
+ * \f$\mathcal{I}_C\f$ and \f$\mathcal{I}_N\f$ of constrained and non-constrained indices, respectively.
+ * The constrained and non-constrained indices should satisfy the following two assumptions:
+ * (a) For any non-constrained index \f$i\in \mathcal{I}_N\f$ we have (using the Kronecker delta)
+ * \f$ L_{ij} = \delta_{ij} \forall j\in\mathcal{I} \f$ and \f$c_i=0\f$.
+ * (b) For any constrained index \f$i\in \mathcal{I}_C\f$ we have
+ * \f$L_{ii}=L_{ij}=0 \forall j\in\mathcal{I}_C\f$.
+ * Then any constrained value \f$x_i\f$ with \f$i\in\mathcal{I}_C\f$
+ * can be interpolated from \f$c_i\f$ and the non-constrained values \f$x_j\f$
+ * according to
  *
- * Consider a linear variational problem
+ * \f[ x_i = \sum_{i\in \mathcal{I}} L_{ij} + c_i = \sum_{i\in \mathcal{I}_N} L_{ij} + c_i. \f]
+ *
+ * Under these assumptions we furthermore have \f$L^2 = L\f$ and \f$L c=0\f$,
+ * i.e. \f$L\f$ is a linear projection and \f$c\f$ is in the kernel of \f$L\f$.
+ * Then the affine map \f$ \Psi:\mathbb{R}^n \to V_{L,c}\f$
+ * given by \f$v \mapsto \Psi(v) = L v +c\f$ is a surjective projection
+ * into \f$V_{L,c}\f$. Hence the set \f$V_{L,c}\f$ of all \f$x\f$ satisfying
+ * the constraint \f$x = Lx+c\f$ is equal to the set of all \f$x\f$ that
+ * can be written as \f$x = Lv+c\f$ for some \f$v\f$, i.e.
+ * \f$ V_{L,c} = \{ x\,|\, x = L x +c \} = \{Lv +c \,|\, v\}\f$.
+ *
+ * Now consider a linear variational problem
  *
  * \f[ u \in V_{L,c}: \qquad \langle A u, v \rangle = \langle b, v \rangle \qquad \forall v \in V_{L,0} \f]
  *
- * in the affine subspace \f$ V_{L,c} = \{ L x +c \,|\, x\in \mathbb{R}^n\}\f$.
- * Using  of the decomposition \f$u= u_0 + c\f$ with \f$ u_0 = L u_0 = L u \in V_{L,0} \f$
- * this can then be written as
+ * in the affine subspace \f$ V_{L,c}\f$.
+ * Using the decomposition \f$u= u_0 + c\f$ with \f$ u_0 = L u_0 = L u \in V_{L,0} \f$
+ * this can be written as
  *
- * \f[ u_0 \in V_{L,c}: \qquad \langle A u_0, v \rangle = \langle b - A c, v \rangle \qquad \forall v \in V_{L,0}. \f]
+ * \f[ u_0 \in V_{L,0}: \qquad \langle A u_0, v \rangle = \langle b - A c, v \rangle \qquad \forall v \in V_{L,0}. \f]
  *
- * which, because \f$L:V_{L,0} \to V_{L,0}\f$ is a linear projection, is equivalent to
+ * Since \f$L:V_{L,0} \to V_{L,0}\f$ is a linear projection, we find that \f$u_0\f$ solves
  *
  * \f[ u_0 \in \mathbb{R}^n: \qquad \langle L^T A L u_0, v \rangle = \langle L^T(b - A c), v \rangle \qquad \forall v \in \mathbb{R}^n \f]
  *
@@ -71,14 +86,16 @@ namespace Dune::Fufem {
  * \f$ker (A) \subset ker(L)\f$ we know
  * that the solution is at least unique in \f$V_{L,0}\f$.
  * To regularize the problem we consider the reduced identity matrix
- * \f$ \hat{I} \in \mathbb{R}^{n \times n}\f$ which is zero for all non-constrained indices.
- * Since \f$LI = L\f$ it is clear that the system
+ * \f$ \hat{I} \in \mathbb{R}^{n \times n}\f$
+ * with \f$\hat{I}_{ij} = \delta_{ij}\f$ if \f$i,j \in \mathcal{I}_C\f$
+ * and \f$\hat{I}_{ij} = 0\f$ if \f$i\in \mathcal{I}_N\f$ or \f$j\in \mathcal{I}_N\f$.
+ * Since \f$L\hat{I} = L\f$ it is clear that the system
  *
  * \f[ \hat{u} \in \mathbb{R}^n: \qquad \langle \hat{A}, v \rangle = \langle \hat{b}, v \rangle \qquad \forall v \in \mathbb{R}^n \f]
  *
  * with \f$ \hat{A} = L^T A L + \hat{I} \f$ and \f$\hat{b} = L^T(b-A c) + c\f$
- * is equivalent in the sense that \f$u_0 = L \hat{u}\f$. Hence we can solve
- * the constrained problem by solving
+ * is equivalent to the constrained problem in the sense that \f$u_0 = L \hat{u}\f$.
+ * Hence we can solve the constrained problem by solving
  *
  * \f[ \hat{A}\hat{u} = \hat{b} \f]
  *
@@ -195,7 +212,7 @@ public:
   /**
    * \brief Const access to i-th row of linear part of affine map
    *
-   * Assuming that the affine constraints is written as \f$ x = L x + c\f$
+   * Assuming that the affine constraints are written as \f$ x = L x + c\f$
    * the result represents the i-th row of the sparse matrix \f$L\f$.
    *
    * This method should only be used if the i-th DOF is constrained.