XLiFE++ website is now hosted here with up-to-date content. The current website, with potentially deprecated content, will be removed soon.
Example 9/18 - Fictitious domain method
A well-known method to avoid some complex meshes consists in using the ficticious domain method where the boundary condition on an inside obstacle is taken into account by a Lagrange multiplier. To illustrate it, we consider the following 2D Laplace problem in the domain with the unit square, the ball with center and radius such that :
Let us introduce the minimization problem:
where
It has a unique solution such that on and there exists such that
where the integrals on are to be
understood as the duality product on .
The key idea of ficticious domain method is to use two different meshes for
and
. As a consequence, the
computation of integrals on
involves shape functions that lie on two different meshes, inducing interpolation operations from one
mesh to the other one. XLiFE++ manages this interpolation operations in an hidden way for the
user:
#include "xlife++.h"
using namespace xlifepp;
Real R;
Real f(const Point& P, Parameters& pa = defaultParameters)
{if(norm(P)<=R) return 0.;
return 2*pi_*pi_*sin(pi_*P(1))*sin(pi_*P(2));}
Real g(const Point& P, Parameters& pa = defaultParameters)
{return sin(pi_*P(1))*sin(pi_*P(2));}
int main(int argc, char** argv)
{
init(argc, argv, _lang=en); // mandatory initialization of xlifepp
// create meshes (rectangle and circle)
R=0.2;
SquareGeo sq(_xmin =0, _xmax = 1, _ymin = 0, _ymax = 1,
_hsteps=0.05,
_domain_name = "Omega", _side_names="Sigma");
Disk circle(_center=Point(0.5, 0.5), _radius =R, _hsteps=0.05,
_domain_name = "Gamma");
Mesh meshS(sq, triangle);
Mesh meshG(circle, segment, 1, gmsh);
Domain omega = meshS.domain("Omega"),
sigmat = meshS.domain("Sigma"),
gammat = meshG.domain("Gamma");
// create spaces and unknowns
Space Vt(_domain=omega, _interpolation=P1, _name="Vt");
Unknown ut(Vt, "ut"); TestFunction vt(ut, "vt");
Space W(_domain=gammat, _interpolation=P0, _name="W");
Unknown p(W, "p"); TestFunction q(p, "q");
// create bilinear forms and Terms
BilinearForm at = intg(omega, grad(ut)|grad(vt))-intg(gammat, p*vt)
- intg(gammat, ut*q);
BoundaryConditions bcst = (ut|sigmat =0);
TermMatrix At(at, bcst, "At");
TermVector Ft(intg(omega, f*vt)-intg(gammat, g*q), "Ft");
// solve linear system and save the solution to a vtu file
TermVector Ut = directSolve(At, Ft);
saveToFile("Ut", Ut, vtu);
return 0;
}
To extract the solution on the real domain ,
a L2 projection is used:
Disk disk(_center=Point(0.5,0.5),_radius =R, _hsteps=0.05,
_domain_name = "Obstacle", _side_names="Gamma");
Mesh mesh(rectangle-disk,_triangle);
Domain omega = mesh.domain("Omega");
Space V(omega,P1,"V"); Unknown u(V,"u"); TestFunction v(u,"v");
//compute L2 projection of ut on omega
TermMatrix M(intg(omega,u*v)), M2(intg(omega,ut*v));
TermVector PUt = directSolve(M,M2*Ut);
saveToFile("PUt",PUt,_vtu);
Note that the computation of the matrix M2 also requires a computation of an integrals involving two meshes!