Version:

Salvus' UnstructuredMesh - Part 1 - Points, connectivity and fields

This tutorial focuses on the data structure of meshes in Salvus -- the UnstructuredMesh objects. It will discuss the internal workings of the UnstructuredMesh class, but it will not focus on creating meshes.
By the end of this tutorial, you will:
  • Understand the basic components of a Salvus unstructured mesh, including points, connectivity, and fields.
  • Understand the distinction between different types of fields (nodal, elemental, and elemental nodal) and how to add or remove them.
Let's look what's in a Salvus mesh!
Copy
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Line3DCollection

import salvus.namespace as sn

Prelude: on spectral elements and GLL shape functions

Before diving into Salvus unstructured meshes, it’s important to understand the role of spectral element modelling in relation to meshes. The spectral element method (SEM) is a high-order finite element method used to solve partial differential equations on unstructured meshes.
The SEM method uses meshes based on finite elements, with given Gauss-Lobatto-Legendre (GLL) points. These GLL points are specific quadrature points used within each element for numerical integration and interpolation. They allow accurate approximation of physical fields within the element. When we speak of a first order (or higher order) mesh, we refer to how many GLL points are used per element. or instance, a first-order mesh uses two GLL points per element edge, first-order (linear) polynomials between the vertices. Higher-order meshes use more, enabling greater accuracy at the cost of increased computational complexity (and higher memory usage).
Let's now explore how Salvus represents these points, elements, and fields defined over them.
In this tutorial, it is assumed you already have a mesh somehow. This can be from an external source (like an 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.
Let's load the mesh from an .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")
Of course, the existence of 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",
)
This will write the mesh again to an HDF5 file. Additionally, it will also write an XDMF file, which can be used to open the mesh in ParaView.
We can visualize this mesh by simply leaving it at the end of a notebook cell without any other statements. This will (for any Python object) try to display it. Salvus meshes have a built-in display method for notebooks that are interactive widgets.
mesh
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7949a2f62590>
That is one boring mesh! Because, for the purposes of this tutorial, the file on disk has no fields attached, the mesh widget is actually not that interesting -- there is no parameter to visualize.
If you are running this tutorial on your own machine (and not reading the website), try click-and-dragging in the widget to rotate the mesh, or zooming to view details. There is a hole in it!
The basic definition of a Spectral Element Method mesh is points and connectivity. In fact, these are all you need to specify a Salvus mesh from scratch.
The points of a mesh describe all the points that make up the elements, and the connectivity describes which points are used for which elements.
Let's investigate the 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      ]])
These are the entries for the different point (node, vertex) locations in our mesh, but they are still unassociated with elements. We can now directly take these points and scatter plot them in 2D:
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)
Finally, we can look at the shape of this object to see how many points are in total in the mesh.
mesh.points.shape
(4476, 3)
This means that our mesh consists of 4476 points in 3 dimensions.
Let's now have a look at the 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]])
What the entries of the connectivity array refer to is which indices from 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.
In 3D, Salvus uses hexahedrals, i.e., deformed 6 sided cubes. For first order shapes (i.e. no curvature), that means one needs 8 nodes to fully define an element. This can be seen directly from the shape of the mesh.connectivity array:
mesh.connectivity.shape
(3641, 8)
This shows that this mesh has 3641 elements, each defined by 8 points, or nodes.
Now, let's have a look at what these entries exactly do.
element_1_point_indices = mesh.connectivity[0, :]

