This tutorial is part of FEMlium.

import os
import urllib

import dolfin
import folium
import matplotlib.pyplot as plt
import meshio
import numpy as np
import numpy.typing
import pyproj

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, the subdomain markers and the boundary markers from file with dolfin.

def gmsh_to_fenics(msh_path: str) -> tuple[
        dolfin.Mesh, dolfin.cpp.mesh.MeshFunctionSizet, dolfin.cpp.mesh.MeshFunctionSizet]:
    """Convert a mesh from gmsh to FEniCS."""
    assert msh_path.endswith(".msh")
    base_path = msh_path[:-4]

    # Read back in the mesh with meshio
    meshio_mesh = meshio.read(msh_path)

    # Save volume mesh in xdmf format
    mesh_xdmf_path = base_path + "_mesh.xdmf"
    if os.path.exists(mesh_xdmf_path):
        os.remove(mesh_xdmf_path)
    if os.path.exists(mesh_xdmf_path.replace(".xdmf", ".h5")):
        os.remove(mesh_xdmf_path.replace(".xdmf", ".h5"))
    points = meshio_mesh.points[:, :2]
    cells = meshio_mesh.cells_dict["triangle"]
    if ("gmsh:physical" in meshio_mesh.cell_data_dict
            and "triangle" in meshio_mesh.cell_data_dict["gmsh:physical"]):
        subdomains_data = meshio_mesh.cell_data_dict["gmsh:physical"]["triangle"]
    else:
        subdomains_data = np.zeros_like(cells)
    meshio.write(
        mesh_xdmf_path,
        meshio.Mesh(
            points=points,
            cells={"triangle": cells},
            cell_data={"subdomains": [subdomains_data]}
        )
    )

    # Save boundary mesh in xdmf format
    boundaries_xdmf_path = base_path + "_boundaries.xdmf"
    if os.path.exists(boundaries_xdmf_path):
        os.remove(boundaries_xdmf_path)
    if os.path.exists(boundaries_xdmf_path.replace(".xdmf", ".h5")):
        os.remove(boundaries_xdmf_path.replace(".xdmf", ".h5"))
    facets = meshio_mesh.cells_dict["line"]
    if ("gmsh:physical" in meshio_mesh.cell_data_dict
            and "line" in meshio_mesh.cell_data_dict["gmsh:physical"]):
        boundaries_data = meshio_mesh.cell_data_dict["gmsh:physical"]["line"]
    else:
        boundaries_data = np.zeros_like(facets)
    meshio.write(
        boundaries_xdmf_path,
        meshio.Mesh(
            points=points,
            cells={"line": facets},
            cell_data={"boundaries": [boundaries_data]}
        )
    )

    # Read back in the mesh with dolfin
    mesh = dolfin.Mesh()
    with dolfin.XDMFFile(mesh_xdmf_path) as infile:
        infile.read(mesh)

    # Read back in subdomains with dolfin
    subdomains_mvc = dolfin.MeshValueCollection("size_t", mesh, mesh.topology().dim())
    with dolfin.XDMFFile(mesh_xdmf_path) as infile:
        infile.read(subdomains_mvc, "subdomains")
    subdomains = dolfin.cpp.mesh.MeshFunctionSizet(mesh, subdomains_mvc)

    # Clean up mesh file
    os.remove(mesh_xdmf_path)
    os.remove(mesh_xdmf_path.replace(".xdmf", ".h5"))

    # Read back in boundaries with dolfin, and explicitly set to 0 any facet
    # which had not been marked by gmsh
    boundaries_mvc = dolfin.MeshValueCollection("size_t", mesh, mesh.topology().dim() - 1)
    with dolfin.XDMFFile(boundaries_xdmf_path) as infile:
        infile.read(boundaries_mvc, "boundaries")
    boundaries_mvc_dict = boundaries_mvc.values()
    for c in dolfin.cells(mesh):
        for f, _ in enumerate(dolfin.facets(c)):
            if (c.index(), f) not in boundaries_mvc_dict:
                boundaries_mvc.set_value(c.index(), f, 0)
    boundaries = dolfin.cpp.mesh.MeshFunctionSizet(mesh, boundaries_mvc)

    # Clean up boundary mesh file
    os.remove(boundaries_xdmf_path)
    os.remove(boundaries_xdmf_path.replace(".xdmf", ".h5"))
    return mesh, subdomains, boundaries

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, subdomains, boundaries = gmsh_to_fenics("data/garda.msh")

