Add description of scheme module: types and functions

doc/userdoc.org

@@ -2463,7 +2463,7 @@ The ~resetRandomSeed~ creates a *new shared* random and displays its value.
 The output is
 #+RESULTS: reset-random-seed
 
-***** ~setRandomSeed: Z -> void~
+***** ~setRandomSeed: Z -> void~ <<set-random-seed>>
 
 In order to reproduce exactly the same results of a calculation, it
 can be interesting to set the random seed to some given value.
@@ -2624,7 +2624,9 @@ The information produced concerns
 
 **** ~mesh~ provided functions
 
-***** ~boundaryName: string -> boundary~
+***** Boundary descriptor functions
+
+****** ~boundaryName: string -> boundary~
 
 Sets a boundary descriptor to a boundary name
 #+BEGIN_SRC pugs :exports both :results none
@@ -2633,7 +2635,7 @@ Sets a boundary descriptor to a boundary name
   let b:boundary, b = boundaryName("boundary_name");
 #+END_SRC
 
-***** ~boundaryTag: Z -> boundary~
+****** ~boundaryTag: Z -> boundary~
 
 Creates a boundary descriptor to a boundary reference
 #+BEGIN_SRC pugs :exports both :results none
@@ -2642,7 +2644,9 @@ Creates a boundary descriptor to a boundary reference
   let b:boundary, b = boundaryTag(12);
 #+END_SRC
 
-***** ~zoneName: string -> zone~
+***** Zone descriptor functions
+
+****** ~zoneName: string -> zone~
 
 Creates a zone descriptor from a ~string~ name
 #+BEGIN_SRC pugs :exports both :results none
@@ -2651,7 +2655,7 @@ Creates a zone descriptor from a ~string~ name
   let z:zone, z = zoneName("zone_name");
 #+END_SRC
 
-***** ~zoneTag: Z -> zone~
+****** ~zoneTag: Z -> zone~
 
 Associates a zone descriptor from zone tag
 #+BEGIN_SRC pugs :exports both :results none
@@ -2936,7 +2940,7 @@ The result of the previous script is given on Figure [[fig:transformed]].
 #+ATTR_HTML: :width 300px;
 #+RESULTS: transformed-img
 
-***** ~relax: mesh*mesh*R -> mesh~
+***** ~relax: mesh*mesh*R -> mesh~ <<relax-function>>
 
 This function is a simple utility that computes a ~mesh~ as the /mean/ of
 two other mesh that share the same connectivity.  The coordinates of
@@ -3006,7 +3010,878 @@ This function is mainly useful to reduce the displacement of nodes
 when using the ~randomizeMesh~ functions (see section [[scheme]]) for
 instance.
 
