XLiFE++ website is now hosted here with up-to-date content. The current website, with potentially deprecated content, will be removed soon.

Example 2/18 - Laplace 2D/3D problem with Neumann condition

We investigate here problems involving laplacian operator in a 2D bounded domain, say Ω : Δu + au = f in Ω(a = k2 for Helmholtz equation)

and various essential conditions (Dirichlet, transmission, quasi periodic, average condition).

First, let us consider the case of the homogeneous Neumann condition on Ω, the boundary of Ω:

u n = 0 on Ω.

The variational formulation we deal with is

| find u V = {v L2(Ω),v (L2(Ω))2}  such that Ωu.v + aΩuv =Ωfvv V.

The following main program corresponds to solving this problem on unity square Ω =]0,1[×]0,1[ with f(x) = cos πxcos πy using P1 Lagrange element (20x20 elements):

 
#include "xlife++.h" 
using namespace xlifepp; 
 
Real cosx2(const Point& P, Parameters& pa = defaultParameters) 
{ 
  Real x=P(1), y=P(2); 
  return cos(pi_ * x) * cos(pi_ * y); 
} 
 
int main(int argc, char** argv) 
{ 
  init(argc, argv, _lang=en); 
 
// mesh square 
  SquareGeo sq(_origin=Point(0., 0.), _length=1, _nnodes=21); 
  Mesh mesh2d(sq, triangle, 1, structured); 
  Domain omega = mesh2d.domain("Omega"); 
 
// build space and unknown 
  Space Vk(_domain=omega, _interpolation=P1, _name="Vk"); 
  Unknown u(Vk, "u"); 
  TestFunction v(u, "v"); 
 
// define variational formulation 
  BilinearForm auv = intg(omega, grad(u) | grad(v)) + intg(omega, u * v); 
  LinearForm fv=intg(omega, cosx2 * v); 
 
// compute matrix and right hand side 
  TermMatrix A(auv, "a(u, v)"); 
  TermVector B(fv, "f(v)"); 
 
// LLt factorize and solve 
  TermMatrix LD; 
  ldltFactorize(A, LD); 
  TermVector U = factSolve(LD, B); 
 
  saveToFile("U_LN", U, vtu); 
  return 0; 
}

pict

Figure 2.1: Solution of the Laplace 2D problem with Neumann condition on the square [0,1]2

Solving this problem with P2 Lagrange interpolation should be the same except the line defining the space:

Space Vh(_domain=omega, _interpolation=P2, _name="Vh", _optimizeNumbering);

Solving this problem in a 3D domain should be the same except the line defining the mesh and the right hand side function. For instance, on the unity cube, the mesh construction command using Gmsh tool is:

Real f(const Point& P, Parameters& pa = defaultParameters) 
{ 
   Real x=P(1), y=P(2), z=P(3); 
   return cos(pi*x) * cos(pi*y) * cos(pi*z); 
} 
... 
Mesh mesh(Cube(_origin=Point(0.,0.,0.), _length=1, _nnodes=10), tetrahedron, 1, _gmsh,"P1 mesh"); 
...

pict

Figure 2.2: Solution of the Laplace 3D problem with Neumann condition on the unit cube [0,1]3