Tutorial/Laplace equation with Neumann boundary condition
Contents
The considered problem
We propose here to solve a first very simple academic example with GMSH and GetDP. We considered the unit square $\Omega$ with boundary $\Gamma$ and unit outward normal $\mathbf{n}$. 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 }\Gamma, \end{cases} \end{equation} where $\Gamma =\partial\Omega$ is the boundary of $\Omega$ and 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). $$ One can easily show that the unique solution of the problem \eqref{eq:problemU} is $$ \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 elements method, we write the weak formulation of the probleme \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. Our solution is composed by 3 different files. - LaplacianNeumann.geo : GMSH file, used to build the domain (the square). 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. - data.geo : this (auxiliary) file contains the index number associated with the geometry. It is used to ensure that GMSH and GetDP have the same numbering of the domains.
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 `getdp/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 `getdp/LaplacianNeumann/GMSH_GETDP/LaplacianNeumann.geo'
