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()
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
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
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
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
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
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
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
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
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
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
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.
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
geo_map = get_garda_geo_map()
dolfinx_plotter.add_vector_field_to(geo_map, vector_field, mode="contour", levels=15, cmap="jet")
geo_map
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