# Intro to FEniCS - Part 3

In the Intro to `FEniCS`

- Part 2 episode we completed the first step in developing a `FEniCS`

model. We characterised the domain. Then, we selected a governing PDE. We then provided suitable boundary conditions. Finally, we expressed the PDE in weak form. We chosen a generalised Helmholtz equation as our PDE. In this episode we will start using the resulting weak form to develop a model.

# Project Files

All the files used for this project are available at the repositories below:

# The Model

In this episode we will build the model of a rectangular room. We already built a series of `Elmer`

models on this topic. This is why this system is well suited to explore a new solver: we have some understanding on it. This will help us figuring out whether the model is correct. In the Intro to `FEniCS`

- Part 2 episode we developed a generalised formalism. We now restrict it to a rectangular room. If $l_{x_0}$, $l_{x_1}$ and $l_{x_2}$ are the dimension of the room along $x_0$, $x_1$ and $x_2$ respectively we have:

$$ l_{x_0} \doteq 4 \space \text{m} \\ l_{x_1} \doteq 5 \space \text{m} \\ l_{x_2} \doteq 3 \space \text{m} $$ $$ \Rightarrow \Omega \doteq \left( 0, l_{x_0} \right) \times \left( 0, l_{x_1} \right) \times \left( 0, l_{x_2} \right) $$

To make things easy, as this is our first `FEniCS`

model, we will do some semplifications.

First, we remove all damping. So:

$$ \eta = 0 \space \frac{\text{rad}}{\text{s}} \text{m}^2 \\ \alpha = 0 \space \frac{\text{rad}}{\text{s}} $$

$$ \Rightarrow \gamma = 1, \space \beta = k^2; \space k = \frac{2 \pi \nu}{c} $$

Then, we only consider Neumann boundary conditions. The room wall identified by the equation $x_1=0$ will have an uniform velocity $w$ of $10$ $\frac{\text{m}}{\text{s}}$. All the other walls will be rigid. Hence, they will have a normal uniform velocity of $0$ $\frac{\text{m}}{\text{s}}$. Hence:

$$ N_{N} = 6 $$

If $\mathbf{x}$ is, again, the generic point in $\overline{\Omega}$:

$$ \mathbf{x} = \begin{bmatrix} x_{0} \\ x_{1} \\ x_{2} \end{bmatrix} \in \overline{\Omega}$$

$$ \begin{gather*} \\ \partial \Omega _{N} ^{0} \doteq \left\{ \mathbf{x} \in \overline{\Omega} : x_0 = 0 \right\} \\ \partial \Omega _{N} ^{1} \doteq \left\{ \mathbf{x} \in \overline{\Omega} : x_0 = l_{x_0} \right\} \\ \partial \Omega _{N} ^{2} \doteq \left\{ \mathbf{x} \in \overline{\Omega} : x_1 = 0 \right\} \\ \partial \Omega _{N} ^{3} \doteq \left\{ \mathbf{x} \in \overline{\Omega} : x_1 = l_{x_1} \right\} \\ \partial \Omega _{N} ^{4} \doteq \left\{ \mathbf{x} \in \overline{\Omega} : x_2 = 0 \right\} \\ \partial \Omega _{N} ^{5} \doteq \left\{ \mathbf{x} \in \overline{\Omega} : x_2 = l_{x_2} \right\} \end{gather*} $$

The normal boundary velocities for each $\partial\Omega_{N}^{l}$, for $l=0,1,2,3,4,5$, are:

$$ w_l \doteq \begin{cases} 10 \space \frac{\text{m}}{\text{s}} & \text{if } \space l = 2 \\ 0 \space \frac{\text{m}}{\text{s}} & \text{otherwise} \end{cases} $$

Finally, we consider the acoustic field as awoken by the Neumann boundary condition only. So, we put the source term $q$ to $0$:

$$ q = 0 $$

With all the choices above the weak form reduces to (after changing signs):

$$ \begin{aligned}

a\left(u,\phi\right) = \int_{\Omega} \nabla u \cdot \nabla \phi^\star d\mathbf{x} - k^2 \int_{\Omega} u \phi^{\star} d\mathbf{x} = \\

= \int_{\partial \Omega_{N}^2} g_{2} \phi^\star ds = L\left(\phi\right)

\end{aligned} $$

With $g$ being the acoustic flux for the boundary $\partial\Omega_{N}^{2}$:

