diff --git a/doc/userdoc.org b/doc/userdoc.org
index e9ba680359c371e480f525378931928127c1177d..54ea305564884ff6363a456990762adfcc378669 100644
--- a/doc/userdoc.org
+++ b/doc/userdoc.org
@@ -2781,6 +2781,20 @@ $\mathbb{R}^{1\times1}$, $\mathbb{R}^{2\times2}$ and $\mathbb{R}^{3\times3}$.
 The output is
 #+RESULTS: trace-examples
 
+The ~doubleDot~ functions compute the double-dot ($A:B = \mathrm{tr}(A^T
+B)$) of two matrices of $\mathbb{R}^{1\times1}$, $\mathbb{R}^{2\times2}$ and
+$\mathbb{R}^{3\times3}$.
+#+NAME: double-dot-examples
+#+BEGIN_SRC pugs :exports both :results output
+   import math;
+   cout << "doubleDot([[1.2]],[[2.3]]) = " << doubleDot([[1.2]],[[2.3]]) << "\n";
+   cout << "doubleDot([[1,2],[3,4]],[[3,6],[2,5]]) = " << doubleDot([[1,2],[3,4]],[[3,6],[2,5]]) << "\n";
+   cout << "doubleDot([[1,2,3],[4,5,6],[7,8,9]], [[5,2,1],[4,2,8],[2,6,2]]) = "
+        << doubleDot([[1,2,3],[4,5,6],[7,8,9]], [[5,2,1],[4,2,8],[2,6,2]]) << "\n";
+#+END_SRC
+The output is
+#+RESULTS: double-dot-examples
+
 Also, one can compute inverses of $\mathbb{R}^{1\times1}$,
 $\mathbb{R}^{2\times2}$ and $\mathbb{R}^{3\times3}$ matrices using the
 ~inverse~ function set.
@@ -2789,7 +2803,7 @@ $\mathbb{R}^{2\times2}$ and $\mathbb{R}^{3\times3}$ matrices using the
    import math;
    cout << "inverse([[1.2]]) = " << inverse([[1.2]]) << "\n";
    cout << "inverse([[1,2],[3,4]]) = " << inverse([[1,2],[3,4]]) << "\n";
-   cout << "inverse([[3,2,1],[5,6,4],[7,8,9]]) = "
+   cout << "inverse([[3,2,1],[5,6,4],[7,8,9]])\n = "
         << inverse([[3,2,1],[5,6,4],[7,8,9]]) << "\n";
 #+END_SRC
 The output is
@@ -3479,6 +3493,13 @@ This function requires that both arguments are defined on the same
 mesh and have the same data type. The result is a
 $\mathbb{P}_0(\mathbb{R})$ function.
 
+Finally the equivalent exists for discrete functions of matrices and
+applies to $\mathbb{P}_0(\mathbb{R}^1x1)$, $\mathbb{P}_0(\mathbb{R}^2x2)$,
+$\mathbb{P}_0(\mathbb{R}^3x3)$
+- ~doubleDot: Vh*Vh -> Vh~
+Both arguments must be defined on the same mesh and have the same data
+type. The result is a $\mathbb{P}_0(\mathbb{R})$ function.
+
 ****** ~R*Vh -> Vh~ and ~Vh*R -> Vh~
 
 These functions are defined for $\mathbb{P}_0(\mathbb{R})$ data and the
@@ -3522,6 +3543,33 @@ The following functions
 - ~dot: Rˆ3*Vh -> Vh~
 - ~dot: Vh*Rˆ3 -> Vh~
 
+****** ~R^1x1*Vh -> Vh~ and ~Vh*R^1x1 -> Vh~
+
+These functions are defined for $\mathbb{P}_0(\mathbb{R}^1x1)$ data and the
+return value is a $\mathbb{P}_0(\mathbb{R})$ function.
+
+The following functions
+- ~doubleDot: Rˆ1x1*Vh -> Vh~
+- ~doubleDot: Vh*Rˆ1x1 -> Vh~
+
+****** ~R^2x2*Vh -> Vh~ and ~Vh*R^2x2 -> Vh~
+
+These functions are defined for $\mathbb{P}_0(\mathbb{R}^2x2)$ data and the
+return value is a $\mathbb{P}_0(\mathbb{R})$ function.
+
+The following functions
+- ~doubleDot: Rˆ2x2*Vh -> Vh~
+- ~doubleDot: Vh*Rˆ2x2 -> Vh~
+
+****** ~R^3x3*Vh -> Vh~ and ~Vh*R^3x3 -> Vh~
+
+These functions are defined for $\mathbb{P}_0(\mathbb{R}^3x3)$ data and the
+return value is a $\mathbb{P}_0(\mathbb{R})$ function.
+
+The following functions
+- ~doubleDot: Rˆ3x3*Vh -> Vh~
+- ~doubleDot: Vh*Rˆ3x3 -> Vh~
+
 ******  ~Vh -> R~
 
 These functions are defined for $\mathbb{P}_0(\mathbb{R})$ data and the