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

From ONELAB
Jump to: navigation, search
(The considered problem)
(TEST)
Line 51: Line 51:
 
{{tutorial|getdp/LaplacianDirichlet/GMSH_GETDP/LaplacianDirichlet.pro|height=22em}}
 
{{tutorial|getdp/LaplacianDirichlet/GMSH_GETDP/LaplacianDirichlet.pro|height=22em}}
  
== TEST ==
+
== Result ==
 
 
[[File:LaplacianNeumannsolutionU.png|thumb|right|alt=Sample magnetostatics demo.|Scalar magnetic potential and magnetic induction.]]
 

Revision as of 11:37, 23 September 2011

The considered problem

In this problem, we consider a Laplace equation, as in that example, except that the boundary condition is here of Dirichlet type. To model this in GetDP, we will introduce a "Constraint". We assume that the reader has already study the previous example.

As in the example with a Neumann boundary condition, we considered the unit square $\Omega = [0,1]\times[0,1]$ and the unknown function $u$ solves the following problem \begin{equation} \begin{cases}\label{eq:problemU} \Delta u + u = f & \text{in } \Omega,\\ \displaystyle{u = 0} & \text{on }\Gamma, \end{cases} \end{equation} where $\Gamma = \partial\Omega$ is the boundary of $\Omega$ and $\displaystyle{\Delta = \frac{\partial^2}{\partial x_1^2} + \frac{\partial^2}{\partial x_2^2} }$ is the Laplacian operator.

