pyCAFE: A Finite Element Framework for Solving Acoustic Problems in Python

Geometry and material selection

pyCAFE operates on externally generated two-dimensional meshes and includes auxiliary scripts that define the expected workflow for geometry creation, mesh generation and material parameter assignment. The geometry module currently supports the interactive definition of simple benchmark domains, namely rectangular and circular cavities, through keyboard input. The user specifies the geometrical parameters, the maximum frequency of interest and the desired element order, which can be either linear or quadratic.

Mesh generation is performed through the Gmsh Python API. In addition to defining the geometrical entities, the script assigns physical groups to the acoustic domain and to each boundary segment. To support reproducibility, the generated mesh is saved in a user-selected location through a file dialog. A companion mesh visualisation utility can then be used to load an existing .msh file, extract the node coordinates, element connectivities and physical boundary groups, and generate a two-dimensional plot for preliminary inspection.

The material properties of the acoustic medium are defined independently from the mesh. A dedicated material utility includes predefined parameter sets for common fluids, such as air and water, and also allows custom media to be defined by the user. The characteristic impedance,

is computed from the selected density and sound speed and returned together with these parameters, so that the same material definition is used consistently in the subsequent analysis steps.

Build matrices

The build matrices folder contains the numerical routines used to construct the acoustic finite element operators. This module implements the element formulations required to assemble the global stiffness, mass and damping matrices from the discretised mesh. Two quadrilateral acoustic elements are currently supported: the four-node element, referred to as CQUAD4, which uses a $2\times 2$ Gauss integration rule, and the eight-node quadratic element, referred to as CQUAD8, which uses a $3\times 3$ Gauss integration rule (20, 21).

Each element formulation is implemented in a dedicated module that contains the corresponding shape functions, Jacobian computation and element matrix construction. The element type is identified automatically from the imported mesh, and the appropriate formulation is then used during assembly. This avoids manual changes when switching between linear and quadratic discretisation. The separation between element formulations also keeps the implementation of each element self-contained and suitable for independent testing.

For regular acoustic fields in simple benchmark domains, CQUAD4 and CQUAD8 elements are expected to show convergence rates of $O(h^2)$ and $O(h^4)$ (22), respectively, according to standard finite element approximation theory (3, 4).

Boundary condition assignment

Boundaries are detected in the mesh, and the corresponding acoustic boundary condition is selected by the user. Multiple boundary condition types can therefore be assigned within the same model, according to the physical groups defined in the mesh. The assignment procedure is separated from both the solver and the matrix assembly routines, so that boundary conditions can be changed without modifying the numerical core.

The boundary condition module currently supports rigid, or so called hard wall boundaries, zero pressure boundaries, constant pressure boundaries, impedance boundaries and prescribed normal velocity boundaries. Rigid boundaries represent the default condition and impose zero normal particle velocity. Zero pressure and constant pressure boundaries prescribe the acoustic pressure along selected boundaries. Impedance boundaries relate pressure to normal particle velocity through a complex acoustic impedance, while prescribed normal velocity boundaries are used to represent vibrating surfaces or velocity-driven acoustic excitation. The treatment of boundary conditions distinguishes between two mathematically distinct types: natural or Neumann conditions and essential or Dirichlet conditions. Rigid walls, prescribed normal velocity and impedance boundaries are treated as natural boundary conditions, since they arise directly from the weak form of the Helmholtz equation. Rigid walls correspond to the homogeneous Neumann condition

In Eq. (2), $p$ is the acoustic pressure and $n$ is the outward normal direction to the boundary. This condition imposes zero normal pressure gradient and corresponds to zero normal particle velocity at a rigid wall.

A prescribed normal velocity $v_n$ contributes to the right-hand side vector through the boundary load term

In Eq. (3), ${\mathbf{\text{f}}}_v$ is the load vector contribution associated with the prescribed normal velocity, i is the imaginary unit, $\omega$ is the angular frequency, ${\rho }_0$ is the density of the acoustic medium, $v_n$ is the prescribed normal particle velocity, ${\mathrm{\Gamma }}_v$ is the boundary where the velocity condition is applied, $\mathbf{\text{N}}$ is the vector of boundary shape functions and $\text{d}s$ is the differential boundary length.

