## Poisson Equation in 2D

In this example we solve the Poisson equation in two space dimensions.

For a domain with boundary , we write the boundary value problem (BVP):

Here, denotes the part of the boundary where we prescribe Dirichlet boundary conditions, and denotes the part of the boundary where we prescribe Neumann boundary conditions. denotes the unit normal of pointing outside .

To obtain the weak form we define the functional spaces and . Then we multiply the strong form by an arbitrary function and integrate over :

Integration by parts of the non-conforming term gives

Recalling that on and that on , the weak form of the BVP is the following.

Find :

To obtain the finite element discretization we then introduce a triangulation (mesh) of the domain and we define a finite dimensional subspace consisting of globally continuous functions that are piecewise polynomial on each element of .

By letting and , the finite element method then reads:

Find such that:

In what follow, we will let be the unit square, be the top boundary, and be the union of the left, bottom, and right boundaries.

The coefficient , , are chosen such that the analytical solution is .

### 1. Imports

We import the following Python packages:

`dolfin`

is the python interface to FEniCS.`matplotlib`

is a plotting library that produces figure similar to the Matlab ones.`math`

is the python built-in library of mathematical functions.

```
from __future__ import print_function, absolute_import, division
# Import FEniCS
import dolfin as dl
import math
# Enable plotting inside the notebook
import matplotlib.pyplot as plt
%matplotlib inline
from hippylib import nb
import logging
logging.getLogger('FFC').setLevel(logging.WARNING)
logging.getLogger('UFL').setLevel(logging.WARNING)
dl.set_log_active(False)
```

### 2. Define the mesh and the finite element space

We define a triangulation (mesh) of the unit square with `n`

elements in each direction. The mesh size is .

We also define the finite element space as the space of globally continuos functions that are piecewise polinomial (of degree ) on the elements of the mesh.

```
n = 64
d = 1
mesh = dl.UnitSquareMesh(n, n)
Vh = dl.FunctionSpace(mesh, "Lagrange", d)
print("Number of dofs", Vh.dim())
nb.plot(mesh, mytitle="Finite Element Mesh", show_axis='on')
plt.show()
```

```
Number of dofs 4225
```

### 3. Define the Dirichlet boundary condition

We define the Dirichlet boundary condition on .

```
def boundary_d(x, on_boundary):
return (x[1] < dl.DOLFIN_EPS or x[0] < dl.DOLFIN_EPS or x[0] > 1.0 - dl.DOLFIN_EPS) and on_boundary
u_d = dl.Expression("sin(DOLFIN_PI*x[0])", degree = d+2)
bcs = [dl.DirichletBC(Vh, u_d, boundary_d)]
```

### 4. Define the variational problem

We write the variational problem . Here, the bilinear form and the linear form are defined as

- .

denotes the trial function and denotes the test function. The coefficients and are also given.

```
uh = dl.TrialFunction(Vh)
vh = dl.TestFunction(Vh)
f = dl.Constant(0.)
g = dl.Expression("DOLFIN_PI*exp(DOLFIN_PI*x[1])*sin(DOLFIN_PI*x[0])", degree=d+2)
a = dl.inner(dl.grad(uh), dl.grad(vh))*dl.dx
L = f*vh*dl.dx + g*vh*dl.ds
```

### 5. Assemble and solve the finite element discrete problem

We now assemble the finite element stiffness matrix and the right hand side vector . Dirichlet boundary conditions are applied at the end of the finite element assembly procedure and before solving the resulting linear system of equations.

```
A, b = dl.assemble_system(a, L, bcs)
uh = dl.Function(Vh)
dl.solve(A, uh.vector(), b)
nb.plot(uh, mytitle="Finite Element Solution")
plt.show()
```

### 6. Compute error norms

We then compute the and the energy norm of the difference between the exact solution and the finite element approximation.