-**** The ~scheme~ module<<scheme>>
+*** The ~scheme~ module<<scheme>>
+
+This module provides a lot of numerical tools.
+#+NAME: get-module-info-scheme
+#+BEGIN_SRC pugs :exports both :results output
+  cout << getModuleInfo("scheme") << "\n";
+#+END_SRC
+#+RESULTS: get-module-info-scheme
+
+#+BEGIN_warning
+This module is very large and will be split in smaller ones in the future
+#+END_warning
+
+This module provides various types and functions. It also provides a
+set of operators overloading.
+
+**** ~scheme~ provided types
+
+***** ~Vh~
+
+The ~Vh~ type is the /abstract/ type of variables that refer to *discrete
+functions*.
+
+****** $\mathbb{P}_0$ scalar functions
+
+These functions are defined on a ~mesh~ and have a constant value in
+each cell.
+
+The type of values in each cells for a $\mathbb{P}_0$ function can be
+$\mathbb{R}$, $\mathbb{R}^1$, $\mathbb{R}^2$, $\mathbb{R}^3$,
+$\mathbb{R}^{1\times1}$, $\mathbb{R}^{2\times2}$ or $\mathbb{R}^{3\times3}$.
+
+For simplicity in this document, we denote these types specific
+$\mathbb{P}_0$ types as $\mathbb{P}_0(\mathbb{R})$,
+$\mathbb{P}_0(\mathbb{R}^1)$, $\mathbb{P}_0(\mathbb{R}^2)$,
+$\mathbb{P}_0(\mathbb{R}^3)$, $\mathbb{P}_0(\mathbb{R}^{1\times1})$,
+$\mathbb{P}_0(\mathbb{R}^{2\times2})$ or
+$\mathbb{P}_0(\mathbb{R}^{3\times3})$.
+
+****** $\mathbb{P}_0$ vector functions
+
+Additionally to scalar values per cell, one can define vectors of real
+values per cell. The size of the vectors is the same for all
+cells. This kind of variables is useful to define mass fractions for
+instance.
+
+Again for convenience, these types are denoted as
+$\vec{\mathbb{P}}_0(\mathbb{R})$.
+
+***** ~Vh_type~
+
+The ~Vh_type~ allows to describe a type of discretization. The available
+types of discretization are
+- ~P0~ for $\mathbb{P}_0(\mathbb{R})$, $\mathbb{P}_0(\mathbb{R}^1)$,
+  $\mathbb{P}_0(\mathbb{R}^2)$, $\mathbb{P}_0(\mathbb{R}^3)$,
+  $\mathbb{P}_0(\mathbb{R}^{1\times1})$, $\mathbb{P}_0(\mathbb{R}^{2\times2})$
+  or $\mathbb{P}_0(\mathbb{R}^{3\times3})$.
+- ~P0Vector~ for $\vec{\mathbb{P}}_0(\mathbb{R})$
+
+***** ~boundary_condition~
+
+This type is used to describe boundary conditions.
+
+***** ~quadrature~
+
+This type is used to describe quadrature types.
+
+**** ~scheme~ provided functions
+
+The ~scheme~ module provides a lot of functions, we categorize their
+description.
+
+***** Mathematical functions
+
+****** ~Vh -> Vh~
+
+These functions are defined for $\mathbb{P}_0(\mathbb{R})$ data and the
+return value is also a $\mathbb{P}_0(\mathbb{R})$ function. The value
+is simply the application of the function to the cell values.
+
+Here is the list of the functions
+- ~abs: Vh -> Vh~
+- ~acos: Vh -> Vh~
+- ~acosh: Vh -> Vh~
+- ~asin: Vh -> Vh~
+- ~asinh: Vh -> Vh~
+- ~atan: Vh -> Vh~
+- ~atanh: Vh -> Vh~
+- ~cos: Vh -> Vh~
+- ~cosh: Vh -> Vh~
+- ~exp: Vh -> Vh~
+- ~log: Vh -> Vh~
+- ~sin: Vh -> Vh~
+- ~sinh: Vh -> Vh~
+- ~sqrt: Vh -> Vh~
+- ~tan: Vh -> Vh~
+- ~tanh: Vh -> Vh~
+
+******  ~Vh*Vh -> Vh~
+
+These functions are defined for $\mathbb{P}_0(\mathbb{R})$ data and the
+return value is also a $\mathbb{P}_0(\mathbb{R})$ function. These
+functions require that the two arguments are defined one the *same
+mesh*. The result is obtained by applying the function cell by cell.
+
+Here is the function list
+- ~atan2: Vh*Vh -> Vh~
+- ~max: Vh*Vh -> Vh~
+- ~min: Vh*Vh -> Vh~
+- ~pow: Vh*Vh -> Vh~
+
+Let us mention another function that applies to
+$\mathbb{P}_0(\mathbb{R}^1)$, $\mathbb{P}_0(\mathbb{R}^2)$,
+$\mathbb{P}_0(\mathbb{R}^3)$ and to $\vec{\mathbb{P}}_0(\mathbb{R})$
+vector functions.
+- ~dot: Vh*Vh -> Vh~
+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.
+
+****** ~R*Vh -> Vh~ and ~Vh*R -> Vh~
+
+These functions are defined for $\mathbb{P}_0(\mathbb{R})$ data and the
+return value is also a $\mathbb{P}_0(\mathbb{R})$ function.
+
+The following functions can be applied using a scalar ~R~ and a ~Vh~
+operand.
+- ~atan2: Vh*R -> Vh~
+- ~atan2: R*Vh -> Vh~
+- ~max: Vh*R -> Vh~
+- ~max: R*Vh -> Vh~
+- ~min: Vh*R -> Vh~
+- ~min: R*Vh -> Vh~
+- ~pow: Vh*R -> Vh~
+- ~pow: R*Vh -> Vh~
+
+****** ~R^1*Vh -> Vh~ and ~Vh*R^1 -> Vh~
+
+These functions are defined for $\mathbb{P}_0(\mathbb{R}^1)$ data and the
+return value is also a $\mathbb{P}_0(\mathbb{R})$ function.
+
+The following functions
+- ~dot: Rˆ1*Vh -> Vh~
+- ~dot: Vh*Rˆ1 -> Vh~
+
+****** ~R^2*Vh -> Vh~ and ~Vh*R^2 -> Vh~
+
+These functions are defined for $\mathbb{P}_0(\mathbb{R}^2)$ data and the
+return value is also a $\mathbb{P}_0(\mathbb{R})$ function.
+
+The following functions
+- ~dot: Rˆ2*Vh -> Vh~
+- ~dot: Vh*Rˆ2 -> Vh~
+
+****** ~R^3*Vh -> Vh~ and ~Vh*R^3 -> Vh~
+
+These functions are defined for $\mathbb{P}_0(\mathbb{R}^3)$ data and the
+return value is also a $\mathbb{P}_0(\mathbb{R})$ function.
+
+The following functions
+- ~dot: Rˆ3*Vh -> Vh~
+- ~dot: Vh*Rˆ3 -> Vh~
+
+******  ~Vh -> R~
+
+These functions are defined for $\mathbb{P}_0(\mathbb{R})$ data and the
+return value is real ($\mathbb{R}$).
+
+The following functions
+- ~min: Vh -> R~\\
+  returns the minimal value of all the cell values
+- ~max: Vh -> R~\\
+  returns the maximal value of all the cell values
+- ~integral_of_R: Vh -> R~\\
+  computes the integral of a $\mathbb{P}_0(\mathbb{R})$ discrete
+  function $f$. If $f_j \in\mathbb{R}$ denotes the value of $f$ in the
+  cell $j$, the return value is $$\sum_{j\in\mathcal{J}}V_j f_j,$$ where $V_j$
+  denotes the volume of the cell $j$.
+- ~sum_of_R: Vh -> R~\\
+  computes the sum of the cell values of a $\mathbb{P}_0(\mathbb{R})$
+  discrete function $f$. If $f_j \in\mathbb{R}$ denotes the value of
+  $f$ in the cell $j$, the return value is $$\sum_{j\in\mathcal{J}} f_j.$$
+
+******  ~Vh -> R^1~
+
+These functions are defined for $\mathbb{P}_0(\mathbb{R}^1)$ data and
+the return value is a $\mathbb{R}^1$ vector.
+
+- ~integral_of_R1: Vh -> R^1~\\
+  computes the integral of a $\mathbb{P}_0(\mathbb{R}^1)$ discrete
+  function $\mathbf{u}$. If $\mathbf{u}_j \in\mathbb{R}^1$ denotes the
+  value of $\mathbf{u}$ in the cell $j$, the return value is
+  $$\sum_{j\in\mathcal{J}}V_j \mathbf{u}_j,$$ where $V_j$ denotes the volume of
+  the cell $j$.
+- ~sum_of_R1: Vh -> R^1~\\
+  computes the sum of the cell values of a $\mathbb{P}_0(\mathbb{R}^1)$
+  discrete function $\mathbf{u}$. If $\mathbf{u}_j \in\mathbb{R}^1$
+  denotes the value of $\mathbf{u}$ in the cell $j$, the return value
+  is $$\sum_{j\in\mathcal{J}} \mathbf{u}_j.$$
+
+******  ~Vh -> R^2~
+
+These functions are defined for $\mathbb{P}_0(\mathbb{R}^2)$ data and
+the return value is a $\mathbb{R}^2$ vector.
+
+- ~integral_of_R2: Vh -> R^2~\\
+  computes the integral of a $\mathbb{P}_0(\mathbb{R}^2)$ discrete
+  function $\mathbf{u}$. If $\mathbf{u}_j \in\mathbb{R}^2$ denotes the
+  value of $\mathbf{u}$ in the cell $j$, the return value is
+  $$\sum_{j\in\mathcal{J}}V_j \mathbf{u}_j,$$ where $V_j$ denotes the volume of
+  the cell $j$.
+- ~sum_of_R2: Vh -> R^2~\\
+  computes the sum of the cell values of a $\mathbb{P}_0(\mathbb{R}^2)$
+  discrete function $\mathbf{u}$. If $\mathbf{u}_j \in\mathbb{R}^2$
+  denotes the value of $\mathbf{u}$ in the cell $j$, the return value
+  is $$\sum_{j\in\mathcal{J}} \mathbf{u}_j.$$
+
+******  ~Vh -> R^3~
+
+These functions are defined for $\mathbb{P}_0(\mathbb{R}^3)$ data and
+the return value is a $\mathbb{R}^3$ vector.
+
+- ~integral_of_R3: Vh -> R^3~\\
+  computes the integral of a $\mathbb{P}_0(\mathbb{R}^3)$ discrete
+  function $\mathbf{u}$. If $\mathbf{u}_j \in\mathbb{R}^3$ denotes the
+  value of $\mathbf{u}$ in the cell $j$, the return value is
+  $$\sum_{j\in\mathcal{J}}V_j \mathbf{u}_j,$$ where $V_j$ denotes the volume of
+  the cell $j$.
+- ~sum_of_R3: Vh -> R^3~\\
+  computes the sum of the cell values of a $\mathbb{P}_0(\mathbb{R}^3)$
+  discrete function $\mathbf{u}$. If $\mathbf{u}_j \in\mathbb{R}^3$
+  denotes the value of $\mathbf{u}$ in the cell $j$, the return value
+  is $$\sum_{j\in\mathcal{J}} \mathbf{u}_j.$$
+
+******  ~Vh -> R^1x1~
+
+These functions are defined for $\mathbb{P}_0(\mathbb{R}^{1\times1})$ data
+and the return value is an $\mathbb{R}^{1\times1}$ matrix.
+
+- ~integral_of_R1x1: Vh -> R^1x1~\\
+  computes the integral of a $\mathbb{P}_0(\mathbb{R}^{1\times1})$
+  discrete function $A$. If $A_j \in\mathbb{R}^{1\times1}$ denotes the
+  value of $A$ in the cell $j$, the return value is
+  $$\sum_{j\in\mathcal{J}}V_j A_j,$$ where $V_j$ denotes the volume of the cell
+  $j$.
+- ~sum_of_R1x1: Vh -> R^1x1~\\
+  computes the sum of the cell values of a
+  $\mathbb{P}_0(\mathbb{R}^{1\times1})$ discrete function $A$. If $A_j
+  \in\mathbb{R}^{1\times1}$ denotes the value of $A$ in the cell $j$, the
+  return value is $$\sum_{j\in\mathcal{J}} A_j.$$
+
+******  ~Vh -> R^2x2~
+
+These functions are defined for $\mathbb{P}_0(\mathbb{R}^{2\times2})$ data
+and the return value is an $\mathbb{R}^{2\times2}$ matrix.
+
+- ~integral_of_R2x2: Vh -> R^2x2~\\
+  computes the integral of a $\mathbb{P}_0(\mathbb{R}^{2\times2})$
+  discrete function $A$. If $A_j \in\mathbb{R}^{2\times2}$ denotes the
+  value of $A$ in the cell $j$, the return value is
+  $$\sum_{j\in\mathcal{J}}V_j A_j,$$ where $V_j$ denotes the volume of the cell
+  $j$.
+- ~sum_of_R2x2: Vh -> R^2x2~\\
+  computes the sum of the cell values of a
+  $\mathbb{P}_0(\mathbb{R}^{2\times2})$ discrete function $A$. If $A_j
+  \in\mathbb{R}^{2\times2}$ denotes the value of $A$ in the cell $j$, the
+  return value is $$\sum_{j\in\mathcal{J}} A_j.$$
+
+******  ~Vh -> R^3x3~
+
+These functions are defined for $\mathbb{P}_0(\mathbb{R}^{3\times3})$ data
+and the return value is an $\mathbb{R}^{3\times3}$ matrix.
+
+- ~integral_of_R3x3: Vh -> R^3x3~\\
+  computes the integral of a $\mathbb{P}_0(\mathbb{R}^{3\times3})$
+  discrete function $A$. If $A_j \in\mathbb{R}^{3\times3}$ denotes the
+  value of $A$ in the cell $j$, the return value is
+  $$\sum_{j\in\mathcal{J}}V_j A_j,$$ where $V_j$ denotes the volume of the cell
+  $j$.
+- ~sum_of_R3x3: Vh -> R^3x3~\\
+  computes the sum of the cell values of a
+  $\mathbb{P}_0(\mathbb{R}^{3\times3})$ discrete function $A$. If $A_j
+  \in\mathbb{R}^{3\times3}$ denotes the value of $A$ in the cell $j$, the
+  return value is $$\sum_{j\in\mathcal{J}} A_j.$$
+
+
+***** Interpolation and integration functions
+
+These functions are ways to define discrete functions from analytic
+data.
+
+****** ~interpolate: mesh*Vh_type*(function) -> Vh~
+
+This functions takes a ~mesh~, a type of discrete function (~Vh_type~) and
+a list of user functions as arguments. It returns a $\mathbb{P}_0$
+function defined at the mesh.
+
+All the user functions of the list must be defined on $\mathbb{R}^d$
+for a mesh of dimension $d$.
+
+There are several situations according to the ~Vh_type~.
+
+******* ~P0~
+
+In that case the list of user functions *must* reduce to a *single* user
+function.
+
+The codomain (or image space) of the user function defines the type of
+the returned discrete function.
+
+Let us give an example.
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+  import scheme;
+
+  let m:mesh, m = cartesianMesh([0,0], [1,1], (10,10));
+  let u:R^2 -> R^3, x -> [x[0], 2*x[1], x[1]-x[0]];
+
+  let uh:Vh, uh = interpolate(m, P0(), u);
+#+END_SRC
+In this exampel the discrete function ~uh~ is a
+$\mathbb{P}_0(\mathbb{R}^3)$ function defined on a 2d ~mesh~. The ~P0()~
+function returns the type of interpolation ($\mathbb{P}_0$).
+
+#+BEGIN_note
+In the case of $\mathbb{P}_0$ interpolation, the function is evaluated
+at the mass center of the mesh cells $$\forall j\in\mathcal{J},\quad
+\mathbf{x}_j=\frac{1}{V_j}\int_j \mathbf{x}.$$
+#+END_note
+
+******* ~P0Vector~
+
+In that case the codomain of each user function in the list must be a
+real function (values in $\mathbb{R}$). The instruction will define a
+$\vec{\mathbb{P}}_0(\mathbb{R})$.
+
+A first example is
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+  import scheme;
+
+  let m:mesh, m = cartesianMesh([0,0], [1,1], (10,10));
+  let f:R^2 -> R, x -> 2*x[0]-x[1];
+
+  let fh:Vh, fh = interpolate(m, P0Vector(), f);
+#+END_SRC
+The obtained ~fh~ is a vector $\vec{\mathbb{P}}_0(\mathbb{R})$ where each
+vector (in each cell) is of dimension 1.
+
+The number of scalar user functions sets the size of the discrete
+vector function.
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+  import scheme;
+  import math;
+
+  let m:mesh, m = cartesianMesh([0,0], [1,1], (10,10));
+  let f0:R^2 -> R, x -> 2*x[0]-x[1];
+  let f1:R^2 -> R, x -> x[0]*x[1]-2;
+  let f2:R^2 -> R, x -> 2*dot(x,x);
+
+  let fh:Vh, fh = interpolate(m, P0Vector(), (f0,f1,f2));
+#+END_SRC
+Here we defined a $\vec{\mathbb{P}}_0(\mathbb{R})$ discrete function of
+dimension 3.
+
+****** ~interpolate: mesh*(zone)*Vh_type*(function) -> Vh~
+
+This function works exactly the same as the previous function. The
+additional parameter, the ~zone~ lists is used to define the cells where
+the user function (or the user function list) is interpolate. For
+cells that are not in the ~zone~ list, the discrete function is set to
+the value $0$.
+
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+  import scheme;
+
+  let m:mesh, m = readGmsh("hybrid-2d.msh");
+  let f:R^2 -> R, x -> 2*x[0]-x[1];
+
+  let fh:Vh, fh = interpolate(m, (zoneName("LEFT"), zoneName("RIGHT")), P0(), f);
+  let fh_l:Vh, fh_l = interpolate(m, zoneName("LEFT"), P0(), f);
+#+END_SRC
+In this example, we define two discrete functions. ~fh~ is defined as
+the interpolation of the function ~f~ in *all* cells of the mesh since the
+mesh is partitioned into two zones labeled ~LEFT~ and ~RIGHT~. The
+discrete function ~fh_l~ has exactly the same values in the ~LEFT~ region,
+but is $0$ in the cells that belong to ~RIGHT~.
+
+For completeness, we give an example in the case of ~P0Vector~.
+The number of scalar user functions sets the size of the discrete
+vector function.
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+  import scheme;
+  import math;
+
+  let m:mesh, m = readGmsh("hybrid-2d.msh");
+  let f0:R^2 -> R, x -> 2*x[0]-x[1];
+  let f1:R^2 -> R, x -> x[0]*x[1]-2;
+  let f2:R^2 -> R, x -> 2*dot(x,x);
+
+  let fh:Vh, fh = interpolate(m, zoneName("RIGHT"), P0Vector(), (f0,f1,f2));
+#+END_SRC
+Here we defined a $\vec{\mathbb{P}}_0(\mathbb{R})$ discrete function of
+dimension 3.
+
+****** ~integrate: mesh*quadrature*function -> Vh~ <<integrate-classic>>
+
+This function integrates the given user function in each cell of a
+~mesh~ using a prescribed ~quadrature~. The result is of type as
+$\mathbb{P}_0(\mathbb{R})$, $\mathbb{P}_0(\mathbb{R}^1)$,
+$\mathbb{P}_0(\mathbb{R}^2)$, $\mathbb{P}_0(\mathbb{R}^3)$,
+$\mathbb{P}_0(\mathbb{R}^{1\times1})$, $\mathbb{P}_0(\mathbb{R}^{2\times2})$
+or $\mathbb{P}_0(\mathbb{R}^{3\times3})$, according to the user function
+codomain.
+
+Let us consider the following example
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+  import scheme;
+  import math;
+
+  let m:mesh, m = readGmsh("hybrid-2d.msh");
+  let u:R^2 -> R^2, x -> [cos(x[0]), sin(x[1])];
+  let U:Vh, U = integrate(m, Gauss(5), u);
+#+END_SRC
+Here, for each cell $j$, the value of the discrete function
+$\mathbf{F}_j$ is computed using a Gauss quadrature formula that is
+exact for polynomials of degree $5$, $\mathbf{F}_j \approx\int_j
+\mathbf{u}$. More details about quadrature formula will be given
+below.
+
+Thus if one wants to project $\mathbf{u}$ on
+$\mathbb{P}_0(\mathbb{R}^2)$ to the sixth order, one can modify the
+previous script by writing
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+  import scheme;
+  import math;
+
+  let m:mesh, m = readGmsh("hybrid-2d.msh");
+  let u:R^2 -> R^2, x -> [cos(x[0]), sin(x[1])];
+  let uh:Vh, uh = (1/cell_volume(m)) * integrate(m, Gauss(5), u);
+#+END_SRC
+The function ~cell_volume: mesh -> Vh~ creates a
+$\mathbb{P}_0(\mathbb{R})$ function whose values are the volume of the
+cells.
+
+****** ~integrate: mesh*quadrature*Vh_type*(function) -> Vh~ <<integrate-P1-vector>>
+
+This function behaves the same, the user function list size defines
+the dimension of the vector value of the produced
+$\vec{\mathbb{P}}_0(\mathbb{R})$ discrete function. Actually the
+~Vh_type~ parameter is there to allow the construction of
+$\vec{\mathbb{P}}_0(\mathbb{R})$ of dimension 1 (since passing a single
+function as a list of user function, the previous function
+[[integrate-classic]], would be used).
+
+****** ~integrate: mesh*(zone)*quadrature*function -> Vh~
+
+This function is an enhancement of the function defined in
+[[integrate-classic]]. It allow to specify a ~zone~ list which defines the
+set of cells where the integration is operated.
+
+#+BEGIN_SRC shell :exports results :results none
+cat << EOF > zones-1d.geo
+//+
+Point(1) = {-1, 0, 0, 0.01};
+//+
+Point(2) = {-0.3, 0, 0, 0.01};
+//+
+Point(3) = {0.3, 0, 0, 0.01};
+//+
+Point(4) = {1, 0, 0, 0.01};
+//+
+Line(1) = {1, 2};
+//+
+Line(2) = {2, 3};
+//+
+Line(3) = {3, 4};
+//+
+Physical Point("XMIN") = {1};
+//+
+Physical Point("XMAX") = {4};
+//+
+Physical Point("INTERFACE1") = {2};
+//+
+Physical Point("INTERFACE2") = {3};
+//+
+Physical Curve("LEFT") = {1};
+//+
+Physical Curve("MIDDLE") = {2};
+//+
+Physical Curve("RIGHT") = {3};
+#+END_SRC
+
+#+BEGIN_SRC shell :exports results :results none
+gmsh -1 zones-1d.geo -format msh2
+#+END_SRC
+
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+  import scheme;
+  import math;
+  import writer;
+
+  let pi:R, pi = acos(-1);
+  let m:mesh, m = readGmsh("zones-1d.msh");
+  let f:R^1 -> R, x -> sin(2*pi*x[0]);
+
+  let fh:Vh,
+      fh = (1/cell_volume(m))
+         * integrate(m, (zoneName("LEFT"), zoneName("RIGHT")), Gauss(5), f);
+
+  write(gnuplot_1d_writer("zone_integrate",0), name_output(fh, "f"), 0);
+#+END_SRC
+In this example, the ~mesh~ provided in the file ~zones-1d.msh~ is a 1d
+~mesh~ of $]-1,1[$ made of $200$ cells that is partitioned into 3
+connected subdomains. The zones corresponding to these 3 subdomains
+are named ~LEFT~ for $]-1,-0.3[$, ~MIDDLE~ for $]-0.3, 0.3[$ and ~RIGHT~ for
+$]0.3,1[$. The result is displayed on Figure [[fig:zone-integrate-1d]]. In
+the ~MIDDLE~ region, cell values are set to 0.
+
+#+NAME: zone-integrate-1d-img
+#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/zone-integrate-1d.png")
+  reset
+  set grid
+  set border
+  unset key
+  set xtics
+  set ytics
+  set terminal png truecolor enhanced size 960,480
+  plot '<(sed "" $PUGS_SOURCE_DIR/doc/zone_integrate.0000.gnu)' lw 2 w l
+#+END_SRC
+
+#+CAPTION: $L^1$ projection of order 6 of the function $\sin(2\pi x)$ on the zones ~LEFT~ and ~RIGHT~. The values in the zone ~MIDDLE~ are set to $0$.
+#+NAME: fig:zone-integrate-1d
+#+ATTR_LATEX: :width 0.38\textwidth
+#+ATTR_HTML: :width 300px;
+#+RESULTS: zone-integrate-1d-img
+
+****** ~integrate: mesh*(zone)*quadrature*Vh_type*(function) -> Vh~
+
+This function behaves essentially as the function described in
+paragraph [[integrate-P1-vector]], it also adds the possibility to choose
+sets of cells where to integrate the list of user functions.
+
+***** Random mesh generators
+
+For numerical it is often useful to create meshes with random vertices
+positions. This is the aim of the functions that are described in this
+section. These function share some properties.
+- The generate mesh is always suitable for calculations in the sense
+  that cells volumes are warrantied to be positive.
+- Generated cells can be non-convex.
+- One has to specify boundary conditions to drive the mesh
+  displacement on boundaries.
+- The obtained mesh does not depend on parallelism: it is exactly the
+  same whichever is the number of ~MPI~ processes. It only depends on
+  the random seed (see paragraph [[set-random-seed]] how to set the random
+  seed to obtain the exact same mesh through different runs).
+
+****** ~randomizeMesh: mesh*(boundary_condition) -> mesh~
+
+This function creates a random mesh by displacing the nodes of a given
+~mesh~ and a list of ~boundary_condition~.
+
+The supported boundary conditions are the following:
+- ~fixed~: the vertices of the ~boundary~ cannot be displaced
+- ~axis~: the vertices are allowed to be displaced in the direction of
+  the ~boundary~
+- ~symmetry~: the vertices are displaced in the plane formed by the
+  ~boundary~
+One should refer to the section [[boundary-condition-descriptor]] for a
+description of the boundary condition descriptors.
+
+#+BEGIN_note
+Let us precise these boundary conditions behavior
+- In dimension 1, ~fixed~, ~axis~ and ~symmetry~ boundary conditions have
+  the same effect.
+- In dimension 2, ~axis~ and ~symmetry~ behave the same. Thus, boundaries
+  supporting this kind of boundary conditions *must* form *straight*
+  lines.
+- In dimension 3, boundaries describing ~axis~ conditions *must* be
+  *straight* lines, and boundaries describing ~symmetry~ conditions *must*
+  be *planar*.
+
+If a boundary does not satisfy geometrical requirements, ~pugs~ produces
+a runtime error.
+#+END_note
+
+Let us consider a simple example
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+  import scheme;
+  import writer;
+
+  setRandomSeed(123456789); // not required
+
+  let m:mesh, m = cartesianMesh([-1,-1],[1,1], (20,20));
+  let bc_list:(boundary_condition),
+      bc_list = ((fixed(boundaryName("XMINYMIN"))),
+                 (fixed(boundaryName("XMINYMAX"))),
+                 (fixed(boundaryName("XMAXYMIN"))),
+                 (fixed(boundaryName("XMAXYMAX"))),
+                 (symmetry(boundaryName("XMIN"))),
+                 (symmetry(boundaryName("XMAX"))),
+                 (symmetry(boundaryName("YMIN"))),
+                 (symmetry(boundaryName("YMAX"))));
+  m = randomizeMesh(m, bc_list);
+
+  write_mesh(gnuplot_writer("random-mesh", 0), m);
+#+END_SRC
+
+Running this script one gets the ~mesh~ displayed on Figure
+[[fig:random-mesh]]. To reduce the vertices displacement, one can use the
+~relax~ function, see section [[relax-function]].
+
+#+NAME: random-mesh-img
+#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/random-mesh.png")
+  reset
+  unset grid
+  unset border
+  unset key
+  unset xtics
+  unset ytics
+  set square
+  set terminal png truecolor enhanced size 640,640
+  plot '<(sed "" $PUGS_SOURCE_DIR/doc/random-mesh.0000.gnu)' w l
+#+END_SRC
+
+#+CAPTION: The obtained random mesh
+#+NAME: fig:random-mesh
+#+ATTR_LATEX: :width 0.38\textwidth
+#+ATTR_HTML: :width 300px;
+#+RESULTS: random-mesh-img
+
+****** ~randomizeMesh: mesh*(boundary_condition)*function -> mesh~
+
+This function is a variation of the previous one. It allows
+additionally to provide a characteristic function that designates the
+vertices that can be displaced.
+
+The characteristic function *must* be a function of $\mathbb{R}^d
+\to\mathbb{B}$ for a ~mesh~ of dimension $d$.
+
+Here is a modification of the previous example, where the random
+displacement is allowed for $x<2y$.
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+  import scheme;
+  import writer;
+  import math;
+
+  setRandomSeed(123456789); // not required
+
+  let m:mesh, m = cartesianMesh([-1,-1],[1,1], (20,20));
+  let bc_list:(boundary_condition),
+      bc_list = ((fixed(boundaryName("XMINYMIN"))),
+                 (fixed(boundaryName("XMINYMAX"))),
+                 (fixed(boundaryName("XMAXYMIN"))),
+                 (fixed(boundaryName("XMAXYMAX"))),
+                 (symmetry(boundaryName("XMIN"))),
+                 (symmetry(boundaryName("XMAX"))),
+                 (symmetry(boundaryName("YMIN"))),
+                 (symmetry(boundaryName("YMAX"))));
+
+  let chi:R^2 -> B, x -> x[0]<2*x[1];
+
+  m = randomizeMesh(m, bc_list, chi);
+
+  write_mesh(gnuplot_writer("random-mesh-chi", 0), m);
+#+END_SRC
+
+Running this script one gets the ~mesh~ displayed on Figure
+[[fig:random-mesh-chi]].
+
+#+NAME: random-mesh-chi-img
+#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/random-mesh-chi.png")
+  reset
+  unset grid
+  unset border
+  unset key
+  unset xtics
+  unset ytics
+  set square
+  set terminal png truecolor enhanced size 640,640
+  plot '<(sed "" $PUGS_SOURCE_DIR/doc/random-mesh-chi.0000.gnu)' w l
+#+END_SRC
+
+#+CAPTION: Random mesh with characteristic function. The original mesh is unchanged for $x \ge 2 y$.
+#+NAME: fig:random-mesh-chi
+#+ATTR_LATEX: :width 0.38\textwidth
+#+ATTR_HTML: :width 300px;
+#+RESULTS: random-mesh-chi-img
+
+#+BEGIN_note
+Since we set the random seed to the same value in both cases (with or
+without using the characteristic function $\chi$), node displacements are
+the same $\forall r\in\mathcal{R}$ such that $\chi(\mathbf{x}_r)$, see Figures
+[[fig:random-mesh]] and [[fig:random-mesh-chi]].
+
+This allows for instance to build a random mesh step-by-step.
+#+END_note
+
+***** Boundary condition descriptors <<boundary-condition-descriptor>>
+
+In this paragraph, we describe the set of boundary condition
+descriptors that are provided by the ~scheme~ module.
+
+#+BEGIN_note
+These functions provide *descriptors*, these are not related to a
+particular implementation. The way they are used in different
+functions dependents of the context or the numerical method itself.
+#+END_note
+
+#+BEGIN_note
+There a three kind of boundaries are supported by ~pugs~, boundaries
+made
+- of sets of nodes,
+- of sets of edges, or
+- of sets of faces.
+
+In dimension 1, nodes, edges and faces denotes the same entities. In
+dimension 2, edges and faces refer the same entities.
+#+END_note
+
+#+BEGIN_note
+~pugs~ integrates some automatic mechanisms to translate boundary types
+in other ones.
+
+For instance, if an algorithm or a method requires a set of nodes to
+set some numerical treatment, it can be deduced from a set of faces.
+
+Obviously, these conversions can be meaningless, for instance, if one
+expects a *line* in 3d, cannot be defined by a set of faces. ~pugs~ will
+forbid this kind of conversion at runtime.
+#+END_note
+
+#+BEGIN_note
+When a method expects a *straight* line or a *planar* surface, ~pugs~ checks
+that the given boundary is actually *straight* or *planar*.
+#+END_note
+
+#+BEGIN_note
+~pugs~ checks that boundaries do not contain /inner/ items.
+#+END_note
+
+We regroup the boundary condition descriptors functions according to
+their arguments
+****** ~boundary -> boundary_condition~
+- ~axis: boundary -> boundary_condition~\\
+  The boundary denotes a *straight* line of the mesh
+- ~fixed: boundary -> boundary_condition~
+- ~symmetry: boundary -> boundary_condition~\\
+  The boundary denotes a *planar* surface of the mesh
+
+****** ~boundary*function -> boundary_condition~
+The provided user function type depends on the numerical method of
+function that utilizes the boundary condition.
+- ~pressure: boundary*function -> boundary_condition~
+- ~velocity: boundary*function -> boundary_condition~
+
+***** Quadrature formulas descriptors
+
+~pugs~ provides a quite complete set of quadrature formulas that can be
+chosen inside the script.
+
+#+BEGIN_warning
+While ~pugs~ is written to deal with general polygonal and polyhedral
+meshes, quadrature formulas are not yet implemented on the general
+elements. Supported elements are triangles, quadrangles, tetrahedra
+pyramids, prisms and hexahedra.
+#+END_warning
+
+****** ~Gauss: N -> quadrature~
+
+This function returns the quadrature descriptor associated to Gauss
+formulas for the given ~N~.
+
+In the following table, we summarize the *maximal degree* quadrature
+that are available in ~pugs~ for various elements.
+| element type           | max. degree |
+|------------------------+-------------|
+| segment                |          23 |
+|------------------------+-------------|
+| triangle               |          20 |
+| square                 |          21 |
+|------------------------+-------------|
+| tetrahedron            |          20 |
+| pyramid (square basis) |          20 |
+| prism (triangle basis) |          20 |
+| cube                   |          21 |
+
+****** ~GaussLegendre: N -> quadrature~
+
+Gets the Gauss-Legendre quadrature descriptor exact for polynomials of
+degree given in argument.
+
+The maximum allowed degree is 23.
+
+For dimension 2 or 3 elements, Gauss-Legendre formulas are defined by
+tensorization. Conform transformations are used to map the cube
+$]-1,1[^d$ to supported elements.
+
+****** ~GaussLobatto: N -> quadrature~
+
+Gets the Gauss-Lobatto quadrature descriptor exact for polynomials of
+degree given in argument.
+
+The maximum allowed degree is "only" 13.
+
+For dimension 2 or 3 elements, Gauss-Lobatto formulas are defined by
+tensorization. Conform transformations are used to map cube $]-1,1[^d$
+to supported elements.
+
+***** ~lagrangian: mesh*Vh -> Vh~
+
+This function is a special function whose purpose is to transport
+lagrangian quantities from one mesh to the other. Obviously, this
+function make lots of sense in the case of lagrangian calculations.
+
+This is a zero cost function, since discrete functions are *constants*
+in ~pugs~, it consists in associating the data of the given discrete
+function to the provided ~mesh~. In other words, the underlying array of
+values is shared by the two discrete functions, which are associated
+to different meshes.
+
+A good example of the use of this kind of function is mass fractions.
+
+***** Numerical methods
+
+We describe rapidly two functions that embed numerical methods. These
+are in some sense /models/ that are used to test evolution of ~pugs~
+itself and can be used as examples.
+
+#+BEGIN_warning
+These functions will become obsolete (soon?), since another interface
+to numerical methods is in preparation.
+#+END_warning
+
+The functions are
+- ~eucclhyd_solver: Vh*Vh*Vh*Vh*Vh*(boundary_condition)*R -> mesh*Vh*Vh*Vh~
+- ~glace_solver: Vh*Vh*Vh*Vh*Vh*(boundary_condition)*R -> mesh*Vh*Vh*Vh~
+Both functions share the same interface. The arguments are provided in
+this order:
+- the mass density $\rho$ of type $\mathbb{P}_0(\mathbb{R})$,
+- the velocity $\mathbf{u}$ of type $\mathbb{P}_0(\mathbb{R}^d)$ in
+  dimension $d$,
+- the total energy density $E$ of type $\mathbb{P}_0(\mathbb{R})$,
+- the sound speed $c$ of type $\mathbb{P}_0(\mathbb{R})$,
+- the pressure $p$ of type $\mathbb{P}_0(\mathbb{R})$,
+- a list of boundary conditions,
+- and a time step of type ~R~.
+Observe that ~pugs~ checks the types of the parameters and that all
+discrete functions are defined on the same mesh.
+
+The functions return a compound type made of
+- the new ~mesh~ $\mathcal{M}$,
+- the new mass density ~\rho~ of type $\mathbb{P}_0(\mathbb{R})$ defined
+  on $\mathcal{M}$,
+- the new velocity $\mathbf{u}$ of type $\mathbb{P}_0(\mathbb{R}^d)$ in
+  dimension $d$, defined on the mesh $\mathcal{M}$,
+- and the new total energy density $E$ of type
+  $\mathbb{P}_0(\mathbb{R})$, defined on the mesh $\mathcal{M}$.
+
+The time step can be calculated through the ~acoustic_dt: Vh -> R~
+function.
+
+**** Operators overloading for ~Vh~ <<Vh-operators>>
 
 
 [fn:pugs-def] ~pugs~: Parallel Unstructured Grid Solvers