Impedance boundaries add an admittance weighted boundary contribution to the system matrix. The corresponding element contribution can be written as

and enters the dynamic stiffness matrix through the term

In Eqs.(4) and (5), ${\mathbf{\text{C}}}_{Z,\mathit{\text{ij}}}$ is the impedance boundary matrix contribution associated with the local shape functions $N_i$ and $N_j$, ${\mathrm{\Gamma }}_Z$ is the impedance boundary, and $Y$ is the specific acoustic admittance. The specific acoustic admittance is defined as

In Eq.(6), $\zeta$ is the normalised acoustic impedance, ${\rho }_0$ is the density of the acoustic medium and $c_0$ is the speed of sound.

Prescribed pressure boundaries are treated as essential boundary conditions and must be imposed explicitly on the algebraic system. In pyCAFE, zero pressure and non zero constant pressure boundaries are handled separately. Zero pressure boundaries are removed from the global system during assembly:

The rows and columns associated with these degrees of freedom are removed from $\mathbf{\text{K}}$, $\mathbf{\text{M}}$ and $\mathbf{\text{C}}$, yielding the reduced operators

Non-zero pressure values are prescribed at selected degrees of freedom (DOF) and enforced through a partitioned formulation (23). The reduced system is split into known degrees of freedom, corresponding to the prescribed pressure values,

and unknown degrees of freedom. The resulting linear system is

In Eq.(10), ${\mathbf{\text{p}}}_i$ is the vector of unknown pressure degrees of freedom, ${\mathbf{\text{p}}}_0$ is the vector of prescribed pressure values, ${\mathbf{\text{f}}}_i$ is the corresponding reduced load vector, ${\mathbf{\text{A}}}_{\mathit{\text{ii}}}$ is the block of the dynamic stiffness matrix associated with the unknown degrees of freedom, and ${\mathbf{\text{A}}}_{\mathit{\text{in}}}$ couples the unknown degrees of freedom with the prescribed pressure degrees of freedom.

The reduced dynamic stiffness matrix is defined as

In Eq.(11), ${\mathbf{\text{K}}}_{\text{red}}$, ${\mathbf{\text{M}}}_{\text{red}}$ and ${\mathbf{\text{C}}}_{\text{red}}$ are the reduced stiffness, mass and damping matrices, respectively.

This treatment follows the standard static condensation procedure commonly used in matrix-based finite element formulations to impose prescribed degrees of freedom. Similar partitioning strategies are adopted in many finite element codes, where known degrees of freedom are separated from the unknowns before solving the reduced linear system.

In addition to boundary conditions defined along mesh boundaries, pyCAFE supports localised acoustic excitation through a point pressure source. The source is defined by its spatial coordinates, and the nearest mesh node is automatically identified. A pressure amplitude is then assigned to this node and included in the acoustic load vector. All boundary condition data are stored in a structured container, which is subsequently passed to the matrix assembly and solver routines.

Solver and post-processing

Once the global acoustic matrices and boundary contributions have been assembled, pyCAFE provides direct and modal solution strategies in the frequency domain. In the direct strategy, the acoustic response is computed independently at each frequency of interest. The final algebraic system is expressed through the dynamic stiffness matrix

where $\mathbf{\text{K}}$, $\mathbf{\text{M}}$ and $\mathbf{\text{C}}$ denote the assembled stiffness, mass and damping matrices after the application of the relevant boundary contributions. At each angular frequency $\omega$, the complex pressure vector is obtained by solving

Sparse direct solvers are employed when available (24). This solution strategy is suitable when the full forced response is required over a prescribed frequency range.

In addition to the direct solution strategy, pyCAFE includes a modal solver for acoustic analysis. The solver computes natural frequencies and corresponding acoustic mode shapes by solving the generalised symmetric eigenvalue problem

