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 :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 :
The variational formulation we deal with is
The following main program corresponds to solving this problem on unity square
with
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;
}
Solving this problem with P2 Lagrange interpolation should be the same except the line defining the
space:
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");
...
{
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");
...