#+SETUPFILE: ../packages/org-themes/src/readtheorg_inline/readtheorg_inline.theme #+STARTUP: org-pretty-entities entitiespretty #+PROPERTY: header-args :comments no #+OPTIONS: timestamp:nil #+OPTIONS: h:3 num:t toc:3 #+TITLE: The pugs user manual #+OPTIONS: author:nil date:nil #+OPTIONS: tex:t #+LANGUAGE: en #+ATTR_LATEX: :width 4cm #+HTML_HEAD_EXTRA: <style> pre.src-pugs:before { content: 'pugs'; } </style> #+HTML_HEAD_EXTRA: <style> pre.src-pugs-error:before { content: 'invalid pugs'; } </style> #+HTML_HEAD_EXTRA: <style> .footref{ color: #2980b9; font-size: 100%; } </style> #+LATEX_CLASS_OPTIONS: [10pt] #+LATEX_HEADER: \usepackage[hmargin=2.5cm,vmargin=1.5cm]{geometry} #+LATEX_HEADER: \usepackage{ae,lmodern} #+LATEX_HEADER: \usepackage[OT1]{fontenc} #+LATEX_HEADER: \usepackage{booktabs} #+LATEX_HEADER: \usepackage{siunitx} #+LATEX_HEADER: \usepackage[table]{xcolor} #+LATEX_HEADER: \usepackage{colortbl} #+LATEX_COMPILER: pdflatex --shell-escape #+LATEX_HEADER_EXTRA: \usepackage{amsmath} #+LATEX_HEADER_EXTRA: \usepackage{amsthm} #+LATEX_HEADER_EXTRA: \usepackage{amssymb} #+LATEX_HEADER_EXTRA: \usepackage{xcolor} #+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!10]} #+LATEX_HEADER_EXTRA: \AfterEndEnvironment{minted}{\end{mdframed}} #+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=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!10]} #+LATEX_HEADER_EXTRA: \AfterEndEnvironment{warning}{\end{mdframed}} #+LATEX_HEADER_EXTRA: \usepackage{mathpazo} * Introduction ~pugs~[fn:pugs-def] is a general purpose solver collection built to approximate solutions of partial differential equations. It is mainly (but not only) designed to deal with finite-volume methods. ~pugs~ is a parallel software that uses ~MPI~[fn:MPI-def] and multi-threading for parallelism. Multi-threading is achieved through an encapsulation of some [[https:github.com/kokkos/kokkos][Kokkos]] mechanisms. The philosophy of ~pugs~ is to provide "simple" numerical tools that are assembled together through a high-level language (a DSL[fn:DSL-def] close to the mathematics) to build more complex solvers. This approach is inspired by the success of [[http://freefem.org][FreeFEM]], which use a similar approach. Before detailing the leading concepts and choices that we have made to develop ~pugs~, we give a simple example. ** An example For instance, the following code builds a uniform Cartesian grid of $]-1,1[^2$ made of $20\times20$ cells and transforms it by displacing the mesh nodes according to a user-defined vector field $T: \mathbb{R}^2 \to \mathbb{R}^2$. #+NAME: introduction-example #+BEGIN_SRC pugs :exports both :results output import mesh; import math; import writer; let pi:R, pi = acos(-1); let theta:R^2 -> R, x -> 0.5*pi*(x[0]*x[0]-1)*(x[1]*x[1]-1); 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 = cartesianMesh([-1,-1], [1,1], (20,20)); m = transform(m, T); write_mesh(gnuplot_writer("transformed"), m); #+END_SRC The example is quite easy to read. - First, some *modules* are loaded: the ~mesh~ module, which contains some mesh manipulation functions. The ~math~ module provides a set of classical mathematical functions ($\sin$, $\cos$, ...). The ~writer~ module is used to save meshes or discrete functions to files using various formats. - The second block of data defines variables of different types - ~pi~ is a real value that is initialized by an approximation of $\pi$. - ~theta~ is the real value function $\theta$ defined by \begin{equation*} \theta: \mathbb{R}^2 \to \mathbb{R},\quad\mathbf{x} \mapsto \frac{\pi}{2} (x_0^2-1)(x_1^2-1) \end{equation*} - ~M~ is the $\mathbb{R}^{2\times2}$ matrix field $M$ defined by \begin{equation*} M: \mathbb{R}^2 \to \mathbb{R}^{2\times2},\quad\mathbf{x} \mapsto \begin{pmatrix} \cos(\theta(\mathbf{x})) & -\sin(\theta(\mathbf{x})) \\ \sin(\theta(\mathbf{x})) & \cos(\theta(\mathbf{x})) \end{pmatrix} \end{equation*} - ~T~ is the vector field $T$ defined by \begin{equation*} T: \mathbb{R}^2 \to \mathbb{R}^2, \quad\mathbf{x} \mapsto (I+M(\mathbf{x}))\mathbf{x} \end{equation*} - Finally ~m~ is defined as the uniform Cartesian mesh grid of $]-1,1[^2$. The last argument: ~(20,20)~ sets the number of cells in each direction: $20$. - The third block is the application of the transformation $T$ to the mesh. Observe that if the resulting mesh is stored in the same variable ~m~, the old one was not modified in the process. It is important to already have in mind that the ~pugs~ language *does not allow* the modifications of values of variables of *non-basic* types. This is discussed in the section [[high-level-types]]. - Finally, the last block consists in saving the obtained mesh in a ~gnuplot~ file. The result is shown on Figure [[fig:intro-example]]. #+NAME: intro-transform-mesh-img #+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/intro-transform-mesh.png") reset unset grid unset border unset key unset xtics unset ytics set terminal png truecolor enhanced set size square plot '<(sed "" $PUGS_SOURCE_DIR/doc/transformed.gnu)' w l #+END_SRC #+CAPTION: Obtained transformed mesh #+NAME: fig:intro-example #+ATTR_LATEX: :width 0.38\textwidth #+ATTR_HTML: :width 300px; #+RESULTS: intro-transform-mesh-img Even if this first example is very simple, some key aspects can already be discussed. - There is no predefined constant in ~pugs~. Here a value is provided for ~pi~. - There are two kinds of variable in ~pugs~: variables of basic types and variable of high-level types. This two kinds of variable behave almost the same but one must know their differences to understand better the underlying mechanisms and choices that we made. See 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~, ~cartesianMesh~,...) are builtin functions and are generally defined when importing a module. These functions behave similarly, one should refer to [[functions]] for details. ** Concepts and design As it was already stated, ~pugs~ can be viewed as a collection of numerical methods or utilities that can be assembled together, using a user friendly language, to build simulation scenarios. Utilities are tools that are often used in numerical simulations. Examples of such utilities are mesh manipulations, definition of initial conditions, post-processing, error calculations,... *** A C++ toolbox driven by a user friendly language Numerical simulation packages are software of a particular kind. Generally, in order to run a calculation, one has to define a set of data and parameters. This can simply be definition of a discretization parameter such as the mesh size. One can also specify boundary conditions, equations of state, source terms for a specific model. Choosing a numerical method or even more, setting the model itself, is common in large codes. In ~pugs~, all these "parameters" are set through a DSL[fn:DSL-def]. Thus, when ~pugs~ is launched, it actually executes a provided script. A ~C++~ function is associated to each instruction of the script. The ~C++~ components of ~pugs~ are completely unaware one of the others. ~pugs~ interpreter is responsible of data flow between the components: it manages the data transfer between those ~C++~ components and ensures that the workflow is properly defined. **** Why? In this section we motivate the choice of a language and not of a more standard approach. ***** Data files are evil There are lots of reasons not to use data files. By data file, we refer to a set of options that describe physical models, numerical methods or their settings. - Data files are not flexible. This implies in then one hand that application scenarios must be known somehow precisely to reflect possible option combinations and in the other hand even defining a specific initial data may require the creation of a new option and the associated code (in ~C++~ for instance). \\ Usually, the last point is addressed by adding a local interpreter to evaluate user functions. - Data files may contain irrelevant information. Actually, it is quite common to allow to define options that are valid but irrelevant to a given scenario. This policy can be changed but it is generally not an easy task and requires more work from the user (which can be a good thing). - Generally data files become rapidly obsolete. An option was not the right one, or its type changed to allow other contexts... This puts pressure on the user. - Even worst, options meaning can depend on other options. Unfortunately, this happens commonly. For instance, a global option can change implicitly the treatment associated to another one. This is dangerous since writing or reading the data file requires an important and up to date knowledge of the code's internals. - Another major difficulty when dealing with data files is to check the compatibility of provided options. ***** Embedded "data files" are not a solution Using directly the general purpose language (~C~, ~C++~, ~Fortran~,...) used to write the code can be tempting. It has the advantage that no particular treatment is necessary to build a parser (to read data files or a script), but it presents several drawbacks. - The first one is probably that it allows to much freedom. While defining the model and numerical options, the user has generally access to the whole code and can change almost everything, even things that should not be changed. - Again, one can easily have access to irrelevant options and it requires a great knowledge of the code to find important ones. - With that regard, defining a simulation properly can be a difficult task. For instance, in the early developments of ~pugs~ (when it was just a raw ~C++~ code) it was tricky to change boundary conditions for coupled physics. - Another difficulty is related to the fact that code's internal API is likely to change permanently in a research code. Thus valid constructions or setting may become rapidly obsolete. In other words keeping up to date embedded "data file" might be difficult. - Finally it requires recompilation of pieces of code (which can be large in some cases) even if one is just changing a simple parameter. ***** Benefits of a DSL Actually, an honest analysis cannot conclude that a DSL is the solution to all problems. However, it offers some advantages. - First, it allows a fine control on what the user can or cannot perform. In some sense, it offers a chosen level of flexibility. - It allows to structure the code in the sense that new developments have to be designed not only focusing on the result but also on the way it should be used (its interactions with the scripting language). - In the same idea, it provides a framework that drives the desired principle of "do simple things and do them well". - There are no hidden dependencies between numerical options: the DSL code is easier to read (than data files) and is less likely to become obsolete (this is not that true in early developments since the language itself and some concepts are still likely to change). - The simulation scenario is *defined* by the script, it is the responsibility of the user to check its meaning (not to the charge of the code). ***** ~pugs~ language purpose ~pugs~ language is used to assemble ~C++~ tools which are well tested. These tools should ideally be small pieces of ~C++~ code that do one single thing and do it well. Another purpose of the language is to allow to perform high-level calculations. In other words, the language defines a data flow and checks that each ~C++~ piece of code is used correctly. Since each piece of code acts as a pure function (arguments are unchanged by calls), the calling context is quite easy to check. Finally it aims at simplifying the definition of new methods since common utilities are available directly in scripts. ***** The framework: divide and conquer ~pugs~ is a research oriented software, thus generally the user is also a developer. If this paragraph is indeed more dedicated to the developer, by reading it, the user will have a better understanding of the development choices and the underlying policy that the code follows. Actually the development framework imposed by the DSL tends to guide writing of new methods. - In the process of writing a *new numerical methods*, one must create *new functions in the language*. Moreover, if a new method is close to an existing one, it is generally *better* to use completely new underlying ~C++~ code than to patch existing methods. Starting from a *copy* of the existing code ~C++~ is *encouraged* for developments. This may sound weird since classical development guidelines encourage inheritance or early redesign. Actually, this policy is the result of the following discussion. - Using this approach, one *separates* clearly the *design* of a numerical method and the *design* of the code. - From the computer science point of view, early design for new numerical methods is generally wrong: usually one cannot anticipate precisely enough eventual problems or method corrections. - It is much more difficult to introduce bugs in existing methods, since previously validated methods are unchanged! - For the same reason, existing methods performances are naturally unchanged by new developments. - Also, when comparing methods, it is better to compare to the original existing code. - If the new method is abandoned or simply not kept, the existing code is not polluted. - Finally when a method is validated and ready to integrate the mainline sources of the code, it is easier to see differences with existing ones and *at this time* one can redesign the ~C++~ code checking that results are unchanged and performances not deteriorated. At this time, it is likely that the numerical method design is finished, thus (re)designing the source code makes more sense. - Another consequence is that utilities are not be developed again and again. - This implies an /a priori/ safer code: utilities are well tested and validated. - It saves development time obviously. - The code of numerical methods is not polluted by environment instructions (data initialization, error calculation, post-processing,...) - The counterpart is somehow classical. In the one hand, the knowledge of existing utilities is required, this document tries to address a part of it. In the other hand, if the developer requires a new utility, a good practice is to discuss with the other ones to check if it could benefit to them. Then one can determine if it should integrate rapidly or not the main development branch. ***** Why not python or any other scripting language? As it was already pointed out above, general purpose languages offer to much freedom: it is not easy to protect data. For instance in the ~pugs~ DSL, non basic variables are constant (see paragraph [[high-level-types]]). It is important since it prevents the user from modifying data in inconsistent ways. Also, one must keep in mind that constraining the expressiveness is actually a strength. As said before, one can warranty coherence of the data, perform calculations without paying attention to the parallelism aspects,... Observe that it is not a limitation: if the DSL's field of application needs to be extended, it is always possible. But these extensions should never break the rule that a DSL must not become a general purpose language. Providing a language close to the application reduces the gap between the application and its expression (code): Domain Specific Languages are made for that. Finally, python is ugly. *** A high-level language Following the previous discussion, the reader should now understand the motivations that drove the design choices that conduct to build ~pugs~ as a ~C++~ toolbox driven by a user friendly language. #+begin_verse Keep it simple and relevant! #+end_verse In the development process of an application, the easy step is generally the implementation itself. It is even more true when the developer feels that changes to the code are natural and that the modifications themselves look easy. Obviously, any experienced programmer knows that writing the code is only the first step of a much longer process which requires rigorous tests and validations (the enrichment of a non-regression database). Generally one also desires to ensure that the development is used within the correct bounds which requires to implement data verification. In a perfect world, an up-to-date documentation of the functionality and its domain of validity should be provided. This is even more true when defining a language (or a DSL). Enriching a language syntax (or grammar) is not something that must be done to answer a specific need. It must not be done /because it is possible to do it/! #+begin_verse When designing a language, the difficulty is not to offer new functionalities,\\ it is generally to decide not to offer them.\\ --- Bjarne Stroustrup, C++ conference 2021. #+end_verse If the grammar of ~pugs~ is still likely to be extended, it should *never* integrate low-level instructions. Low-level instructions give too much freedom and thus are a source of errors. Several things are already done to forbid this kind of evolution. The constness of high-level data is a good illustration. For instance, meshes or discrete functions *cannot* be modified. This is not only a security to protect the user from doing "dangerous" manipulations, but it also permits to define high-level optimizations. - Since meshes are constant objects, one can for instance compute geometric data on demand. These data are kept into memory as long as the mesh lives. - Forbidding the modification of values of a discrete function ensures that parallel communication instructions should never appear in a ~pugs~ script. Another benefit of not providing low-level instructions is that the scripts are more easy to write and read, and it is more difficult to write errors. * Language ** Variables In order to simplify the presentation, before going further, we introduce the syntax that allows to print data to the terminal. It follows ~C++~ streams construction for convenience. For instance #+NAME: cout-preamble-example #+BEGIN_SRC pugs :exports both :results output cout << "2+3 = " << 2+3 << "\n"; #+END_SRC produces the following output #+results: cout-preamble-example The code is quite obvious for ~C++~ users, note that ~"\n"~ is the linefeed string (there is no character type in ~pugs~, just strings). Actually, ~cout~ is itself a variable, we will come to this later. ~pugs~ is a strongly typed language. It means that a variable *cannot* change of type in its lifetime. *** Declaration and affectation syntax **** Declaration of simple variables To declare a variable ~v~ of a given type ~V~, one writes #+BEGIN_SRC pugs :exports source let v:V; #+END_SRC This instruction is read as #+begin_verse Let $v\in V$. #+end_verse Actually, - ~let~ is the declaration keyword, - ~v~ is the variable name, - ~:~ is a separation token (it can be read as "/in/" in this context), - ~V~ is the identifier of the type, and - ~;~ marks the end of the instruction. For instance to declare a real variable ~x~, one writes #+NAME: declare-R #+BEGIN_SRC pugs :exports both :results none let x:R; #+END_SRC The instructions that follow the declaration of a variable can use it while defined in the same scope, but it cannot be used before. Also, after its declaration, one cannot declare another variable with the same name in the same scope (see [[blocks-and-life-time]]). #+NAME: redeclare-variable #+BEGIN_SRC pugs-error :exports both :results output let x:R; // first declaration let x:R; // second declaration #+END_SRC produces the following error #+results: redeclare-variable **** Affectation of simple variables To affect the value of an expression ~expression~ to variable ~v~ one simply uses the ~=~ operator. Thus assuming that a variable ~v~ has already been /declared/, one writes simply #+BEGIN_SRC pugs :exports source v = expression; #+END_SRC There is not to much to comment, reading is quite natural - ~v~ is the variable name, - ~=~ is the affectation operator, - ~expression~ is some code that provides a value of the same type as ~v~ (it can be another variable, an arithmetic expression, the result of a function,...), and - ~;~ marks the end of the instruction. For instance, one can write #+NAME: simple-affectations-example #+BEGIN_SRC pugs :exports both :results output import math; // to load the sin function let a:N; let b:N; let x:R; a=3; b=2+a; x=sin(b)+a; cout << "a = " << a << " b = " << b << " x = " << x << "\n"; #+END_SRC In this example, we import the ~math~ module which provides the ~sin~ function and we use the ~N~ data type of natural integers ($\mathbb{N}\equiv\mathbb{Z}_{\ge0}$). Running the example gives the following result. #+results: simple-affectations-example #+BEGIN_warning Actually, *separating* the *declaration* from the *initialization* of a variable is quite *dangerous*. This is prone to errors and can lead to the use of undefined values. During the compilation of scripts, ~pugs~ detects uninitialized values. For instance, #+NAME: uninitialized-variable #+BEGIN_SRC pugs-error :exports both :results output let x:R; let y:R; y=x+1; #+END_SRC produces the following compilation error #+results: uninitialized-variable For more complex constructions, it can be very difficult to detect it at /compile time/, this is why it is *encouraged* to use variable definition (see [[definition-simple-variables]]). Observe nonetheless that ~pugs~ checks at /run time/ that used variables are correctly defined. If not, a /runtime/ error is produced. #+END_warning **** Definition of simple variables<<definition-simple-variables>> The best way to define variables is the following. To define a variable ~v~ of a given type ~V~, from an expression one writes #+BEGIN_SRC pugs :exports source let v:V, v = expression; #+END_SRC - ~let~ is the declaration keyword, - ~v~ is the variable name, - ~:~ is a separation token (can be read as "/in/" in this context), - ~V~ is the identifier of the type, - ~,~ is the separator that indicates the beginning of the affectation, - ~expression~ is some code that provides a value of type ~V~ (it can be another variable, an expression, the result of a function,...), and - ~;~ marks the end of the instruction (here the definition of ~v~). A practical example is #+NAME: simple-definition-example #+BEGIN_SRC pugs :exports both :results output import math; // to load the sin function let a:N, a=3; let b:N, b=2+a; let x:R, x=sin(b)+a; cout << "a = " << a << " b = " << b << " x = " << x << "\n"; #+END_SRC which produces the result #+results: simple-definition-example *** Blocks and lifetime of variables<<blocks-and-life-time>> In pugs scripts, variables have a precise lifetime. They are defined within scopes. The main scope is implicitly defined and sub-scopes are enclosed between curly bracket pairs: ~{~ and ~}~. Following ~C++~, a variable exists (can be used) as soon as it has been declared and until the end of the scope (where it has been declared). At this point we give a few examples. **** A variable cannot be used before its declaration #+NAME: undeclare-variable #+BEGIN_SRC pugs-error :exports both :results output n = 3; let n:N; #+END_SRC #+results: undeclare-variable **** A variable cannot be used after its definition scope #+NAME: out-of-scope-variable-use #+BEGIN_SRC pugs-error :exports both :results output { let n:N, n = 3; n = n+2; } cout << n << "\n"; #+END_SRC #+results: out-of-scope-variable-use **** Variable name can be reused in an enclosed scope #+NAME: nested-scope-variable-example #+BEGIN_SRC pugs :exports both :results output let n:N, n = 0; // global variable { cout << "global scope n = " << n << "\n"; let n:N, n = 1; // scope level 1 variable { cout << "scope level 1 n = " << n << "\n"; let n:N, n = 2; // scope level 2 variable cout << "scope level 2 n = " << n << "\n"; } { cout << "scope level 1 n = " << n << "\n"; let n:N, n = 4; // scope level 2.2 variable cout << "scope level 2.2 n = " << n << "\n"; } cout << "scope level 1 n = " << n << "\n"; } cout << "global scope n = " << n << "\n"; #+END_SRC #+results: nested-scope-variable-example This example is self explanatory. Obviously such constructions are generally *bad ideas*. This kind of constructions can appear in loops where the variables defined in blocks follow the same lifetime rules. *** Basic types<<basic-types>> Basic types in ~pugs~ are boolean (~B~), natural integers (~N~), integers (~Z~), real (~R~), small vectors (~R^1~, ~R^2~ and ~R^3~), small matrices (~R^1x1~, ~R^2x2~ and ~R^3x3~) and strings (~string~). #+BEGIN_note Observe that while mathematically, obviously $\mathbb{R} = \mathbb{R}^1 = \mathbb{R}^{1\times1}$, the data types ~R~, ~R^1~ and ~R^1x1~ are different in ~pugs~ and are *not implicitly* convertible from one to the other! This may sound strange but there are few reasons for that. - First, these are the reflect of internal ~pugs~ ~C++~-data types that are used to write algorithms. In its core design pugs aim at writing numerical methods generically with regard to the dimension. One of the ingredients to achieve this purpose is to use dimension $1$ vectors and matrices when some algorithms reduce to dimension $1$ instead of ~double~ values. To avoid ambiguity that may arise in some situations (this can lead to very tricky code), we decided to forbid automatic conversions of these types with ~double~. When designing the language, we adopted the same rule to avoid ambiguity. - A second reason is connected to the first one. Since ~pugs~ aims at providing numerical methods for problems in dimension $1$, $2$ or $3$, this allows to distinguish the nature of the underlying objects. - It is natural to consider that the coordinates of the vertices defining a mesh in dimension $d$ are elements of $\mathbb{R}^d$, - or that a velocity or a displacement are also defined as $\mathbb{R}^d$ values. Thus using ~R^1~ in dimension $1$ for this kind data precise their nature in some sense . #+END_note **** Expression types The language associates statically some types to several special expressions. - Special boolean (type ~B~) expressions. Two *keywords* allow to define boolean values: ~true~ and ~false~. - Integers (type ~Z~) are defined as a contiguous list of digits. - For instance, the code ~0123~ is interpreted as the integer $123$. - However the sequence ~123 456~ is *not interpreted* as $123456$ but as the two integers $123$ and $456$. - Real values (type ~R~) use the same syntax as in ~C++~. For instance, the following expressions are accepted to define the number $1.23$. #+BEGIN_SRC pugs :exports both 1.23; 0.123E1; .123e+1; 123e-2; 12.3E-1; #+END_SRC - Small vectors values (types ~R^1~, ~R^2~ and ~R^3~) are written through square brackets. Each component of vector values must be a scalar expression. #+BEGIN_SRC pugs :exports both [1]; // R^1 value [1,2]; // R^2 value [1,2,3]; // R^3 value #+END_SRC - Small matrices values (types ~R^1x1~, ~R^2x2~ and ~R^3x3~) are written by lines through square brackets. Each line is enclosed between square brackets. #+BEGIN_SRC pugs :exports both [[1]]; // R^1x1 value [[1,2],[3,4]]; // R^2x2 value [[1,2,3],[4,5,6],[7,8,9]]; // R^3x3 value #+END_SRC - ~string~ values are defined as the set of characters enclosed between two double quotes ( ~"~ ). The string /Hello world!/ would be simply written as ~"Hello world!"~. Strings support the following escape sequences (similarly to ~C++~): #+LATEX: \definecolor{contiYellow}{RGB}{255,165,0} #+LATEX: \rowcolors[]{2}{contiYellow!5}{contiYellow!20} | escape | meaning | |--------+-----------------| | ~\'~ | single quote | | ~\"~ | double quote | | ~\?~ | question mark | | ~\\~ | backslash | | ~\a~ | audible bell | | ~\b~ | backspace | | ~\f~ | new page | | ~\n~ | new line | | ~\r~ | carriage return | | ~\t~ | horizontal tab | | ~\v~ | vertical tab | These special characters are not interpreted by ~pugs~ itself but interpreted by the system when an output is created. They are just allowed in the definition of a ~string~. **** Variables of basic types are stored by value This means that each variable of basic type allocates its own memory to store its data. This is the natural behavior for variables. It is illustrated in the following example #+NAME: basic-type-value-storage #+BEGIN_SRC pugs :exports both :results output let a:N, a = 3; let b:N, b = a; a = 1; cout << "a = " << a << " b = " << b << "\n"; #+END_SRC which produces the expected result #+results: basic-type-value-storage When defining ~b~, the *value* contained in ~a~ is copied to set the value of ~b~. Thus changing ~a~'s value does not impact the variable ~b~. **** Variables of basic types are mutable In ~pugs~ the only variables that are mutable (their value can be *modified*) are of basic types. Executing the following code #+NAME: basic-type-mutable-value #+BEGIN_SRC pugs :exports both :results output let a:N, a = 2; a += 3; cout << "a = " << a << "\n"; #+END_SRC gives #+results: basic-type-mutable-value which is not a surprise. However, the use of the ~+=~ operator results in the modification of the stored value. There is no copy. Actually, this is not really important from the user point of view. One just have to keep in mind that, as it will be depicted after, high-level variables *are not mutable*: their values can be *replaced* by a new ones but *cannot be modified*. *** Implicit type conversions<<implicit-conversion>> In order to avoid ambiguities, in ~pugs~, there is *no implicit* conversion in general. #+BEGIN_note Actually, there are only three situations for which implicit type conversion can occur: when the value is - given as a parameter of a function, - used as the returned value of a function, or - used to define a tuple This will be illustrated in section [[functions]] and [[tuples]] #+END_note This means that all affectations, unary operators and binary operators are defined explicitly for supported types. Here is a table of implicit type conversions *when allowed*. | expected type | convertible types | |---------------+------------------------------------------------| | ~N~ | ~B~, ~Z~ | | ~Z~ | ~B~, ~N~ | | ~R~ | ~B~, ~N~, ~Z~ | | ~R^1~ | ~B~, ~N~, ~Z~, ~R~ | | ~R^2~ | ~0~ (special value) | | ~R^3~ | ~0~ (special value) | | ~R^1x1~ | ~B~, ~N~, ~Z~, ~R~ | | ~R^2x2~ | ~0~ (special value) | | ~R^3x3~ | ~0~ (special value) | | ~string~ | ~B~, ~N~, ~Z~, ~R~, ~R^1~, ~R^2~, ~R^3~, ~R^1x1~, ~R^2x2~, ~R^3x3~ | *** Operators **** Affectation operators In the ~pugs~ language, the affectation operators are the following. | operator | description | |----------+----------------------| | ~=~ | affectation operator | |----------+----------------------| | ~+=~ | increment operator | | ~-=~ | decrement operator | | ~*=~ | multiply operator | | ~/=~ | divide operator | #+BEGIN_note It is important to note that in ~pugs~ language, affectation operators have *no* return value. This is a notable difference with ~C~ or ~C++~. Again, this is done to avoid common mistakes that can be difficult to address. For instance, the following ~C++~ code is valid but does not produce the expected result. #+BEGIN_SRC C++ :exports source bool b = true; // do things ... if (b=false) { // do other things ... } #+END_SRC Obviously the mistake is that the test should have been ~(b==false)~, since otherwise, the conditional block is never executed (in ~C++~, the result of an affectation is the value affected to the variable, which is always false in that case. This cannot happen with ~pugs~. A similar example #+NAME: no-affectation-result #+BEGIN_SRC pugs-error :exports both :results output let b:B, b=true; // do things if (b=false) { // do other things } #+END_SRC produces the following error #+results: no-affectation-result Actually, affectations are /expressions/ in ~C++~. In ~pugs~, affectations are /instructions/. #+END_note ***** List of defined operator ~=~ for basic types. As already mentioned, operator ~=~ is defined for *all* types in ~pugs~ if the right hand side expression has the *same* type as the left hand side variable. This is true for basic types as well as for high-level types. We now give the complete list of supported ~=~ affectations. The lists are sorted by type of left hand side variable. - ~B~: boolean left hand side variable. One is only allowed to affect boolean values. | ~B =~ allowed expression type | |-----------------------------| | ~B~ | - ~N~: natural integer ($\mathbb{N}$ or $\mathbb{Z}_{\ge0}$) left hand side variable. | ~N =~ allowed expression type | |-----------------------------| | ~B~ | | ~N~ | | ~Z~ (for convenience) | - ~Z~: integer ($\mathbb{Z}$) left hand side variable. | ~Z =~ allowed expression type | |-----------------------------| | ~B~ | | ~N~ | | ~Z~ | - ~R~: real ($\mathbb{R}$) left hand side variable. | ~R =~ allowed expression type | |-----------------------------| | ~B~ | | ~N~ | | ~Z~ | | ~R~ | - ~R^1~: vector of dimension 1 ($\mathbb{R}^1$) left hand side variable. | ~R^1 =~ allowed expression type | |-------------------------------| | ~B~ | | ~N~ | | ~Z~ | | ~R~ | | ~R^1~ | - ~R^2~: vector of dimension 2 ($\mathbb{R}^2$) left hand side variable. | ~R^2 =~ allowed expression type | |---------------------------------------------| | ~R^2~ | | ~0~ (special value) | | list of 2 scalar (~B~, ~N~, ~Z~ or ~R~) expressions | An example of initialization using an $\mathbb{R}^2$ value or the special value ~0~ is #+NAME: affectation-to-R2-from-list #+BEGIN_SRC pugs :exports both :results output let u:R^2, u = [-3, 2.5]; let z:R^2, z = 0; cout << "u = " << u << "\n"; cout << "z = " << z << "\n"; #+END_SRC which produces #+RESULTS: affectation-to-R2-from-list Observe that ~0~ is a special value and *treated as a keyword*. There is no conversion from integer values. For instance: #+NAME: R2-invalid-integer-affectation #+BEGIN_SRC pugs-error :exports both :results output let z:R^2, z = 1-1; #+END_SRC produces the compilation error #+results: R2-invalid-integer-affectation - ~R^3~: vector of dimension 3 ($\mathbb{R}^3$) left hand side variable. | ~R^3 =~ allowed expression type | |---------------------------------------------| | ~R^3~ | | ~0~ (special value) | | list of 3 scalar (~B~, ~N~, ~Z~ or ~R~) expressions | An example of initialization is #+NAME: affectation-to-R3-from-list #+BEGIN_SRC pugs :exports both :results output let u:R^3, u = [-3, 2.5, 1E-2]; let z:R^3, z = 0; cout << "u = " << u << "\n"; cout << "z = " << z << "\n"; #+END_SRC the output is #+RESULTS: affectation-to-R3-from-list - ~R^1x1~: matrix of dimensions $1\times1$ ($\mathbb{R}^{1\times1}$) left hand side variable. | ~R^1x1 =~ allowed expression type | |---------------------------------| | ~B~ | | ~N~ | | ~Z~ | | ~R~ | | ~R^1x1~ | - ~R^2x2~: matrix of dimension $2\times2$ ($\mathbb{R}^{2\times2}$) left hand side variable. | ~R^2x2 =~ allowed expression type | |---------------------------------------------| | ~R^2x2~ | | ~0~ (special value) | | list of 4 scalar (~B~, ~N~, ~Z~ or ~R~) expressions | An example of initialization using an $\mathbb{R}^{2\times2}$ value or the special value ~0~ is #+NAME: affectation-to-R2x2-from-list #+BEGIN_SRC pugs :exports both :results output let u:R^2x2, u = [[-3, 2.5], [ 4, 1.2]]; let z:R^2x2, z = 0; cout << "u = " << u << "\n"; cout << "z = " << z << "\n"; #+END_SRC which produces #+RESULTS: affectation-to-R2x2-from-list - ~R^3x3~: matrix of dimension $3\times3$ ($\mathbb{R}^{3\times3}$) left hand side variable. | ~R^3x3 =~ allowed expression type | |---------------------------------------------| | ~R^3x3~ | | ~0~ (special value) | | list of 9 scalar (~B~, ~N~, ~Z~ or ~R~) expressions | An example of initialization is #+NAME: affectation-to-R3x3-from-list #+BEGIN_SRC pugs :exports both :results output let u:R^3x3, u = [[ -3, 2.5, 1E-2], [ 2, 1.7, -2], [1.2, 4, 2.3]]; let z:R^3x3, z = 0; cout << "u = " << u << "\n"; cout << "z = " << z << "\n"; #+END_SRC the output is #+RESULTS: affectation-to-R3x3-from-list - ~string~ left hand side variable. Expressions of any basic type can be used as the right hand side. | ~string =~ allowed expression type | |----------------------------------| | ~B~ | | ~N~ | | ~Z~ | | ~R~ | | ~R^1~ | | ~R^2~ | | ~R^3~ | | ~R^1x1~ | | ~R^2x2~ | | ~R^3x3~ | | ~string~ | ***** List of defined operator ~+=~ for basic types. - ~B~: the ~+=~ operator is not defined for left hand side boolean variables. - ~N~: natural integer ($\mathbb{N}$ or $\mathbb{Z}_{\ge0}$) left hand side variable. | ~N +=~ allowed expression type | |------------------------------| | ~B~ | | ~N~ | | ~Z~ (for convenience) | - ~Z~: integer ($\mathbb{Z}$) left hand side variable. | ~Z +=~ allowed expression type | |------------------------------| | ~B~ | | ~N~ | | ~Z~ | - ~R~: real ($\mathbb{R}$) left hand side variable. | ~R +=~ allowed expression type | |------------------------------| | ~B~ | | ~N~ | | ~Z~ | | ~R~ | - ~R^1~: vector of dimension 1 ($\mathbb{R}^1$) left hand side variable. | ~R^1 +=~ allowed expression type | |--------------------------------| | ~R^1~ | - ~R^2~: vector of dimension 2 ($\mathbb{R}^2$) left hand side variable. | ~R^2 +=~ allowed expression type | |--------------------------------| | ~R^2~ | - ~R^3~: vector of dimension 3 ($\mathbb{R}^3$) left hand side variable. | ~R^3 +=~ allowed expression type | |--------------------------------| | ~R^3~ | - ~R^1x1~: matrix of dimensions $1\times1$ ($\mathbb{R}^{1\times1}$) left hand side variable. | ~R^1x1 +=~ allowed expression type | |----------------------------------| | ~R^1x1~ | - ~R^2x2~: matrix of dimension $2\times2$ ($\mathbb{R}^{2\times2}$) left hand side variable. | ~R^2x2 +=~ allowed expression type | |----------------------------------| | ~R^2x2~ | - ~R^3x3~: matrix of dimension $3\times3$ ($\mathbb{R}^{3\times3}$) left hand side variable. | ~R^3x3 +=~ allowed expression type | |----------------------------------| | ~R^3x3~ | - ~string~ left hand side variable. Expressions of any basic type can be used as the right hand side. | ~string +=~ allowed expression type | |-----------------------------------| | ~B~ | | ~N~ | | ~Z~ | | ~R~ | | ~R^1~ | | ~R^2~ | | ~R^3~ | | ~R^1x1~ | | ~R^2x2~ | | ~R^3x3~ | | ~string~ | ***** List of defined operator ~-=~ for basic types. - ~B~: the ~-=~ operator is not defined for left hand side boolean variables. - ~N~: natural integer ($\mathbb{N}$ or $\mathbb{Z}_{\ge0}$) left hand side variable. | ~N -=~ allowed expression type | |------------------------------| | ~B~ | | ~N~ | | ~Z~ (for convenience) | - ~Z~: integer ($\mathbb{Z}$) left hand side variable. | ~Z -=~ allowed expression type | |------------------------------| | ~B~ | | ~N~ | | ~Z~ | - ~R~: real ($\mathbb{R}$) left hand side variable. | ~R -=~ allowed expression type | |------------------------------| | ~B~ | | ~N~ | | ~Z~ | | ~R~ | - ~R^1~: vector of dimension 1 ($\mathbb{R}^1$) left hand side variable. | ~R^1 -=~ allowed expression type | |--------------------------------| | ~R^1~ | - ~R^2~: vector of dimension 2 ($\mathbb{R}^2$) left hand side variable. | ~R^2 -=~ allowed expression type | |--------------------------------| | ~R^2~ | - ~R^3~: vector of dimension 3 ($\mathbb{R}^3$) left hand side variable. | ~R^3 -=~ allowed expression type | |--------------------------------| | ~R^3~ | - ~R^1x1~: matrix of dimensions $1\times1$ ($\mathbb{R}^{1\times1}$) left hand side variable. | ~R^1x1 -=~ allowed expression type | |----------------------------------| | ~R^1x1~ | - ~R^2x2~: matrix of dimension $2\times2$ ($\mathbb{R}^{2\times2}$) left hand side variable. | ~R^2x2 -=~ allowed expression type | |----------------------------------| | ~R^2x2~ | - ~R^3x3~: matrix of dimension $3\times3$ ($\mathbb{R}^{3\times3}$) left hand side variable. | ~R^3x3 -=~ allowed expression type | |----------------------------------| | ~R^3x3~ | - ~string~: the ~-=~ operator is not defined for left hand side string variables. ***** List of defined operator ~*=~ for basic types. - ~B~: the ~*=~ operator is not defined for left hand side boolean variables. - ~N~: natural integer ($\mathbb{N}$ or $\mathbb{Z}_{\ge0}$) left hand side variable. | ~N *=~ allowed expression type | |------------------------------| | ~B~ | | ~N~ | | ~Z~ (for convenience) | - ~Z~: integer ($\mathbb{Z}$) left hand side variable. | ~Z *=~ allowed expression type | |------------------------------| | ~B~ | | ~N~ | | ~Z~ | - ~R~: real ($\mathbb{R}$) left hand side variable. | ~R *=~ allowed expression type | |------------------------------| | ~B~ | | ~N~ | | ~Z~ | | ~R~ | - ~R^1~: vector of dimension 1 ($\mathbb{R}^1$) left hand side variable. | ~R^1 *=~ allowed expression type | |--------------------------------| | ~B~ | | ~N~ | | ~Z~ | | ~R~ | - ~R^2~: vector of dimension 2 ($\mathbb{R}^2$) left hand side variable. | ~R^2 *=~ allowed expression type | |--------------------------------| | ~B~ | | ~N~ | | ~Z~ | | ~R~ | - ~R^3~: vector of dimension 3 ($\mathbb{R}^3$) left hand side variable. | ~R^3 *=~ allowed expression type | |--------------------------------| | ~B~ | | ~N~ | | ~Z~ | | ~R~ | - ~R^1x1~: matrix of dimensions $1\times1$ ($\mathbb{R}^{1\times1}$) left hand side variable. | ~R^1x1 *=~ allowed expression type | |----------------------------------| | ~B~ | | ~N~ | | ~Z~ | | ~R~ | - ~R^2x2~: matrix of dimension $2\times2$ ($\mathbb{R}^{2\times2}$) left hand side variable. | ~R^2x2 *=~ allowed expression type | |----------------------------------| | ~B~ | | ~N~ | | ~Z~ | | ~R~ | - ~R^3x3~: matrix of dimension $3\times3$ ($\mathbb{R}^{3\times3}$) left hand side variable. | ~R^3x3 *=~ allowed expression type | |----------------------------------| | ~B~ | | ~N~ | | ~Z~ | | ~R~ | #+BEGIN_note Observe that for these small matrix types ($\mathbb{R}^{d\times d}$) the construction ~A *= B;~ where ~B~ is a matrix of the same type as ~A~ is not allowed. The main reason for that is that for $d>1$ this operation has no interests since it requires a temporary. One will see bellow that it is possible to write ~A = A*B;~ if needed. #+END_note - ~string~: the ~*=~ operator is not defined for left hand side string variables. ***** List of defined operator ~/=~ for basic types. - ~B~: the ~/=~ operator is not defined for left hand side boolean variables. - ~N~: natural integer ($\mathbb{N}$ or $\mathbb{Z}_{\ge0}$) left hand side variable. | ~N /=~ allowed expression type | |------------------------------| | ~N~ | | ~Z~ (for convenience) | - ~Z~: integer ($\mathbb{Z}$) left hand side variable. | ~Z /=~ allowed expression type | |------------------------------| | ~N~ | | ~Z~ | - ~R~: real ($\mathbb{R}$) left hand side variable. | ~R /=~ allowed expression type | |------------------------------| | ~N~ | | ~Z~ | | ~R~ | - ~R^d~: the ~/=~ operator is not defined for left hand side vector (of dimension $d\in\{1,2,3\}$) variables. - ~R^dxd~: the ~/=~ operator is not defined for left hand side matrix (of dimension $d\times d$ with $d\in\{1,2,3\}$) variables. **** Unary operators The ~pugs~ language allows the following tokens as unary operators | operator | description | |----------+----------------------| | ~not~ | not operator | | ~+~ | plus unary operator | | ~-~ | minus unary operator | |----------+----------------------| | ~++~ | increment operator | | ~--~ | decrement operator | |----------+----------------------| | ~[]~ | access operator | The ~not~, ~+~ and ~-~ operators apply to the *expression* on their right. ~++~ and ~--~ operators apply only to a *variable* that can be positioned before (pre increment/decrement) or after the token (post increment/decrement). These operators are also inspired from their ~C++~ counterparts for commodity. The ~+~ unary operator is a convenient operator that is *elided* when parsing the script. For basic types, when the operators ~not~, ~+~, ~-~, ~++~ or ~--~ are defined, they return a value of the same type as the argument (except for the operator ~-~ if the argument is a ~B~ or a ~N~, then the result is a ~Z~). These operators can be defined for high-level types. - The ~not~ operator is only defined for boolean values (~B~). - The ~-~ unary operator is defined for numeric basic types: ~B~, ~N~, ~Z~, ~R~, ~R^1~, ~R^2~, ~R^3~, ~R^1x1~, ~R^2x2~ and ~R^3x3~. It is not defined for ~string~ variables. - Pre and post increment operators, ~--~ and ~++~, are defined for integer types: ~N~ and ~Z~. They are not defined for ~B~, ~R~, ~R^1~, ~R^2~, ~R^3~, ~R^1x1~, ~R^2x2~, ~R^3x3~ and ~string~ variables. Note that the pre increment/decrement operators behave slightly differently than their ~C++~ counterparts since they are not allowed to be chained. In ~C++~, the following code is allowed #+BEGIN_SRC C++ :exports source int i = 0; int j = ++ ++i; #+END_SRC In ~pugs~, it is forbidden: #+NAME: double-pre-incr-result #+BEGIN_SRC pugs-error :exports both :results output let i:N, i=0; let j:N, j = ++ ++i; #+END_SRC produces the compilation error #+results: double-pre-incr-result Again, this is done to simplify the syntax and to avoid weird constructions. - Access operators are only defined for small vectors ~R^d~ and small matrices ~R^dxd~. To avoid use of uninitialized variables (or partially uninitialized variables), these are ~read-only~ access operators. Their syntax is the following. #+NAME: Rd-Rdxd-access-operator #+BEGIN_SRC pugs :exports both :results output let x:R^2, x = [1,2]; let A:R^3x3, A = [[1,2,3],[4,5,6],[7,8,9]]; cout << "x[0] = " << x[0] << "\nx[1] = " << x[1] << "\n"; cout << "A[0,0] = " << A[0,0] << "\nA[2,1] = " << A[2,1] << "\n"; #+END_SRC This code produces #+results: Rd-Rdxd-access-operator **** Binary operators Syntax for binary operators follows again a classical structure: if ~exp1~ and ~exp2~ denotes two expressions and if ~op~ denotes a binary operator, one simply writes ~exp1 op exp2~. Here is the list of binary operators | keyword | operator | |---------+-----------------------| | ~and~ | logic and | | ~or~ | logic or | | ~xor~ | logic exclusive or | |---------+-----------------------| | ~==~ | equality | | ~!=~ | non-equality | | ~<~ | lower than | | ~<=~ | lower or equal than | | ~>~ | greater than | | ~>=~ | greater or equal than | |---------+-----------------------| | ~<<~ | shift left | | ~>>~ | shift right | |---------+-----------------------| | ~+~ | sum | | ~-~ | difference | | ~*~ | product | | ~/~ | division | Binary operators can be defined for high-level types. For basic types, they follow a few rules. - Logical operators ~and~, ~or~ and ~xor~ are defined for boolean operands (type is ~B~) only. The result of the expression is a boolean. #+begin_src latex :results drawer :exports results \begin{equation*} \left| \begin{array}{rl} \mathtt{and}:&\quad \mathbb{B} \times \mathbb{B} \to \mathbb{B}\\ \mathtt{or}:& \quad\mathbb{B} \times \mathbb{B} \to \mathbb{B}\\ \mathtt{xor}:& \quad \mathbb{B} \times \mathbb{B} \to \mathbb{B} \end{array} \right. \end{equation*} #+end_src - Comparison operators ~==~, ~!=~, ~<~, ~<=~, ~>~ and ~>=~ are defined for all basic scalar type and return a boolean value. #+begin_src latex :results drawer :exports results \begin{equation*} \forall \mathbb{S}_1, \mathbb{S}_2 \in \{\mathbb{B},\mathbb{N},\mathbb{Z},\mathbb{R}\}, \quad \left| \begin{array}{rl} \mathtt{==}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{B}\\ \mathtt{!=}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{B}\\ \mathtt{<}: & \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{B}\\ \mathtt{<=}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{B}\\ \mathtt{>}: & \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{B}\\ \mathtt{>=}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{B} \end{array} \right. \end{equation*} #+end_src When comparing a boolean value (type ~B~) with another scalar value type (~N~, ~Z~ or ~R~), the value ~true~ is interpreted as $1$ and the value ~false~ as $0$. \\ For vector and matrix basic types, the only allowed operators are ~==~ and ~!=~. #+begin_src latex :results drawer :exports results \begin{equation*} \forall d \in \{1,2,3\},\quad \left| \begin{array}{rl} \mathtt{==}:& \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{B}\\ \mathtt{==}:& \mathbb{R}^{d \times d} \times \mathbb{R}^{d \times d} \to \mathbb{B}\\ \mathtt{!=}:& \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{B}\\ \mathtt{!=}:& \mathbb{R}^{d \times d} \times \mathbb{R}^{d \times d} \to \mathbb{B} \end{array} \right. \end{equation*} #+end_src \\ This is also the case for ~string~ values: only allowed operators are ~==~ and ~!=~. #+begin_src latex :results drawer :exports results \begin{equation*} \left| \begin{array}{rl} \mathtt{==}:& \mathtt{string} \times \mathtt{string} \to \mathbb{B}\\ \mathtt{!=}:& \mathtt{string} \times \mathtt{string} \to \mathbb{B} \end{array} \right. \end{equation*} #+end_src - Shift operators ~<<~ and ~>>~ are not used to define binary operators between two basic types. - Arithmetic operators (~+~, ~-~, ~*~ and ~/~) are defined for a set of combinations of basic types. We classify them by their returned types. - Their is no arithmetic operation that returns a boolean ~B~. - Operators that return a natural integer ~N~. #+begin_src latex :results drawer :exports results \begin{equation*} \forall \mathbb{S}_1, \mathbb{S}_2 \in \{\mathbb{B},\mathbb{N}\}, \quad \left| \begin{array}{rl} \mathtt{+}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{N}\\ \mathtt{*}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{N} \end{array} \right. \end{equation*} #+end_src Boolean values (type ~B~) are not allowed as the right operand of a division. #+begin_src latex :results drawer :exports results \begin{equation*} \forall \mathbb{S} \in \{\mathbb{B},\mathbb{N}\}, \mathtt{/}: \mathbb{S} \times \mathbb{N} \to \mathbb{N} \end{equation*} #+end_src Observe that ~-~ is *not* in the list. - Operators that return an integer ~Z~. #+begin_src latex :results drawer :exports results \begin{equation*} \forall \mathbb{S}_1, \mathbb{S}_2 \in \{\mathbb{B},\mathbb{N},\mathbb{Z}\}, \mbox{ such that }\mathbb{S}_1 = \mathbb{Z}\mbox{ or }\mathbb{S}_2 = \mathbb{Z}, \quad \left| \begin{array}{rl} \mathtt{+}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{Z}\\ \mathtt{-}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{Z}\\ \mathtt{*}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{Z}\\ \mathtt{/}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{Z} \end{array} \right. \end{equation*} #+end_src The divide operator does not allow boolean values as a right operand. #+begin_src latex :results drawer :exports results \begin{align*} \forall \mathbb{S}_1 \in \{\mathbb{B},\mathbb{N},\mathbb{Z}\}&\mbox{ and } \mathbb{S}_2 \in \{\mathbb{N},\mathbb{Z}\}, \mbox{ such that }\mathbb{S}_1 = \mathbb{Z}\mbox{ or }\mathbb{S}_2 = \mathbb{Z},\\ &\mathtt{/}: \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{Z}. \end{align*} #+end_src Finally the following operator is also defined #+begin_src latex :results drawer :exports results \begin{equation*} \mathtt{-}: \mathbb{N} \times \mathbb{N} \to \mathbb{Z} \end{equation*} #+end_src The difference of two natural integers returns a /signed/ integer. - Operators that return a real ~R~. #+begin_src latex :results drawer :exports results \begin{equation*} \forall \mathbb{S}_1, \mathbb{S}_2 \in \{\mathbb{B},\mathbb{N},\mathbb{Z},\mathbb{R}\}, \mbox{ such that }\mathbb{S}_1 = \mathbb{R}\mbox{ or }\mathbb{S}_2 = \mathbb{R}, \quad \left| \begin{array}{rl} \mathtt{+}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{R}\\ \mathtt{-}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{R}\\ \mathtt{*}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{R} \end{array} \right. \end{equation*} #+end_src The divide operator does not allow boolean values as a right operand. #+begin_src latex :results drawer :exports results \begin{align*} \forall \mathbb{S}_1 \in \{\mathbb{B},\mathbb{N},\mathbb{Z},\mathbb{R}\}&\mbox{ and } \mathbb{S}_2 \in \{\mathbb{N},\mathbb{Z},\mathbb{R}\}, \mbox{ such that }\mathbb{S}_1 = \mathbb{R}\mbox{ or }\mathbb{S}_2 = \mathbb{R},\\ &\mathtt{/}: \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{Z}. \end{align*} #+end_src - Operators that return a small vector ~R^d~. #+begin_src latex :results drawer :exports results \begin{equation*} \forall d \in\{1,2,3\}, \quad \left| \begin{array}{rl} \mathtt{+}:& \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}^d\\ \mathtt{-}:& \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}^d\\ \mathtt{*}:& \mathbb{B} \times \mathbb{R}^d \to \mathbb{R}^d\\ \mathtt{*}:& \mathbb{N} \times \mathbb{R}^d \to \mathbb{R}^d\\ \mathtt{*}:& \mathbb{Z} \times \mathbb{R}^d \to \mathbb{R}^d\\ \mathtt{*}:& \mathbb{R} \times \mathbb{R}^d \to \mathbb{R}^d\\ \mathtt{*}:& \mathbb{R}^{d\times d} \times \mathbb{R}^d \to \mathbb{R}^d \end{array} \right. \end{equation*} #+end_src - Operators that return a small matrix ~R^dxd~. #+begin_src latex :results drawer :exports results \begin{equation*} \forall d \in\{1,2,3\}, \quad \left| \begin{array}{rl} \mathtt{+}:& \mathbb{R}^{d\times d} \times \mathbb{R}^{d\times d} \to \mathbb{R}^{d\times d}\\ \mathtt{-}:& \mathbb{R}^{d\times d} \times \mathbb{R}^{d\times d} \to \mathbb{R}^{d\times d}\\ \mathtt{*}:& \mathbb{B} \times \mathbb{R}^{d\times d} \to \mathbb{R}^{d\times d}\\ \mathtt{*}:& \mathbb{N} \times \mathbb{R}^{d\times d} \to \mathbb{R}^{d\times d}\\ \mathtt{*}:& \mathbb{Z} \times \mathbb{R}^{d\times d} \to \mathbb{R}^{d\times d}\\ \mathtt{*}:& \mathbb{R} \times \mathbb{R}^{d\times d} \to \mathbb{R}^{d\times d}\\ \mathtt{*}:& \mathbb{R}^{d\times d} \times \mathbb{R}^{d\times d} \to \mathbb{R}^{d\times d} \end{array} \right. \end{equation*} #+end_src - Operators that return a ~string~. #+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 \left| \begin{array}{rl} \mathtt{+}:& \mathbb{S} \times \mathtt{string} \to \mathtt{string}\\ \mathtt{+}:& \mathtt{string} \times \mathbb{S} \to \mathtt{string}\\ \mathtt{+}:& \mathtt{string} \times \mathtt{string} \to \mathtt{string} \end{array} \right. \end{equation*} #+end_src **** Operator precedence and associativity To avoid confusions, the operators precedence in ~pugs~ language follows the same rules as in ~C++~. This is summarized in the following table, where ~a~ and ~b~ denotes two expressions. | Precedence | Operator | |------------+----------| | 1 | ~a++~ | | | ~a--~ | | | ~a[]~ | |------------+----------| | 2 | ~++a~ | | | ~--a~ | | | ~not a~ | |------------+----------| | 3 | ~a*b~ | | | ~a/b~ | |------------+----------| | 4 | ~a+b~ | | | ~a-b~ | |------------+----------| | 5 | ~a<<b~ | | | ~a>>b~ | |------------+----------| | 6 | ~a<b~ | | | ~a<=b~ | | | ~a>b~ | | | ~a>=b~ | |------------+----------| | 7 | ~a==b~ | | | ~a!=b~ | |------------+----------| | 8 | ~a=b~ | | | ~a+=b~ | | | ~a-=b~ | | | ~a*=b~ | | | ~a/=b~ | As already said, we forbid some constructions to try to avoid spurious constructions. By construction associativity of many operators makes no sense in ~pugs~ language (affectations and increment/decrement operators for instance). For all other operators, the associativity rule is *left to right*. Thus the following code #+BEGIN_SRC pugs :exports source - 1 + 3 - 4 + 2; #+END_SRC is equivalent to #+BEGIN_SRC pugs :exports source (((- 1) + 3) - 4) + 2; #+END_SRC Obviously ~pugs~ allows the use of parenthesis to write expressions. It is enough to give a simple example. #+NAME: parenthesis-precedence #+BEGIN_SRC pugs :exports both :results output cout << " 2 + 3 * 4 = " << 2 + 3 * 4 << "\n"; cout << "(2 + 3) * 4 = " << (2 + 3) * 4 << "\n"; #+END_SRC the output is #+RESULTS: parenthesis-precedence *** High-level types<<high-level-types>> Aside from the basic types described in the previous section, ~pugs~ also deals with "high-level" types. This term is more to understand as "non-basic types". The ~pugs~ language is not object oriented to keep it simple. To fix ideas, let us give examples of a few high-level data types. These can be meshes (the ~mesh~ type), output streams (the ~ostream~ type), boundary descriptors (~boundary~), quadrature formula (~quadrature~), discrete functions (~Vh~),... One can already see that the complexity of these types may vary a lot. The main difference between these types and basic types is that, high-level types are not available in directly the language but they are loaded on demand (using the ~import~ keyword in the preamble of the script). The second difference is that data of these types are *constant*. More precisely, the content of high-level variables can be replaced by a new one *but* it cannot be modified. For this reason, the following operators can never be applied to variables of this kind | forbidden operators | description | |---------------------+-----------------------------| | ~++~ | pre/post increment operator | | ~--~ | pre/post decrement operator | | ~+=~ | assignment by sum | | ~-=~ | assignment by difference | | ~*=~ | assignment by product | | ~/=~ | assignment by quotient | We conclude by stating that if access operator ~[]~ can eventually be overloaded for high-level types, it should be done with care. It is not recommended. Finally, the last difference lies in the fact that high-level types use shallow copies and not value copies as it is the case for basic types. This is transparent to the user and provides the intuitive (similar) behavior since data of high-level variables are constant. To illustrate this, let us consider the following example. #+BEGIN_SRC pugs :exports both import mesh; let m1:mesh, m1 = cartesianMesh(0, [1,1], (10,10)); let m2:mesh, m2 = m1; let m3:mesh, m3 = cartesianMesh(0, [1,1], (10,10)); m3 = m1; #+END_SRC In this example, we are dealing with 3 ~mesh~ variables. - First, ~m1~ is defined as a uniform cartesian mesh in dimension 2 of $]0,1[^2$ made of $10\times10$ identical square cells. Let us call this resident mesh $\mathcal{M}_a$ for clarity. - Next, the variable ~m2~ is defined as a copy of the variable ~m1~. It is a copy of the variable, but its content is *the same*. The variable ~m2~ is also referring to the mesh $\mathcal{M}_a$: there is no duplication of the mesh in memory. - Then, the variable ~m3~ is defined to refer to a *new* mesh (called $\mathcal{M}_b$. It is /similar/ to $\mathcal{M}_a$, but it is *not* the same one! When ~m3~ is defined two meshes $\mathcal{M}_a$ and $\mathcal{M}_b$ (and two distinct connectivities) are resident in memory. - Finally, the last instruction (~m3 = m1;~) sets ~m3~ to also refer $\mathcal{M}_a$. Since no other variable refers to $\mathcal{M}_b$, this mesh is destroyed (memory is freed). At the end of the program, all the variables ~m1~, ~m2~ and ~m3~ are referring to $\mathcal{M}_a$ which is the only mesh that resides in memory. *** Compound types The ~pugs~ language allows to deal with compound types. The idea is to define a list of variables as a member of a product space (each variable belongs to one of them, each variable has a simple type: basic or high-level). **** Compound declaration Let us provide an example to fix ideas. #+NAME: compound-declaration #+BEGIN_SRC pugs :exports both :results output let (A,x,n): R^2x2*R^3*N; A = [[1,2],[3,4]]; x = [2,4,6]; n = 2; cout << "A = " << A << "\nx = " << x << "\nn = " << n << "\n"; #+END_SRC the output is #+RESULTS: compound-declaration This is completely equivalent to declaring the variables one after the other. **** Compound definition One can also use the following definition instruction #+NAME: compound-definition #+BEGIN_SRC pugs :exports both :results output let (A,x,n): R^2x2*R^3*N, (A,x,n) = ([[1,2],[3,4]], [2,4,6], 2); cout << "A = " << A << "\nx = " << x << "\nn = " << n << "\n"; #+END_SRC which gives the same result #+RESULTS: compound-definition A potential mistake with this construction, is that variables must appear in the same order in the declaration part and in the affectation part of a definition. For instance, the following code is invalid. #+NAME: invalid-compound-definition #+BEGIN_SRC pugs-error :exports both :results output let (x,y):R*R, (y,x) = (0,1); #+END_SRC It produces the following error #+results: invalid-compound-definition which is easy to fix. Another potential mistake is that all variables of the list are marked as defined at the same time. Thus one cannot use construction like this: #+NAME: undeclared-compound-definition #+BEGIN_SRC pugs-error :exports both :results output let (x,y):R*R, (x,y) = (0,2+x); #+END_SRC It produces the following error #+results: undeclared-compound-definition While the variable ~x~ is defined *before* ~y~, this kind of construction is forbidden. From a technical point of view, this behavior would be easy to change (allow to use the fresh value of ~x~ in the definition of ~y~), but this make the code unclear and this is not the purpose of compound types. #+BEGIN_note Observe that there is no implicit conversion when dealing with compound types. The ~=~ operators are used sequentially to set the different variables. #+END_note **** Compound affectation The last way to use compound types is illustrated by the following example. #+NAME: compound-affectation #+BEGIN_SRC pugs :exports both :results output let A: R^2x2; let x: R^3; let n: N; (A,x,n) = ([[1,2],[3,4]], [2,4,6], 2); cout << "A = " << A << "\nx = " << x << "\nn = " << n << "\n"; #+END_SRC It produces again #+results: compound-affectation #+BEGIN_note Observe that the only operator allowed for this kind of construction is the operator ~=~. #+END_note **** Use of compound types Actually using compound types the way it is presented in this paragraph is not recommend. The purpose of compounds types and compound affectations is related to functions. As one will see in section [[functions]], functions can return compound values, thus compound affectations (or definitions) are needed to get returned values in that case. *** Tuple types<<tuples>> The last special type construction is the ability to deal with tuples in the ~pugs~ language. The tuples we are describing here are lists of data of a *unique* and *simple* type (one cannot use compound types to define tuples). Tuples in the ~pugs~ language have *no relationship* with ~C++~'s ~std::tuple~ these are two completely different notions. The list of values given to the tuple must be *implicitly* convertible to the type of tuple elements (see the conversion table in section [[implicit-conversion]]). There is no ambiguity, the implicit conversions follow the rules of operator ~=~. **** Tuple declaration and affectation For instance, one may write #+NAME: tuple-declaration-affectation #+BEGIN_SRC pugs :exports both :results output let x:(R); x = (1,2,3.4,2); cout << x << "\n"; #+END_SRC Executing this code, one gets #+results: tuple-declaration-affectation **** Tuple definition The definition syntax is also possible. #+NAME: tuple-definition #+BEGIN_SRC pugs :exports both :results output let x:(R^2), x = ([1,2],[3.4,2], 0, [2,3]); cout << x << "\n"; #+END_SRC This code gives #+results: tuple-definition Observe that the special value ~0~ is used there. **** Tuple purpose Tuples variables are just lists of data of the same type. In the ~pugs~ language, one cannot access to a specific value of the list nor alter one of them. This is not something that should ever change. Tuples are not arrays! The ~pugs~ language is not meant to allow low-level instructions. The use case of tuples is to provide lists of data to the ~C++~ underlying methods. A classical example is to provide a set of boundary conditions to a method. ** Statements The ~pugs~ language supports classical statements to control the data flow. For simplicity, these statements syntax follow their ~C++~ counterpart. The only statement that is not implemented in ~pugs~ is the ~switch...case~. This may change but in the one hand, up to now it was never necessary (up to now, we did not encountered the need to chain ~if...else~ statements), and on the other hand, keeping the language as simple as possible remains the policy in ~pugs~ development. *** ~if...else~ statement The simplest statement is the conditional statement ~if...else~. If a condition is satisfied, the ~truestatement~ (a single instruction or block of instructions) is executed. In the other case, if the /optional/ ~else~ keyword is used, the ~falsestatement~ (a single instruction or block of instructions) is executed. #+BEGIN_SRC pugs :exports code if (condition) truestatement if (condition) truestatement else falsestatement #+END_SRC The condition itself *must be* a boolean value (of type ~B~). #+NAME: if-instruction #+BEGIN_SRC pugs :exports both :results output if (2>1) cout << "yes: 2 > 1\n"; if (not (2>1)) cout << "hmm... 2 <= 1 ?\n"; #+END_SRC generates the following expected output. The second output is hopefully not executed. #+results: if-instruction It is probably better to use the ~else~ keyword to obtain the same result. An example of the use of the ~else~ keyword follows. #+NAME: if-else-instruction #+BEGIN_SRC pugs :exports both :results output if (2<=1) cout << "hmm... 2 <= 1 ?\n"; else cout << "ok: 2 > 1\n"; #+END_SRC #+results: if-else-instruction However, it is recommended even in the case of single instructions in the conditional statements to use instruction blocks: #+NAME: if-else-block-of-one #+BEGIN_SRC pugs :exports both :results output import math; let sin100_positive:B; if (sin(100) > 0){ sin100_positive = true; } else { sin100_positive = false; } cout << "sin(100)>0: " << sin100_positive << "\n"; #+END_SRC #+results: if-else-block-of-one We give a final illustration #+NAME: if-block #+BEGIN_SRC pugs :exports both :results output let m:R, m = 2.5; let M:R, M = 2; let has_swapped:B, has_swapped = false; if (m>M){ let tmp:R, tmp = m; m = M; M = tmp; has_swapped = true; } cout << "min = " << m << " max = " << M << " has_swapped = " << has_swapped << "\n"; #+END_SRC #+results: if-block *** ~for~ loops ~pugs~ allows to write ~for~ loops. It follows the ~C++~ syntax #+BEGIN_SRC pugs :exports code for (declarationinstruction ; condition ; postinstruction) statement #+END_SRC Following ~C++~, ~declarationinstruction~, ~condition~ and ~postinstruction~ are optional. The ~condition~ argument, if it is present, *must* be a boolean value (type ~B~). If it is absent, the default value ~true~ is used. - The ~declarationinstruction~ is execute only *once* /before/ the beginning of the loop. The lifetime of the variable declared here is defined by the ~for~ instruction itself. - The ~condition~ is evaluate /before/ each loop. - The ~postinstruction~ is executed /after/ each loop. - The ~statement~ is either a single instruction or a block of instructions. The ~statement~ is executed if the ~condition~ has the value ~true~. For instance, one can write #+NAME: for-block #+BEGIN_SRC pugs :exports both :results output let n:N, n = 10; let sum:N, sum = 0; for (let i:N, i=1; i<=n; ++i) { sum += i; } cout << "sum = " << sum << "\n"; #+END_SRC which gives as expected #+results: for-block The lifetime of the declared variable (in the ~declarationinstruction~ statement) is illustrated by the following example #+NAME: for-scope-error #+BEGIN_SRC pugs-error :exports both :results output for (let i:N, i=0; i<2; ++i) { cout << "i = " << i << "\n"; } cout << "i = " << i << "\n"; #+END_SRC Running this example produces the following error #+results: for-scope-error To fix the previous code, one can write #+NAME: for-no-decl #+BEGIN_SRC pugs :exports both :results output let i:N, i=0; for (; i<2; ++i) { cout << "i = " << i << "\n"; } cout << "i = " << i << "\n"; #+END_SRC One gets #+results: for-no-decl *** ~do...while~ loops The second kind of loops that is supported is the ~do...while~ construction which executes /at least/ one time the enclosed ~statement~. #+BEGIN_SRC pugs :exports code do statement while (condition); #+END_SRC The ~statement~ is either a single instruction or a block of instructions. The ~condition~ is an expression of boolean value (type ~B~). #+NAME: do-while-block #+BEGIN_SRC pugs :exports both :results output let sum:N, sum = 0; let i:N, i = 0; do { sum += i; ++i; } while (i<=10); cout << "sum = " << sum << "\n"; #+END_SRC It gives #+results: do-while-block *** ~while~ loops The last kind of loops that is allowed in ~pugs~ language is the ~while~ loop. #+BEGIN_SRC pugs :exports code while (condition) statement #+END_SRC The ~statement~ is either a single instruction or a block of instructions. The ~condition~ is an expression of boolean value (type ~B~). This time is the ~condition~ is not satisfied (~false~ when reaching the ~while~ instruction), the ~statment~ is never executed. An example of the ~while~ loop is the following. #+NAME: while-block #+BEGIN_SRC pugs :exports both :results output let sum:N, sum = 0; let i:N, i = 1; while (i<=10) { sum += i; ++i; } cout << "sum = " << sum << "\n"; #+END_SRC The result is #+results: while-block *** ~break~ and ~continue~ keywords. These *keywords* are used to control loops behavior from enclosed loop statements. They follow their ~C++~ counterparts. - The ~continue~ keyword is used to skip the instructions that follow in the loop statement and continue the loop. - The ~break~ keyword leaves the current loop. An example of use of the ~continue~ keyword is #+NAME: nested-continue #+BEGIN_SRC pugs :exports both :results output for (let i:N, i = 0; i < 5; ++i) { cout << i << ": "; for (let j:N, j=0; j < 5; ++j) { if (j<i) continue; cout << j << " "; } cout << "\n"; } #+END_SRC The result is #+results: nested-continue Replacing the ~continue~ keyword by a ~break~ (and changing the test) #+NAME: nested-break #+BEGIN_SRC pugs :exports both :results output for (let i:N, i = 0; i < 5; ++i) { cout << i << ": "; for (let j:N, j=0; j < 5; ++j) { if (j>i) break; cout << j << " "; } cout << "\n"; } #+END_SRC changes the output to #+results: nested-break Obviously the behavior is the same using ~do...while~ or ~while~ loops. ** Functions<<functions>> The ~pugs~ language allows the manipulation and definition of functions. These are *mathematical functions* and are not functions as functions in ~C++~: functions in the ~pugs~ language are *not subroutines*. To be more precise, a function $f$ follows the following structure \begin{align*} \mbox{let }f:& X_1\times \cdots \times X_n \to Y_1\times \cdots \times Y_m,\\ & (x_1,\ldots,x_n)\mapsto (y_1,\ldots,y_m)=f(x_1,\ldots,x_n), \end{align*} where $n,m\in\mathbb{N}_{\ge1}$, and where ${(X_i)}_{1\le i\le n}$ and ${(Y_i)}_{1\le i\le m}$ are /simple/ types. Actually $X_1\times \cdots \times X_n$ and $Y_1\times \cdots \times Y_m$ are /compound/ types. Thus assuming that the function $f$ is defined in the language as ~f~, one can use the following syntax to evaluate it. This is a pseudo-code, real examples will follow. Assuming that ~(x1...,xn)~ has been properly defined in ~X1*...*Xn~ one can write #+BEGIN_SRC pugs :exports code let (y1,...,ym):Y1*...*Ym, (y1,...,ym) = f(x1,...,xn); #+END_SRC or if ~(y1,...,ym)~ has already been declared in ~Y1*...*Ym~ #+BEGIN_SRC pugs :exports code (y1,...,ym) = f(x1,...,xn); #+END_SRC The ~pugs~ language handles two kinds of functions. User-defined functions (see [[user-defined-functions]]) and builtin functions (see [[builtin-functions]]). They behave essentially the same and are both constant. In this context, it means that one cannot just declare or modify a function. *** Pure functions In the ~pugs~ language, functions are *pure functions* in the sense that arguments given to a function are *never* modified by the function. They act as operators. #+BEGIN_note Actually these functions are not strictly /pure functions/ in the computer science sense. The reason is that they can eventually have side effects. As an example, it is possible to modify the random seed used by the code. In that case, the modified value is not a variable of the language itself but the internal random seed. #+END_note *** Implicit type conversion for parameters and returned values As already said, in general, there is no implicit type conversion in the ~pugs~ language. The only exceptions are when initializing tuples (see [[tuples]]), passing arguments to functions and dealing with returned values. See the conversion table in section [[implicit-conversion]]. *** User-defined functions<<user-defined-functions>> To define user functions, the syntax mimics mathematics. The definition of the function itself is a *single* expression. This means that one cannot use ~if...else~, loops and there is not such keyword as ~return~ in the ~C++~. The syntax of the definition of functions is #+BEGIN_SRC pugs :exports code let f1 : X -> Y, x -> e; let f2 : X -> Y1*...*Ym, x -> (e1,...,em); let f3 : X1*...*Xn -> Y, (x1,...,xn) -> e; let f4 : X1*...*Xn -> Y1*...*Ym, (x1,...,xn) -> (e1,...,em); #+END_SRC ~X~, ~Y~, ~X1~, ..., ~Xn~ and ~Y1~, ..., ~Ym~ are /simple/ data types (not compound). The type of ~x~ is ~X~ and for lists, each argument ~xi~ belongs to ~Xi~. The function themselves are defined by the expressions ~e~ (or ~e1~, ..., ~em~) which types are set by ~Y~ (or ~Y1~, ..., ~Ym~). Let us give a few examples. #+NAME: R-to-R-function #+BEGIN_SRC pugs :exports both :results output let f: R -> R, x -> -x; cout << "f(2.3) = " << f(2.3) << "\nf(1) = " << f(1) << "\nf(true) = " << f(true) << "\n"; #+END_SRC One observes the implicit conversion of the values passed by argument. #+results: R-to-R-function The following example illustrates the implicit conversion to the returned type #+NAME: R-to-R1-function #+BEGIN_SRC pugs :exports both :results output let f: R -> R^1, x -> 2*x; cout << "f(3.2) = " << f(3.2) << "\n"; #+END_SRC #+results: R-to-R1-function Using compound types as input and output, one can write #+NAME: R22-R-string-to-R-string-function #+BEGIN_SRC pugs :exports both :results output let f : R^2x2*R*string -> R*string*R^3x3, (A,x,s) -> (x*A[0,0]*A[1,1]-A[1,0], A+","+x+","+s, [[A[0,0], A[0,1], 0], [A[1,0], A[1,1], 0], [ 0, 0, x]]); let x : R, x=0; let s : string, s= ""; let A : R^3x3, A = 0; (x,s,A) = f([[1,2],[3,4]], 2.3, "foo"); cout << "x = " << x << "\ns = " << s << "\nA = " << A << "\n"; let (y,t,A2):R*string*R^3x3, (y,t,A2) = f([[3,1],[4,2]], -5.2, "bar"); cout << "y = " << y << "\nt = " << t << "\nA2 = " << A2 << "\n"; #+END_SRC This meaningless example produces the following result. #+results: R22-R-string-to-R-string-function **** Lifetime of functions arguments The arguments used to define a function are *local* variables that exists only during the evaluation of the function. Let us consider the following example #+NAME: lifetime-of-function-args #+BEGIN_SRC pugs :exports both :results output let a:R, a = 1.4; let plus: R -> R, a -> (a>0)*a; cout << "plus(1) = " << plus(1) << "\n"; cout << "plus(-3.2) = " << plus(-3.2) << "\n"; cout << "a = " << a << "\n"; #+END_SRC This gives the expected result: the value of the variable ~a~ is unchanged. #+results: lifetime-of-function-args **** Non-arguments variables in function expressions Here we discuss rapidly of using variables (which are not arguments) in function expressions. #+NAME: non-arg-variables-in-functions #+BEGIN_SRC pugs :exports both :results output let a:R, a = 1.4; let plus: R -> R, x -> (x>0)*a; cout << "a = " << a << "\n"; cout << "plus(2.2) = " << plus(2.2) << "\n"; cout << "plus(-3.2) = " << plus(-3.2) << "\n"; a = 2; cout << "a = " << a << "\n"; cout << "plus(2.2) = " << plus(2.2) << "\n"; cout << "plus(-3.2) = " << plus(-3.2) << "\n"; #+END_SRC Running the example, one gets #+results: non-arg-variables-in-functions While the function itself is a constant object, one sees that since the value of ~a~ is changed, the value function is implicitly modified. /This is a dangerous feature and should be avoid!/ Since functions themselves are variables one can use functions in function expressions. #+NAME: functions-in-functions #+BEGIN_SRC pugs :exports both :results output let plus: R -> R, x -> (x>0)*x; let minus: R -> R, x -> -(x<0)*x; let pm : R -> R*R, x -> (plus(x), minus(x)); let toR2: R*R -> R^2, (x,y) -> [x,y]; cout << "pm(2) = " << toR2(pm(2)) << " pm(-3) = " << toR2(pm(-3)) << "\n"; #+END_SRC One observes the utility function ~toR2~ that is used to perform the output since ~cout~ does not handle compound types output. One gets #+results: functions-in-functions **** Lifetime of user-defined functions Since functions are somehow variables, the lifetime of functions follows the similar rules. Let us give an example #+NAME: functions-lifetime #+BEGIN_SRC pugs :exports both :results output let f: R -> R, x -> 2.3*x+2; { let f: R -> R, x -> 1.25*x+3; cout << "block 1: f(2) = " << f(2) << "\n"; { let f: R -> R, x -> 3.25*x-0.3; cout << "block 2: f(2) = " << f(2) << "\n"; } } for (let i:N, i = 1; i<=2; ++i) { let f: R -> R, x -> i*x+2; cout << "for(i=" << i << "): f(2) = " << f(2) << "\n"; } cout << "global: f(2) = " << f(2) << "\n"; #+END_SRC Running this example produces #+results: functions-lifetime **** Recursion Since functions cannot be simply declared but must be defined, double recursion is not possible. Thus, there is no equivalent construction in ~pugs~ of the following ~C++~ code #+BEGIN_SRC C++ :exports source int f(int i); int g(int i) { if (i < 0) return 0; else return f(i-1); } int f(int i) { return g(i); } #+END_SRC Moreover simple recursion is forbidden for similar reasons, since functions are made of a single expression (one cannot use statements), there would be no way to end the recursion. Thus, the code #+NAME: no-recursion #+BEGIN_SRC pugs-error :exports both :results output let f: N->N, n -> f(2*n); #+END_SRC produces the following compilation time error #+results: no-recursion *** Builtin functions<<builtin-functions>> In ~pugs~ language, builtin functions are ~C++~ pieces of code that can be called in scripts. There usage is very similar to user-defined functions. They differ from user-defined functions in three points. - Builtin functions can have no parameter or no returned value. - Builtin functions are polymorphic. More precisely, this means that the signature of a builtin function is also defined by its expected arguments types. - Builtin functions can take user-defined functions as parameters. - user-defined functions cannot take functions as parameters - builtin functions cannot take builtin functions as parameters (actually, this is not a limitation since it is trivial to embed a builtin function into a user-defined one). Here is a simple example of builtin function embedding in a user function. #+NAME: builtin-function-embedding #+BEGIN_SRC pugs :exports both :results output import math; let cosinus: R -> R, x -> cos(x); cout << "cosinus(2) = " << cosinus(2) << "\n"; #+END_SRC Running this example produces #+results: builtin-function-embedding Builtin functions are defined when modules are imported, see [[modules]]. ** Modules<<modules>> Modules are sets of *builtin functions* and *high-level types* that are not provided by default by ~pugs~. To load a module, one has to use an ~import~ instruction in the preamble of the script. For instance to load the ~math~ module which contains builtin mathematical functions, one writes in the preamble of the script #+BEGIN_SRC pugs :exports source import math; #+END_SRC #+BEGIN_warning A work in progress - At the time of writing this documentation, one should note that module inter-dependencies is still not implemented. - Also, (and especially with regard to the ~scheme~ module), module contents are likely to change and to be reorganized. - Finally it is almost sure that modules will be equipped with a /namespace/-like functionality to avoid conflicts. This actually should make the scripts cleaner since, the naming of functions would be more natural. #+END_warning One can access to the list of available modules inside the language. #+NAME: get-available-modules #+BEGIN_SRC pugs :exports both :results output cout << getAvailableModules() << "\n"; #+END_SRC The output lists all available modules #+RESULTS: get-available-modules Let us comment a bit this output. One notices that there are two kind of modules. Modules that are automatically imported (tagged with a ~*~) and the other ones. In this section we will not describe exhaustively the whole module contents but will give the basic information that should allow the user to find his way. To do so, it is important to examine carefully the content of the ~core~ module, since it contains some helper functions. *** The ~core~ module As already said, the ~core~ module cannot be imported manually since it is imported automatically. **** ~core~ provided types ***** ~ostream~ The ~core~ module provides an important type: ~ostream~. This type has already been encountered in this documentation. This is the type of some specific predefined objects: ~cout~, ~cerr~ and ~clog~. These objects are very similar to their ~C++~ counterparts, respectively ~std::cout~, ~std::cerr~ and ~std::clog~. Objects of type ~ostream~ can be used as the left operand for the binary operator ~<<~ for basic types and their associated tuples. | ~ostream <<~ allowed expression type | |------------------------------------| | ~B~ | | ~N~ | | ~Z~ | | ~R~ | | ~R^1~ | | ~R^2~ | | ~R^3~ | | ~R^1x1~ | | ~R^2x2~ | | ~R^3x3~ | | ~string~ | As one should have noticed, when the left operand of the binary operator ~<<~ is an ~ostream~, the result of the operation is also an ~ostream~, which allows to chain output. One can overload the ~ostream <<~ construction for high-level types. Other variables of the type ~ostream~ can be created (in order to write to files for instance) as one should see below. **** ~core~ provided functions ***** ~getAvailableModules: void -> string~ This function that is used in the preamble of this section is a function that returns a ~string~ that contains the list of modules that are available in the current version of pugs. ***** ~getModuleInfo: string -> string~ This is a *very* helpful function. It lists *all* the available *builtin functions* of a module and *all* the provided *types*. It is important to remark that the name of the module is given as a ~string~ and not as is. #+BEGIN_note This function provides a very simple documentation of the modules. Hopefully in future developments it would be really useful to obtain a complete documentation of each function to offer a precise online documentation. #+END_note In the following of the documentation we will run this function on all modules. #+NAME: get-module-info-core #+BEGIN_SRC pugs :exports both :results output cout << getModuleInfo("core") << "\n"; #+END_SRC For the ~core~ module, it gives #+RESULTS: get-module-info-core This output is quite rustic, but it contains the main information: the name of the function and its input and output sets. #+BEGIN_warning Observe that this function does not provide the list of operators that are defined in the module (eventually associated to the defined types). #+END_warning ***** ~getPugsBuildInfo: void -> string~ Sets the building information of the executable into a ~string~. #+NAME: get-pugs-build-info #+BEGIN_SRC pugs :exports both :results output cout << getPugsBuildInfo() << "\n"; #+END_SRC It gives for instance #+RESULTS: get-pugs-build-info ***** ~getPugsVersion: void -> string~ Sets the ~pugs~ version of the executable into a ~string~. #+NAME: get-pugs-version #+BEGIN_SRC pugs :exports both :results output cout << getPugsVersion() << "\n"; #+END_SRC The output contains also ~git~ information #+RESULTS: get-pugs-version ***** ~resetRandomSeed: void -> void~ At the start of ~pugs~, a *common* random seed is defined for all ~MPI~ processes. The seed is written automatically written in the preamble of the listing (in the documentation the ~--no-preamble~ option is used this is why it is not displayed in examples). The ~resetRandomSeed~ creates a *new shared* random and displays its value. #+NAME: reset-random-seed #+BEGIN_SRC pugs :exports both :results output resetRandomSeed(); #+END_SRC The output is #+RESULTS: reset-random-seed ***** ~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. #+NAME: set-random-seed #+BEGIN_SRC pugs :exports both :results output setRandomSeed(123456789); #+END_SRC Running this example gives #+RESULTS: set-random-seed ***** ~ofstream: string -> ostream~ This function is used to create an ~ostream~ that actually write to the file which name is given by the ~string~ argument. One should notice that the file is *created* at the function call. If a file with the same name existed, it is *erased*. #+NAME: ofstream-example #+BEGIN_SRC pugs :exports both :results output let fout:ostream, fout = ofstream("filename.txt"); fout << [1,2] << " is a vector of R^2\n"; #+END_SRC Running this example produces no output #+RESULTS: ofstream-example But a file is created (in the execution directory), with the name ~"filename.txt"~. Its content is #+NAME: cat-filename-txt #+BEGIN_SRC shell :exports results :results output cat filename.txt #+END_SRC #+RESULTS: cat-filename.txt Following ~C++~, the file is closed when the variable is destroyed. *** The ~dev_utils~ module This module contains utilities for ~pugs~ developers. **** ~dev_utils~ provided functions #+NAME: get-module-info-dev-utils #+BEGIN_SRC pugs :exports both :results output cout << getModuleInfo("dev_utils") << "\n"; #+END_SRC #+RESULTS: get-module-info-dev-utils ***** ~getAST: void -> string~ Gets the abstract syntax tree (AST) that will be executed into a ~string~. ***** ~getFunctionAST: function -> string~ Gets the abstract syntax tree associated to a user function into a ~string~. ***** ~saveASTDot: string -> void~ Saves the AST of the script into the file (whose name is given as argument) using the dot format. *** The ~math~ module The ~math~ module is a small utility module that provides a set of standard mathematical functions. #+NAME: get-module-info-math #+BEGIN_SRC pugs :exports both :results output cout << getModuleInfo("math") << "\n"; #+END_SRC For the ~math~ module, it gives #+RESULTS: get-module-info-math There is not much to say. One can see that this module does not provide new types of data. All the classical mathematical functions follow their ~C++~ counterparts. One can see that there are integer versions (type ~Z~) of ~abs~, ~min~ and ~max~ functions. The ~dot~ function family provides the dot product for vectors of $\mathbb{R}^1$, $\mathbb{R}^2$ and $\mathbb{R}^3$. #+NAME: dot-examples #+BEGIN_SRC pugs :exports both :results output import math; cout << "([1],[2]) = " << dot([1],[2]) << "\n"; cout << "([1,2],[3,4]) = " << dot([1,2],[3,4]) << "\n"; cout << "([1,2,3],[4,5,6]) = " << dot([1,2,3],[4,5,6]) << "\n"; #+END_SRC The output is #+RESULTS: dot-examples #+BEGIN_note Observe that the use of a proper rounding or truncation function is the right way to convert a real value to an integer one. Available 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 provides 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 ***** Boundary descriptor 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 ***** Zone descriptor functions ****** ~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"),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.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"), 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.gnu)' lt rgb "green" w l, '<(sed "" $PUGS_SOURCE_DIR/doc/diamond.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 variables ~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"), 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.gnu)' lt rgb "green" w l, '<(sed "" $PUGS_SOURCE_DIR/doc/median.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"), 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.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~ <<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 the vertices 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*} #+BEGIN_warning The connectivity must be the *same in memory*, this means that constructing two identical meshes with /equivalent/ connectivity is not allowed. #+END_warning 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"), 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"), m1); let m2: mesh, m2 = relax(m0, m1, 0.3); write_mesh(gnuplot_writer("relax_example_m2"), m2); #+END_SRC In this example, the relaxation parameter is set to $\theta=0.3$. 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.gnu)' lt rgb "green" w l, '<(sed "" $PUGS_SOURCE_DIR/doc/relax_example_m1.gnu)' lt rgb "blue" w l, '<(sed "" $PUGS_SOURCE_DIR/doc/relax_example_m2.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]]) for instance. *** 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"), name_output(fh, "f")); #+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.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"), 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.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"), 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.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>> 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 *** The ~linear_solver~ module This module provides the following functions #+NAME: get-module-info-linear-solver #+BEGIN_SRC pugs :exports both :results output cout << getModuleInfo("linear_solver") << "\n"; #+END_SRC #+RESULTS: get-module-info-linear-solver **** ~linear_solver~ provided functions The set of provided functions is used to define the global behavior of ~pugs~ linear system solver. ***** Utility functions An important function is ****** ~getLSAvailable: void -> string~ This shows the available options or libraries that are available. It depends on the compilation options of the code. #+NAME: get-ls-available-example #+BEGIN_SRC pugs :exports both :results output import linear_solver; cout << getLSAvailable() << "\n"; #+END_SRC #+results: get-ls-available-example ****** ~getLSOptions: void -> string~ This function show the current tuning of the global linear solver #+NAME: get-ls-options-example #+BEGIN_SRC pugs :exports both :results output import linear_solver; cout << getLSOptions() << "\n"; #+END_SRC #+results: get-ls-options-example ***** Tuning functions ****** ~setLSEpsilon: R -> void~ Sets the relative epsilon criterion for iterative methods. ****** ~setLSLibrary: string -> void~ Selects the library to use. ****** ~setLSMaxIter: N -> void~ Sets the maximum number of iterations. ****** ~setLSMethod: string -> void~ Sets the method to solve linear systems. ****** ~setLSPrecond: string -> void~ Sets the preconditioner type. ****** ~setLSVerbosity: B -> void~ Sets verbosity mode to ~true~ or ~false~ for linear solvers. *** The ~writer~ module This module provides functionalities to numerical results to files for post processing. It provides the following functions and types. #+NAME: get-module-info-writer #+BEGIN_SRC pugs :exports both :results output cout << getModuleInfo("writer") << "\n"; #+END_SRC #+RESULTS: get-module-info-writer [fn:pugs-def] ~pugs~: Parallel Unstructured Grid Solvers [fn:MPI-def] ~MPI~: Message Passing Interface [fn:DSL-def] ~DSL~: Domain Specific Language~