$$ g_2 \doteq j \omega \rho w_2 = j \omega \rho 10 $$

As we noticed at the end of the Intro to `FEniCS`

- Part 2 episode we need to expand the integrands in real and imaginary parts. This will make most sense when we start writing the code. For the time being, let’s just prepare the ingredients:

$$ \begin{gather*}

u = \Re \left( u \right) + j \Im \left( u \right) \\
\phi = \Re \left( \phi \right) + j \Im \left( \phi \right) \\
g_2 = \Re \left( g_2 \right) + j \Im \left( g_2 \right) \\
\Rightarrow \nabla u = \Re\left(\nabla u \right) + j \Im\left(\nabla u \right) = \nabla \Re \left(u\right) + j \nabla \Im \left(u\right) \\
\Rightarrow \nabla \phi = \Re\left(\nabla \phi \right) + j \Im\left(\nabla \phi \right) = \nabla \Re \left(\phi\right) + j \nabla \Im \left(\phi\right)

\end{gather*}$$

Hence, we are ready to rewrite the integrals. The first inegral for $a\left(u,\phi\right)$ reads:

$$ \begin{aligned}

\int_{\Omega} \nabla u \cdot \nabla \phi^\star d\mathbf{x} = \\

= \int_{\Omega} \nabla\Re\left( u \right) \cdot \nabla \Re\left( \phi \right) + \nabla\Im \left( u \right) \cdot \nabla\Im \left( \phi \right) d\mathbf{x} + \\

+ j \int_{\Omega} \nabla\Im \left( u \right) \cdot \nabla\Re\left( \phi \right) - \nabla\Re\left( u \right) \cdot \nabla\Im \left( \phi \right) d\mathbf{x} = \\

= \int_{\Omega} \nabla\Re\left( u \right) \cdot \nabla \Re\left( \phi \right) d\mathbf{x} + \int_{\Omega} \nabla\Im \left( u \right) \cdot \nabla\Im \left( \phi \right) d\mathbf{x} + \\

+ j\left( \int_{\Omega} \nabla\Im \left( u \right) \cdot \nabla\Re\left( \phi \right) d\mathbf{x} - \int_{\Omega} \nabla\Re\left( u \right) \cdot \nabla\Im \left( \phi \right) d\mathbf{x} \right)

\end{aligned} $$

The second integral for $a\left(u,\phi\right)$ reads:

$$ \begin{aligned}

\int_{\Omega} u \phi^{\star} d\mathbf{x} = \\

= \int_{\Omega} \Re\left(u\right) \Re\left(\phi\right) + \Im\left(u\right) \Im\left(\phi\right) d\mathbf{x} + \\

+ j \int_{\Omega} \Im\left(u\right) \Re\left(\phi\right) - \Re\left(u\right) \Im\left(\phi\right) d\mathbf{x} = \\

= \int_{\Omega} \Re\left(u\right) \Re\left(\phi\right) d\mathbf{x} + \int_{\Omega} \Im\left(u\right) \Im\left(\phi\right) d\mathbf{x} + \\

+ j \left( \int_{\Omega} \Im\left(u\right) \Re\left(\phi\right) d\mathbf{x} - \int_{\Omega} \Re\left(u\right) \Im\left(\phi\right) d\mathbf{x} \right)

\end{aligned} $$

While $L\left(\phi\right)$ reads:

$$ \begin{aligned}

\int_{\partial \Omega_{N}^2} g_{2} \phi^\star ds = \\

= \int_{\partial \Omega_{N}^2} \Re\left(g_{2}\right) \Re\left(\phi\right) + \Im\left(g_{2}\right) \Im\left(\phi\right) ds + \\

+ j \int_{\partial \Omega_{N}^2} \Im\left(g_{2}\right) \Re\left(\phi\right) - \Re\left(g_{2}\right) \Im\left(\phi\right) ds = \\

= \int_{\partial \Omega_{N}^2} \Re\left(g_{2}\right) \Re\left(\phi\right) ds + \int_{\partial \Omega_{N}^2} \Im\left(g_{2}\right) \Im\left(\phi\right) ds + \\

+ j \left( \int_{\partial \Omega_{N}^2} \Im\left(g_{2}\right) \Re\left(\phi\right) ds - \int_{\partial \Omega_{N}^2} \Re\left(g_{2}\right) \Im\left(\phi\right) ds \right)

\end{aligned} $$

