NOTE: You can use the code that was presented in lectures.
Let us consider a problem of de-noising a color image, by solving the following PDE: $$ \begin{align*} \frac{\partial}{\partial t} u _i &= {\rm div} (\nabla u _i) + \lambda (u^{(0)}-u _i) \\ u _i &= u^{(0)} _i \\ u _i| _{\partial\Omega} &= u^{(0)} _i| _{\partial\Omega} \\ \end{align*} $$ where $u_i$, $i=1,2,3$, are the intensity of the red, green, and blue components of the image; $u^{(0)}_i$ represent the original image. The parameter $\lambda$ should be inversely proportional to the level of the noise, and can be taken, for a start, to be 1/100 or 1/200. The time of integration should be taken so that we can approximate a steady-state solution.
Take a noisy image from here or generate it yourself. (If you generate yourself then make the noise level to be 64 out of 256.)
(10pt) Write a code that solves the system
Let us consider a more advanced approach: $$ \begin{align*} \frac{\partial}{\partial t} u _i &= {\rm div} (k \nabla u _i) + \lambda (u^{(0)}-u _i) \\ u _i &= u^{(0)} _i \\ u _i| _{\partial\Omega} &= u^{(0)} _i| _{\partial\Omega} \\ \end{align*} $$ where now $k=\big(1 + \max_i |\nabla u_i|^2\big)^{-1/2}$.
Consider a domain $\Omega$ with the boundary consisting of lines between the following 12 points (given in ${\rm mm}$ in this very order):
(250, -2 S - H), (100 - H, -2 S - H), (100 - H, -S), (150 - H, -S), (150 - H, -H), (0, -H), (0, 0), (150, 0), (150, -S - H), (100, -S - H), (100, -2 S), (250, -2 S)for some values of $S$ and $H$. The domain is as show here:
Next, we assume a piece of elastic material (similar to a beam) of this shape. It is fixed at the leftmost end and the force $P$ is applied to its rightmost end (the illustration below is for different values of $H$ and $S$):
To solve this problem, we consider the linear elasticity equations. The elasticity equations describe the displacement $(u,v)$ of the material in the $x$ and $y$ coordiate, respectively. We split the boudary of $\Omega$ into three parts: $\Gamma_1 = \{(x,y)\in\partial\Omega : x = 0\}$ (leftmost end), $\Gamma_2 = \{(x,y)\in\partial\Omega : 0 < x < 250\}$ (middle part), and $\Gamma_3 = \{(x,y)\in\partial\Omega : x = 250\}$ (rightmost part). It is then easy to write the boundary conditions on $\Gamma_1$: $$ u=v=0\qquad\text{on }\Gamma_1 $$
To write down the equations and the rest of the boundary conditions, we need the elasticity tensors which are denoted as $$ \begin{align*} C_{11} &= \begin{pmatrix} \frac{E (1-\nu)}{(1-2\nu)(1+\nu)} & 0 \\ 0 & \frac{E}{2(1+\nu)} \end{pmatrix}, &\qquad C_{12} &= \begin{pmatrix} 0 & \frac{E \nu}{(1-2\nu)(1+\nu)} \\ \frac{E}{2(1+\nu)} & 0 \end{pmatrix}, \\ C_{21} &= C_{12}^{\rm T}, &\qquad C_{22} &= \begin{pmatrix} \frac{E}{2(1+\nu)} & 0 \\ 0 & \frac{E (1-\nu)}{(1-2\nu)(1+\nu)} \end{pmatrix} \end{align*} $$ Here $E$ is called the modulus of elasticity and $\nu$ is called the Poisson ratio.
The equations are then written as $$ \begin{align*} -{\rm div} (C_{11} \nabla u + C_{12} \nabla v) &= 0 \qquad\text{in }\Omega \\ -{\rm div} (C_{21} \nabla u + C_{22} \nabla v) &= 0 \qquad\text{in }\Omega \\ u = v &= 0 \qquad\text{on }\Gamma_1 \\ C_{11} u_n + C_{12} v_n = C_{21} u_n + C_{22} v_n &= 0 \qquad\text{on }\Gamma_2 \\ C_{11} u_n + C_{12} v_n &= 0 \qquad\text{on }\Gamma_3 \\ C_{21} u_n + C_{22} v_n &= \frac{P}{H W} \qquad\text{on }\Gamma_3, \end{align*} $$ where $u_n$ and $v_n$ are normal derivatives of the solution. Note that in the last equation $C_{21} \nabla u + C_{22} \nabla v$ is the stress on the boundary and we set it equal to the stress $\frac{P}{H W}$ of the external force $P$. (The stress has the same dimension as the pressure.)
Let us now convert this to the variational formulation. We denote the test functions by $p=p(x,y)$ and $q = q(x,y)$. We require that the test functions satisfy the Dirichlet boundary conditions (i.e., $p(0,y)=q(0,y)=0$) and then the weak formulation becomes $$ \int_{\Omega} \big[ (C_{11} \nabla u + C_{12} \nabla v)\cdot \nabla p + (C_{21} \nabla u + C_{22} \nabla v)\cdot \nabla q \big] = \frac{P}{H W} \int_{\Gamma_3} q \qquad\forall p, q $$
(15pt) Write a code that generates a mesh for this domain (you may, but do not need to, use mesh generation packages)
(35pt) Write a code that solves this problem by the finite element method and outputs the $y$-component of the deflection of this piece of material, to be precise, $v(250,2 S)$
(10pt) For three sets of parameters, $S=20$, $W=3$, $H\in\{7,8,9\}$ (all in ${\rm mm}$) compute the displacement $v(250,2 S)$, $E = 3.1\cdot 10^9 {\rm Pa}$, $\nu=0.35$, $P = m\cdot g$, $m=200 {\rm grams}$, $g=-9810 {\rm mm}/{\rm s}^2$.
(Optional) Plot the three solutions found. For that you need to change the coordinates of each of the mesh nodes from $(x,y)$ by $(x,y) + (u(x,y),v(x,y))$ and plot the displaced mesh
In [2]:
from IPython.core.display import HTML
def css_styling():
styles = open("./styles/alex.css", "r").read()
return HTML(styles)
css_styling()
Out[2]:
In [ ]: