2D/3D Laplace problem with Neumann condition#

Keywords: PDE:Laplace Method:FEM FE:Lagrange Solver:Direct

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

\frac{\partial u}{\partial n} = 0 on \partial Ω .

The variational formulation we deal with is

\begin{array}{r} | \begin{array}{c} find u \in V = {v \in L^{2} (Ω), \nabla v \in (L^{2} (Ω))^{2}} such that \\ \int_{Ω} \nabla u . \nabla v + a \int_{Ω} u v = \int_{Ω} f v \forall v \in V . \end{array} \end{array}

The following main program corresponds to solving this problem on unity square $Ω =] 0, 1 [\times] 0, 1 [$ with $f (x) = \cos π x \cos π 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, _shape=triangle, _generator=structured);
  Domain omega = mesh2d.domain("Omega");

  // build space and unknown
  Space Vk(_domain=omega, _interpolation=P1, _name="Vk");
  Unknown u(Vk, _name="u");
  TestFunction v(u, _name="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, _name="a(u, v)");
  TermVector B(fv, _name="f(v)");

  // LLt factorize and solve
  TermMatrix LD;
  ldltFactorize(A, LD);
  TermVector U = factSolve(LD, B);

  saveToFile("U_LN", U, _format=vtu);
  return 0;
}

../_images/ex_2d.png — Fig. 25 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), _generator=gmsh, _name="P1 mesh");
...

../_images/ex_3d.png — Fig. 26 Solution of the Laplace 3D problem with Neumann condition on the unit cube $[0, 1]^{3}$ #