2D Laplace problem with transmission condition#

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

We turn to the Laplace problem with transmission condition:

\begin{array}{r} {\begin{cases} - Δ u^{-} = f & in Ω^{-} \\ - Δ u^{+} = f & in Ω^{+} \\ u_{| Σ^{-}} = 1 and u_{| Σ^{+}} = 1 \\ u_{| Γ}^{-} = u_{| Γ}^{+} and \partial_{x} u_{| Γ}^{-} = \partial_{x} u_{| Γ}^{+} \end{cases} \end{array}

Its variational formulation in

V = {(v^{-}, v^{+}) \in H^{1} (Ω^{-}) \times H^{1} (Ω^{+}), v_{| Σ^{-}}^{-} = 0, v_{| Σ^{-}}^{+} = 0, v_{| Γ}^{-} = v_{| Γ}^{+}}

is

\begin{array}{r} | \begin{array}{c} find (u^{-}, u^{+}) \in H^{1} (Ω^{-}) \times H^{1} (Ω^{+}), u_{| Σ^{-}}^{-} = 1, u_{| Σ^{-}}^{+} = 1, u_{| Γ}^{-} = u_{| Γ}^{+} such that \\ \int_{Ω^{-}} \nabla u^{-} . \nabla v^{-} + \int_{Ω^{+}} \nabla u^{+} . \nabla v^{+} = \int_{Ω^{-}} f v^{-} + \int_{Ω^{+}} f v^{+} \forall v \in V . \end{array} \end{array}

Note that derivatives matching is taken into account in a weak sense. The implementation in XLiFE++, using $P^{2}$ Lagrange finite element, looks like:\

#include "xlife++.h"
using namespace xlifepp;

Real f(const Point& P, Parameters& pa = defaultParameters)
{
  return -8.;
}

int main(int argc, char** argv)
{
  init(argc, argv, _lang=en); // mandatory initialization of xlifepp

  // mesh domain
  Strings sn(4);
  sn[1] = "x=1/2-"; sn[3] = "x=0";
  Rectangle r1(_xmin=0, _xmax=0.5, _ymin=0, _ymax=1, _nnodes=Numbers(20, 40), _domain_name="Omega-", _side_names=sn);
  Mesh mesh2d(r1, triangle, 1, _structured);
  sn[1] = "x=1"; sn[3] = "x=1/2+";
  Rectangle r2(_xmin=0.5, _xmax=1, _ymin=0, _ymax=1, _nnodes=Numbers(20, 40), _domain_name="Omega+", _side_names=sn);
  Mesh mesh2d_p(r2, _shape=triangle, _generator=structured);
  mesh2d.merge(mesh2d_p);
  Domain omegaM=mesh2d.domain("Omega-"), omegaP=mesh2d.domain("Omega+");
  Domain sigmaM=mesh2d.domain("x=0"), sigmaP=mesh2d.domain("x=1");
  Domain gamma=mesh2d.domain("x=1/2- or x=1/2+");

  // create P2 interpolation
  Space VM(_domain=omegaM, _interpolation=P2, _name="VM");
  Unknown uM(VM, _name="u-");
  TestFunction vM(uM, _name="v-");
  Space VP(_domain=omegaP, _interpolation=P2, _name="VP");
  Unknown uP(VP, _name="u+");
  TestFunction vP(uP, _name="v+");

  // create bilinear form and linear form
  BilinearForm auv=intg(omegaM, grad(uM)|grad(vM))+intg(omegaP, grad(uP)|grad(vP));
  LinearForm fv=intg(omegaM, f*vM)+intg(omegaP, f*vP);
  EssentialConditions ecs= (uM|sigmaM = 1) & (uP|sigmaP = 1) & ((uM|gamma) - (uP|gamma) = 0);
  TermMatrix A(auv, ecs, _name="A");
  TermVector B(fv, _name="B");

  // solve linear system AX=B using LU factorization
  TermVector U=directSolve(A, B);
  saveToFile("U_LT.vtu", U);
  return 0;
}

Here, a tool merging mesh is used to create a two domains mesh. gmsh should be used also. The picture below shows that the solution is continuous across the boundary $Γ$ .

../_images/ex_2d_transmission.png — Fig. 16 Solution of the Laplace 2D problem with transmission condition on the unit square $[0, 1]^{2}$ #