Difference between revisions of "Tutorial/Laplace equation with Neumann boundary condition"

From ONELAB
Jump to: navigation, search
(param.geo: the auxiliary file)
(LaplacianNeumann.geo: creation of the geometry with GMSH)
Line 39: Line 39:
  
 
== LaplacianNeumann.geo: creation of the geometry with GMSH ==
 
== LaplacianNeumann.geo: creation of the geometry with GMSH ==
 +
 +
{{tutorial|getdp/LaplacianNeumann/GMSH_GETDP/LaplacianNeumann.geo|height=22em}}
  
 
== LaplacianNeumann.pro: weak formulation ==
 
== LaplacianNeumann.pro: weak formulation ==

Revision as of 10:07, 5 September 2011

The considered problem

File:Magnetostatics.png
Scalar magnetic potential and magnetic induction.

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'


LaplacianNeumann.pro: weak formulation