Mondaic DocumentationIn Part 1 we explored how points, connectivity and the various field types
(nodal, elemental, elemental nodal) are represented on Salvus meshes.
This second part focuses on a very common practical task:
Converting an elemental nodal field (ENF) (values on every GLL point per
element) to a single representative elemental value (EF) per element.
Walking through this process helps to understand the way the spectral element
method computes integrals, and how to correctly handle fields on meshes.
We will:
- Construct a small 2D spectral element mesh.
- Define a smoothly varying field as an elemental nodal field.
- Compute the mean per element in various ways.
- Show why the naive approach can be wrong.
- Show how to do it correctly using the mass matrix.
Just like in Part 1 the emphasis is on understanding shapes and assumptions.
import salvus.namespace as sn
import sympy as sp
import numpy as np
import numpy.typing as nptWe'll use the simple homogeneous Cartesian mesh generation helper. A higher
polynomial order (here 7) ensures non-uniform GLL node spacing inside each
element, so that a naive mean is visibly different from a quadrature-based
average. This might not be as obvious for zero-order (linear) or first-order
(quadratic) elements.
mesh = sn.simple_mesh.CartesianHomogeneousAcoustic2D(
vp=1, rho=1, x_max=1, y_max=1, max_frequency=1
).create_mesh()
mesh.change_tensor_order(7)
mesh
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x78b317576450>
We create a scalar field that varies only with x. Evaluating it at the global
point coordinates gives a nodal array. Indexing with
mesh.connectivity
reshapes it to an elemental nodal field: (nelem, nodes_per_element).points = mesh.points
def v(x, y):
# y is unused
return np.sin(2 * np.pi * x)
nodal_velocity = v(*points.T)
elemental_nodal_velocity = nodal_velocity[mesh.connectivity]
mesh.elemental_fields.clear()
mesh.attach_field("VP", elemental_nodal_velocity)
mesh
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x78b317576450>
First attempt: take the plain mean over the nodal axis. Implicit assumptions:
- GLL nodes are uniformly spaced (they are not).
- All quadrature weights are equal (they are not).
- Geometry mapping is uniform (can fail for curved / deformed elements).
For higher order or deformed elements this mean introduces bias.
naive_average = elemental_nodal_velocity.mean(axis=-1)
print("Naive average (all elements):", naive_average)Naive average (all elements): [ 0.43620842 -0.43620842 0.43620842 -0.43620842]
To anchor the results, we can evaluate the element-average analytically for
the first element. Since we have a function of only x and uniform mesh
elements, this is doable.
x = sp.Symbol("x")
v_sp = sp.sin(2 * sp.pi * x)
element_length = 1 / (mesh.nelem**0.5)
area = element_length**2
average_analytical = (
sp.integrate(v_sp, (x, 0, element_length)) * element_length
) / area
print(f"Analytical average (element 0): {float(average_analytical):.6f}")
print(
f"Element 0 naive vs analytical: {naive_average[0]:.6f} vs "
f"{float(average_analytical):.6f}"
)Analytical average (element 0): 0.636620 Element 0 naive vs analytical: 0.436208 vs 0.636620
The mismatch already illustrates that the plain mean does not correspond
to integrating the field over the element and dividing by the element area.
Spectral elements use Gauss–Lobatto–Legendre (GLL) quadrature. For each local
node in element we have a weight from the quadrature. Salvus
exposes this weight through the mass matrix , which has shape
(nelem, nodes_per_element):A consistent average per element is then:
This collapses to a simple mean only if all
mm weights are identical.Let's first have a look at the mass matrix:
mm = mesh.get_mass_matrix()
mesh.elemental_fields.clear()
mesh.attach_field("mass_matrix", mm)
mesh
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x78b317576450>
We can see how for seventh-order elements the GLL nodes weights are larger
towards the center of the element, and smaller towards the edges and
vertices. This has a strong influence on how the individual elemental nodal
values contribute to the average.
print(
"Mass matrix shape:",
mm.shape,
"| ENF shape:",
elemental_nodal_velocity.shape,
)Mass matrix shape: (4, 64) | ENF shape: (4, 64)
Now we can compute the quadrature weighted average.
VP_elemental = (mm * elemental_nodal_velocity).sum(axis=-1) / mm.sum(axis=-1)
print("Quadrature (SEM) average (all elements):", VP_elemental)
print(
f"Element 0 weighted vs analytical: {VP_elemental[0]:.6f} vs {float(average_analytical):.6f}"
)Quadrature (SEM) average (all elements): [ 0.63661977 -0.63661977 0.63661977 -0.63661977] Element 0 weighted vs analytical: 0.636620 vs 0.636620
Great! The weighted average matches the analytical result very closely.
Below is a small helper mirroring the pattern above for reuse. It
deliberately pulls the mass matrix inside (cheap vs clarity) and performs the
weighted reduction.
def average_enf_to_ef(
enf: npt.NDArray, mesh: sn.UnstructuredMesh
) -> npt.NDArray:
"""
Convert an elemental nodal field (Part 2) into an elemental field using
the lumped mass matrix weights.
"""
mm = mesh.get_mass_matrix()
return (mm * enf).sum(axis=-1) / mm.sum(axis=-1)
average_enf_to_ef(elemental_nodal_velocity, mesh)array([ 0.63661977, -0.63661977, 0.63661977, -0.63661977])
PAGE CONTENTS