Fast PDE/IE course, Skoltech, Spring 2015

Pre-Problem Set 2:

A strongly suggested problem to solve

Problem 1 (not graded)

Take a Poisson equation $$ \begin{align*} -\Delta u &= 1\qquad\text{on $\Omega$} \\ u &= 0\qquad\text{on $\Gamma$} \end{align*} $$ where $\Omega = B_1(0,0)\setminus B_{0.3}(0.7,0)$, $B_r(x,y)$ is the ball with radius $r$ and center at $(x,y)$.

Use finite elements.

Your goal is to compute $\int_{\Omega} u$

Many mesh generators will fail to construct a mesh in that region, so consider two regions, $\Omega_1 = B_1(0,0)\setminus B_{0.3}(0.7-\epsilon,0)$ and $\Omega_2 = B_1(0,0)\setminus B_{0.3}(0.7+\epsilon,0)$ and show that the computed solutions (or integrals of the solutions) are close to each other for small values of $\epsilon$.

Step 1. Build a mesh

You can do it either by using MeshPy or by downloading the prepared files.

MeshPy short guide and functions that read the prepared files can be found here.

Step 2. Assemble the stiffness matrix

• Write a function that assembes $3\times 3$ matrix $M$ for a given triangle $$M = \frac{|T|}{2} G G^T, \qquad\text{where}\qquad G = \left( \begin{array}{ccc} 1 & 1 & 1 \\ x_1 & x_2 & x_3 \\ y_1 & y_2 & y_3 \end{array} \right)^{-1} \left( \begin{array}{ccc} 0 & 0 \\ 1 & 0 \\ 0 & 1 \end{array} \right),$$ (see lecture 4 for more details).

• Run a loop for all triangles and assemble the stiffness matrix

• Correct rows and columns that correspond to Dirichlet boundary points

Note: use $\verb|scipy.sparse|$ to work with sparse matrices. Note that in the $\verb|scipy.sparse|$ the only sparse format that allows to change elements is lil.

Step 3. Assemble the right-hand side

The RHS for the $i$-th node is $$f_i = \int_\Omega f \eta_i d\Omega .$$ To assemble it run a loop for all triangles and use quadrature rule from the lecture 4.

Note: do not forget to put zeros in the rows that correspond to Dirichlet boundary points.

Step 4. Verify your code

Check your code using the fact that for $\Omega = B_1(0,0)$ the solution of $$ \begin{align*} -\Delta u &= 1\qquad\text{on $\Omega$} \\ u &= 0\qquad\text{on $\Gamma$} \end{align*} $$ is $$ u = \frac{1-x^2 - y^2}{4} $$ Note: you can plot the solution using the following code

%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
xi = np.linspace(-1, 1, 200) 
yi = np.linspace(-1, 1, 200) 
z = mlab.griddata(points[:, 0], points[:, 1], sol, xi, yi)
plt.contourf(xi, yi, z)
plt.colorbar()

Step 5. Calculate fluid flux

• Calculate fluid fluxes $\int_{\Omega} u\, d \Omega$ in $B_1(0,0)\setminus B_{0.3}(0.7,0)$ and in $B_1(0,0)\setminus B_{0.3}(0.5,0)$. Which one is bigger? Do you have intuition behind this fact?

Note: $\int_{\Omega} u\, d \Omega$ is not just a sum over all elemets of the solution vector as $u(x,y) \approx \sum с_i \eta_i(x,y)$.