1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
|
# ---
# jupyter:
# jupytext:
# text_representation:
# extension: .py
# format_name: light
# format_version: '1.5'
# jupytext_version: 1.13.6
# ---
# # Interpolation and IO
#
# Copyright (C) 2022 Garth N. Wells
#
#
# ```{admonition} Download sources
# :class: download
# * {download}`Python script <./demo_interpolation-io.py>`
# * {download}`Jupyter notebook <./demo_interpolation-io.ipynb>`
# ```
# This demo shows how to:
# - Interpolate functions into vector-element $H(\mathrm{curl})$
# finite element spaces
# - Interpolate these special finite elements into discontinuous Lagrange
# spaces for artifact-free visualisation.
# $H(\mathrm{curl})$ finite element spaces, and the interpolation of
# these special finite elements in discontinuous Lagrange spaces for
# artifact-free visualisation.
# +
from pathlib import Path
from mpi4py import MPI
import numpy as np
from sys import byteorder
from dolfinx import default_scalar_type, has_adios2, plot
from dolfinx.fem import Function, functionspace
from dolfinx.mesh import CellType, create_rectangle, locate_entities
# -
# Create a mesh. For later in the demo we need to ensure that a boundary
# between cells is located at $x_0=0.5$.
msh = create_rectangle(MPI.COMM_WORLD, ((0.0, 0.0), (1.0, 1.0)), (16, 16), CellType.triangle)
# Create a Nédélec function space and finite element Function
V = functionspace(msh, ("Nedelec 1st kind H(curl)", 1))
u = Function(V, dtype=default_scalar_type)
# Find cells with *all* vertices (0) $x_0 <= 0.5$ or (1) $x_0 >= 0.5$:
tdim = msh.topology.dim
cells0 = locate_entities(msh, tdim, lambda x: x[0] <= 0.5)
cells1 = locate_entities(msh, tdim, lambda x: x[0] >= 0.5)
# Interpolate in the Nédélec/H(curl) space a vector-valued expression
# `f`, where $f \cdot n$ is discontinuous at $x_0 = 0.5$ and $f \cdot
# e$ is continuous.
u.interpolate(lambda x: np.vstack((x[0], x[1])), cells0)
u.interpolate(lambda x: np.vstack((x[0] + 1, x[1])), cells1)
# Create a vector-valued discontinuous Lagrange space and function, and
# interpolate the $H({\rm curl})$ function `u`
gdim = msh.geometry.dim
V0 = functionspace(msh, ("Discontinuous Lagrange", 1, (gdim,)))
u0 = Function(V0, dtype=default_scalar_type)
u0.interpolate(u)
# We save the interpolated function `u0` in VTX format. When visualising
# the field, at $x_0 = 0.5$ the $x_0$-component should appear
# discontinuous and the $x_1$-component should appear continuous.
# We use the {py:class}`dolfinx.io.VTXWriter` to store the data.
out_folder = Path("output_nedelec")
out_folder.mkdir(parents=True, exist_ok=True)
if has_adios2:
from dolfinx.io import VTXWriter
with VTXWriter(msh.comm, out_folder / "output_nedelec.bp", u0) as f:
f.write(0.0)
else:
print("ADIOS2 required for VTX output")
# Plot the functions
try:
import pyvista
cells, types, x = plot.vtk_mesh(V0)
grid = pyvista.UnstructuredGrid(cells, types, x)
values = np.zeros((x.shape[0], 3), dtype=np.float64)
values[:, : msh.topology.dim] = u0.x.array.reshape(x.shape[0], msh.topology.dim).real
grid.point_data["u"] = values
pl = pyvista.Plotter(shape=(2, 2))
pl.subplot(0, 0)
pl.add_text("magnitude", font_size=12, color="black", position="upper_edge")
pl.add_mesh(grid.copy(), show_edges=True)
pl.subplot(0, 1)
glyphs = grid.glyph(orient="u", factor=0.08)
pl.add_text("vector glyphs", font_size=12, color="black", position="upper_edge")
pl.add_mesh(glyphs, show_scalar_bar=False)
pl.add_mesh(grid.copy(), style="wireframe", line_width=2, color="black")
pl.subplot(1, 0)
pl.add_text("x-component", font_size=12, color="black", position="upper_edge")
pl.add_mesh(grid.copy(), component=0, show_edges=True)
pl.subplot(1, 1)
pl.add_text("y-component", font_size=12, color="black", position="upper_edge")
pl.add_mesh(grid.copy(), component=1, show_edges=True)
pl.view_xy()
pl.link_views()
# If pyvista environment variable is set to off-screen (static)
# plotting save png
if pyvista.OFF_SCREEN:
pl.screenshot(out_folder / "uh.png")
else:
pl.show()
except ModuleNotFoundError:
print("pyvista is required to visualise the solution")
|
