```
%matplotlib inline
import matplotlib.pyplot as plt
import firedrake
import icepack, icepack.plot
```

# Meshes, functions¶

This demo will show some of the basic features of the finite element modeling package firedrake, which icepack is built on. These are the kinds of things you’d need to know before you begin to model anything. First, we’ll show to create a simple mesh, and how to define a function on that mesh symbolically using handles for the x- and y-coordinates of the points of the domain. Then we’ll demonstrate some of the plotting routines included in icepack, which are meant to mimic as closely as possible the corresponding routines in matplotlib. Finally, we’ll show how to analyze a function by evaluating some area and contour integrals.

## Meshes¶

First, we need to make a mesh. Firedrake includes routines for generating meshes of simple domains; to see all of the meshes that are available, you can type

```
help(firedrake.utility_meshes)
```

from the python interpreter. In this demo we’ll use a mesh of the unit square.

```
nx, ny = 16, 16
mesh = firedrake.UnitSquareMesh(nx, ny)
```

Mesh plots are shown with a different color for each segments of the boundary and a legend. This is useful to when applying different boundary conditions to different segments, for example, Dirichlet boundary conditions on part of the domain and Neumann boundary conditions on the other part.

```
fig, axes = icepack.plot.subplots()
icepack.plot.triplot(mesh, axes=axes)
plt.show(fig)
```

## Functions¶

We can make scalar and vector fields defined on this mesh by interpolating the values of an algebraic expression to some space of functions.

First, we have to construct a function space `Q`

. We already have a
mesh, and we need to decide what element family and polynomial degree.
In almost all the cases you’ll encounter, the element family consists of
continuous piecewise-polynomial functions in each triangle, which is
abbreviated to `CG`

for “continuous Galerkin”.

Next, we have to make an expression for the function and interpolate it
to the function space `Q`

. The function
`firedrake.SpatialCoordinate`

returns two symbolic objects `x`

,
`y`

for the coordinates of each point of the mesh. We can then use
these symbols to define expressions for the function we’d like to
analyze. In this case, I’ve chosen the Rosenbrock
function.
Firedrake has built-in functions for evaluating various transcendental
functions of the coordinates, for example the sine, cosine, exponential,
logarithm, etc. To see all of the available functions, you can check the
namespace `ufl.operators`

.

Finally, the function `firedrake.interpolate`

takes in an expression
and a function space, and returns a field from that function space.

```
Q = firedrake.FunctionSpace(mesh, family='CG', degree=2)
x, y = firedrake.SpatialCoordinate(mesh)
a, b = 0.5, 10.0
expr = (a - x)**2 + b*(y - x**2)**2
q = firedrake.interpolate(expr, Q)
fig, axes = icepack.plot.subplots()
contours = icepack.plot.tricontourf(q, 36, axes=axes)
fig.colorbar(contours)
plt.show(fig)
```

We can make vector fields in much the same way as we make scalar fields, the only difference being the type of function space. There are many ways to plot a vector field. Here we show a plot of the streamlines of the vector field, colored according to the magnitude of the vector. This method is expensive, especially if you set the resolution or spacing of the streamlines to be too small, but nonetheless it produces really nice graphics.

```
V = firedrake.VectorFunctionSpace(mesh, family='CG', degree=2)
v = firedrake.interpolate(firedrake.grad(expr), V)
fig, axes = icepack.plot.subplots()
opts = {'precision': 1/128, 'density': 1/16, 'max_num_points': 128}
streamlines = icepack.plot.streamplot(v, axes=axes, **opts)
fig.colorbar(streamlines)
plt.show(fig)
```

## Post-processing¶

We’ve already seen about how to visualize a scalar or a vector field;
this is just a call to the function `icepack.plot`

. To analyze the
results of a simulation, you might also want to evaluate scalar or
vector fields at points in the domain:

```
z = (a, a**2)
print(q(z))
print(v(z))
```

```
-2.6731072308911782e-17
[ 3.46901282e-15 -2.52175273e-15]
```

Firedrake also provides a rich set of operations for evaluating integral
expressions of scalar and vector fields. Much like how the function
`SpatialCoordinate`

gives you two symbols `x`

, `y`

that represent
the coordinates of each point of the mesh, firedrake also provides an
object `dx`

that represents the differential area element.

To define an integral, we multiply an expression by `dx`

, and then
call the function `firedrake.assemble`

to evaluate it.

```
from firedrake import inner, grad, dx, assemble
print(assemble(x*y*dx))
print(assemble(q*dx))
print(assemble(inner(v, v)*dx))
```

```
0.25000000000000006
2.0833346048990866
172.0160966979131
```

We can also evaluate integrals over the boundary of the mesh. Boundary
integrals are specified using the differential surface element `ds`

instead of `dx`

.

Sometimes we need to integrate over only part of the boundary. We can
specify which part by passing the right numeric ID to `ds`

; these are
the same IDs that are color-coded in the mesh plot above.

Finally, you might also want to evaluate the flux of a vector field
across the boundary. To evaluate a flux, we need to get the unit outward
normal vector to the boundary. You can get a symbol representing the
unit outward normal by calling the function `firedrake.FacetNormal`

.

```
from firedrake import ds
print(assemble(q*ds))
print(assemble(q*ds(2)))
n = firedrake.FacetNormal(mesh)
print(assemble(inner(v, n) * ds))
```

```
14.666669209798183
3.5833333333333344
42.000000000000014
```

## Conclusion¶

Icepack uses the functions in firedrake to implement solvers for the various physics problems that show up in ice sheet modeling. The firedrake routines that we’ve demonstrated here can be used for analyzing simulation results, either through visualization, evaluating fields at points, or evaluating integrals of fields. They can also be used to define input fields of simulations, if those fields have a simple analytical expression. In the next demo, we’ll show how to use icepack to solve for the velocity and thickness of a floating ice shelf with a synthetic geometry. In later demos, we’ll show how to use real observational data sets, as well as more complicated geometries imported from mesh generators.

To read more about firedrake, you can visit their documentation or check out some of the demos.