```
u_ex = dl.Expression("exp(DOLFIN_PI*x[1])*sin(DOLFIN_PI*x[0])", degree = d+2, domain=mesh)
grad_u_ex = dl.Expression( ("DOLFIN_PI*exp(DOLFIN_PI*x[1])*cos(DOLFIN_PI*x[0])",
"DOLFIN_PI*exp(DOLFIN_PI*x[1])*sin(DOLFIN_PI*x[0])"), degree = d+2, domain=mesh )
norm_u_ex = math.sqrt(dl.assemble(u_ex**2*dl.dx))
norm_grad_ex = math.sqrt(dl.assemble(dl.inner(grad_u_ex, grad_u_ex)*dl.dx))
err_L2 = math.sqrt(dl.assemble((uh - u_ex)**2*dl.dx))
err_grad = math.sqrt(dl.assemble(dl.inner(dl.grad(uh) - grad_u_ex, dl.grad(uh) - grad_u_ex)*dl.dx))
print ("|| u_ex - u_h ||_L2 / || u_ex ||_L2 = ", err_L2/norm_u_ex)
print ("|| grad(u_ex - u_h)||_L2 / = || grad(u_ex)||_L2 ", err_grad/norm_grad_ex)
```

```
|| u_ex - u_h ||_L2 / || u_ex ||_L2 = 0.0004267539774457824
|| grad(u_ex - u_h)||_L2 / = || grad(u_ex)||_L2 0.02453704048509311
```

## Hands on: Advection-Diffusion-Reaction PDEs

For a bounded domain with boundary , consider the boundary value problem

Here, denotes the part of the boundary where we prescribe Dirichlet boundary conditions, and denotes the part of the boundary where we prescribe Neumann boundary conditions. denotes the unit normal of pointing outside .

### Question 1

Derive the weak formulation of the above boundary value problem.

To derive the weak form, we first define the function spaces and subspaces of associated with the inhomogeneuos and homogeneuos form of the Dirichlet boundary condition:

Then we multiply the strong form by an arbitrary function and integrate over :

Integration by parts of the non-conforming term gives

Recalling that on and that on , the weak form of the BVP is the following.

Find :

### Question 2

Compute the finite element solution of the above problem using FEniCS with , , , , , . Choose , , and as in the example above.

```
n = 64
d = 1
mesh = dl.UnitSquareMesh(n, n)
Vh = dl.FunctionSpace(mesh, "Lagrange", d)
print("Number of dofs", Vh.dim())
def boundary_d(x, on_boundary):
return (x[1] < dl.DOLFIN_EPS or x[0] < dl.DOLFIN_EPS or x[0] > 1.0 - dl.DOLFIN_EPS) and on_boundary
u_d = dl.Constant(0.)
bcs = [dl.DirichletBC(Vh, u_d, boundary_d)]
uh = dl.TrialFunction(Vh)
vh = dl.TestFunction(Vh)
k = dl.Expression("exp(-x[0]*x[0] -x[1]*x[1])", degree = d+2)
w = dl.Constant( (1.,0.))
c = dl.Constant(1.)
f = dl.Constant(1.)
g = dl.Constant(0.)
a = k*dl.inner(dl.grad(uh), dl.grad(vh))*dl.dx + dl.dot(w, dl.grad(uh))*vh*dl.dx + c*uh*vh*dl.dx
L = f*vh*dl.dx + g*vh*dl.ds
A, b = dl.assemble_system(a, L, bcs)
uh = dl.Function(Vh)
dl.solve(A, uh.vector(), b)
nb.plot(uh, mytitle="Finite Element Solution")
plt.show()
```

```
Number of dofs 4225
```

Copyright © 2018, The University of Texas at Austin & University of California, Merced. All Rights reserved. See file COPYRIGHT for details.

This file is part of the hIPPYlib library. For more information and source code availability see https://hippylib.github.io.

hIPPYlib is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License (as published by the Free Software Foundation) version 2.0 dated June 1991.