diff --git a/doc/userdoc.org b/doc/userdoc.org
index a1ec61d21629c2c6ec34a17ce338eb03e3b43943..97c6c7ee7b47c3ac08e2544c0ddaa06e63cbec01 100644
--- a/doc/userdoc.org
+++ b/doc/userdoc.org
@@ -3883,6 +3883,159 @@ function.
 
 **** Operators overloading for ~Vh~ <<Vh-operators>>
 
+The ~scheme~ module provides overload of binary operators. Since ~Vh~ is
+an abstract type, some operators may not be defined (or allowed) for
+concrete types.
+
+***** Unary operators
+
+The only supported unary operators for ~Vh~ are given in the following
+table
+| operator | description          |
+|----------+----------------------|
+| ~+~        | plus unary operator  |
+| ~-~        | minus unary operator |
+We recall that the unary ~+~ operator is a convenience operator that has
+no effect.
+
+***** Binary operators
+
+The supported binary operators for ~vh~ data types are arithmetic
+operators.
+
+#+begin_src latex :results drawer :exports results
+  \begin{equation*}
+    \left|
+      \begin{array}{rl}
+        \mathtt{+}:& \mathbb{V}_h \times \mathbb{V}_h \to \mathbb{V}_h\\
+        \mathtt{-}:& \mathbb{V}_h \times \mathbb{V}_h \to \mathbb{V}_h\\
+        \mathtt{*}:& \mathbb{V}_h \times \mathbb{V}_h \to \mathbb{V}_h\\
+        \mathtt{/}:& \mathbb{V}_h \times \mathbb{V}_h \to \mathbb{V}_h
+      \end{array}
+    \right.
+  \end{equation*}
+#+end_src
+
+Observe that in the case of $\vec{\mathbb{P}}_0(\mathbb{R})$, the only
+available operators are ~+~ and ~-~.
+
+#+BEGIN_note
+In the case of $\mathbb{P}_0$ functions, an operator is defined as soon
+as it is defined for the value type.
+#+END_note
+For instance, one can multiply a $\mathbb{P}_0(\mathbb{R}^{2\times2})$
+discrete function by a $\mathbb{P}_0(\mathbb{R}^2)$. The result is then
+a $\mathbb{P}_0(\mathbb{R}^2)$ discrete function.
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+  import scheme;
+
+  let m:mesh, m = cartesianMesh(0,[1,1],(10,10));
+
+  let A:R^2 -> R^2x2, x -> [[2*x[0], x[1]],[-x[0], 3*x[1]]];
+  let u:R^2 -> R^2,   x -> 2*[x[0], x[1]*x[0]];
+
+  let Ah:Vh, Ah = interpolate(m, P0(), A);
+  let uh:Vh, uh = interpolate(m, P0(), u);
+
+  let Auh:Vh, Auh = Ah*uh;
+#+END_SRC
+
+Another illustration is: trying to add $\mathbb{P}_0(\mathbb{R})$ and
+$\mathbb{P}_0(\mathbb{R}^1)$
+#+NAME: invalid-Vh-sum-type
+#+BEGIN_SRC pugs-error :exports both :results output
+import mesh;
+import scheme;
+
+let m:mesh, m = cartesianMesh([0],[1],10);
+
+let f:R^1  -> R,   x -> 2*x[0];
+let f1:R^1 -> R^1, x -> 2*x;
+
+let fh:Vh,  fh  = interpolate(m, P0(), f);
+let f1h:Vh, f1h = interpolate(m, P0(), f1);
+
+fh+f1h;
+#+END_SRC
+produces the following error
+#+results: invalid-Vh-sum-type
+
+****** Additional ~+~ and ~-~ operators
+
+#+begin_src latex :results drawer :exports results
+  \begin{equation*}
+    \forall \mathbb{S}, \mathbb{S}_2 \in \{\mathbb{B},\mathbb{N},\mathbb{Z},\mathbb{R},\mathbb{R}^1,\mathbb{R}^2,\mathbb{R}^3,\mathbb{R}^{1\times1},\mathbb{R}^{2\times2},\mathbb{R}^{3\times3}\},
+    \quad
+    \left|
+      \begin{array}{rl}
+        \mathtt{+}:& \mathbb{S} \times \mathbb{V}_h \to \mathbb{V}_h\\
+        \mathtt{-}:& \mathbb{S} \times \mathbb{V}_h \to \mathbb{V}_h\\
+        \mathtt{+}:& \mathbb{V}_h \times \mathbb{S} \to \mathbb{V}_h\\
+        \mathtt{-}:& \mathbb{V}_h \times \mathbb{S} \to \mathbb{V}_h
+      \end{array}
+    \right.
+  \end{equation*}
+#+end_src
+
+Let us consider the following example
+#+NAME: substract-mean-value-to-Vh
+#+BEGIN_SRC pugs :exports both :results output
+  import mesh;
+  import scheme;
+  import math;
+
+  let m:mesh, m = cartesianMesh(0,[0.3,1.1],(5,15));
+
+  let mesh_volume:R, mesh_volume = sum_of_R(cell_volume(m));
+
+  let u:R^2 -> R, x -> 2*dot(x,x);
+
+  let uh:Vh, uh = interpolate(m, P0(), u);
+  let uh0:Vh, uh0 = uh - (integral_of_R(uh) / mesh_volume);
+
+  cout << "integral(uh)  = " << integral_of_R(uh) << "\n";
+  cout << "integral(uh0) = " << integral_of_R(uh0) << "\n";
+#+END_SRC
+Here we substract the mean value of a discrete function.
+#+results: substract-mean-value-to-Vh
+
+****** Additional ~*~ operators
+
+The following constructions are allowed for ~*~ operator.
+#+begin_src latex :results drawer :exports results
+  \begin{equation*}
+    \forall \mathbb{S} \in \{\mathbb{B},\mathbb{N},\mathbb{Z},\mathbb{R},\mathbb{R}^{1\times1},\mathbb{R}^{2\times2},\mathbb{R}^{3\times3}\},
+    \quad
+    \mathtt{*}: \mathbb{S} \times \mathbb{V}_h \to \mathbb{V}_h.
+  \end{equation*}
+#+end_src
+Obviously, if $\mathbb{S}=\mathbb{R}^{d\times d}$, for $d\in\{1,2,3\}$,
+the right operand must be have a compatible type, for instance
+$\mathbb{P}_0(\mathbb{R}^{d\times d})$ or $\mathbb{P}_0(\mathbb{R}^d)$.
+
+Additionally these operators are defined
+#+begin_src latex :results drawer :exports results
+  \begin{equation*}
+    \forall \mathbb{S} \in \{\mathbb{B},\mathbb{N},\mathbb{Z},\mathbb{R},\mathbb{R}^1,\mathbb{R}^2,\mathbb{R}^3,\mathbb{R}^{1\times1},\mathbb{R}^{2\times2},\mathbb{R}^{3\times3}\},
+    \quad
+    \mathtt{*}: \mathbb{V}_h \times  \mathbb{S}\to \mathbb{V}_h.
+  \end{equation*}
+#+end_src
+Following logic, if for instance the right operand is an ~R^2~
+expression, the left operand ~Vh~ must be of type
+$\mathbb{P}_0(\mathbb{R})$ or $\mathbb{P}_0(\mathbb{R}^{2\times 2})$.
+
+****** Additional ~/~ operators
+
+The following operators are defined
+#+begin_src latex :results drawer :exports results
+  \begin{equation*}
+    \forall \mathbb{S} \in \{\mathbb{N},\mathbb{Z},\mathbb{R}\},
+    \quad
+    \mathtt{/}: \mathbb{S} \times \mathbb{V}_h \to \mathbb{V}_h.
+  \end{equation*}
+#+end_src
 
 [fn:pugs-def] ~pugs~: Parallel Unstructured Grid Solvers
 [fn:MPI-def] ~MPI~: Message Passing Interface