UnstructuredMesh
class, but it will not focus on creating meshes.import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Line3DCollection
import salvus.namespace as sn
exodus
or .h5
file) or from a Salvus
meshing routine. For transferring meshes between scripts, to other machines,
or to Mondaic for debugging, we always prefer the .h5
format..h5
file. This is an HDF5 file, but that only
specifies in what format the data is stored. As HDF5 files are only data
containers, the actual contents of HDF5 files is program specific, and no
other program would be expected to understand these meshes without
additionaly information.mesh = sn.UnstructuredMesh.from_h5("data/block_with_hole_small.h5")
from_h5
implies the existence of its
counterpart, and indeed, this is the Salvus-native way to write your mesh to
a file:mesh.write_h5( "data/block_with_hole_v2.h5", )
XDMF
file, which can be used to open the mesh in ParaView.mesh
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7949a2f62590>
points
array on the mesh we just loaded.mesh.points
array([[ 83.52724086, 93.75 , -2.52800323], [ 83.52724086, 112.5 , -2.52800323], [ 83.52724086, 131.25 , -2.52800323], ..., [500. , 300. , 110.4375 ], [500. , 300. , 128.84375 ], [500. , 300. , 147.25 ]])
fig = plt.figure(figsize=(4, 6))
ax1, ax2, ax3 = fig.subplots(3, 1)
ax1.scatter(mesh.points[:, 0], mesh.points[:, 1], s=1)
ax1.set_xlabel("x [m]")
ax1.set_ylabel("y [m]")
ax1.set_title("XY")
ax1.set_aspect("equal")
ax2.scatter(mesh.points[:, 0], mesh.points[:, -1], s=1)
ax2.set_xlabel("x [m]")
ax2.set_ylabel("z [m]")
ax2.set_title("XZ")
ax2.set_aspect("equal")
ax3.scatter(mesh.points[:, 1], mesh.points[:, -1], s=1)
ax3.set_xlabel("Y [m]")
ax3.set_ylabel("z [m]")
ax3.set_title("YZ")
ax3.set_aspect("equal")
fig.suptitle("Mesh points")
fig.set_tight_layout(True)
mesh.points.shape
(4476, 3)
connectivity
array.mesh.connectivity
array([[2842, 2662, 3034, ..., 2663, 3035, 2819], [3010, 2842, 3190, ..., 2843, 3191, 3035], [2866, 2794, 3010, ..., 2795, 3011, 2843], ..., [4286, 4285, 4128, ..., 4302, 4129, 4105], [4287, 4286, 4140, ..., 4303, 4141, 4129], [4288, 4287, 4153, ..., 4304, 4155, 4141]])
mesh.points
are used per element. The amount of elements in
mesh.connectivity
will be smaller than the amount of points defining each
element. Additionally, connected elements share indices from the
mesh.points
.mesh.connectivity
array:mesh.connectivity.shape
(3641, 8)
element_1_point_indices = mesh.connectivity[0, :]
element_1_point_indices
array([2842, 2662, 3034, 2818, 2843, 2663, 3035, 2819])
[2842, 2662, 3034, 2818, 2843, 2663, 3035, 2819]
mesh.points
array:element_1_points = mesh.points[element_1_point_indices]
element_1_points
array([[361.85480312, 281.25 , 8.46142574], [345.94590417, 281.25 , 17.39697932], [378.30937558, 281.25 , 26.90014588], [356.55802305, 281.25 , 34.9196765 ], [361.85480312, 300. , 8.46142574], [345.94590417, 300. , 17.39697932], [378.30937558, 300. , 26.90014588], [356.55802305, 300. , 34.9196765 ]])
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, projection="3d")
# Scatter the nodes of the element
scatter = ax.scatter(
element_1_points[:, 0],
element_1_points[:, 1],
element_1_points[:, 2],
s=50,
c="b",
)
# This whole facet-edge plotting is simply to give an impression of depth.
facet_0 = [0, 1, 3, 2, 0]
facet_1 = [4, 5, 7, 6, 4]
facet_2 = [0, 1, 5, 4, 0]
facet_3 = [2, 3, 7, 6, 2]
# Two SEM facets are omitted as they are not needed
# for the plot.
facets = [facet_0, facet_1, facet_2, facet_3]
lines = [
(element_1_points[facet[iv], ...], element_1_points[facet[iv + 1], ...])
for facet in facets
for iv in range(len(facet) - 1)
]
line_segments = Line3DCollection(lines, colors="r", linewidths=1.0, alpha=0.4)
ax.add_collection3d(line_segments)
ax.set_xlabel("X [m]")
ax.set_ylabel("Y [m]")
ax.set_zlabel("Z [m]")
ax.set_title("Element 1 nodes")
ax.set_box_aspect([1, 1, 1])
ax.view_init(elev=35, azim=-30)
fig.set_tight_layout(True)
mesh.points
array to obtain all points per
element.mesh.points[mesh.connectivity]
array([[[361.85480312, 281.25 , 8.46142574], [345.94590417, 281.25 , 17.39697932], [378.30937558, 281.25 , 26.90014588], ..., [345.94590417, 300. , 17.39697932], [378.30937558, 300. , 26.90014588], [356.55802305, 300. , 34.9196765 ]], [[374.34573945, 281.25 , -1.9986235 ], [361.85480312, 281.25 , 8.46142574], [388.37532567, 281.25 , 7.58267686], ..., [361.85480312, 300. , 8.46142574], [388.37532567, 300. , 7.58267686], [378.30937558, 300. , 26.90014588]], [[365.80311735, 281.25 , -15.51511688], [351.08043863, 281.25 , -7.93131701], [374.34573945, 281.25 , -1.9986235 ], ..., [351.08043863, 300. , -7.93131701], [374.34573945, 300. , -1.9986235 ], [361.85480312, 300. , 8.46142574]], ..., [[500. , 93.75 , 110.4375 ], [500. , 93.75 , 92.03125 ], [481.05034833, 93.75 , 110.22090727], ..., [500. , 112.5 , 92.03125 ], [481.05034833, 112.5 , 110.22090727], [480.92216434, 112.5 , 91.75211413]], [[500. , 93.75 , 128.84375 ], [500. , 93.75 , 110.4375 ], [481.11587253, 93.75 , 128.71455525], ..., [500. , 112.5 , 110.4375 ], [481.11587253, 112.5 , 128.71455525], [481.05034833, 112.5 , 110.22090727]], [[500. , 93.75 , 147.25 ], [500. , 93.75 , 128.84375 ], [481.13207547, 93.75 , 147.25 ], ..., [500. , 112.5 , 128.84375 ], [481.13207547, 112.5 , 147.25 ], [481.11587253, 112.5 , 128.71455525]]])
mesh.points[mesh.connectivity].shape
(3641, 8, 3)
mesh.connectivity.shape
(3641, 8)
point
array, it's range is [0, number_of_points - 1]
(as Python and we at Mondaic use zero based indexing):mesh.connectivity.min(), mesh.connectivity.max(), mesh.points.shape[0]
(0, 4475, 4476)
# The same value, multiple names for convenience
mesh.npoint, mesh.number_of_nodes
(4476, 4476)
# The same value, multiple names for convenience
mesh.nelem, mesh.number_of_elements
(3641, 3641)
mesh.ndim
3
mesh.nodes_per_element
8
# Intimatily related with mesh.nodes_per_element
mesh.shape_order
1
# This is the relationship between shape order and number of elements.
# This is defined by the way GLL integration points work.
assert (mesh.shape_order + 1) ** mesh.ndim == mesh.nodes_per_element
fields
. These fields are NumPy arrays where every entry
is understood to refer to a specific position in the mesh. There are various
types of fields one can add, that correspond to different ways of describing
spatial variation:mesh.points
). The value of a field between points is understood to be
smoothly varying.mesh.get_element_centroid()
). As this
completely spans the domain of a mesh, elemental fields fully define field
values everywhere. Since elements have discrete boundaries, these fields
can be discontinuous in space.mesh.points[mesh.connectivity]
). These
fields are understood to be smoothly varying using the shape functions of
the elements within elements, but are allowed to be discontinuous across
element boundaries. This also means that a single spatial point on element
boundaries can have multiple values for the same field, as every element
allows its own value to be set.mesh.attach_field(name, data)
.
All fields that can be added have distinct shapes. Since the number of
elements and the number of points in a mesh can never be the same (every
element has multiple points, so n_elem < n_points
, always), the
attach_field
method can autoderive which field a user wants to add. Let's
have a look at how to define these fields on a mesh, and attach them.mesh.points
. Below is demonstrated how the
shape of these fields is determined:node_locations = mesh.points
print(f"Shape of node locations array: {node_locations.shape}")
nodal_field_shape = (mesh.points.shape[0],)
print(f"Shape of a nodal field: {nodal_field_shape}")
Shape of node locations array: (4476, 3) Shape of a nodal field: (4476,)
# The last axis is coordinate, so let's compute distance from zero as a nice
# field to visualize
nodal_field = np.linalg.norm(node_locations, axis=-1)
# normalize
nodal_field -= nodal_field.min()
nodal_field /= nodal_field.max()
# Add some variation
nodal_field += np.random.rand(*nodal_field.shape)
# Check the shape is correct
assert nodal_field.shape == nodal_field_shape
mesh.attach_field("nodal_field_A", nodal_field)
mesh
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7949a2f62590>
nodal_fields
dictionary.mesh.nodal_fields
{'nodal_field_A': array([0.40645513, 0.38656187, 0.50480609, ..., 1.06504736, 1.10565353, 1.50676217])}
mesh.nodal_fields
dictionary is a fully
featured Python dictionary, so use any method you are comfortable with, e.g.
dict.pop
.mesh.nodal_fields.pop("nodal_field_A")
array([0.40645513, 0.38656187, 0.50480609, ..., 1.06504736, 1.10565353, 1.50676217])
element_locations = mesh.get_element_centroid()
print(f"Shape of element locations array: {element_locations.shape}")
elemental_field_shape = (mesh.nelem,)
print(f"Shape of a elemental field: {elemental_field_shape}")
Shape of element locations array: (3641, 3) Shape of a elemental field: (3641,)
# The last axis is coordinate, so let's compute distance from zero as a nice
# field to visualize
elemental_field = np.linalg.norm(element_locations, axis=-1)
# normalize
elemental_field -= elemental_field.min()
elemental_field /= elemental_field.max()
# Add some variation
elemental_field += np.random.rand(*elemental_field.shape)
# Check the shape is correct
assert elemental_field.shape == elemental_field_shape
mesh.attach_field("elemental_field_A", elemental_field)
mesh
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7949a2f62590>
UnstructuredMesh
objects, namely at mesh.elemental_fields
:mesh.elemental_fields
{'elemental_field_A': array([1.40685787, 1.22147655, 0.993209 , ..., 1.57602486, 0.90437845, 1.37639816])}
del mesh.elemental_fields["elemental_field_A"]
element_node_locations = mesh.points[mesh.connectivity]
print(f"Shape of element node locations array: {element_node_locations.shape}")
element_nodal_field_shape = mesh.element_nodes_shape
print(f"Shape of a element nodal field: {element_nodal_field_shape}")
Shape of element node locations array: (3641, 8, 3) Shape of a element nodal field: (3641, 8)
nelem, nodes_per_element, 3
) to a scalar
field over those positions (i.e. shape nelem, nodes_per_element
), and
perturb it a little with randomness.# The last axis is coordinate, so let's compute distance from zero as a nice
# field to visualize
element_nodal_field = np.linalg.norm(element_node_locations, axis=-1)
# normalize
element_nodal_field -= element_nodal_field.min()
element_nodal_field /= element_nodal_field.max()
# Add some variation
element_nodal_field += np.random.rand(*element_nodal_field.shape)
# Check the shape is correct
assert element_nodal_field.shape == element_nodal_field_shape
mesh.attach_field("element_nodal_field_A", element_nodal_field)
mesh
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7949a2f62590>
mesh.attach_field("elemental_field_A", elemental_field)
list(mesh.elemental_fields)
['element_nodal_field_A', 'elemental_field_A']
elemental_fields
dictionary. However, they are not of equal shape:{name: field.shape for name, field in mesh.elemental_fields.items()}
{'element_nodal_field_A': (3641, 8), 'elemental_field_A': (3641,)}
element_nodal_fields
:list(mesh.element_nodal_fields)
['element_nodal_field_A']
# This won't work!
mesh.element_nodal_fields.pop("element_nodal_field_A")
print(list(mesh.element_nodal_fields), list(mesh.elemental_fields))
['element_nodal_field_A'] ['element_nodal_field_A', 'elemental_field_A']
elemental_fields
.mesh.elemental_fields.pop("element_nodal_field_A")
print(list(mesh.element_nodal_fields), list(mesh.elemental_fields))
mesh.elemental_fields.pop("elemental_field_A")
print(list(mesh.element_nodal_fields), list(mesh.elemental_fields))
[] ['elemental_field_A'] [] []
def human_readable_size(byte_size):
"""
Convert byte size into a human-readable format (KB, MB, GB, etc.).
"""
for unit in ["bytes", "KB", "MB", "GB", "TB"]:
if byte_size < 1024.0:
return f"{byte_size:.2f} {unit}"
byte_size /= 1024.0
return f"{byte_size:.2f} PB"
print("Memory size of numpy arrays in bytes:")
for name, array in zip(
["nodal\t\t", "elemental\t", "elemental nodal\t"],
[nodal_field, elemental_field, element_nodal_field],
):
print(f"\t{name}\t{human_readable_size(array.itemsize * array.size)}.")
Memory size of numpy arrays in bytes: nodal 34.97 KB. elemental 28.45 KB. elemental nodal 227.56 KB.
# Periodic in X
point_periodic = np.sin(mesh.points[:, 0] / 25)
# Boundary at z coordinate
elemental_discrete = mesh.get_element_centroid()[:, -1] < 0
# To use a nodal field in a elemental nodal field, we need to use the
# connectivity
A = point_periodic[mesh.connectivity]
# To use an elemental field in a elemental nodal field, it needs to be
# broadcast.
B = np.broadcast_to(elemental_discrete[:, None], (mesh.nelem, 8))
total_field = np.zeros(mesh.element_nodes_shape)
# Create a mask based on elements
mask_A = mesh.get_element_centroid()[:, 1] > 200
mask_B = mesh.get_element_centroid()[:, 1] < 150
total_field[mask_A, ...] = A[mask_A, ...]
total_field[mask_B, ...] = B[mask_B, ...]
mesh.attach_field("element_nodal_field_B", total_field)
mesh
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7949a2f62590>
mesh.elemental_fields.pop("element_nodal_field_B")
array([[0.94374256, 0.95553404, 0.54432835, ..., 0.95553404, 0.54432835, 0.99217846], [0.66994425, 0.94374256, 0.17208932, ..., 0.94374256, 0.17208932, 0.54432835], [0.87998878, 0.99559 , 0.66994425, ..., 0.99559 , 0.66994425, 0.94374256], ..., [0. , 0. , 0. , ..., 0. , 0. , 0. ], [0. , 0. , 0. , ..., 0. , 0. , 0. ], [0. , 0. , 0. , ..., 0. , 0. , 0. ]])
mesh.attach_field("elemental_field_B", total_field.mean(axis=-1))
mesh
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7949a2f62590>
del mesh.elemental_fields["elemental_field_B"]
mesh_higher_order = mesh.copy()
# In-place operation
mesh_higher_order.change_tensor_order(4)
mesh_higher_order
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7949a05aa190>
# Compute would-be array sizes (total number of elements)
elemental_nodal_field_size = np.product(mesh_higher_order.element_nodes_shape)
elemental_field_size = int(mesh_higher_order.nelem)
nodal_field_size = int(mesh_higher_order.npoint)
print(
f"Dtype used: {nodal_field.dtype}, dtype size: {nodal_field.itemsize}.\n"
)
itemsize = int(nodal_field.itemsize)
print("Memory size of numpy arrays in bytes:")
for name, array_size in zip(
["nodal\t\t", "elemental\t", "elemental nodal\t"],
[nodal_field_size, elemental_field_size, elemental_nodal_field_size],
):
print(f"\t{name}\t{human_readable_size(itemsize * array_size)}.")
Dtype used: float64, dtype size: 8. Memory size of numpy arrays in bytes: nodal 1.88 MB. elemental 28.45 KB. elemental nodal 3.47 MB.