In [1]:

```
# This magic makes plots appear in the browser
%matplotlib notebook
import matplotlib.pyplot as plt
```

$$ -\nabla^2 u + u = f \text{ on } \Omega = [0, 1] \times [0, 1], \\
\nabla u \cdot \vec{n} = 0 \text{ on } \Gamma, $$

where $\vec{n}$ is the outward-pointing unit normal to $\Gamma$ and $f \in V$ is a known source function. Note that the Laplacian and mass terms have opposite signs, which makes this equation much more benign than the symmetric, indefinite operator $(\nabla^2 + I)$.

Now, we write the problem in its variational form by multiplying by a test function $v \in V$ and integrating by parts over the domain $\Omega$:

$$ \int_\Omega \nabla u\cdot\nabla v + uv\ \mathrm{d}x = \int_\Omega
vf\ \mathrm{d}x + \underbrace{\int_\Gamma v \nabla u \cdot \vec{n}\, \mathrm{d} s}_{=0}. $$

$$ \int_\Omega \nabla u\cdot\nabla v + uv\ \mathrm{d}x = \int_\Omega
vf\ \mathrm{d}x, \text{ for all } v \in V. $$

We can choose the source function $f$, so lets take it to be:

$$ f(x, y) = (1.0 + 8.0\pi^2)\cos(2\pi x)\cos(2\pi y). $$

This conveniently yields the analytical solution:

$$ u(x, y) = \cos(2\pi x)\cos(2\pi y). $$

`UnitSquareMesh`

.

`firedrake`

. To save on typing, we will import all of the public API into the current namespace

In [2]:

```
from firedrake import *
```

`UnitSquareMesh`

will give us an appropriate mesh of the domain $\Omega$, but what are the arguments to the constructor? The user-facing API is usually well-documented, and we can access this information via Python using the builtin `help`

command:

In [3]:

```
help(UnitSquareMesh)
```

In [4]:

```
# In the notebook, we can also use a special "?" magic to pop up a help box
?UnitSquareMesh
```

In [5]:

```
mesh = UnitSquareMesh(10, 10)
```

In [6]:

```
# NBVAL_IGNORE_OUTPUT
# We test the output of the notebooks with out continuous integration system to
# make sure they run.
# Unfortunately we can't compare the plots from run to run, so the above line
# indicates that the output of this cell should be ignored for the purposes
# of testing.
fig, axes = plt.subplots()
triplot(mesh, axes=axes)
axes.legend();
```

In [7]:

```
V = FunctionSpace(mesh, "Lagrange", 1)
```

*symbolic* objects that correspond to test and trial spaces, and the linear and bilinear forms in our problem.

In [8]:

```
u = TrialFunction(V)
v = TestFunction(V)
```

In [9]:

```
x, y = SpatialCoordinate(mesh)
f = (1 + 8*pi*pi)*cos(2*pi*x)*cos(2*pi*y)
```

In [10]:

```
a = (dot(grad(v), grad(u)) + v * u) * dx
L = f * v * dx
```

In [11]:

```
uh = Function(V)
```

In [12]:

```
solve(a == L, uh, solver_parameters={'ksp_type': 'cg', 'pc_type': 'none'})
```

In [13]:

```
# NBVAL_IGNORE_OUTPUT
fig, axes = plt.subplots()
collection = tripcolor(uh, axes=axes, cmap='coolwarm')
fig.colorbar(collection);
```

In [14]:

```
u_exact = cos(2*pi*x)*cos(2*pi*y)
```

In [15]:

```
# NBVAL_IGNORE_OUTPUT
difference = assemble(interpolate(u_exact, V) - uh)
fig, axes = plt.subplots()
collection = tripcolor(difference, axes=axes, cmap='coolwarm')
fig.colorbar(collection);
```

Solve the same problem, only this time, use a piecewise quadratic approximation space.

- Hint: check the help for
`FunctionSpace`

to see how to specify the degree.

Solve the same problem, but using quadrilateral, rather than triangular, cells.

- Hint 1: check the help for
`UnitSquareMesh`

to see how to make a quadrilateral mesh - Hint 2: To specify a piecewise continuous space on quadrilaterals, use the family name
`"Q"`

.

In [ ]:

```
```