Domain Omega.

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). $$ As a consequence, the unique solution $u$ of the problem \eqref{eq:problemU} is clearly given by $$ \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 problem \eqref{eq:problemU}: \begin{equation}\label{eq:WeakFormulation} \left\{\begin{array}{l} \text{Find } u\in H^1_0(\Omega) \text{ such that, }\\ \displaystyle{\forall v\in H^1_0(\Omega), \qquad \int_{\Omega} \nabla u\cdot\nabla v \;{\rm d}\Omega + \int_{\Omega}uv \;{\rm d}\Omega = 0}, \end{array}\right. \end{equation} where the functions $v$ are the test functions and $H^1_0(\Omega) = \left\{ v \in H^1(\Omega) \text{ such that } v|_{\Gamma} = 0\right\}$.

Outline of the program

The proposed solution is composed by 3 different files:

param.geo
Auxiliary file that contains the index number associated with the geometry.
LaplacianNeumann.geo
GMSH file, used to build both the square $\Omega$ and its boundary $\Gamma$. Except for the building of $\Gamma$, this file is the same as in the previous example.
LaplacianNeumann.pro
GetDP file, contains the weak formulation \eqref{eq:WeakFormulation} of the problem \eqref{eq:problemU}. The main difference between this example and the Neumann case is the addition of the Constraint part, which is located in this GetDP file.

param.geo: the auxiliary file

// File "param.geo"

//Numbers that caracterise the interior of the square (Omega) and its boundary (Gama):
Omega = 1000;
Gama = 1001;

// 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/LaplacianDirichlet/GMSH_GETDP/param.geo'


LaplacianDirichlet.geo: creation of the geometry with GMSH

// File "LaplacianDirichlet.geo".

// This file is exactly the same as "LaplacianNeumann.geo" EXCEPT
// for the creation of the Physical Entity "Gama" (the boundary), at the end of the file

// We include the file containing the numbering of the geometry.
Include "param.geo";

//Caracteristic length of the finite elements (reffinement is also possible after the mesh is built):
lc = 0.05;

// 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};

//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".
// Both "Omega" and "Gama" are imported from the file "param.geo".
Physical Surface(Omega) = {SurfaceOmega};
Physical Line(Gama) = {L1,L2,L3,L4}; // the Physical Line Gama is composed of the four lines of te boundary (not of the Line Loop !)
// Do not forget to let a blank line at the end, this could make GMSH crash...

Direct link to file `getdp/LaplacianDirichlet/GMSH_GETDP/LaplacianDirichlet.geo'


LaplacianDirichlet.pro: weak formulation

// File LaplacianDirichlet.pro
// This file is similar to LaplacianDirichlet.pro
// Only changes are detailed.

// we import the indices of Gama and Omega
Include "param.geo";

// Group
//======
// We have 2 Groups : Omega and Gama (boundary)
Group{
Omega = Region[{Omega}];
Gama = Region[{Gama}];
}
//Again note that "Gamma" is a syntaxic word of GetDP
// That is why we write here "Gama" and not "Gamma".
// (not a mistake from the author :-))

// Function
//=========
Function{
// Pi is a special value in GetDp (=3.1415...)
Coef = Pi;
// Definition of the function f(x,y)
f[] = -(1+2*Coef*Coef)*Sin[Coef*X[]]*Sin[Coef*Y[]];
// Here f exist on Omega AND Gama. As a consequence, we have to let "f" lives everywhere,
// and we let the bracket [] empty.
}



//Jacobian
//========
Jacobian {
  { Name JVol ;
    Case {
      { Region All ; Jacobian Vol ; }
    }
  }
  { Name JSur ;
    Case {
      { Region All ; Jacobian Sur ; }
    }
  }
  { Name JLin ;
    Case {
      { Region All ; Jacobian Lin ; }
    }
  }
}

//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 ; }
        }
      }
    }
  }
}

//Contraint
//=============
// This part is new, compare to LaplacianDirichlet.pro.
Constraint{
	{Name DirichletBC; Type Assign;
	Case{
        {Region Gama;  Value 0; }
		}
	}
}
/*
The Dirichlet Boundary Condition is here called "DirichletBC.
"Assign": type means that it is a Dirichlet type condition
"Region": domain where the condition is applied.
"Value": value assigned (here 0).
*/



//FunctionSpace
//=============
FunctionSpace{
  { Name Vh; Type Form0;
    BasisFunction{
      {Name wn; NameOfCoef vn; Function BF_Node;
    Support Omega; Entity NodesOf[All];}
    }
    // Contrary to Neumann case, we had here the constraint :
	Constraint{
	{NameOfCoef vn; EntityType NodesOf;
	 NameOfConstraint DirichletBC;}
		}
  }
}

/* On the constraint: 
- "NameOfCoef" : a constraint is applied to specific coefficients, here "vn" (in our example, there are the only coefficients
because of functions are P1)
- "Entitytype" : The constraint is applied on the nodes (e.g : not on the edge)
- "NameOfConstraint" : Last (but not least), we have to specify which Constraint is called.

With that definition, every function u (even the test-functions) of V_h satisfy the Dirichlet Boundary Condition:
		u = 0 on Gama.
In other words, this spaces approximate the Sobolev space H^1_0(Omega).
*/


//(Weak) Formulation
//==================

Formulation{
  {Name LaplacianDirichlet; Type FemEquation;
    Quantity{
      {Name u; Type Local; NameOfSpace Vh;}
    }
    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 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  = 0

The Dirichlet Boundary Condition is contained in the function space V_h !
*/


// Resolution (of the problem)
//============================

Resolution{
// nothing change ... (except the name :-)
  {Name LaplacianDirichlet;
    System{
      {Name Syst; NameOfFormulation LaplacianDirichlet;}
    }
    Operation{
      Generate[Syst]; Solve[Syst]; SaveSolution[Syst];
    }
  }
}


//Post Processing
//===============

PostProcessing{
  {Name LaplacianDirichlet; NameOfFormulation LaplacianDirichlet;
    Quantity{
      {Name u; Value {Local{[{u}];In Omega;Jacobian JVol;}}}
      // We also ask the computation of the absolute value of "u" and the function "f".
      //These new variable are called respectively "AbsU" and "f"
      {Name AbsU; Value {Local{[Norm[{u}]];In Omega;Jacobian JVol;}}}
      {Name f; Value {Local{[f[]];In Omega;Jacobian JVol;}}}
    }
  }
}

//Post Operation
//==============

PostOperation{
  {Name Map_u; NameOfPostProcessing LaplacianDirichlet;
    Operation{
      Print[u, OnElementsOf Omega, File "u_Dirichlet.pos"];
      // Plot of Abs(u) and f
      Print[AbsU, OnElementsOf Omega, File "u_Abs_Dirichlet.pos"];
      Print[f, OnElementsOf Omega, File "f.pos"];
    }
  }
}

Direct link to file `getdp/LaplacianDirichlet/GMSH_GETDP/LaplacianDirichlet.pro'


Result