From 14c1acab8b1aeae8fbb8425f7860ecdb2bedc87c Mon Sep 17 00:00:00 2001 From: Stephane Del Pino <stephane.delpino44@gmail.com> Date: Fri, 15 Jul 2022 16:36:16 +0200 Subject: [PATCH] Add description of scheme module: types and functions --- doc/userdoc.org | 889 +++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 882 insertions(+), 7 deletions(-) diff --git a/doc/userdoc.org b/doc/userdoc.org index e34c7c8db..a1ec61d21 100644 --- a/doc/userdoc.org +++ b/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 -- GitLab