Now that all integrals are split in real and imaginary parts we have no trouble putting together the real and imaginary parts for $a \left(u,\phi\right) $ and $L\left(\phi\right)$. We will do this last step right in the code. For now note how each real part and imaginary part is the addition of two terms.

# Implementation

We will refer to the implementation as in the repository.

In the following we will assume you followed the installation instructions in the Intro to `FEniCS`

- Part 1 episode

## Preparation

First, let’s create a folder in your favorite location. Let’s `cd`

into it and create a `Python`

virtual environment:

```
mkdir project && cd project
python3 -m venv --system-site-packages venv
```

the `--system-site-packages`

ensures our new virtual environment has `FEniCS`

in it.

Then, you can open `PyCharm`

. You want to open an existing project and select our `project`

folder created above. `PyCharm`

could launch in different ways. If it opens on the editor select *File* > *Open*. For more information, refer to the `PyCharm`

Documentation.

However you end up opening the `project`

folder with `PyCharm`

you should have it loaded correctly with the `Python`

interpreter from our new virtual environment. If not, check the settings under *File* > *Settings* > *Project: project* > *Python Interpreter* and *File* > *Settings* > *Build, Execution, Deployment* > *Console* > *Python Console*. As an example, see the pictures below.

Note that it is entirely possible to create e new project all from the `PyCharm`

GUI. Make sure you tick *Inherit global site-packages* when you do so. Creating the environment beforehand is typically easier when cloning projects, although also this can be done with the `PyCharm`

GUI.

In the `PyCharm`

window you will see, on the left, a browser. Right click *project* and select *New* > *Python Package* as shown below.

We can call this package `simulation`

. You will see this crates a directory under `project`

called `simulation`

. In it you will find an empty `__init__.py`

file. Right click the `simulation`

folder in the browser and choose `New`

> `Python File`

. Let’s call the new file `helmholtz`

. We will implement the solver in this file.

Note that it is not necessary to create a `Python`

package to use `FEniCS`

. We can simply code a script. However, it is useful to make a package as it makes things nice and modular and simple to reuse.

## Code

The entire code is available here

First, let’s import the packages we need:

```
import numpy as np
import fenics
import pathlib
```

Then, we wrap all our simulation code in a function:

```
def simulate(
real_part_output_path: pathlib.Path,
imag_part_output_path: pathlib.Path
) -> None:
```

This function takes as arguments the paths to the real part and imaginary part files. That is, this code will produce two files, one for the real part of the field and one for the imaginary part. It will save them at the specified locations. These locations are provided as `Path`

objects from `pathlib`

. This is not necessary, but `pathlib`

actually makes working with files very easy and convenient, and it is then recommended. We will use `numpy`

for a few mathematical functions.

Now, let’s put all of our project parameters in a few variables:

```
nu = 57.17 # Frequency of the Simulation, Hz
c = 343 # Speed of sound in air, m/s
rho = 1.205 # Density of air, kg/m^3
w = 10 # Velocity normal to boundary, for inflow Neumann condition, m/s
l_x0 = 4 # Room size along x0, m
l_x1 = 5 # Room size along x1, m
l_x2 = 3 # Room size along x2, m
tol = 1e-10 # Tolerance for boundary condition definitions
```

The `tol`

parameter will become clearer later. Note how these paramters could all be arguments of the function. That is, you can make a function to simulate any room. In this code we kept them in the body for simplicity.

We now compute the mesh size. It should be smaller than one tenth of the wavelenght. Once we have the size we can easily compute the number of element per side of our room.

```
# Computing some useful stuff
omega = 2 * np.pi * nu # Angular frequency, rad/s
k = omega / c # Wave number, rad/m
s = c / (10 * nu) # Element Size
n_x0 = np.intp(np.ceil(l_x0 / s)) # Number of elements for each direction
n_x1 = np.intp(np.ceil(l_x1 / s))
n_x2 = np.intp(np.ceil(l_x2 / s))
```

We use `np.ceil`

to make sure that the number is the smallest integer bigger or equal to the exact size. `np.intp`

is used to convert the type to integer (`np.ceil`

will still return a floating point value).

Now it comes an important step. We are using `FEniCS`

. But `FEniCS`

does not support complex PDEs. But our PDE is complex. Turns out that a complex PDE is actually two coupled PDEs, one for the real part and one for the imaginary part. This is why we expanded the weak form in real and imaginary parts, so that we can easily see them. Solving coupled PDEs can be done in `FEniCS`

