This tutorial is part of FEMlium.

import os
import urllib

import dolfinx.fem
import dolfinx.io
import dolfinx.plot
import folium
import mpi4py.MPI
import numpy as np
import numpy.typing as npt
import pyproj
import viskex

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 dolfinx.

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, *_ = dolfinx.io.gmshio.read_from_msh(
    "data/garda.msh", comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2)
assert subdomains is not None
assert boundaries is not None

Info    : Reading 'data/garda.msh'...
Info    : 2390 entities
Info    : 3937 nodes
Info    : 7872 elements
Info    : Done reading 'data/garda.msh'

Plot the mesh using viskex.

viskex.dolfinx.plot_mesh(mesh)

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
error: XDG_RUNTIME_DIR is invalid or not set in the environment.
error: XDG_RUNTIME_DIR is invalid or not set in the environment.

Plot the subdomain markers using viskex.

viskex.dolfinx.plot_mesh_tags(mesh, subdomains, "subdomains")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
error: XDG_RUNTIME_DIR is invalid or not set in the environment.

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

dolfinx_plotter = femlium.DolfinxPlotter(transformer)

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

geo_map = get_garda_geo_map()
dolfinx_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()
dolfinx_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
}
dolfinx_plotter.add_mesh_to(
    geo_map, mesh, face_mesh_tags=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()
dolfinx_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"
}
dolfinx_plotter.add_mesh_to(geo_map, mesh, cell_mesh_tags=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)
dolfinx_plotter.add_mesh_to(
    geo_map, mesh,
    cell_mesh_tags=subdomains, face_mesh_tags=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.geometry.x[:, :2], axis=0)

We may plot the centroid on top of the mesh.

geo_map = get_garda_geo_map()
dolfinx_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 = dolfinx.fem.functionspace(mesh, ("CG", 2))

def scalar_field_eval(x: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
    """Evaluate the scalar field."""
    rho = np.sqrt((x[0] - centroid[0])**2 + (x[1] - centroid[1])**2)
    theta = np.arctan2(x[1] - centroid[1], x[0] - centroid[0])
    return rho / np.sqrt(1 - 0.5 * np.cos(theta)**2)

scalar_field = dolfinx.fem.Function(scalar_function_space)
scalar_field.interpolate(scalar_field_eval)

We next show a filled contour plot using viskex.

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

error: XDG_RUNTIME_DIR is invalid or not set in the environment.

error: XDG_RUNTIME_DIR is invalid or not set in the environment.

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

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

geo_map = get_garda_geo_map()
dolfinx_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()
dolfinx_plotter.add_mesh_to(geo_map, mesh, face_colors="grey")
dolfinx_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 = dolfinx.fem.functionspace(mesh, ("CG", 2, (mesh.geometry.dim, )))

def vector_field_eval(x: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
    """Evaluate the vector field."""
    rho = np.sqrt((x[0] - centroid[0])**2 + (x[1] - centroid[1])**2)
    theta = np.arctan2(x[1] - centroid[1], x[0] - centroid[0])
    values = np.zeros((2, x.shape[1]))
    values[0] = - rho * np.sin(theta)
    values[1] = rho * np.cos(theta)
    return values

vector_field = dolfinx.fem.Function(vector_function_space)
vector_field.interpolate(vector_field_eval)

We first see a plot obtained with viskex, which shows either the magnitude of the vector field or its representation using glyphs.

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

error: XDG_RUNTIME_DIR is invalid or not set in the environment.

error: XDG_RUNTIME_DIR is invalid or not set in the environment.

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

error: XDG_RUNTIME_DIR is invalid or not set in the environment.

error: XDG_RUNTIME_DIR is invalid or not set in the environment.

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

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

Also a quiver plot can rendered on top of the geographic map.

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