FreeFEM Documentation on GitHub

stars - forks

# Transmission problem

Consider an elastic plate whose displacement change vertically, which is made up of three plates of different materials, welded on each other.

Let $$\Omega_i$$, $$i=1,2,3$$ be the domain occupied by $$i$$-th material with tension $$\mu_i$$ (see Soap film).

The computational domain $$\Omega$$ is the interior of $$\overline{\Omega_1}\cup \overline{\Omega_2}\cup \overline{\Omega_3}$$. The vertical displacement $$u(x,y)$$ is obtained from:

(59) $\begin{split}\begin{array}{rcll} -\mu_i\Delta u &=& f & \textrm{ in }\Omega_i\\ \mu_i\p_n u|_{\Gamma_{i}} &=& -\mu_j\p_n u|_{\Gamma_{j}} & \textrm{ on }\overline{\Omega_{i}}\cap\overline{\Omega_{j}} \textrm{ if }1\le i< j\le 3 \end{array}\end{split}$

where $$\p_n u|_{\Gamma_{i}}$$ denotes the value of the normal derivative $$\p_n u$$ on the boundary $$\Gamma_i$$ of the domain $$\Omega_i$$.

By introducing the characteristic function $$\chi_i$$ of $$\Omega_i$$, that is:

$\chi_i(x)=1\ \textrm{ if }x\in\Omega_i;\ \chi_i(x)=0\ \textrm{ if }x\not\in\Omega_i$

we can easily rewrite (59) to the weak form. Here we assume that $$u=0$$ on $$\Gamma=\p\Omega$$.

Transmission problem: For a given function $$f$$, find $$u$$ such that:

$a(u,v) = \ell(f,v) \textrm{ for all }v\in H^1_0(\Omega)$
$\begin{split}\begin{array}{rcl} a(u,v) &=& \int_{\Omega}\mu \nabla u\cdot \nabla v\nonumber\\ \ell(f,v) &=& \int_{\Omega}fv\nonumber \end{array}\end{split}$

where $$\mu=\mu_1\chi_1+\mu_2\chi_2+\mu_3\chi_3$$. Here we notice that $$\mu$$ become the discontinuous function.

This example explains the definition and manipulation of region, i.e. sub-domains of the whole domain. Consider this L-shaped domain with 3 diagonals as internal boundaries, defining 4 sub-domains:

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 // Mesh border a(t=0, 1){x=t; y=0;}; border b(t=0, 0.5){x=1; y=t;}; border c(t=0, 0.5){x=1-t; y=0.5;}; border d(t=0.5, 1){x=0.5; y=t;}; border e(t=0.5, 1){x=1-t; y=1;}; border f(t=0, 1){x=0; y=1-t;}; border i1(t=0, 0.5){x=t; y=1-t;}; border i2(t=0, 0.5){x=t; y=t;}; border i3(t=0, 0.5){x=1-t; y=t;}; mesh th = buildmesh(a(6) + b(4) + c(4) +d(4) + e(4) + f(6) + i1(6) + i2(6) + i3(6)); // Fespace fespace Ph(th, P0); //constant discontinuous functions / element Ph reg=region; //defined the P0 function associated to region number // Plot plot(reg, fill=true, wait=true, value=true); 

region is a keyword of FreeFEM which is in fact a variable depending of the current position (is not a function today, use Ph reg=region; to set a function). This variable value returned is the number of the sub-domain of the current position. This number is defined by buildmesh which scans while building the mesh all its connected component.

So to get the number of a region containing a particular point one does:

 1 2 3 4 5 6 7 8 // Characteristic function int nupper = reg(0.4, 0.9); //get the region number of point (0.4,0.9) int nlower = reg(0.9, 0.1); //get the region number of point (0.4,0.1) cout << "nlower = " << nlower << ", nupper = " << nupper<< endl; Ph nu = 1 + 5*(region==nlower) + 10*(region==nupper); // Plot plot(nu, fill=true,wait=true); 

This is particularly useful to define discontinuous functions such as might occur when one part of the domain is copper and the other one is iron, for example.

We this in mind we proceed to solve a Laplace equation with discontinuous coefficients ($$\nu$$ is 1, 6 and 11 below).

  1 2 3 4 5 6 7 8 9 10 11 12 13 // Problem solve lap (u, v) = int2d(th)( nu*(dx(u)*dx(v) + dy(u)*dy(v)) ) + int2d(th)( - 1*v ) + on(a, b, c, d, e, f, u=0) ; // Plot plot(u); 
Table of content