Plot the mesh using dolfin.plot.

fig = plt.figure(figsize=(12, 12))
dolfin.plot(mesh)
fig.gca().axis("equal")

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

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 dolfin format, which is implemented in femlium.DolfinPlotter.

dolfin_plotter = femlium.DolfinPlotter(transformer)

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

geo_map = get_garda_geo_map()
dolfin_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()
dolfin_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
}
dolfin_plotter.add_mesh_to(
    geo_map, mesh, face_mesh_function=boundaries, 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()
dolfin_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"
}
dolfin_plotter.add_mesh_to(geo_map, mesh, cell_mesh_function=subdomains, 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)
dolfin_plotter.add_mesh_to(
    geo_map, mesh,
    cell_mesh_function=subdomains, face_mesh_function=boundaries,
    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(), axis=0)

We may plot the centroid on top of the mesh.

geo_map = get_garda_geo_map()
dolfin_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 = dolfin.FunctionSpace(mesh, "CG", 2)

Calling FFC just-in-time (JIT) compiler, this may take some time.

class ScalarField(dolfin.UserExpression):
    """Expression of the scalar field."""

    def eval_cell(
        self, value: np.typing.NDArray[np.float64], x: np.typing.NDArray[np.float64], cell: int
    ) -> None:
        """Evaulate the expression."""
        rho = np.sqrt((x[0] - centroid[0])**2 + (x[1] - centroid[1])**2)
        theta = np.arctan2(x[1] - centroid[1], x[0] - centroid[0])
        value[0] = rho / np.sqrt(1 - 0.5 * np.cos(theta)**2)

    def value_shape(self) -> tuple[int]:
        """Shape of a scalar expression."""
        return ()

scalar_field = dolfin.interpolate(ScalarField(), scalar_function_space)

We next show a filled contour plot with 15 levels using dolfin.plot.

fig = plt.figure(figsize=(12, 12))
trif = dolfin.plot(scalar_field, mode="contourf", 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))

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

geo_map = get_garda_geo_map()
dolfin_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 = dolfin.plot(scalar_field, mode="contour", 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()
dolfin_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()
dolfin_plotter.add_mesh_to(geo_map, mesh, face_colors="grey")
dolfin_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 = dolfin.VectorFunctionSpace(mesh, "CG", 2)

Calling FFC just-in-time (JIT) compiler, this may take some time.

class VectorField(dolfin.UserExpression):
    """Expression of the vector field."""

    def eval_cell(
        self, value: np.typing.NDArray[np.float64], x: np.typing.NDArray[np.float64], cell: int
    ) -> None:
        """Evaulate the expression."""
        rho = np.sqrt((x[0] - centroid[0])**2 + (x[1] - centroid[1])**2)
        theta = np.arctan2(x[1] - centroid[1], x[0] - centroid[0])
        value[0] = - rho * np.sin(theta)
        value[1] = rho * np.cos(theta)

    def value_shape(self) -> tuple[int]:
        """Shape of a vector expression."""
        return (2, )

vector_field = dolfin.interpolate(VectorField(), vector_function_space)

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

geo_map = get_garda_geo_map()
dolfin_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()
dolfin_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 with dolfin.plot.

fig = plt.figure(figsize=(12, 12))
quiv = dolfin.plot(vector_field, mode="glyphs", 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))

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

geo_map = get_garda_geo_map()
dolfin_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