The eigenproblem is solved using ARPACK through scipy.sparse.linalg.eigsh in shift invert mode. This formulation allows only the first $k$ eigenvalues and eigenvectors to be requested, which is appropriate for low- and mid-frequency studies where only a limited number of modes is required. The modal solver can be used for standalone modal studies or as a basis for reduced order models (25).

Post-processing utilities are provided to facilitate result interpretation. The acoustic pressure field can be reconstructed and visualised over the entire computational domain, while pressure values can also be extracted at specific spatial locations to obtain pointwise frequency response functions. Both amplitude and phase spectra are available and can be exported as plain text files containing the real part, imaginary part, modulus and phase of the complex pressure at the selected point.

An animated visualisation of the pressure field evolution across the analysed frequency range can also be generated as a video file. For modal analyses, individual acoustic mode shapes are reconstructed on the full mesh and displayed using a diverging colour map, which highlights nodal lines and pressure extrema. Finally, all field results, including frequency domain pressure distributions and acoustic mode shapes, can be exported in VTK Unstructured Grid (VTK = Visualization Toolkit) format together with a ParaView Data PVD collection file. This enables direct import into ParaView for interactive inspection and frequency by frequency visualisation.

Automated test suite

The software is tested using pytest. The test suite includes approximately 70 tests organised into four categories, each implemented in a dedicated module:

• Element matrices (test_element_matrices.py): Unit tests verifying symmetry of K_e and M_e; the null-space condition $K_e𝟏=𝟎$ (constant pressure is a free mode under Neumann boundary conditions); mass conservation ($\sum M_e=A_e/c_0^2$, where A_e is element area); partition of unity for shape functions (${\sum }_i{N}_i=1$); and the Kronecker-delta interpolation property at element nodes. Tests are implemented for both CQUAD4 and CQUAD8 elements.

• Boundary conditions (test_boundary_conditions.py): Tests verifying correct Dirichlet elimination via the partitioned system, correct identification of free-Degree of Freedom (DOF) index sets, and preservation of matrix symmetry after boundary condition application.

• Modal solver (test_solver_modal.py): Validation of natural frequencies against the analytical closed- form solution for a rigid rectangular cavity, with mesh-dependent tolerances. Tests also verify monotone convergence of eigenvalues under mesh refinement and $\mathbf{\text{M}}$-orthogonality of the computed mode shapes.

• Helmholtz solver (test_solver_helmholtz.py): Verification of the direct frequency domain solver on a reference problem with known pressure distribution.

Rectangular rigid-wall cavity validation

The test domain is a two-dimensional rectangular acoustic cavity,

with dimensions $L_x=1.0 \text{m}$ and $L_y=0.5 \text{m}$. The cavity is filled with air at $20\mathrm{°}\text{C}$. All four boundaries are modelled as rigid walls, corresponding to the homogeneous Neumann condition

This condition is naturally satisfied by the weak formulation and therefore does not require any modification of the assembled system matrices.

The discrete acoustic eigenvalue problem is written as

where $\mathbf{\text{K}}$ and $\mathbf{\text{M}}$ are the acoustic stiffness and mass matrices assembled element by element from the quadrilateral interpolation functions.

For a rigid rectangular cavity, the analytical eigenfrequencies are given by (26):

A mesh size of $h=0.070 \text{m}$ is used for the eigenfrequency comparison and for the histogram representation. This value lies in the middle of the convergence range, which will be discussed in the next section and provides a good compromise between numerical accuracy and computational cost. At this resolution, CQUAD8 already achieves a mean relative error below $0.016\%$, whereas CQUAD4 requires a finer mesh to reach a comparable accuracy level.

Table 1 reports the first 10 eigenfrequencies computed by pyCAFE using CQUAD4 and CQUAD8 elements, and compares them with the analytical solution. The balanced mesh size is $h=0.070 \text{m}$, corresponding to 120 nodes for CQUAD4 and 337 nodes for CQUAD8. The relative error is defined as

