# Tutorial/Laplace equation with Neumann boundary condition

## Contents

## The considered problem

We propose here to solve a first very simple academic problem with the help of GMSH and GetDP. We considered the unit square $\Omega = [0,1]\times[0,1]$ and we seek $u$, solution of the following problem \begin{equation} \begin{cases}\label{eq:problemU} -\Delta u + u = f & \text{in } \Omega,\\ \displaystyle{\frac{\partial u}{\partial \mathbf{n}} = 0} & \text{on }\partial\Omega, \end{cases} \end{equation} where $\mathbf{n}$ is the unit outward normal of $\Omega$ and $\displaystyle{\Delta = \frac{\partial^2}{\partial x_1^2} + \frac{\partial^2}{\partial x_2^2} }$ is the Laplace operator.

To simplify, we suppose that the function $f$ is defined by
$$
\forall (x,y)\in [0,1]^2,\qquad f(x,y) = (1+2\pi^2)\cos(\pi x)\cos(\pi y).
$$
With this choice of $f$, one can easily show that the unique solution of the problem \eqref{eq:problemU} reads as
$$
\forall (x,y)\in[0,1]^2, \qquad u(x,y) = \cos(\pi x)\cos(\pi y).
$$
In order to solve problem \eqref{eq:problemU} with the finite element method, we write the weak formulation of problem \eqref{eq:problemU}:
\begin{equation}\label{eq:WeakFormulation}
\left\{\begin{array}{l}
\text{Find } u\in H^1(\Omega) \text{ such that, }\\
\displaystyle{\forall v\in H^1(\Omega), \qquad \int_{\Omega} \nabla u\cdot\nabla v \;{\rm d}\Omega + \int_{\Omega}uv \;{\rm d}\Omega - \int_{\Omega}fv\;{\rm d}\Omega = 0},
\end{array}\right.
\end{equation}
where $H^1(\Omega)$ is the classical Sobolev space and the functions $v$ are the *test functions*.

## Outline of the program

We give here a (very) detailled solution which is composed by 3 different files:

- param.geo
- This (auxiliary) file contains the index number associated with the geometry. This ensures that GMSH and GetDP use the same numbering of the domains.
- LaplacianNeumann.geo
- GMSH file, used to build the domain (the square $\Omega$). The extension ".geo" is mainly used to design a GMSH file
- LaplacianNeumann.pro
- GetDP file, contains the weak formulation \eqref{eq:WeakFormulation} of the problem \eqref{eq:problemU}. The extension ".pro" is associated with GetDP files.

## param.geo: the auxiliary file

// File "param.geo" //Numbers that caracterise the interior of the square (Omega) and its boundary (Gama): Omega = 1000; // Three remarks on these numbers : // - They are arbitrary choosen. // - They are placed in a separated file to be readable by both GMSH and GetDP. // - "Gamma" is a special word used by GMSH/GetDP, that is why the boundary is named "Gama", with one "m"... // Do not forget to let a blank line at the end, this could make GMSH crash...

Direct link to file `LaplacianNeumann/GMSH_GETDP/param.geo'

## LaplacianNeumann.geo: creation of the geometry with GMSH

// File "LaplacianNeumann.geo". // We include the file containing the numbering of the geometry. // This is usefull at the end of this file, and used to "synchronise" GMSH and GetDP Include "param.geo"; //Caracteristic length of the finite elements (reffinement is also possible after the mesh is built): lc = 0.05; // This parameter could be placed for instance in "param.geo", to separate more easyly the geometry // and the discretization parameters. // The parameters of the border of the domain : x_max = 1; x_min = 0; y_max= 1; y_min = 0; //Creation of the 4 angle points of the domain Omega (=square) p1 = newp; Point(p1) = {x_min,y_min,0,lc}; p2 = newp; Point(p2) = {x_min,y_max,0,lc}; p3 = newp; Point(p3) = {x_max,y_max,0,lc}; p4 = newp; Point(p4) = {x_max,y_min,0,lc}; // Remarks: // -"newp" is a GMSH function that give the first available number for describing a point. // For any other entity, like Line, Surface, etc. We recommand the use of "newreg" (see below). // - By default, GMSH create a 3D domain. The z-coordinate must always be precised. //The four edges of the square L1 = newreg; Line(L1) = {p1,p2}; L2 = newreg; Line(L2) = {p2,p3}; L3 = newreg; Line(L3) = {p3,p4}; L4 = newreg; Line(L4) = {p4,p1}; // Line Loop (= boundary of the square) Bound = newreg; Line Loop(Bound) = {L1,L2,L3,L4}; //Surface of the square SurfaceOmega = newreg; Plane Surface(SurfaceOmega) = {Bound}; // To conclude, we define the physical entities, that is "what GetDP could see/use". // "Omega" is a number imported from the file "param.geo". Physical Surface(Omega) = {SurfaceOmega}; // Do not forget to let a blank line at the end, this could make GMSH crash...

Direct link to file `LaplacianNeumann/GMSH_GETDP/LaplacianNeumann.geo'

## LaplacianNeumann.pro: weak formulation

// File LaplacianNeumann.pro // As in the .geo file, we include the file containing the numbering of the geometry. Include "param.geo"; // Group //====== //We now build the "Groups", that is, the geometrical entities, the different domains of computation. //Here we only need the interior of the open domain Omega Group{ Omega = Region[{Omega}]; } // Function //========= Function{ // Pi is a special value in GetDp (=3.1415...) Coef = Pi; // Definition of the source function f(x,y) f[] = (1+2*Coef*Coef)*Cos[Coef*X[]]*Cos[Coef*Y[]]; } /* Remark: - the argument (as "x" and "y") are not writen between the bracket []. Indeed, between the bracket is writen the domain of definition of the function. - The argument "x" and "y" are here designed by the GetDP inner functions X[] and Y[], which give respectively the x-coordinate and the y-coordinate of a considered point. - To define a function globaly (i.e: not only on a subdomain), we write: f[] = ... - In our example, we could also have written f[Omega] = (1+2*Coef*Coef)*Cos[Coef*X[]]*Cos[Coef*Y[]]; */ //Jacobian //======== Jacobian { { Name JVol ; Case { { Region All ; Jacobian Vol ; } } } { Name JSur ; Case { { Region All ; Jacobian Sur ; } } } { Name JLin ; Case { { Region All ; Jacobian Lin ; } } } } /* Remark: roughly speacking, we make use of...: - Jacobian "Vol" as far as the integration domain is of same dimension than the problem (e.g: calculation in a 3D domain, in a 3D problem, a 2D domain (surface) in a 2 problem, a 1D domain (line) in a 1D problem) - Jacobian "Sur" when the domain of integration has one dimension less than the global problem (e.g: calculation on a surface (2D) in a 3D problem, calculation on a line (1D) in a 2D problem). - Jacobian "Lin" when the domain of integration has 2 dimensions less than the problem. That is, for example, a calculation on a line (1D) in a 3D problem. - Here, we just have define some Jacobian, you will see later that we just use the Jacobian "JVol" */ //Integration (parameters) //======================= Integration { { Name I1 ; Case { { Type Gauss ; Case { { GeoElement Point ; NumberOfPoints 1 ; } { GeoElement Line ; NumberOfPoints 4 ; } { GeoElement Triangle ; NumberOfPoints 6 ; } { GeoElement Quadrangle ; NumberOfPoints 7 ; } { GeoElement Tetrahedron ; NumberOfPoints 15 ; } { GeoElement Hexahedron ; NumberOfPoints 34 ; } } } } } } // There is no Constrain because of the Neumann boundary condition // We go directly to the FuntionSpace //FunctionSpace //============= FunctionSpace{ { Name Vh; Type Form0; BasisFunction{ {Name wn; NameOfCoef vn; Function BF_Node; Support Omega; Entity NodesOf[All];} } } } /* Explanation: The space of approximation is called "Vh". Its type is "Form0", that means "scalar". The basis functions are called "wn", and the associated coefficients "vn". In other words, a function "v" of "Vh" can be written as v(x,y) = sum_{n} vn.wn(x,y) The functions "wn" are nodale functions P1 ("BF_Node"), supported on the domain Omega ("Support Omega"). */ //(Weak) Formulation //================== Formulation{ {Name LaplacianNeumann; Type FemEquation; // We decided to call the formulation "LaplacianNeumann". // Its type is "FemEquation" ("Finite Element Method") Quantity{ {Name u; Type Local; NameOfSpace Vh;} } // Here, we introduce a quantity "u", which belongs to the function space Vh, defined above. Equation{ Galerkin{ [Dof{Grad u}, {Grad u}]; In Omega; Jacobian JVol; Integration I1;} Galerkin{ [ Dof{u}, {u}]; In Omega; Jacobian JVol; Integration I1;} Galerkin{ [-f[], {u}]; In Omega; Jacobian JVol; Integration I1;} } } } /* The variationnal formulation is written between the acolades {}, after "Equation". Let us first give some vocabulary: - "Galerkin" : GetDP syntaxic word. This could be translated mathematicaly by "integration" (see below) - Dof{u}: "Degree Of Freedom". This is used to specify that the quantity is the unknown. If "Dof" is not written, then "u" is seen as a test function and not as the unknown. (Be carreful, the unknown and the test functions has the same name in GetDP ! The "Dof" is here to distinguish the unknown to the test-functions) Now, let us detail the program written between the acolades {}, : Galerkin{ [Dof{Grad u}, {Grad u}]; In Omega; Jacobian JVol; Integration Int;} This could be translated mathematicaly by : Integration on "Omega" of Grad(u).Grad(v) where "u" is the unknown and "v" a test function. Note the use of the Jacobian JVol(2D problem and integration on a 2D domain). Moreover the number of integrations point is given in "I1" (see above "Integration"). The total variationnal formulation can be read as : search "Dof{u}" in V_h such that, for every "{u}" in V_h, \int_{\Omega} Grad( Dof{u}) . Grad( {u}) d\Omega + \int_{\Omega} Dof{u}.{u} d\Omega - ... ... - \int_{\Omega} f.{u} d\Omega = 0 (this is exactly the weak formulation of our problem !) Remarks: - Between two "Galerkin", a positive sign "+" is implicitly written - The sum of all integrals "Galerkin" is equal to 0 (do not forget the "minus" sign of the right hand side) - Why do we use a volumic jacobian even in a 2D problem? The problem is a 2 dimensional problem and the integral is defined on a 2D domain (Omega). If the integral was writen on, e.g, the boundary, then the Jacobian "Jsur" would have been used. */ // Resolution (of the problem) //============================ Resolution{ {Name LaplacianNeumann; // We chose the name LaplacianNeumann for the resolution // Remark: in GetDP, every entity has a name. The same same can be used for different entities but of different kind, of course // Here we chose the same name as the formulation, but this is just a choice, no obligation ! System{ {Name Syst; NameOfFormulation LaplacianNeumann;} } // A system is linked to a weak formulation // Here, we only have one weak formulation, which is "LaplacianNeumann" Operation{ Generate[Syst]; Solve[Syst]; SaveSolution[Syst]; } /* When we launch GetDP, the program will ask the user to chose a Resolution. When calling the Resolution "LaplacianNeumann", GetDP will... - Generate the system - Solve the system - Save the solution Note that GetDP respects the order of the operation ! */ } } //Post Processing //=============== PostProcessing{ {Name LaplacianNeumann; NameOfFormulation LaplacianNeumann; Quantity{ {Name u; Value {Local{[{u}];In Omega;Jacobian JVol;}}} } } } /* The name of the PostProcessing is LaplacianNeumann. It call the weak formulation LaplacianNeumann, then take the solution "u" and compute it on all the domain Omega, by the operation: Value {Local{[{u}];In Carre; Jacobian JVol;}}} This means that "u" is the interpolate of the solution "Dof{u}" computed in the weak formulation "LaplacianNeumann" Remark: Again, we chose the same name, but this is just a choice, this has no influence on GetDP resolution. */ //Post Operation //============== PostOperation{ {Name Map_u; NameOfPostProcessing LaplacianNeumann; Operation{ Print[u, OnElementsOf Omega, File "u_Neumann.pos"]; } } } /* The only PostOperation we write is to display "u" introduced in the PostProcessing. This PostOperation is called "Map_u". When launching GetDP, it will ask the user to chose : - the Resolution - the PostOperation */

Direct link to file `LaplacianNeumann/GMSH_GETDP/LaplacianNeumann.pro'

## How to launch

All the files (.geo and .pro) must be located in the same directory.

- Meshing the domain

- Graphical way: launch GMSH and open "LaplacianNeumann.geo". Then chose the "mesh" menu and click on "2D". Finally, do not forget to save the mesh by clicking on the "Save" button
- With a terminal: go to the directory and then type:

- gmsh -LaplacianNeumann.geo -algo iso -2

- Note that here we chose the isobare algorithm ("iso"). One can also chose the following algorithms:
- del2d
- frontal

After the mesh is built, a file LaplacianNeumann.msh should have been created in the directory.

- Solving the problem with GetDP

In a Terminal, type (in the right directory)

- getdp LaplacianNeumann.pro -solve -pos

GetDP will then propose the Resolution ("-solve") and the PostOperation ("-pos"). This should build a file called "u_Neumann.pos".

Remarks :

- the choice can be pre-selected by typing:
- getdp LaplacianNeumann.pro -solve#1 -pos#1

- There exists other options, like the choice of the linear solver, \ldots

- Showing the result

Finally, open the file "u_Neumann.pos" with GMSH, with the graphical way ("open \ldots') of with a terminal (by typing "gmsh u_Neumann.pos" in a terminal in the right directory).

## Result

Here is the result obtained with GMSH/GetDP. On the left is an example of the mesh of $\Omega$ and on the right, the modulus of the solution $u$.

## All files

A zipfile containing all files (.geo, .pro, .mesh and .pos) can be downloaded here.