by using *Mixed Spaces*. The Solving PDEs in Python book covers this in details in the chapter 3.5. For more information refer to said chapter. In a nutshell, we will be solving two coupled PDEs, one for the real part and one for the imaginary one. This means that $j$, the imaginary unit, will disappear.

First, we crate the mesh:

```
# Making a mesh, then an element. Then, Mixed Space so that we can solve for real and imaginary parts.
mesh = fenics.BoxMesh(
fenics.Point(0, 0, 0),
fenics.Point(l_x0, l_x1, l_x2),
n_x0,
n_x1,
n_x2
)
```

The `BoxMesh`

function creates a box with one corner at the origin, the opposite corner at the chosen lengths. It also meshes it with the number of elements we previously computed. The mesh is tetrahedral.

Then, we define the type of element. We use a second order Lagrange element:

```
P = fenics.FiniteElement("Lagrange", mesh.ufl_cell(), 2)
```

This means that each cell in the mesh will be a second order element.

Finally, we can define our vector space `V`

as a mixed space over the mesh:

```
V = fenics.FunctionSpace(mesh, P * P)
```

The mixed aspect is given by `P * P`

, which signifies a mixed element for two fields, each of the `P`

kind. Note that `V`

is not the vector space $V$ we defined in the Intro to `FEniCS`

- Part 2 episode. The space `V`

is the space of functions represented as linear combinations of basis function (second order polynomials) over the mesh. This space is finite dimensional (finite number of basis functions). This is what allows to compute an approximate solution. That is: FEM computes a numerical solution by simplifying $V$. This step is automated by `FEniCS`

, as you see. `V`

largely works as $V$ in the sense that it is a linear space in which we search for a solution to the (now) simplified problem. So we still refer to it with the same letter. For more information refer to the FEM in a Nutshell episode.

In the `V`

space the unknown function we are searching for is called `u`

. This also goes by the name of *Trial Function*. The test function (which in $V$ was called $\phi$) will be called `phi`

. Hence:

```
u_re, u_im = fenics.TrialFunctions(V)
phi_re, phi_im = fenics.TestFunctions(V)
```

Since we are using a mixed space these functions are returned already split in real and imaginary parts, as we did in the maths above. For convenience we also define $k^2$ as a `Constant`

:

```
k_sq = fenics.Constant(k**2)
```

Now, we have all the ingredients ready for $a\left(u,\phi\right)$. Since we are using mixed space we simply sum the terms. That is, in the maths above we have split all terms that constitute $a\left(u,\phi\right)$ in real and imaginary parts. In mixed space formulation is like we had one $a\left(u,\phi\right)$ for the real part and one $a\left(u,\phi\right)$ for the imaginary part. These two are constituted of the real part and imaginary part of the original complex $a\left(u,\phi\right)$. The mixed space `a`

is then simply the sum of them. We report it below:

```
a = \
fenics.inner(fenics.nabla_grad(u_re), fenics.nabla_grad(phi_re)) * fenics.dx - k_sq * phi_re * u_re * fenics.dx + \
fenics.inner(fenics.nabla_grad(u_im), fenics.nabla_grad(phi_im)) * fenics.dx - k_sq * phi_im * u_im * fenics.dx + \
fenics.inner(fenics.nabla_grad(u_im), fenics.nabla_grad(phi_re)) * fenics.dx - k_sq * u_im * phi_re * fenics.dx - \
fenics.inner(fenics.nabla_grad(u_re), fenics.nabla_grad(phi_im)) * fenics.dx + k_sq * u_re * phi_im * fenics.dx
```

You will remember how each real and imaginary part of the integrals for $a\left(u,\phi\right)$ is composed of two terms. The first line for `a`

above puts together the first terms, the second one the second terms… and so on. `a`

does not need to be constructed in this specific order, but this form is clearer. It shows that a mixed space `a`

has a term with real parts only, one with an imaginary part only, and 2 mixed terms. The mixed terms are important, as they dictate the relationship the imaginary and real parts have: they are not independent. Note that `a`

is defined without any symbol for integral. `FEniCS`

already know it will need to take the integral of this thing, so we do not need it. We have to specify with respect what we have to integrate. Since these are volume integrals we use the default volume measure `fenics.dx`

