This tutorial is part of FEMlium.

import os
import urllib

import folium
import matplotlib.pyplot as plt
import meshio
import numpy as np
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 from file with meshio.

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 = meshio.read("data/garda.msh")

Plot the mesh using matplotlib.

fig = plt.figure(figsize=(12, 12))
fig.gca().triplot(mesh.points[:, 0], mesh.points[:, 1], mesh.cells_dict["triangle"])
fig.gca().axis("equal")

(616658.039149, 647061.960851, 5029877.225841, 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 meshio format, which is implemented in femlium.MeshioPlotter.

mesh_plotter = femlium.MeshioPlotter(transformer)

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

geo_map = get_garda_geo_map()
mesh_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()
mesh_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
}
mesh_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()
mesh_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"
}
mesh_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)
mesh_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.points[:, 0:2], axis=0)

We may plot the centroid on top of the mesh.

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

We define polar coordinates centered at the centroid.

rho = np.sqrt((mesh.points[:, 0] - centroid[0])**2 + (mesh.points[:, 1] - centroid[1])**2)
theta = np.arctan2(mesh.points[:, 1] - centroid[1], mesh.points[:, 0] - centroid[0])

The scalar field is defined as $s(\rho, \theta) = \frac{\rho}{\sqrt{1 - 0.5 \cos^2 \theta}}$.

scalar_field = rho / np.sqrt(1 - 0.5 * np.cos(theta)**2)

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

fig = plt.figure(figsize=(12, 12))
trif = fig.gca().tricontourf(
    mesh.points[:, 0], mesh.points[:, 1], mesh.cells_dict["triangle"], scalar_field, levels=15, cmap="jet")
fig.colorbar(trif)
fig.gca().axis("equal")

(618040.03559, 645679.96441, 5032400.20531, 5082859.79469)

In order to plot a field on a geographic map, we use a femlium.NumpyPlotter.

solution_plotter = femlium.NumpyPlotter(transformer)

We plot the same filled contour plot on the geographic map.

geo_map = get_garda_geo_map()
solution_plotter.add_scalar_field_to(
    geo_map, mesh.points[:, 0:2], mesh.cells_dict["triangle"], scalar_field,
    mode="contourf", levels=15, cmap="jet")
geo_map

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

If the optional arguments are not provided then the contourf mode will be used, with 10 levels and jet colormap.

geo_map = get_garda_geo_map()
solution_plotter.add_scalar_field_to(
    geo_map, mesh.points[:, 0:2], mesh.cells_dict["triangle"], scalar_field)
geo_map

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

Instead of using an integer levels input to generate an equispaced array of levels between the minimum and the maximum of the scalar field, the user may pass levels as a numpy array.

geo_map = get_garda_geo_map()
solution_plotter.add_scalar_field_to(
    geo_map, mesh.points[:, 0:2], mesh.cells_dict["triangle"], scalar_field,
    mode="contourf", levels=np.linspace(scalar_field.min(), scalar_field.max(), 15), cmap="jet")
geo_map

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

levels may also be passed in as a list of numbers, rather than a numpy array.

geo_map = get_garda_geo_map()
solution_plotter.add_scalar_field_to(
    geo_map, mesh.points[:, 0:2], mesh.cells_dict["triangle"], scalar_field,
    mode="contourf", levels=np.linspace(scalar_field.min(), scalar_field.max(), 15).tolist(), 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 = fig.gca().tricontour(
    mesh.points[:, 0], mesh.points[:, 1], mesh.cells_dict["triangle"], scalar_field, levels=15, cmap="jet")
fig.colorbar(tri)
fig.gca().axis("equal")

(618040.03559, 645679.96441, 5032400.20531, 5082859.79469)

geo_map = get_garda_geo_map()
solution_plotter.add_scalar_field_to(
    geo_map, mesh.points[:, 0:2], mesh.cells_dict["triangle"], 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()
mesh_plotter.add_mesh_to(geo_map, mesh, face_colors="grey")
solution_plotter.add_scalar_field_to(
    geo_map, mesh.points[:, 0:2], mesh.cells_dict["triangle"], 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_field = np.array([- rho * np.sin(theta), rho * np.cos(theta)]).T

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

geo_map = get_garda_geo_map()
solution_plotter.add_vector_field_to(
    geo_map, mesh.points[:, 0:2], mesh.cells_dict["triangle"], vector_field,
    mode="contourf", levels=15, cmap="jet")
geo_map

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

If the optional arguments are not provided then the contourf mode will be used, with 10 levels and jet colormap.

geo_map = get_garda_geo_map()
solution_plotter.add_vector_field_to(
    geo_map, mesh.points[:, 0:2], mesh.cells_dict["triangle"], vector_field)
geo_map

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

Instead of using an integer levels input to generate an equispaced array of levels between the minimum and the maximum of the scalar field, the user may pass levels as a numpy array.

geo_map = get_garda_geo_map()
vector_field_magnitude = np.linalg.norm(vector_field, axis=1)
solution_plotter.add_vector_field_to(
    geo_map, mesh.points[:, 0:2], mesh.cells_dict["triangle"], vector_field,
    mode="contourf", levels=np.linspace(vector_field_magnitude.min(), vector_field_magnitude.max(), 15),
    cmap="jet")
geo_map

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

levels may also be passed in as a list of numbers, rather than a numpy array.

geo_map = get_garda_geo_map()
solution_plotter.add_vector_field_to(
    geo_map, mesh.points[:, 0:2], mesh.cells_dict["triangle"], vector_field,
    mode="contourf", levels=np.linspace(vector_field_magnitude.min(), vector_field_magnitude.max(), 15).tolist(),
    cmap="jet")
geo_map

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

The mode can be changed to contour.

geo_map = get_garda_geo_map()
solution_plotter.add_vector_field_to(
    geo_map, mesh.points[:, 0:2], mesh.cells_dict["triangle"], vector_field,
    mode="contour", levels=15, cmap="jet")
geo_map

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

Vector fields can also be plotted using a quiver. We first see the quiver plot in matplotlib.

fig = plt.figure(figsize=(12, 12))
quiv = fig.gca().quiver(
    mesh.points[:, 0], mesh.points[:, 1],
    vector_field[:, 0], vector_field[:, 1], np.linalg.norm(vector_field, axis=1),
    cmap="jet")
fig.colorbar(quiv)
fig.gca().axis("equal")

(616658.039149, 647061.960851, 5029877.225841, 5085382.774159)

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

geo_map = get_garda_geo_map()
solution_plotter.add_vector_field_to(
    geo_map, mesh.points[:, 0:2], mesh.cells_dict["triangle"], vector_field,
    mode="quiver", scale=1e-1, cmap="jet")
geo_map

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