For solutions with sufficient smoothness (like the choice we have here), this method with a piecewise linear approximation space should converge in the $L_2$ error with rate $\mathcal{O}(h^{-2})$, where $h$ is the typical mesh spacing. Confirm this for the example in question by computing the $L_2$ error in the solution for a sequence of finer and finer meshes.

- Hint 1: You can compute errors using errornorm
- Hint 2: If the error is $\mathcal{O}(h^{-2})$ then $\log_2 (e_H/e_h) \approx 2$.
The Python module
`math`

contains an implementation of`log`

.

Instead of (or as well as!) refining the mesh, we can increase the degree of the approximating polynomial space.

To help, here's the complete problem setup.

In [16]:

```
mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, "Lagrange", 1)
u = TrialFunction(V)
v = TestFunction(V)
x, y = SpatialCoordinate(mesh)
f = (1 + 8*pi*pi)*cos(2*pi*x)*cos(2*pi*y)
a = (dot(grad(v), grad(u)) + v * u) * dx
L = f * v * dx
uh = Function(V)
solve(a == L, uh, solver_parameters={'ksp_type': 'cg',
'pc_type': 'none'})
u_exact = cos(2*pi*x)*cos(2*pi*y)
```

In [ ]:

```
```

Let us recall again the statement of our problem. We seek $u \in V$ satisfying

$$ -\nabla^2 u + u = f \text{ on } \Omega = [0, 1] \times [0, 1], \\ \nabla u \cdot \vec{n} = 0 \text{ on } \Gamma. $$However, we will not always want to constrain the normal derivative on the boundary to be zero. Let's consider the modification:

$$ -\nabla^2 u + u = f \text{ on } \Omega = [0, 1] \times [0, 1], \\ \nabla u \cdot \vec{n} = g = 1 \text{ on } \Gamma_1, \\ \nabla u \cdot \vec{n} = 0 \text{ on } \Gamma \setminus \Gamma_1,$$where $\Gamma_1$ is the boundary $x = 0$.

As previously, we introduce a test function $v \in V$, multiply and integrate. After integrating by parts, we obtain

$$ \int_\Omega \nabla u\cdot\nabla v + uv\ \mathrm{d}x = \int_\Omega vf\ \mathrm{d}x + \underbrace{\int_{\Gamma\setminus \Gamma_1} v \nabla u \cdot \vec{n}\, \mathrm{d} s}_{=0} + \int_{\Gamma_1} v \nabla u \cdot \vec{n}\,\mathrm{d} s.$$The first surface integral over $\Gamma \setminus \Gamma_1$ vanishes, since $\nabla u \cdot \vec{n} = 0$. The second, however, does not. Substituting the boundary value we obtain

$$ \int_\Omega \nabla u\cdot\nabla v + uv\ \mathrm{d}x = \int_\Omega vf\ \mathrm{d}x + \int_{\Gamma_1} v g\,\mathrm{d} s.$$We see that the linear form for the right hand side has gained an integral over part of the boundary (some of its *exterior* facets). We've already seen how to express integrals over the cells of the domain, using `dx`

. Unsurprisingly, we can also write integrals over the boundary, for which we use the `ds`

measure.

There's one final wrinkle, just as `dx`

integrates over all the cells in the mesh, `ds`

integrates over all the exterior facets of the mesh. Firedrake uses *mesh markers* to select which parts of the mesh the integral measure should integrate over. These are integers (or tuples thereof) that have some meaning for the mesh. Most external mesh generators will have some way of marking meshes in this way. For the utility meshes, the docstring indicates which markers correspond to which parts of the mesh. So let's look:

In [17]:

```
?UnitSquareMesh
```

`1`

. Our full problem specification is now

In [18]:

```
mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, "Lagrange", 1)
u = TrialFunction(V)
v = TestFunction(V)
x, y = SpatialCoordinate(mesh)
f = (1 + 8*pi*pi)*cos(2*pi*x)*cos(2*pi*y)
a = (dot(grad(v), grad(u)) + v * u) * dx
```

`ds`

would integrate over all exterior facets, we select the facets corresponding to $x = 0$ by specifying the appropriate mesh marker.

In [19]:

```
g = Constant(1)
L = f*v*dx + g*v*ds(1)
```

Now to solve the problem

In [20]:

```
uh = Function(V)
solve(a == L, uh, solver_parameters={'ksp_type': 'cg', 'pc_type': 'none'})
```

Now plot your solution.

In [ ]:

```
```