. One can define custom ones, for example if some integral are within a sub-volume. We will do something similar right now for the boundaries. Before leaving this section, note the use of `inner`

and `nabla_grad`

. `inner`

is used to implement the scalar product between the gradients. As these are 1D tensors (vectors) then the result is the same as `dot`

. Similarly, since `u_re`

and `u_im`

are scalar fields, using `nabla`

or `nabla_grad`

is not really different. However, for vector fields the choice between `inner`

and `dot`

, `nabla`

and `nabla_grad`

becomes important. These differences are explained at page 25 and page 58 of Solving PDEs in Python.

All is left to do is to specify the boundary conditions. One very useful resource (together with Solving PDEs in Python) is Douglas N Arnold page. What we want to do is to assign a different Neumann boundary condition to each wall (although only one wall has a non-zero normal velocity). So, we need to split the boundary. Since we are already dealing with the mesh we need to mark the nodes that belong to the each boundary. Then, on each boundary, we define a custom measure so that when we define $L\left(\phi\right)$ each integral is restricted to the correct boundary. In other words: since `FEniCS`

does not have a symbol for the integral, but one for the measure (`dx`

, `ds`

…) we need the measure itself to be associated to the subset to which it is applied. Let’s go through.

First, we make a `MeshFunction`

that returns integers. We initialise it to return `0`

everywhere.

```
mf = fenics.MeshFunction("size_t", mesh, 2)
mf.set_all(0)
```

This function returns an unsigned integer of type `size_t`

(the largest unsigned integer type on your machine). Note how `size_t`

is a C\C++ type. This because all of these functions act on the C\C++ backend. Of course, this mesh function is defined over our mesh `mesh`

which we defined above. Since we are using `mf`

to mark boundaries, which are 2D entities, its dimensionality is `2`

. We will now edit it to return a different number for each boundary. This way, when our mesh function takes in a mesh node, it will return a number that tells `FEniCS`

on which boundary that is (if any).

To do so, we need to define the boundaries. A boundary is a domain. So we use `SubDomain`

. Defining a `SubDomain`

is easy. We simply create a class which has `SubDomain`

as a parent and override the `inside`

method. This is the method that `FEniCS`

will use to figure out if a mesh point is inside our subdomain. Below are all the subdomain for the various room walls:

```
class BX0(fenics.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and fenics.near(x[0], 0, tol)
class BXL(fenics.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and fenics.near(x[0], l_x0, tol)
class BY0(fenics.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and fenics.near(x[1], 0, tol)
class BYL(fenics.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and fenics.near(x[1], l_x1, tol)
class BZ0(fenics.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and fenics.near(x[2], 0, tol)
class BZL(fenics.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and fenics.near(x[2], l_x2, tol)
```

Note how the `inside`

method takes the mesh point `x`

, a 3D vector with the Cartesian coordinates of the point, and a boolean variable `on_boundary`

as input. That is, in the method `inside`

we will already know if `x`

is on the boundary or not based on the value of `on_boundary`

. We simply need to figure out which of our sub-boundaries that is. To do so, we use the function `near`

. Since everything is in floating point, we cannot check for equality, so we need a tolerance. This is what our `tol`

variable is for. For more information about this, check the excellent LEARN C++ tutorial about floating point comparison. As you can see, the `SubDomain`

definitions directly mirror those for $\partial \Omega _{n}^l$ that we provided above.

Now that we have the sub-domains we just need to instruct our mesh function to return a different integer on each of them. We first instantiate one object for each `SubDomain`

type we crated and use `mf`

to mark the nodes with an integer of our choice.

```
bx0 = BX0()
bx0.mark(mf, 1)
bxl = BXL()
bxl.mark(mf, 2)
by0 = BY0()
by0.mark(mf, 3)
byl = BYL()
byl.mark(mf, 4)
bz0 = BZ0()
bz0.mark(mf, 5)
bzl = BZL()
bzl.mark(mf, 6)
```

Now, we define the fluxes $g_l$. Or better, we define their “digital” counterpart. We only have dead boundaries (zero velocity) and alive boundaries (non-zero velocity). So:

```
# The flux is 0 for all rigid walls
flux_rig = 0
g_rig_re = fenics.Constant(np.real(flux_rig))
g_rig_im = fenics.Constant(np.imag(flux_rig))
# The flux is this for the active walls
flux_in = 1j * omega * rho * w
g_in_re = fenics.Constant(np.real(flux_in))
g_in_im = fenics.Constant(np.imag(flux_in))
```

