Project 9
A finite elements solver for the Poisson problem in 2D
with mixed boundary conditions

1. Solve the boundary value problem shown in the following figure.
   

   This models the steady-state temperature in a square where the boundary temperature is maintained at $0$ all around except for the lower-half of the left edge where heat flux is specified instead of temperature.

   1. Modify square-spec.c by inserting an extra point (point #4 in the figure) and segment #4. Now your domain will be defined through 5 points and 5 segments.
   2. Set all boundary conditions to FEM_BC_DIRICHLET except for segment #4 which would be FEM_BC_NEUMANN.
   3. Your square-spec.c needs to define the functions $f$, $g$, $h$, $\eta$. We have $f\equiv0$, $g\equiv0$, $h\equiv2$, and $\eta\equiv1$. We have no $u_\text{exact}$, so set it to NULL in struct problem_spec.
   4. Plot the solution in python. Here is what mine looks like:

      ./demo-square 4 0.001
      domain is a square
      nodes = 835, edges = 2402, elems = 1568
      system stiffness matrix is 835x835 (=697225) has 5101 nonzero entries
      pyplot output written to file demo-square.py
      

      
2. Let $\Omega$ be the annulus with inner and outer radii $a=0.325$ and $b=2a$ as before. Solve the PDE $\Delta u = 0$ on $\Omega$ subject to $u=u_0$ on the inner boundary and $u=0$ on the outer boundary. [Take $u_0 = 0.5$ to get a good-looking picture.] The exact solution, $u_\text{exact} = u_0 \frac{\ln (b/r)}{\ln(b/a)}$, represents the steady-state temperature in the annulus where the inner boundary is heated to temperature $u_0$ while the outer boundary in maintained at zero temperature.
   1. Your program needs to define the functions $f$, $g$, $h$, $\eta$, and $u_\text{exact}$. We have $f≡0$, $h≡0$, $\eta≡1$.
   2. Defining $g$ is a bit tricky, since it needs to distinguish between the inner and outer boundaries. If $n$, the number of edges on the inner and outer boundaries, is sufficiently large—in our case $n \ge 5$ will do— then a circle of radius $(a+b)/2$ separates the inner and outer boundaries. [Draw a picture to see that.] So the function $g$, which receives the coordinates $(x,y)$ of a vertex, can compare $\sqrt{x^2+y^2}$ with $(a+b)/2$ and return $0$ or $u_0$, as the case may be.
   3. The natural logarithm function, $\ln$, is named log in the Standard C Library. Remember to put #include<math.h> in your annulus-spec.c.
   4. Plot the solution in python. Here is what mine looks like, but the "energy norm = 0" looks suspicious. I will investigate and update later.

      ./demo-annulus 5 0.005 24
      domain is an annulus (really a 24-gon) of radii 0.325 and 0.65
      nodes = 200, edges = 535, elems = 335
      system stiffness matrix is 200x200 (=40000) has 860 nonzero entries
      errors: L^infty = 0.00619778, L^2 = 0.00268236, energy norm = 0
      pyplot output written to file demo-annulus.py
      

      

      Update: Okay, I see why we get zero for the error in the energy norm. That's because my code incorrectly applies the formula (25.17) from page 347. That formula was derived for boundary value problems with zero Dirichlet data. Here we have nonzero data on the inner boundary and therefore we need a generalization of (25.17) to nonzero data. I don't know how to do that at the moment, so for now ignore that zero value, or better yet, just don't print it.