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 Ω = CB(O,R)with C =]0,1[×]0,1[ the unit square, B(O,R) the ball with center C and radius R such that B C:

{ Δu = f in Ω u = 0  on C = Σ u = g  on  B = Γ

pict

Let us introduce the minimization problem:

min (g)1 2C||2 Cf~

where

(g) = { H1(C), |Σ = 0 and |Γ = g} and f~ = { f in Ω 0 in B .

It has a unique solution ũ (g) such that ũ = u on Ω and there exists p H1 2 (Γ) such that

{ Cũ. Γp =Cf~ H01(C) Γũq =Γgq q H1 2 (Γ),

where the integrals on Γ are to be understood as the duality product on H1 2 (Γ) × H1 2 (Γ).

The key idea of ficticious domain method is to use two different meshes for C 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:

// create a mesh for Omega 
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!


pict

Figure 5.1: Solution of the Laplace 2D problem with the ficticious domain method