How to impose interface condition ∂u⁺/∂n = k ∂u⁻/∂n in FreeFem++

General Discussion
1

Hello everyone,

I am solving a PDE on a domain Ω split into two sub-domains Ω⁺ and Ω⁻ separated by an interface Γ.

On Γ I need to impose the interface (flux) condition
[
\frac{\partial u^{+}}{\partial n} = k,\frac{\partial u^{-}}{\partial n}
\quad \text{on } \Gamma,
]
where k is a given constant (or function defined on Γ).

I can mesh the two sub-domains separately (or use a single mesh with a label for the interface), and I can write the weak form of the PDE in each sub-domain.
However, I am not sure how to correctly enforce this interface condition in FreeFEM++.

Could someone provide guidance or a small example showing:

  • how to write the weak formulation that incorporates this interface condition,

  • or how to implement it using Lagrange multipliers or jump terms in FreeFEM++.

A minimal code example would be greatly appreciated.

Thank you in advance for your help!

2

This is a classical problem, if you do heat equation, this two domain two thermal diffusivity k1,k2

the the boundary condition say the thermal flux is continuous , and the u is continuous and in this case the

$$ k_1 \Delta u_1 = f_1 $$ on domain \\Omaga_1 + BC on external boundary

$$ k_2 \Delta u_2 = f_2 $$ on domain \\Omaga_2 + BC on external boundary

and k_1 \\frac{\\partial u_1}{\\partial n} = k_2 \\frac{\\partial u_2}{\\partial n} on \\Gamma

if you call u the solution on \\Omega = \\Omega1 \\cup \\Omega2

k the discontinuous diffusivity of \\Omega

then the Varf is

find u (u given on \\partial \\Omega = g$)

$$ \forall v \H^1_0(\Omega), int_\Omega k \nabla u . \nabla v = int_ f v$$

3

Thank you for your answer.
I’m working on a problem where
[
\Delta u^+ = 0 \quad \text{in } \Omega_1
]
and
[
\Delta u^- = 0 \quad \text{in } \Omega_2
]
with (\Omega_1) being a square domain and (\Omega_2) an ellipse inside (\Omega_1).

On the interface (\Gamma = \partial\Omega_2), I impose
[
u^+ = u^- \quad \text{on } \Gamma , \qquad
k_1 \frac{\partial u^+}{\partial n} = k_2 \frac{\partial u^-}{\partial n}.
]

So the solution is continuous across the interface, and the flux is continuous as well.

I tried to reproduce the discontinuity of (k) as suggested in the FreeFem++ documentation, using:

solve Poisson(u,v,solver=LU) =
    int2d(tout)(dx(u)*dx(v) + dy(u)*dy(v))
  - int1d(tout)(k*(nx*dx(u) + ny*dy(u))*v)
  + int2d(tout)(f*v)
  + on(1, u=0)
  + on(2, u=g);

and I defined (k) as follows:

// We create the variable k = 1 or -1
fespace V(tout,P0);
V pl = region;
int inside = pl((x1+x0)/2, (y1+y0)/2);
V chi = (inside == region);
V k = (1 - chi) * (-1) + 1 * chi;

I think this should be sufficient to achieve a good modeling of the problem.

4

kamjd.edp (571 Bytes)

I have build a small example

first a mesh Th of ]-2,2[^2 square with a ellipse inside label 2 on ellipse and label 1 on square.

solve Poisson(u,v,solver=LU) =
    int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v))
  - int1d(Th,2)(k*v)
  + int2d(Th)(f*v)
  + on(1, u=0);

The variationnal formulation is 

$$ \forall v \in H^1(\Omega_1), \int_\Omega  \nabla u.\nabla v = \int_\Gamma (k_1-k_2) v$$
and $k = k_1-k_2$;