Figure 2 presents the same comparison in graphical form. The upper panel compares the absolute eigenfrequencies, while the lower panel reports the relative error with respect to the analytical solution.

Figure 2

hp-refinement convergence study

The convergence behaviour is assessed through an $h$-refinement study. For a polynomial element of order $p$, the eigenfrequency error is expected to scale as

Therefore, CQUAD4 elements, based on linear interpolation, are expected to show an $\mathcal{O}(h^2)$ convergence rate, while CQUAD8 elements, based on quadratic interpolation, are expected to show an $\mathcal{O}(h^4)$ convergence rate (27).

Table 2 reports the mean relative error over the first 10 modes for a sequence of mesh sizes. The balanced configuration, $h=0.070 \text{m}$, is the mesh used for the eigenfrequency comparisons reported in the previous section. Figure 3 shows the corresponding convergence curves, together with the theoretical reference slopes and the solver time associated with each mesh.

Figure 3

Modal analysis with hard walls

Table 3 compares the first 10 eigenfrequencies. For the balanced configuration, pyCAFE and FEniCSx achieve the same mean relative error to four significant digits, namely 0.0163%, in both cases. However, pyCAFE reaches this accuracy with 337 pressure unknowns instead of 435 for FEniCSx, corresponding to approximately 22.5% fewer DOFs. In addition, the reduced eigensolve stage is faster in pyCAFE, with a computational time of 0.004 s compared with 0.018 s for FEniCSx. This corresponds to a speed ratio of 4.48 in favour of pyCAFE for this test case. Figure 4 summarises the corresponding normalised error, DOF count, and eigensolve time metrics.

Figure 4

To complement the eigenfrequency comparison, Figure 5 compares the corresponding pressure mode shapes for modes 1, 4, 6 and 7 computed with the two solvers. Each field is normalised to the interval $[–1,1]$, so that the visual comparison focuses on the nodal pattern and phase distribution rather than on the absolute scaling.

Figure 5

Modal analysis with zero pressure walls

As a second modal benchmark, all four walls are assigned a homogeneous Dirichlet condition, $p=0$, turning the domain into a pressure-release cavity. In pyCAFE, this boundary condition is handled by the built-in Guyan reduction (23), where the boundary nodes are eliminated from the system before the eigensolve. In FEniCSx, the same reduction is applied by extracting the free degrees of freedom before calling the sparse eigenvalue solver.

Table 4 reports the first 10 eigenfrequencies. Both solvers achieve very low errors, with a mean error of 0.0502% for pyCAFE and 0.0495% for FEniCSx. pyCAFE requires 257 reduced DOFs, while FEniCSx requires 351 free DOFs. In this pressure-release configuration, the reduced eigensolve stage is slightly slower in pyCAFE than in FEniCSx. However, the accuracy remains essentially the same for the two solvers, confirming that the Dirichlet reduction adopted in pyCAFE preserves the expected modal behaviour. Figure 6 reports the corresponding summary chart, following the same format used for the rigid-wall case. To further assess the spatial structure of the numerical solution, Figure 7 also compares selected pressure mode shapes obtained with the two solvers.

Figure 6

Figure 7

To further support reproducibility, additional direct-frequency-response comparisons between pyCAFE and FEniCSx are provided in the examples and validation folders of the GitHub repository. These examples are intended for users interested in inspecting the direct Helmholtz sweep cases and the corresponding boundary condition implementations.

Comparison with commercial software

To implement the pyCAFE validation, the same cavity is also compared against COMSOL Multiphysics 6.2. Two modal datasets are considered: rigid walls, corresponding to a hard-wall Neumann condition, and zero pressure walls on all four sides. As reported in Tables 5 and 6, pyCAFE provides slightly lower frequency errors than COMSOL for both boundary condition configurations. This result confirms that the matrix-based implementation adopted in pyCAFE can reproduce the analytical modal solution with accuracy at least comparable to, and in this benchmark slightly higher than, that obtained with a commercial finite element package.