This tutorial is not for the latest stable Salvus version (viewing 2025.1.3).

Salvus’ UnstructuredMesh - Part 2 - Averaging and integrating mesh parameters

In 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 npt

1. Creating a small higher-order mesh

We’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

<unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7143514d2810>

2. Defining an elemental nodal field

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

<unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7143514d2810>

3. Naive arithmetic mean

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]

4. Analytical reference for one element

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.

5. Mass-matrix (quadrature) weighted averaging

Spectral elements use Gauss–Lobatto–Legendre (GLL) quadrature. For each local node $i$ in element $e$ we have a weight $w_i$ from the quadrature. Salvus exposes this weight through the mass matrix $mm$ , which has shape (nelem, nodes_per_element):

A consistent average per element is then:

\text{average}_e = \frac{\sum_i \mathrm{mm}_{e,i} \cdot \mathrm{value}_{e,i}}{\sum_i \mathrm{mm}_{e,i}}

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

<unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7143514d2810>

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.

6. Utility function

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])