Finally, we define our custom measure `ds`

over our `mesh`

. We use our mesh function `mf`

to ensure that `ds`

has the various measures for each of the sub-domains:

```
ds = fenics.Measure('ds', domain=mesh, subdomain_data=mf)
```

Now, if we use `ds(1)`

we will access the measure for boundary marked with `1`

by `mf`

, with `ds(2)`

that marked with `2`

and so forth. Hence, we can now put together `L`

, the FEM analogue of $L\left(\phi\right)$:

```
L = \
(g_rig_re * v_re + g_rig_im * v_im + g_rig_im * v_re - g_rig_re * v_im) * ds(1) + \
(g_rig_re * v_re + g_rig_im * v_im + g_rig_im * v_re - g_rig_re * v_im) * ds(2) + \
(g_rig_re * v_re + g_rig_im * v_im + g_rig_im * v_re - g_rig_re * v_im) * ds(6) + \
(g_rig_re * v_re + g_rig_im * v_im + g_rig_im * v_re - g_rig_re * v_im) * ds(4) + \
(g_rig_re * v_re + g_rig_im * v_im + g_rig_im * v_re - g_rig_re * v_im) * ds(5) + \
(g_in_re * v_re + g_in_im * v_im + g_in_im * v_re - g_in_re * v_im) * ds(3)
```

Since `g_rig_re`

and `g_rig_im`

are simply `0`

we could simply omit the terms with them. In that case we would be left with the last term only, as our equations from above tells. However, it is instructive to show the fully expanded form for `L`

, with all the different terms for each boundary.

All we have left to do is solve and save.

To solve, we define a new `Function`

over the vector space `V`

and invoke the `solve`

method:

```
u = fenics.Function(V)
fenics.solve(a == L, u)
```

This will solve the equation and put the result in `u`

. This is the easiest way to solve, but more advanced ways are also available. Control over the solver algorithm is provided by the more advanced methods. For now, let’s stick to the simple stuff.

Since `u`

is a mixed space solution we need to split it to access real and imaginary parts. We do so and rename as well for easy access in the final `VTK`

files. `VTK`

files are simple to create with the `File`

class. The code is below:

```
u_1, u_2 = u.split()
u_1.rename('re', 're')
u_2.rename('im', 'im')
vtk_file_1 = fenics.File(real_part_output_path.absolute().as_posix())
vtk_file_1 << u_1
vtk_file_2 = fenics.File(imag_part_output_path.absolute().as_posix())
vtk_file_2 << u_2
```

Running the code is simple. From your python environment simply run:

```
from simulation import helmholtz
import pathlib
path_re = pathlib.Path().joinpath('re.pvd')
path_im = pathlib.Path().joinpath('im.pvd')
helmholtz.run(path_re, path_im)
```

This will crate the files `re.pvd`

and `im.pvd`

in your current directory. Or, you can create a `run.py`

file in your `project`

folder with these contents. Then, simply run the file from your `project`

folder:

```
venv/bin/python run.py
```

# Results

We computed the results and solved the same problem with Elmer as well. We will not discuss the Elmer project setup as there are plenty of tutorials on this website about that. You can refer to the Elmer project files for more details. The results are comapred below.

Results are overall very similar. The real part is negligible, as expected from a damping free material with pure imaginary boundary flux. In fact, the weak form could be much simplified in this special case to yield a faster solution. However, we shown the entirety of it as it is more instructive. Interestingly the polarity of the real part is flipped in `FEniCS`

. The imaginary parts agree, overall, fairly well. The agreement is not exact. However, since we validated Elmer before, we know that our solution with `FEniCS`

is already a good starting point.

# Conclusion

We solved a simplified equation for a damping free material in a rectangular room. We derived the equations to do so and written the `FEniCS`

code to compute the solution. This solution was computed also with Elmer for benchmark. Results are very similar, a sign that our `FEniCS`

solution was setup properly.

Getting to solve PDEs in `FEniCS`

is more work with respect doing so with Elmer. However, we can setup to solve any PDE in a programmatic way. This means that we can write more complex simulation and analysis code all within `Python`

. Additionally, getting intimate with the weak form provides more insight about the underlying formalism at the base of PDEs.

In the next episodes we will explore solving PDEs with `FEniCS`

further.

# License Information

This work is licensed under a Creative Commons Attribution 4.0 International License.