This tutorial is part of FEMlium.

import os
import urllib

import firedrake
import firedrake.pyplot
import folium
import matplotlib.pyplot as plt
import numpy as np
import pyproj
import ufl
import viskex

/home/firedrake/firedrake/lib/python3.12/site-packages/pytools/persistent_dict.py:52: RecommendedHashNotFoundWarning: Unable to import recommended hash 'siphash24.siphash13', falling back to 'hashlib.sha256'. Run 'python3 -m pip install siphash24' to install the recommended hash.
  warn("Unable to import recommended hash 'siphash24.siphash13', "

import femlium

Auxiliary function to get a folium Map close to Lake Garda.

def get_garda_geo_map(boundary_icons: bool = False) -> folium.Map:
    """Get a map close to Lake Garda, and possibly add some boundary markers."""
    # Add map close to Lake Garda
    geo_map = folium.Map(location=[45.6389113, 10.7521368], zoom_start=10.3)

    # Add markers
    if boundary_icons:
        location_markers = {
            "Sarca": [45.87395405, 10.87087005],
            "Mincio": [45.43259035, 10.7007715]
        }
        location_colors = {
            "Sarca": "red",
            "Mincio": "green"
        }

        for key in location_markers.keys():
            folium.Marker(
                location=location_markers[key],
                tooltip=key,
                icon=folium.Icon(color=location_colors[key])
            ).add_to(geo_map)

    # Return folium map
    return geo_map

get_garda_geo_map()

Make this Notebook Trusted to load map: File -> Trust Notebook

Read the mesh from file with firedrake.

msh_filename = "data/garda.msh"
if not os.path.isfile(msh_filename):
    os.makedirs("data", exist_ok=True)
    msh_url = (
        "https://raw.githubusercontent.com/FEMlium/FEMlium/main/"
        "tutorials/01_introduction/data/garda.msh")
    with urllib.request.urlopen(msh_url) as response, open(msh_filename, "wb") as msh_file:
        msh_file.write(response.read())

mesh = firedrake.Mesh("data/garda.msh")

Plot the mesh using firedrake.pyplot.triplot.

fig = plt.figure(figsize=(12, 12))
firedrake.pyplot.triplot(mesh, axes=fig.gca())
fig.gca().axis("equal")

(np.float64(616658.039149),
 np.float64(647061.960851),
 np.float64(5029877.225841),
 np.float64(5085382.774159))

Plot the mesh and its markers using viskex.

viskex.firedrake.plot_mesh(mesh)

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

viskex.firedrake.plot_mesh_sets(mesh, 2, "subdomains")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

viskex.firedrake.plot_mesh_sets(mesh, 1, "boundaries")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

Define a pyproj Transformer to map between different reference systems, because the points read from file are stored a $(x, y)$ pairs in the EPSG32632 reference system, while the map produced by folium is based on (latitude, longitude) pairs in the EPSG4326 reference system.

transformer = pyproj.Transformer.from_crs("epsg:32632", "epsg:4326", always_xy=True)

We define a mesh plotter for meshes in firedrake format, which is implemented in femlium.FiredrakePlotter.

firedrake_plotter = femlium.FiredrakePlotter(transformer)

We use the firedrake_plotter to draw the mesh on top of the geographic map.

geo_map = get_garda_geo_map()
firedrake_plotter.add_mesh_to(geo_map, mesh)
geo_map

Make this Notebook Trusted to load map: File -> Trust Notebook

We may change the color and the weight of the line.

geo_map = get_garda_geo_map()
firedrake_plotter.add_mesh_to(geo_map, mesh, face_colors="red", face_weights=2)
geo_map

Make this Notebook Trusted to load map: File -> Trust Notebook

Furthermore, we may set the colors and the weights of the face representation to depend on the markers associated to each segment.

geo_map = get_garda_geo_map(boundary_icons=True)
face_colors = {
    0: "gray",
    1: "blue",
    2: "red",
    3: "green"
}
face_weights = {
    0: 1,
    1: 2,
    2: 5,
    3: 5
}
firedrake_plotter.add_mesh_to(geo_map, mesh, face_colors=face_colors, face_weights=face_weights)
geo_map

Make this Notebook Trusted to load map: File -> Trust Notebook

Cells can be colored as well, with a uniform color or depending on the cell markers. We start from a uniform color.

geo_map = get_garda_geo_map()
firedrake_plotter.add_mesh_to(geo_map, mesh, cell_colors="orange")
geo_map

Make this Notebook Trusted to load map: File -> Trust Notebook

We also show the case of colors being set from cell markers. There are two cell markers in this mesh, equal to 1 for the region close to the shoreline (colored in purple) and 2 for the rest of the domain (colored in yellow).

geo_map = get_garda_geo_map()
cell_colors = {
    1: "purple",
    2: "yellow"
}
firedrake_plotter.add_mesh_to(geo_map, mesh, cell_colors=cell_colors)
geo_map

Make this Notebook Trusted to load map: File -> Trust Notebook

Once can use colors associated to both cell and face markers on the same plot.

geo_map = get_garda_geo_map(boundary_icons=True)
firedrake_plotter.add_mesh_to(
    geo_map, mesh,
    cell_colors=cell_colors, face_colors=face_colors, face_weights=face_weights)
geo_map

Make this Notebook Trusted to load map: File -> Trust Notebook

In order to define a simple scalar field, we compute the centroid of the domain.

centroid = np.mean(mesh.coordinates.dat.data_ro, axis=0)

We may plot the centroid on top of the mesh.

geo_map = get_garda_geo_map()
firedrake_plotter.add_mesh_to(geo_map, mesh)
folium.Marker(location=transformer.transform(*centroid)[::-1], tooltip="Centroid").add_to(geo_map)
geo_map

Make this Notebook Trusted to load map: File -> Trust Notebook

The scalar field is defined as $s(\rho, \theta) = \frac{\rho}{\sqrt{1 - 0.5 \cos^2 \theta}}$, and is interpolated on a $\mathbb{P}^2$ finite element space. Here $(\rho, \theta)$ are the polar coordinates centered at the centroid.

scalar_function_space = firedrake.FunctionSpace(mesh, "CG", 2)

x = ufl.SpatialCoordinate(mesh)
rho = ufl.sqrt((x[0] - centroid[0])**2 + (x[1] - centroid[1])**2)
theta = ufl.atan2(x[1] - centroid[1], x[0] - centroid[0])
scalar_field = firedrake.interpolate(rho / ufl.sqrt(1 - 0.5 * ufl.cos(theta)**2), scalar_function_space)

/home/firedrake/firedrake/src/loopy/loopy/schedule/__init__.py:2255: LoopyWarning: in kernel expression_kernel: Falling back to a slow scheduler implementation due to: cannot schedule kernels with priority dependencies between sibling loop nests (add 'v1_scheduler_fallback' to silenced_warnings kernel attribute to disable)
  warn_with_kernel(

We next show a filled contour plot with 15 levels using firedrake.pyplot.tricontourf.

fig = plt.figure(figsize=(12, 12))
trif = firedrake.pyplot.tricontourf(scalar_field, axes=fig.gca(), levels=15, cmap="jet")
fig.colorbar(trif)
fig.gca().axis("equal")

(np.float64(618040.03559),
 np.float64(645679.96441),
 np.float64(5032400.20531),
 np.float64(5082859.79469))

We also show the same plot using viskex.

viskex.firedrake.plot_scalar_field(scalar_field, "scalar")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

In order to plot a field on a geographic map, we use again the firedrake_plotter. We may plot the same filled contour plot on the geographic map.

geo_map = get_garda_geo_map()
firedrake_plotter.add_scalar_field_to(geo_map, scalar_field, mode="contourf", levels=15, cmap="jet")
geo_map

Make this Notebook Trusted to load map: File -> Trust Notebook

Similarly, we can also use (unfilled) contour plots.

fig = plt.figure(figsize=(12, 12))
tri = firedrake.pyplot.tricontour(scalar_field, axes=fig.gca(), levels=15, cmap="jet")
fig.colorbar(tri)
fig.gca().axis("equal")

(np.float64(618040.03559),
 np.float64(645679.96441),
 np.float64(5032400.20531),
 np.float64(5082859.79469))

geo_map = get_garda_geo_map()
firedrake_plotter.add_scalar_field_to(geo_map, scalar_field, mode="contour", levels=15, cmap="jet")
geo_map

Make this Notebook Trusted to load map: File -> Trust Notebook

One may also combine mesh plots and solution plots.

geo_map = get_garda_geo_map()
firedrake_plotter.add_mesh_to(geo_map, mesh, face_colors="grey")
firedrake_plotter.add_scalar_field_to(geo_map, scalar_field, mode="contour", levels=15, cmap="jet")
geo_map

Make this Notebook Trusted to load map: File -> Trust Notebook

We next define a vector field $\mathbf{v}(\rho, \theta) = \begin{bmatrix}-\rho \sin \theta\\\rho \cos\theta \end{bmatrix}$.

vector_function_space = firedrake.VectorFunctionSpace(mesh, "CG", 2)

vector_field = firedrake.interpolate(
    ufl.as_vector((- rho * ufl.sin(theta), rho * ufl.cos(theta))), vector_function_space)

/home/firedrake/firedrake/src/loopy/loopy/schedule/__init__.py:2255: LoopyWarning: in kernel expression_kernel: Falling back to a slow scheduler implementation due to: cannot schedule kernels with priority dependencies between sibling loop nests (add 'v1_scheduler_fallback' to silenced_warnings kernel attribute to disable)
  warn_with_kernel(

We may obtain contourf or contour plots of the magnitude of the vector field.

viskex.firedrake.plot_vector_field(vector_field, "vector")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

geo_map = get_garda_geo_map()
firedrake_plotter.add_vector_field_to(geo_map, vector_field, mode="contourf", levels=15, cmap="jet")
geo_map

Make this Notebook Trusted to load map: File -> Trust Notebook

geo_map = get_garda_geo_map()
firedrake_plotter.add_vector_field_to(geo_map, vector_field, mode="contour", levels=15, cmap="jet")
geo_map

Make this Notebook Trusted to load map: File -> Trust Notebook

Vector field can also be plotted using a quiver. We first see the quiver plot obtained either with firedrake.quiver or with viskex.

fig = plt.figure(figsize=(12, 12))
quiv = firedrake.pyplot.quiver(vector_field, axes=fig.gca(), cmap="jet")
fig.colorbar(quiv)
fig.gca().axis("equal")

(np.float64(616658.039149),
 np.float64(647061.960851),
 np.float64(5029877.225841),
 np.float64(5085382.774159))

viskex.firedrake.plot_vector_field(vector_field, "vector", glyph_factor=1e-1)

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

A similar plot can rendered on top of the geographic map.

geo_map = get_garda_geo_map()
firedrake_plotter.add_vector_field_to(geo_map, vector_field, mode="quiver", scale=1e-1, cmap="jet")
geo_map

Make this Notebook Trusted to load map: File -> Trust Notebook