element_1_point_indices
array([2842, 2662, 3034, 2818, 2843, 2663, 3035, 2819])
In other words, the very first element of our mesh, consists of the points
[2842, 2662, 3034, 2818, 2843, 2663, 3035, 2819]
These indices refer to the 3D coordinates of the points in the 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 ]])
Let's look at these points.
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)
We can reshape the entire 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]]])
That is a bit overwhelming, so instead let's look at its shape.
mesh.points[mesh.connectivity].shape
(3641, 8, 3)
In this case, for every element (of the 3641 elements), we now have 8 points (more if using higher order shapes), with each 3 coordinates.
The shape of the connectivity array is unrelated to the shape of the point array.
Many points can be used to define few elements, or the other way around. It is even possible to pass on more points to Salvus meshes than are used in the connectivity array -- these simply won't be used to construct elements!
To understand this better, let's again look at the connectivity shape.
mesh.connectivity.shape
(3641, 8)
As this array indexes into the 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)
That matches! We have 4476 points, and the maximum index is exactly 4475 (note again we use zero-based indexing).
Finally, there are a few attributes on meshes that are a lot faster than constantly getting the shapes of points and connectivity:
# 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
Salvus meshes allow us to define quantities that vary with location in data structures we call 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:
  • Nodal fields. These are fields defined directly on the points of a mesh (mesh.points). The value of a field between points is understood to be smoothly varying.
  • Elemental fields. These are fields defined per element of the mesh (around the centroids found by 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.
  • Elemental nodal fields. These are fields defined on every GLL point of elements in the mesh (found by 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.
Fields are added to mesh using the method 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.

Nodal Fields

Nodal fields are fields that are defined once per point (node, vertex) in the mesh, before those points are grouped into elements. As such, they are defined exactly once per row of 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,)
Creating a nodal field can be easily demonstrate by applying some scalar (specifically, vector to scalar) function to the coordinates of the points. By adding some randomness to the field as well, one can get a good feeling how spatial variations manifest in nodal fields:
# 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>
This fields shows smoothly varying transitions between the vertices, indicating that nodal fields are appropriate for storing smoothly varying data.
Since it is a nodal field, it can be found in the nodal_fields dictionary.
mesh.nodal_fields
{'nodal_field_A': array([0.40645513, 0.38656187, 0.50480609, ..., 1.06504736, 1.10565353,
        1.50676217])}
Let's remove thise field again. The 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])
Elemental fields are fields that are constant value over an element, and thus require one entry per element:
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,)
As elemental fields are constant over an element, they are lower resolution than nodal fields.
In the same way, we can use the locations array to define some method that maps it to a scalar value, and add it to the mesh. Because this field is now a different shape, Salvus automatically detects what type of field you are adding.
# 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>
The nodal field as visualized shows many discontinuities, indicating that elemental fields are very appropriate when dealing with discontinuous fields.
Elemental fields are stored in a different location on 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])}
Let's clean up again, now using another approach to delete entries from dictionaries.
del mesh.elemental_fields["elemental_field_A"]
Finally, we can define fields per node of elements. That have the same shape as the first two dimensions of the reconstructed points array for elements, as was used earlier:
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)
These are fields that allow the largest flexibility in Salvus: between nodes they are smoothly varying, while between elements they are allowed to have discontinuities. Discontinuities are possible because neighboring elements are allowed to have different field values for the same point.
Let's again create a function, this time one that transforms the node locations per element (i.e. shape 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>
One can see above that the elemental nodal fields are very flexible: they are smoothly varying, and discrete varying.
It is important to make a clear distinction between elemental fields and element nodal fields. Let's see what happens when both are present on a mesh.
mesh.attach_field("elemental_field_A", elemental_field)

list(mesh.elemental_fields)
['element_nodal_field_A', 'elemental_field_A']
Salvus defines elemental fields as any element that varies per element. Thus, strictly speaking, element nodal fields are a subset of elemental fields. Regardless of what we call both of these fields, both are stored in the 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,)}
Salvus meshes provide a convenient way to query only those elemental fields that are defined per node, via element_nodal_fields:
list(mesh.element_nodal_fields)
['element_nodal_field_A']
It is important, however, to understand that this is a read-only property. Modifying it doesn't write to the original elemental fields.
# 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']
To modify either of these fields, one has to go through 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']
[] []
Finally, let's have a look at what it means in terms of memory usage when defining fields of different types.
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.
Although nodal and element fields are quite similar in size, elemental nodal fields are almost an order of magnitude larger, for first order meshes. This is further compounded when using higher order meshes.
Sometimes, one might want to combine different kinds of fields. Typically, it is easiest to upcast something to element nodal fields, as this, coming from elemental fields or point fields, adds degrees of freedom.
# 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.        ]])
We can relatively easily average the element nodal fields back down to elemental fields as well, using NumPy's mean. Note, however, that for strongly deformed or strongly spatially varying fields, the simple mean of these fields over the nodal values will not be an accurate spatial average: for that one would need to account for the GLL integration weights.
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"]
To go from a first order mesh to a higher order mesh allows one to add more spatial detail to a mesh, as well as to more smoothly conform to curved surfaces. Ideally, these higher order meshes are created by whatever meshing routine creates the mesh in the first place, but any lower order mesh can be raised to a higher order (in-place) as shown below. Note however, that this linearly interpolates the GLL locations, thus possibly forsaking some of the benefits higher order meshes would've brought in the first place.
Please be aware that using this method will delete any potentially set element nodal fields and nodal fields. Elemental fields and side sets will be retained.
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>
And of course, using higher order meshes also significantly raises the memory usage of these fields, both on disk and in-memory.
# 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.
Those (nodal, elemental nodal) are some large fields!
PAGE CONTENTS