Add the mesh module section

.gitignore

@@ -17,3 +17,10 @@ GTAGS
 /doc/*.png
 /doc/*.gnu
 /doc/filename.txt
+/doc/*pvd
+/doc/*.geo
+/doc/*.py
+/doc/*.pvtu
+/doc/lisp/.fltk/
+/doc/*.msh
+/doc/*.vtu

doc/userdoc.org

@@ -30,15 +30,15 @@
 #+LATEX_HEADER_EXTRA: \usepackage{mdframed}
 
 #+LATEX_HEADER_EXTRA: \BeforeBeginEnvironment{tabular}{\rowcolors[]{2}{orange!5}{orange!10}}
-#+LATEX_HEADER_EXTRA: \BeforeBeginEnvironment{minted}{\begin{mdframed}[linecolor=blue,backgroundcolor=blue!5]}
+#+LATEX_HEADER_EXTRA: \BeforeBeginEnvironment{minted}{\begin{mdframed}[linecolor=blue,backgroundcolor=blue!10]}
 #+LATEX_HEADER_EXTRA: \AfterEndEnvironment{minted}{\end{mdframed}}
-#+LATEX_HEADER_EXTRA: \BeforeBeginEnvironment{verbatim}{\begin{mdframed}[linecolor=gray,backgroundcolor=green!5]}
+#+LATEX_HEADER_EXTRA: \BeforeBeginEnvironment{verbatim}{\begin{mdframed}[linecolor=gray,backgroundcolor=gray!5]}
 #+LATEX_HEADER_EXTRA: \AfterEndEnvironment{verbatim}{\end{mdframed}}
 #+LATEX_HEADER_EXTRA: \newtheorem{note}{Note}
-#+LATEX_HEADER_EXTRA: \BeforeBeginEnvironment{note}{\begin{mdframed}[linecolor=orange,backgroundcolor=orange!5]}
+#+LATEX_HEADER_EXTRA: \BeforeBeginEnvironment{note}{\begin{mdframed}[linecolor=green,backgroundcolor=green!10]}
 #+LATEX_HEADER_EXTRA: \AfterEndEnvironment{note}{\end{mdframed}}
 #+LATEX_HEADER_EXTRA: \newtheorem{warning}{Warning}
-#+LATEX_HEADER_EXTRA: \BeforeBeginEnvironment{warning}{\begin{mdframed}[linecolor=red,backgroundcolor=red!5]}
+#+LATEX_HEADER_EXTRA: \BeforeBeginEnvironment{warning}{\begin{mdframed}[linecolor=red,backgroundcolor=red!10]}
 #+LATEX_HEADER_EXTRA: \AfterEndEnvironment{warning}{\end{mdframed}}
 
 #+LATEX_HEADER_EXTRA: \usepackage{mathpazo}
@@ -77,7 +77,7 @@ mesh nodes according to a user-defined vector field $T: \mathbb{R}^2
   let M:R^2 -> R^2x2, x -> [[cos(theta(x)), -sin(theta(x))],
                             [sin(theta(x)),  cos(theta(x))]];
   let T: R^2 -> R^2, x -> x + M(x)*x;
-  let m:mesh, m = cartesian2dMesh([-1,-1], [1,1], (20,20));
+  let m:mesh, m = cartesianMesh([-1,-1], [1,1], (20,20));
 
   m = transform(m, T);
 
@@ -150,7 +150,7 @@ already be discussed.
   section [[basic-types]] and [[high-level-types]] for details.
 - Also, there are two types of functions: *user-defined* functions and
   *builtin* functions. In this example, ~theta~, ~M~ and ~T~ are user-defined
-  functions. All other functions (~cos~, ~cartesian2dMesh~,...) are
+  functions. All other functions (~cos~, ~cartesianMesh~,...) are
   builtin functions and are generally defined when importing a
   module. These functions behave similarly, one should refer to
   [[functions]] for details.
@@ -1632,9 +1632,9 @@ illustrate this, let us consider the following example.
 #+BEGIN_SRC pugs :exports both
   import mesh;
 
-  let m1:mesh, m1 = cartesian2dMesh(0, [1,1], (10,10));
+  let m1:mesh, m1 = cartesianMesh(0, [1,1], (10,10));
   let m2:mesh, m2 = m1;
-  let m3:mesh, m3 = cartesian2dMesh(0, [1,1], (10,10));
+  let m3:mesh, m3 = cartesianMesh(0, [1,1], (10,10));
 
   m3 = m1;
 #+END_SRC
@@ -2490,7 +2490,7 @@ Running this example produces no output
 But a file is created (in the execution directory), with the name
 ~"filename.txt"~. Its content is
 #+NAME: cat-filename-txt
-#+BEGIN_SRC bash :exports results :results output
+#+BEGIN_SRC shell :exports results :results output
   cat filename.txt
 #+END_SRC
 #+RESULTS: cat-filename.txt
@@ -2558,7 +2558,453 @@ rounding or truncation functions are ~ceil~, ~floor~, ~round~ and ~trunc~. See
 their ~C++~ man pages for details.
 #+END_note
 
+#+BEGIN_note
+Let us comment the use of the ~pow~ function. Actually one can wonder
+why we did not use a syntax like ~x^y~? The reason is that if
+mathematically ${x^y}^z = x^{(y^z)}$, many software treat it (by mistake)
+as ${(x^y)}^z$. Thus, using the ~pow~ function avoids any confusion.
+#+END_note
+
+*** The ~mesh~ module
+
+This is an important module. It provide mesh utilities tools.
+#+NAME: get-module-info-mesh
+#+BEGIN_SRC pugs :exports both :results output
+  cout << getModuleInfo("mesh") << "\n";
+#+END_SRC
+#+RESULTS: get-module-info-mesh
+
+**** ~mesh~ provided types
+
+***** ~boundary~
+
+The ~boundary~ type is a boundary descriptor: it refers to a boundary of
+a ~mesh~ that is either designated by an integer or by a ~string~.
+
+A ~boundary~ can be used refer to an interface, it can designate a set
+of nodes, edges or faces. The ~boundary~ (descriptor) itself is not
+related to any ~mesh~, thus the nature of the ~boundary~ is precised when
+it is used with a particular ~mesh~.
+
+***** ~zone~
+
+Following the same idea, a ~zone~ is a descriptor of set of cells. It
+can be either defined by an integer or by a ~string~. Its meaning is
+precised when it is associated with a ~mesh~.
+
+***** ~mesh~
+
+The type ~mesh~ is an *abstract* type that is used to store meshes. A
+variable of that type can refer for instance unstructured meshes of
+dimension 1, 2 or 3.
+
+The following binary operator is provided.
+| ~ostream <<~ allowed expression type |
+|------------------------------------|
+| ~mesh~                               |
+
+It enables to write mesh information to an ~ostream~, here is a simple
+example.
+
+#+NAME: ostream-mesh-example
+#+BEGIN_SRC pugs :exports both :results output
+  import mesh;
+
+  let m:mesh, m = cartesianMesh(0,[1,1,2], (2,3,2));
+  cout << m << "\n";
+#+END_SRC
+Running this script generates the following output.
+#+RESULTS: ostream-mesh-example
+The information produced concerns
+- the dimension of the connectivity,
+- the number of cells and their references,
+- the number of faces and their references,
+- the number of edges and their references,
+- and the number of nodes and their references.
+
+**** ~mesh~ provided functions
+
+***** ~boundaryName: string -> boundary~
+
+Sets a boundary descriptor to a boundary name
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+
+  let b:boundary, b = boundaryName("boundary_name");
+#+END_SRC
+
+***** ~boundaryTag: Z -> boundary~
+
+Creates a boundary descriptor to a boundary reference
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+
+  let b:boundary, b = boundaryTag(12);
+#+END_SRC
+
+***** ~zoneName: string -> zone~
+
+Creates a zone descriptor from a ~string~ name
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+
+  let z:zone, z = zoneName("zone_name");
+#+END_SRC
+
+***** ~zoneTag: Z -> zone~
+
+Associates a zone descriptor from zone tag
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+
+  let z:zone, z = zoneTag(5);
+#+END_SRC
+
+***** ~cartesianMesh: Rˆd*Rˆd*(N) -> mesh~ (with $d\in\{1, 2, 3\}$)
+
+Creates a cartesian mesh of dimension $d\in\{1, 2, 3\}$. The produced
+cartesian grid is aligned with the axis and made of identical cells.
+
+The two first arguments are two opposite corners of the box (or
+segment in 1d) and the list of natural integers (type ~(N)~) sets the
+number of *cells* in each direction. Thus size of the list of ~N~ is $d$.
+
+For instance one can write:
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+
+  let m1d:mesh, m1d = cartesianMesh([0], [1], 10);
+  let m2d:mesh, m2d = cartesianMesh([0,0], [1,1], (3,6));
+  let m3d:mesh, m3d = cartesianMesh([-1,-1,-1], [1,1,1], (3,2,4));
+#+END_SRC
+- The ~m1d~ variable contains a uniform grid of $]0,1[$ made of 10 cells.
+- The ~m2d~ variable refers to a uniform grid of $]0,1[^2$, made of 3
+  cells in the $x$ direction, and of 6 cells along the $y$ axis.
+- The ~m3d~ variable designates a mesh of $]-1,1[^3$ made of $3\times2\times4$
+  cells.
+
+
+***** ~readGmsh: string -> mesh~
+
+Reads a ~mesh~ from a file.
+#+BEGIN_warning
+The file must conform to the mesh format ~msh2~.
+#+END_warning
+
+#+BEGIN_SRC shell :exports results :results none
+cat << EOF > hybrid-2d.geo
+//+
+Point(1) = {0, 0, 0, 1.0};
+//+
+Point(2) = {1, 0, 0, 1.0};
+//+
+Point(3) = {1, 1, 0, 1.0};
+//+
+Point(4) = {0, 1, 0, 1.0};
+//+
+Point(5) = {2, 0, 0, 1.0};
+//+
+Point(6) = {2, 1, 0, 1.0};
+//+
+Line(1) = {4, 1};
+//+
+Line(2) = {2, 3};
+//+
+Line(3) = {5, 6};
+//+
+Line(4) = {3, 4};
+//+
+Line(5) = {6, 3};
+//+
+Line(6) = {1, 2};
+//+
+Line(7) = {2, 5};
+//+
+Curve Loop(1) = {1, 6, 2, 4};
+//+
+Surface(1) = {1};
+//+
+Curve Loop(2) = {-2, 7, 3, 5};
+//+
+Surface(2) = {2};
+//+
+Characteristic Length {4, 3, 6, 5, 2, 1} = 0.3;
+//+
+Recombine Surface {1};
+//+
+Physical Curve("XMIN") = {1};
+//+
+Physical Curve("XMAX") = {3};
+//+
+Physical Curve("YMAX") = {4, 5};
+//+
+Physical Curve("YMIN") = {6, 7};
+//+
+Physical Curve("INTERFACE") = {2};
+//+
+Physical Surface("LEFT") = {1};
+//+
+Physical Surface("RIGHT") = {2};
+//+
+Physical Point("XMINYMIN") = {1};
+//+
+Physical Point("XMINYMAX") = {4};
+//+
+Physical Point("XMAXYMIN") = {5};
+//+
+Physical Point("XMAXYMAX") = {6};
+#+END_SRC
+
+#+BEGIN_SRC shell :exports results :results none
+gmsh -2 hybrid-2d.geo -format msh2
+#+END_SRC
+
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+  import writer;
+
+  let m:mesh, m = readGmsh("hybrid-2d.msh");
+  write_mesh(gnuplot_writer("hybrid-2d", 0), m);
+#+END_SRC
+
+The ~mesh~ is represented on Figure [[fig:gmsh-hybrid-2d]].
+
+#+NAME: gmsh-hybrid-2d-img
+#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/gmsh-hybrid-2d.png")
+  reset
+  unset grid
+  unset border
+  unset key
+  unset xtics
+  unset ytics
+  set terminal png truecolor enhanced size 960,480
+  plot '<(sed "" $PUGS_SOURCE_DIR/doc/hybrid-2d.0000.gnu)' w l
+#+END_SRC
+
+#+CAPTION: The mesh that was read from the file ~hydrid-2d.msh~ and then saved to the ~gnuplot~ format
+#+NAME: fig:gmsh-hybrid-2d
+#+ATTR_LATEX: :width 0.38\textwidth
+#+ATTR_HTML: :width 300px;
+#+RESULTS: gmsh-hybrid-2d-img
+
+***** ~diamondDual: mesh -> mesh~
+
+This function creates the diamond dual ~mesh~ of a primal ~mesh~.
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+  import writer;
+
+  let primal_mesh:mesh, primal_mesh = readGmsh("hybrid-2d.msh");
+  let dual_mesh:mesh, dual_mesh = diamondDual(primal_mesh);
+  write_mesh(gnuplot_writer("diamond", 0), dual_mesh);
+#+END_SRC
+
+The diamond dual mesh is defined by joining the nodes of the faces to
+the center of the adjacent cells of the primal mesh.
+
+The ~mesh~ is represented on Figure [[fig:gmsh-hybrid-2d]].
+
+#+NAME: diamond-dual-img
+#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/diamond-dual.png")
+  reset
+  unset grid
+  unset border
+  unset key
+  unset xtics
+  unset ytics
+  set terminal png truecolor enhanced size 960,480
+  plot '<(sed "" $PUGS_SOURCE_DIR/doc/hybrid-2d.0000.gnu)' lt rgb "green" w l, '<(sed "" $PUGS_SOURCE_DIR/doc/diamond.0000.gnu)'  lt rgb "black" w l
+#+END_SRC
+
+#+CAPTION: The primal mesh in green and the diamond dual mesh in black
+#+NAME: fig:diamond-dual
+#+ATTR_LATEX: :width 0.38\textwidth
+#+ATTR_HTML: :width 300px;
+#+RESULTS: diamond-dual-img
+
+#+BEGIN_note
+The mesh storage mechanisms in ~pugs~ is such that the diamond dual mesh
+is built only once. This means that is one writes for instance
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+
+  let primal_mesh:mesh, primal_mesh = readGmsh("hybrid-2d.msh");
+  let dual_mesh_1:mesh, dual_mesh_1 = diamondDual(primal_mesh);
+  let dual_mesh_2:mesh, dual_mesh_2 = diamondDual(primal_mesh);
+#+END_SRC
+then the meshes ~dual_mesh_1~ and ~dual_mesh_2~ will refer to the same
+~mesh~ in memory.
+#+END_note
+
+#+BEGIN_warning
+The diamond dual mesh construction is not yet implemented in parallel
+#+END_warning
+
+***** ~medianDual: mesh -> mesh~
+
+This function creates the median dual ~mesh~ of a primal ~mesh~.
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+  import writer;
+
+  let primal_mesh:mesh, primal_mesh = readGmsh("hybrid-2d.msh");
+  let dual_mesh:mesh, dual_mesh = medianDual(primal_mesh);
+  write_mesh(gnuplot_writer("median", 0), dual_mesh);
+#+END_SRC
+
+The median dual mesh is defined by joining the centers of the faces to
+the centers of the adjacent cells of the primal mesh.
+
+The ~mesh~ is represented on Figure [[fig:gmsh-hybrid-2d]].
+
+#+NAME: median-dual-img
+#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/median-dual.png")
+  reset
+  unset grid
+  unset border
+  unset key
+  unset xtics
+  unset ytics
+  set terminal png truecolor enhanced size 960,480
+  plot '<(sed "" $PUGS_SOURCE_DIR/doc/hybrid-2d.0000.gnu)' lt rgb "green" w l, '<(sed "" $PUGS_SOURCE_DIR/doc/median.0000.gnu)'  lt rgb "black" w l
+#+END_SRC
+
+#+CAPTION: The primal mesh in green and the median dual mesh in black
+#+NAME: fig:median-dual
+#+ATTR_LATEX: :width 0.38\textwidth
+#+ATTR_HTML: :width 300px;
+#+RESULTS: median-dual-img
+
+#+BEGIN_note
+In ~pugs~, the storage mechanisms of median dual meshes follows the same
+rules as the diamond dual meshes. As long as the primary mesh lives
+and as long as the median dual mesh is referred, it is kept in memory,
+thus constructed only once.
+#+END_note
+
+#+BEGIN_warning
+The median dual mesh construction is not yet implemented in 3d and not
+available in parallel
+#+END_warning
+
+***** ~transform: mesh*function -> mesh~
+
+This function allows to compute a new mesh as the transformation a
+given mesh by displacing its nodes through a user defined function.
+
+For a mesh of dimension $d$, the mesh must be a function $\mathbb{R}^d
+\to\mathbb{R}^d$.
+
+#+BEGIN_SRC pugs :exports both :results none
+import mesh;
+import math;
+import writer;
+
+let m0:mesh, m0 = cartesianMesh([0,0], [1,1], (10,10));
+
+let pi_over_3: R, pi_over_3 = acos(-1)/3;
+let t: R^2 -> R^2, x -> [(0.2+0.8*x[0])*cos(pi_over_3*x[1]), (0.2+0.8*x[0])*sin(pi_over_3*x[1])];
+
+let m1: mesh, m1 = transform(m0, t);
+write_mesh(gnuplot_writer("transformed", 0), m1);
+#+END_SRC
+
+#+BEGIN_note
+One should keep in mind that the mesh produced the ~transform~ function
+*shares* the same connectivity than the given mesh. This means that in
+~pugs~ internals, there is only one connectivity object.
+#+END_note
+
+The result of the previous script is given on Figure [[fig:transformed]].
+
+#+NAME: transformed-img
+#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/transformed.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/transformed.0000.gnu)' w l
+#+END_SRC
+
+#+CAPTION: The transformed unit square
+#+NAME: fig:transformed
+#+ATTR_LATEX: :width 0.38\textwidth
+#+ATTR_HTML: :width 300px;
+#+RESULTS: transformed-img
+
+***** ~relax: mesh*mesh*R -> mesh~
+
+This function is a simple utility that computes a ~mesh~ as the /mean/ of
+two other mesh that share the same connectivity. The connectivity must
+be the *same in memory*, this means that constructing two identical
+meshes with /equivalent/ connectivity is not allowed.
+
+Thus for instance, the following code
+#+NAME: relax-with-similar-connecticities
+#+BEGIN_SRC pugs-error :exports both :results output
+  import mesh;
+  let m1:mesh, m1 = cartesianMesh(0, [1,1], (10,10));
+  let m2:mesh, m2 = cartesianMesh(0, [1,1], (10,10));
+  let m3:mesh, m3 = relax(m1,m2,0.3);
+#+END_SRC
+produces the runtime error
+#+results: relax-with-similar-connecticities
+
+A proper use is
+#+BEGIN_SRC pugs :exports both :results none
+import mesh;
+import math;
+import writer;
+
+let m0:mesh, m0 = cartesianMesh([0,0], [1,1], (10,10));
+write_mesh(gnuplot_writer("relax_example_m0", 0), m0);
+
+let t: R^2 -> R^2, x -> [x[0]+0.25*x[1]*x[1], x[1]+0.3*sin(x[0])];
+let m1: mesh, m1 = transform(m0, t);
+write_mesh(gnuplot_writer("relax_example_m1", 0), m1);
+
+let m2: mesh, m2 = relax(m0, m1, 0.3);
+write_mesh(gnuplot_writer("relax_example_m2", 0), m2);
+#+END_SRC
+
+In this example, the relaxation parameter is set to $\theta=0.3$, this
+means that the coordinates of any vertex of the relaxed mesh
+$\mathcal{M}_2$, are given by
+\begin{equation*}
+\forall r\in\mathcal{R},\quad\mathbf{x}_r^{\mathcal{M}_2} =  (1-\theta) \mathbf{x}_r^{\mathcal{M}_0} + \theta \mathbf{x}_r^{\mathcal{M}_1}
+\end{equation*}
+
+The different meshes produced in this example are displayed on Figure [[fig:relax]].
+
+#+NAME: relax-img
+#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/relax.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/relax_example_m0.0000.gnu)' lt rgb "green" w l, '<(sed "" $PUGS_SOURCE_DIR/doc/relax_example_m1.0000.gnu)' lt rgb "blue"  w l, '<(sed "" $PUGS_SOURCE_DIR/doc/relax_example_m2.0000.gnu)' lt rgb "black" w l
+#+END_SRC
+
+#+CAPTION: Example of meshes relaxation. The relaxed mesh $\mathcal{M}_2$ (black) and the original meshes in green ($\mathcal{M}_0$) and blue ($\mathcal{M}_1$).
+#+NAME: fig:relax
+#+ATTR_LATEX: :width 0.38\textwidth
+#+ATTR_HTML: :width 300px;
+#+RESULTS: relax-img
+
+This function is mainly useful to reduce the displacement of nodes
+when using the ~randomizeMesh~ functions (see section [[scheme]]).
+
+**** The ~scheme~ module<<scheme>>
+
 
 [fn:pugs-def] ~pugs~: Parallel Unstructured Grid Solvers
 [fn:MPI-def] ~MPI~: Message Passing Interface
-[fn:DSL-def] ~DSL~: Domain Specific Language
+[fn:DSL-def] ~DSL~: Domain Specific Language~