NHSJS Reports

Surrogate Model for Prediction of Displacement Field of Violin Plate Driven by Material Parameters and Point Cloud

November 5, 2025

119

Abstract

Traditional physics-based modeling of instruments often relies on the Finite Element Method (FEM), which is known for its accuracy. However, it is computationally expensive and not suitable for real-time prediction or flexible applications. Recently, neural operators have proven effective for solving partial differential equations (PDEs). However, their application in musical instrument acoustics remains largely unexplored. This research proposes a neural network-based surrogate model for predicting the full-field displacement distribution and eigenfrequencies of violin plates. Unlike traditional surrogate models, the proposed method does not require high-end GPUs and exhibits strong generalization capabilities and computational efficiency. Its lightweight structure and deployment flexibility make it particularly well-suited for acoustic evaluation during the instrument-making process.

Keywords: Surrogate Model, Neural Operator, Violin Acoustic Analysis

Introduction

Acoustic modeling of musical instruments traditionally relies on the finite element method for solving complex partial differential equations. Such methods have extremely high accuracy, but are often accompanied by considerable computational cost. In order to improve computational efficiency, surrogate models have been studied as an alternatives to the traditional finite element methods.

In recent years, deep learning applied to physics modeling, extending to Physics-Informed Neural Networks (PINNs), a neural operator, and mesh-based surrogate models, has been successfully applied in solid mechanics, fluid mechanics, and aerodynamics.

However, as far as currently known, the application of such surrogate models to instrument acoustic simulation remains largely unexplored. Some existing surrogate models can handle irregular 3D data and complex geometries, but they typically require region-based inputs and do not flexibly support fixed geometries with varying material parameters.

In contrast, the proposed model is designed for higher flexibility: it takes as input the geometry data together with global material parameters, and predicts the displacement at each point as well as the eigenfrequency of the plate. The model can operate directly on point cloud data or incorporate mesh connectivity information. The model provides a lightweight and efficient alternative, achieving a speedup of approximately 9.4x compared to FEM simulations.

Physical of Violin Making

Many luthiers continue to rely on traditional craftsmanship and subjective judgment rather than physics-based models or experimental data. Although modern luthiers often attempt to replicate the shapes of Old Italian violins using historical data, the acoustic outcome is affected by many variable physical factors, such as wood density, aging effects, and material properties. As a result, copying the geometry and thickness of a historical instrument is insufficient to reproduce its timbre. To thoroughly understand and approximate the timbre of a famous musical instrument, a physics-based model and scientific evaluation approach are essential.

One widely adopted scientific approach is the analysis of vibration modes and eigenfrequencies of violin plates. Eigenfrequencies the natural frequencies at which an object vibrates freely are influenced by the mass, shape, and material properties of the plate. In practice, luthiers often measure the eigenfrequencies of the top and back plates before assembly, using them as indicators of acoustic quality. Curtin, J¹. has recorded the vibration modes, eigenfrequencies, and weights of numerous Old Italian violins, providing valuable benchmarks for modern luthiers. Hutchins, C. M². observed that specific frequency ranges of certain vibrational modes tend to be favored by performers.

By adjusting the thickness and weight of the plates, luthiers can shift the eigenfrequencies into preferred ranges, thus enhancing the tonal qualities of the finished instrument. Specifically, luthiers hold the violin plate at a designated location with one hand and tap a selected area with the other to identify the corresponding modal frequency based on the produced sound.

The choice of holding position and tapping point is guided by Chladni patterns—nodal shapes formed in regions of minimal vibration during acoustic resonance.

These patterns are visualized by sprinkling fine sand onto the plate and applying vibrations; the sand accumulates along the nodal lines where the amplitude is nearly zero. Luthiers typically hold the plate at these nodal points and tap elsewhere to excite the corresponding vibrational mode and observe its characteristic frequency.

In practice, achieving the desired eigenfrequencies typically involves an iterative process of measurement, thickness adjustment, and re-measurement. However, the eigenfrequency of a violin plate is influenced by more than just thickness. Factors such as material properties (e.g., density, Young’s modulus, poisson ratio), plate outline, and arching geometry also play critical roles. Therefore, the primary contribution of this research is the development of a surrogate model capable of predicting the eigenfrequencies of a violin plate in real time based on its material and geometric parameters. This model offers a assist luthiers in achieving target eigenfrequencies more efficiently and accurately.

Moreover, while the eigenfrequencies of isolated plates are correlated with the resonant behavior of the assembled instrument, the final acoustic response also depends on complex structural interactions and coupling effects. To address this, the proposed model is extended beyond single frequency prediction: it also provides visualizations of pressure distributions and full-field displacement across the plate, thereby offering a more comprehensive foundation for structural-acoustic analysis.

Methodology

When parametric partial differential equation modeling, in general:

(1) $\begin{eqnarray*}\mathrm{R}(u, a)=0, \quad \text{In } \Omega \subset \mathbb{R}^d.\end{eqnarray*}$

Here, $R$ denotes a differential operator, $a: \mathbb{R}^{d} \rightarrow \mathbb{R}^{d_a}$ represents the input parameter function (e.g., material properties),

and $u: \mathbb{R}^{d} \rightarrow \mathbb{R}^{d_u}$ is the solution of the PDE, i.e., the target physical field to be predicted (such as displacement or pressure). $\Omega$ denotes the domain, which in this context refers to the geometric region of the violin plate. The input parameters $a$ include material properties: density $\rho$ , Young’s modulus $E$ , and Poisson ratio $\nu$ . The output field $u$ corresponds to the three-dimensional displacement field $(u, v, w)$ on the violin plate. Additionally, the global physical quantity of interest is the eigenfrequency.

To build a surrogate model for simulating the acoustic response of a violin, first define a solution mapping function $g^+$ that maps material parameters and geometric information to the displacement field in the solution space:

(2) $\begin{eqnarray*}g^+:(a, \Omega) \mapsto u.\end{eqnarray*}$

Here, $a = [\rho, E, \nu] \in \mathbb{R}^{du}$ denotes the global material parameter vector defined over the entire structure, $\Omega$ represents the geometric domain of the violin plate, and $u(x): \Omega \rightarrow \mathbb{R}^{3}$ is the three-dimensional displacement field $(u, v, w)$ defined on $\Omega$ . This mapping $g^+$ represents a correspondence from global features and geometry to a field solution. Due to the high geometric complexity of $\Omega$ , the model must output a continuous full-field solution $u(x)$ . This is a mapping problem in high-dimensional function space where the degree of freedom can reach tens of thousands, and geometric information cannot be easily represented with standard grids.

To address these challenges, neural operator are designed to directly approximate the solution mapping $g^+$ , with the ability to handle arbitrary geometries and output fields.

Literature Review

Convolutional Neural Networks (CNNs) are typically designed for structured, regular grid data. However, they face significant limitations when applied to unstructured meshes and irregular geometries commonly found in Finite Element Method (FEM) discretizations. Directly mapping such geometries onto a regular grid often introduces interpolation errors and compromises geometric fidelity. Moreover, representing connected three-dimensional domains using high-resolution regular grids leads to high-dimensional inputs, resulting in reduced computational efficiency.

To overcome these limitations, PointNet³. was proposed as a method for directly processing point cloud data, enabling the capture of complex 3D geometries without the need for structured meshing PointNet treats each point independently and employs multilayer perceptrons (MLPs), followed by symmetric aggregation function (e.g., max pooling), to extract global features. However, PointNet lacks the ability to capture local spatial relationships and hierarchical structures, which are essential for many physical modeling tasks.

To address this shortcoming, PointNet++⁴. was introduced as an extension of PointNet. It incorporates a hierarchical feature learning framework by dividing the point cloud into multiple local regions and applying PointNet within each region to extract localized features. This multi-scale architecture enables the network to capture both fine-grained local patterns and global structural information, making it more suitable for modeling complex physical systems.

Traditional Graph Neural Networks (GNN)⁵. based surrogate models typically require a reference mesh $\Omega$ and a deformation mapping $\Psi$ , where mesh nodes are treated as graph vertices and mesh elements or edges are treated as connections. This formulation allows the model to effectively capture both geometric and topological information, and has demonstrated strong performance on domains with varying shapes. However, as the number of mesh nodes increases, the computational efficiency of GNN-based models tends to decline significantly.

(3) $\begin{eqnarray*}\Psi:\Omega_{ref} \mapsto \Omega,\end{eqnarray*}$

Deep Operator Network (DeepONet)⁶. is a neural operator architecture designed to learn mappings between functions, formulated as:

(4) $\begin{eqnarray*}g: u(x) \rightarrow g(u)(y).\end{eqnarray*}$

Here, $u(x)$ is the input function, and $g(u)(y)$ is the corresponding output field. The objective is to approximate the operator mapping $\mathcal{g}$ itself. DeepONet consists of two subnetworks: the Branch Net receives discrete observations of the function $u$ , and the Trunk Net takes the query point $y$ as input. The outputs of the two networks are combined via inner product to generate the final prediction. This design has strong generalization capabilities and can efficiently learn high-dimensional function mappings.

However, the standard DeepONet assumes a fixed domain and cannot handle complex geometries or variable boundary conditions. To address this, Geom-DeepONet⁷, an extended version, introduces the Signed Distance Function (SDF) as a representation of geometry, enabling the model to work with irregular shapes. In Geom-DeepONet, the SDF is input into the Branch Net as a continuous field, while the Trunk Net can receive material parameters and other control variables. This allows stable operator learning under varying geometric conditions.

The geometry of 3D objects can be represented either explicitly (e.g., point clouds or meshes) or implicitly (e.g., SDF). The SDF describes the shortest distance from each point to the geometric boundary. Its continuity and differentiability make it especially suitable for constructing geometry-aware latent spaces in neural operator architectures, offering robust generalization to geometric variations.

Fourier Neural Operator (FNO)⁸ is an integral operator learning method that enables convergence after discretization, using a set of Fourier kernel functions defined as $K_{k}(x,y) = e^{2 \pi i k \cdot (x - y)}$ . One limitation of this approach is its inefficiency on irregular grids.

Geometry-Informed Neural Operator (GINO)⁹ is a recently proposed framework that combines geometric field representation with neural operator learning. It integrates the efficiency of FNO in spectral learning with the representational capacity of Graph Neural Networks (GNNs) for unstructured geometric data. The key idea is to jointly describe geometry using SDF and point clouds, then map the input fields to a regular latent space (typically a uniform grid), allowing efficient global feature learning via Fourier convolution. The latent representation is then decoded back to make predictions on the original irregular geometry.

GINO and Point Cloud Neural Operator (PCNO)¹⁰ are similar in terms of training data and application scenarios, both capable of handling arbitrarily shaped domains. Compared to PCNO, which performs computation directly on point clouds, GINO emphasizes aligning irregular geometric fields to a structure suitable for spectral learning. This improves computational efficiency and model convergence.

Physics-Informed Neural Network (PINN) is a modeling framework that integrates deep learning with physical laws. Unlike purely data-driven neural networks, PINNs are trained not only with observational data but also under constraints from known physical equations, such as partial differential equations (PDEs), to enhance generalization and physical consistency.

For example, ResoNet¹¹ applies a PINN framework to simulate resonance in a one-dimensional acoustic tube in the time domain. Its underlying physical equation is the wave equation:

(5) $\begin{eqnarray*}\frac{\partial^2 u(x,t)}{\partial t^2} = c^2\frac{\partial^2u (x,t)}{\partial x^2}.\end{eqnarray*}$

In the PINN framework, the neural network is trained to satisfy this PDE. The total loss function consists of two parts:

(6) $\begin{eqnarray*}\mathcal{L}total = \mathcal{L}_{data} + \lambda \mathcal{L}_{physics}.\end{eqnarray*}$

where $\mathcal{L}{\text{data}}$ is the prediction error with respect to the observed data, and $\mathcal{L}{\text{physics}}$ is the physical error computed from the PDE residual.

The parameter $\lambda$ controls the weight of the physical term.

PINN focuses on learning the solution of a specific PDE, rather than learning an operator that maps functions to functions as in neural operator frameworks. This makes PINNs less suitable for handling complex geometries or high-dimensional PDEs. Furthermore, when boundary conditions, initial values, or physical parameters change, the model must be retrained.

To date, there are only two known studies applying neural networks to violin physics¹²’¹³, both authored by Gonzalez, S., Salvi, D., and collaborators. These studies use data-driven approaches, applying a single hidden layer feed-forward neural network with a sigmoid activation and linear output layer. The input consists of 20 geometry-related parameters, and the model predicts the first 10 eigenfrequencies. However, this method only estimates global eigenfrequencies based on coarsely described geometry, thickness, and material parameters, without resolving or visualizing the actual vibration patterns across the violin plate. Therefore, deep learning in musical acoustics remains an underdeveloped research area.

Training Data

The training data used in this study were generated using finite element simulations performed in COMSOL Multiphysics. The model archive created by Linus Andersson is available on the COMSOL Blog website. Andersson reproduced the experiments in COMSOL from the paper《The evolution of air resonance power efficiency in the violin and its ancestors》, a study conducted by graduate students at MIT and violin manufactures on North Bennett Street in Boston. The study investigated how variations in the shape of the violin’s f-holes affect the acoustic performance of the instrument body.

Fig 1. Acoustic-Structure Interaction (ASI) analysis result.

The original COMSOL model displays the sound pressure distribution and eigenfrequencies of a violin body. A spherical air domain was added around the violin body to allow sound modes to decay freely and avoid pressure truncation above the f-holes. In the mode shape dataset, the violin ribs originally defined as sound hard boundaries were removed. The model also included a soundpost connecting the top and back plates, showing how internal body pressure causes the plates to deform outward.

Building upon this setup, the present study modifies the original mode shape database by isolating only the top plate of the violin. Under fixed geometric dimensions, acoustic-structure interaction (ASI) simulations were performed using COMSOL Multiphysics. The top plate material is spruce, referenced from studies by Ulrike G. K. Wegst¹⁴ and Gonzalez et al¹³., have identified various material parameters that influence acoustic performance. Among these, density, Young’s modulus, and Poisson ratio were selected in this study as the primary variables for parametric analysis. The plate thickness varied from 2.0 mm to 3.5 mm.

The simulated top plate mesh consists of 28,309 nodes, each corresponding to a 3D spatial coordinate (x, y, z) and associated with displacement outputs, including magnitude and directional components (u, v, w). Internal pressure applied within the body cavity generates deformation on the plate surface, resulting in a spatially varying displacement field and sound pressure distribution. To focus on vibrational characteristics and material response, only the top plate was analyzed. The resulting dataset captures changes in eigenfrequencies and displacement fields across different combinations of material parameters and thickness values.

Fig 2. The sound pressure distribution of violin plate.

Only the top plate is selected for analysis, and consider the variations of the displacement field and eigenfrequencies under different material parameters and top plate thickness.

8 sets of material parameters have been simulated, each set with 16 sets of thickness, a total of 128 sets of sample data, each material parameter distribution:

Density	350 kg / m³	500 kg / m³
Young’s modulus	8 GPa	12 GPa
Poisson ratio	0.2	0.4

Table 1: Material parameter ranges.

Density $\rho$ Mean: $\mu = 400kg/m^{3}$ , Standard deviation $\sigma = 50kg/m^{3}$
Young’s modulus $E$ Mean: $\mu = 9.875 GPa$ , Standard deviation $\sigma = 1.366 GPa$
Poisson ratio $\nu$ Mean: $\mu = 0.2625$ , Standard deviation $\sigma = 0.0737$

Data Preprocessing

In the COMSOL Multiphysics model, output plots were selected for domain and edge. Domain is the entire area of the violin top plate, and edge is the contour lines of the boundary. exported as CSV files to obtain node data containing displacement information. The CSV filename stores global features: thickness, eigenfrequency, density, Young’s modulus and Poisson’s ratio.

The original domain contains 28309 nodes, and the edges contain 4934 nodes. Each raw CSV file includes the following columns: $(x, y, z, \text{displacement magnitude}, u, v, w)$ , which through downsampling, reconstruction nodes connectivity are converted into NumPy binary format during the preprocessing stage.

Downsampling

To reduce the computational cost associated with overly dense point distributions in certain regions, applied farthest point sampling to downsample the point cloud, which promotes a more uniform distribution of points while still preserving higher density in complexity regions. The algorithm begins by randomly selecting an initial point as the centroid $C$ and finding the minimum Euclidean distance point $P$ to the existing centroid. At each iteration, the next centroid is chosen as the point with the largest minimum distance to the already selected centroids.

(7) $\begin{eqnarray*}c_{k+1}=\arg \max\limits_{p \in P} \max\limits_{c \in C} ||p - c||^2.\end{eqnarray*}$

Connectivity

After combining domain and edge points and downsampling to 5000, the edge points account for 20% of the total, ensuring the f-hole shape and outline geometry are preserved. Since the data are downsampled, the original mesh connectivity cannot be used, so the node connectivity must be reconstructed. An element connected with three nodes forms a triangle mesh. To preserve the outline, edge points are connected to their nearest neighbors to create a smooth boundary.

Synthetic Data

In practice, generating additional samples from the COMSOL Multiphysics eigenvalue solver is computationally demanding, as each simulate requires large-scale PDE system on fine meshes. To evaluate how model performance scale with training data size without higher performance required, an alternative synthetic dataset was generated. The output value at each node was defined using a superposition of trigonometric basis function:

(8) $\begin{eqnarray*} u(x, y; 0) = A_1(0)\sin(w_{1x}(0)x+w_{1y}(0)y), \\ + A_2(0)\cos(w_{2x}(0)x-w_{2y}(0)y), \\ + A_3(0)\sin(w_{3x}(0)x)\sin(w_{3y}(0)y).\end{eqnarray*}$

$0$ is a constant, represent a material parameters. $A_i$ denotes amplitudes and $w_i$ is frequencies. The trigonometric formulation provides controllable variability and ensures smoothness consistent with PDE-based solutions. By varying these parameters, synthetic dataset of 256 samples were generated using same geometry and material parameter ranges as in the original setup.

In addition to the 2D synthetic dataset, a 3D violin back plate dataset was also generated from STL geometry using the same synthetic formulation. This allowed verify the model on higher-dimensional data without incurring the prohibitive computational cost of full 3D FEM simulations.

Methods

This work introduces a model of the Material Aware Spectral Operator (MASO).By directly input point cloud data with material parameters to predict displacement magnitude in each point and a plate eigenfrequency, mesh connectivity information is optional. displacement magnitude predict is through spectral for global branch, and comparison for different local branch process method, FNO⁸ use convolution to capture local effect, PCNO¹⁰ use differential operator, while MASO method using propose diffusion GNN.

The MASO model consists of four main modules:

Input encoders: received point clouds (x, y, z) and corresponding displacements (u, v, w).
Global features encoders: Process material properties, thickness, density, Young modulus, Poisson ratio, and eigenfrequency.
Spectral Layers: based on the Fourier Neural Operator architecture, perform frequency domain conversion and learn input features.
Local Branch Layers: preserves local interactions and the identity mapping.
Output decoders: combined local and global features, predict the displacement of each point and global eigenfrequency.

The nonlinear activation function $\sigma$ used in the network is the Gaussian Error Linear Unit (GeLU). GeLU is a smooth blend of Sigmoid and ReLU, featuring continuous derivatives that avoid the non-differentiability at $x=0$ encountered in ReLU. It also allows the retention of small negative values, helping to mitigate the vanishing gradient problem.

Spectral Layers

In numerical simulations, continuous geometries $\Omega$ must be discretized, this refers to the process of converting continuous mathematical entities into a finite set of discrete elements. Through discretization, an originally infinite-dimensional continuous structure is approximated by a composition of a finite number of discrete elements (points, lines, surfaces). The manifold $\Omega \subset \mathbb{R}^{3}$ is discretized into a point cloud $X = {x^{(i)}}_{i=1}^N$ , where $x^{(i)} \in \mathbb{R}^3$ denotes the coordinates of each sampled point.

The MASO architecture employs spectral layers inspired by Fourier neural operators, where the kernel function is parameterized using Fourier modes rather then direct point-to-point interactions. The spectral decomposition of the input function into frequency components and learnable weights is applied to the retained modes. For each spectral layer, the transformation is defined as:

(9) $\begin{eqnarray*}f_{\text{out}}(x) = \mathcal{F}^{-1}\left(\sum_{k \in \mathcal{K}} W_k \cdot \mathcal{F}(f_{\text{in}})(k)\right)(x) + W_0 f_{\text{in}}(x).\end{eqnarray*}$

where $\mathcal{F}$ and $\mathcal{F}^{-1}$ denote the Fourier transform and inverse Fourier transform, $\mathcal{K}$ is the set of retained frequency modes, $W_k$ are learnable weights for each mode $k$ , and $W_0$ is a linear transformation for the residual connection.

The operation can be expressed in the spatial domain as an integral operator:

(10) $\begin{eqnarray*} f_{out}(x)=\int_{\Omega} k(x,y,f_{in}) \cdot W_{v}f_{in}(y)dy.\end{eqnarray*}$

where $k(x,y,f_{\text{in}})$ is a learnable kernel function and $W_v$ is a linear transformation matrix. Since the exact integration over the continuous domain $\Omega$ is not tractable, it is approximated by discrete sampling:

(11) $\begin{eqnarray*} f_{\text{out}}(x) \approx \frac{1}{|\Omega|}\sum_{i=1}^{N} k(x,x^{(i)},f_{\text{in}})\cdot W_u f_{\text{in}}(x^{(i)}) d\Omega_i.\end{eqnarray*}$

Input Parameters

Each node scaled into $[-0.5, 0.5]$ and concatenated with the global features, forming the input representation $x=[coords, g] \in {\mathbb{R}^{Nxd}}$ . A set of $N$ -dimensional frequency vectors ${k_i}_{i=1}^M$ is then precomputed, where only the lowest $M$ modes are retained, and truncation of high-frequency components thus reduce computational cost, obtain frequency indices $k_x \in [-8,8]$ and $k_y \in [0, 8]$ due to the Fourier transform being symmetric, only half the frequencies is unique. producing a symmetric set of Fourier basis vectors that capture the most significant frequencies of the solution field.

Projection

Project each node $x_n \in \mathbb{R}^{C_{in}}$ onto the Fourier basis $\cos\langle x \cdot k_x + y \cdot k_y \rangle$ and $\sin\langle x \cdot k_x + y \cdot k_y \rangle$ . This projection step allows the model to represent irregular points in the Fourier domain and learn mode-specific coefficients $W$ for the input and output channels. The spectral components are computed as follows, producing an $(N \times M)$ matrix:

(12) $\begin{eqnarray*} y_{spectral} = \sum_{i=1}^{C_{in}} \sum_{k=1}^{K} X (\sin{(k \cdot x)} \cos{(k \cdot x)}) \cdot W.\end{eqnarray*}$

Each input feature is projected into a Fourier feature space and then linearly combined across frequency modes $k$ , with an additional $y_{dc}$ term included to represent the linear projection without frequency variation.

The spectral process follows the same idea as the FNO⁸. However, since the data lies on irregular points rather than a uniform grid, the Fast Fourier Transform (FFT) cannot be directly applied. Instead, each point is projected onto a truncated frequency basis, similar to the FNO formulation but adapted for irregular nodes.

In 2D space, the frequency vector $k$ can be expressed as:

(13) $\begin{eqnarray*} k=\bigg(\frac{2\pi}{L_x}k_x, \frac{2\pi}{L_y}k_y\bigg).\end{eqnarray*}$

Represents the $k$ -th Fourier mode, and a total of $K$ modes are used. $L$ represents the length of the domain (e.g., Fourier transform on $[0, L]^d$ ), which is used for frequency scaling.

$k (x - y)$ denotes Fourier kernel function, In the FNO⁸, the core of learning the integral operator is to capture the spatial correlation in the input function $f_{\text{in}}(x)$ through Fourier series. The kernel function $k(x - y)$ used here can be approximated as a translation-invariance and positive semi-definite form, i.e. it depends only on the relative distance between points, not on their absolute position in space. This feature allows the model to be generalized and consistent in different spatial regions and strengthens its ability to learn the spatial structure.

$k(x - y)$ is considered as a periodic function and can therefore be expanded into a Fourier series:

(14) $\begin{eqnarray*} k(x-y)=\sum_k(W^{c}_{k} \cos (k \cdot (x-y)) + W^{s}_{k} \sin(k \cdot (x-y))).\end{eqnarray*}$

$W_k^c$ and $W_k^s$ are learnable weights in the model, representing the response intensity to different frequency components.

This kernel function can be considered to be an optimal kernel learning from a set of smooth orthogonal bases. Unlike traditional CNN, the core is global and does not depend on the area of regional reception. It is suitable for processing field functions with non-local dependency, such as PDE solutions.

Local Branch Layers

The local branch is designed to capture short-range interactions that cannot be fully represented by spectral operators. In FNO⁸, this is achieved through point-wise convolution in the spatial domain, which preserves both local interactions and the identity mapping. In contrast, PCNO¹⁰ introduces a differential operator approach, where point connectivity is precomputed: an adjacency matrix $A(x)$ encodes the neighbors $x^{(1)}, x^{(2)},...,x^{(n)}$ ,

and local gradients $\nabla f_{in}(x)$ are estimated by solving a least-squares $A(x) \nabla f_{in}(x)$ . Meanwhile, MASO incorporates a diffusion-based graph neural network as its local branch, modeling neighborhood effects through message passing on irregular graphs.

In this section, compare four local-branch designs: the convolutional scheme of FNO, the differential operator method of PCNO, and the diffusion-based approach of MASO. While the spectral branch remains consistent across models.

Diffusion Operators

Implement a diffusion-based local operator that updates node features to capture spatial correlations on irregular meshes. This operator aggregates information from neighboring nodes according to the mesh connectivity. Directed edges are defined by the mesh structure, and the graph Laplacian is formulated as a diagonal matrix minus the adjacency matrix, where the diagonal entries $D_{ii} = |\mathcal{N}(i)|$ denote the number of neighbors of node $I$ . The adjacency matrix $A$ encodes the directed edges. This formulation is equivalent to the random-walk normalized Laplacian $L=I-D^{-1}A$ , which updates each node by moving it toward the average of its neighbors. The local diffusion operator measures the difference between a node and its neighbors, and can be expressed as:

(15) $\begin{eqnarray*} f^{new}_{n} = f_n + \alpha \biggl(\frac{1}{|N(n)|}\sum_{j \in N(n)} f_j - f_n\biggr).\end{eqnarray*}$

where $\alpha \in (0,1)$ is a diffusion coefficient that controls the strength of neighbor aggregation. Dividing by the number of neighbors normalizes the update, ensuring stability across irregular meshes with varying connectivity. Here, $n$ denotes the current node, and $j$ runs over its neighbors.

Diffusion Operator with Gradient Weights

In the weighted formulation, each directed edge $(i, j)$ is associated with a gradient-based weight. A pseudo-inverse of local neighbor differences is used to compute an edge weight vector $W_{ij}$ that approximates local derivatives. The update rule becomes:

(16) $\begin{eqnarray*} f^{new}_n = f_n + \eta \sum_{j \in N(n)} W_{nj}(f_j - f_n).\end{eqnarray*}$

where $\eta$ is a diffusion coefficient, and $W_{nj}$ encodes geometry-aware weights derived from local gradients. This operator generalizes the simple averaging scheme by incorporating directional information from the mesh geometry. Each node feature is represented as $f_n \in \mathbb{R}^C$ , where $C$ is the feature dimension. In addition, a global input vector $p \in \mathbb{R}^{g}$ , containing material parameters and target eigenfrequencies, is processed through 1D convolutional layers and then fused with the local node features to produce the final output.

Eigenfrequency

The global features are processed through two fully connected layers with a width of 128 to predict the eigenfrequency; after 20 epochs, the test relative L2 error is at 0.0142.

Results

Using 8 Fourier modes with 1.22M parameters achieves good performance, while increasing to 16 modes with 4.49M parameters slightly reduces overfitting after 80 epochs. 8 Fourier modes have a test loss of 0.088, while 16 Fourier modes have a test loss of 0.082. However, increasing the Fourier modes comes at a significant computational cost; 8 Fourier modes in this case are more efficient overall.

Training was performed using 80\% for training data and 20\% for testing data. Adam optimizer with a OneCycleLR scheduler, employing a weight decay of $1 \times 10^{-5}$ and a base learning rate of $5 \times 10^{-5}$ , while keeping the default hyperparameters $(\beta_1 = 0.9, \beta_2 = 0.999)$ .

Fig 5: 128 samples with 5000 point cloud data in 100 epochs.

Data Size	Loss	Per Epoch
128	0.1320	42s
256	0.0883	77s

Table 2. Spectral + Conv1d (1.22M Parameters).

2D violin plates only required spectral convolution layers and point cloud input without mesh connectivity information, which is the most efficient way, and the lowest loss is close to adding a local branch, such as convolution 1D, diffusion GNN, and the differential operator by PCNO. Adding to convolution 1D and local branch layers, similar to PCNO architecture, may have the risk of under-fitting.

Model	Loss	Per Epoch
Spectral+PCNO	0.1864	1231s
Spectral+Conv1D	0.1617	1147s
Spectral+Diffusion	0.1543	1143s
Spectral+GNN	0.1452	1289s

Table 3. 3D Violin Spectral + Local (20.17M Parameters).

For 3D violin plates using 256 samples of synthetic data. The Spectral+Diffusion variant applies a fixed local diffusion operator on the graph. The Spectral+GNN variant extends this by directed edges and edge gradient weights, making the diffusion operator trainable.

All training and inference experiments were conducted on a 2021 Apple MacBook Pro with the M1 Pro chip. For comparison, a single COMSOL Multiphysics finite element simulation involved approximately 455,000 degrees of freedom and required around 10 hours 59 minutes and 40 seconds to complete one simulation and data export.

Only the 2D violin plate dataset (128 samples) was generated using COMSOL. For larger datasets (256 samples) and the 3D violin plate, synthetic data based on STL models was used, since direct COMSOL simulation of such scales would exceed the computational resources.

In contrast, the proposed Material-Aware Spectral Operator model completed training after 100 epochs, with a total training time of only 1 hour 10 minutes.

Fig 6: Visualization of predictions; Left: Ground truth data; Middle: Predictions; Right: Absolute error.

Discussion

As far as currently known, this is the first neural operator model specifically designed to predict the full-field displacement and regress global material parameters for violin plates. The architecture provides high flexibility: it can operate directly on point cloud data or incorporate mesh connectivity information. The input includes both point-wise displacement field data and global material parameters, enabling real-time prediction under varying material and geometric configurations. Compared with traditional finite element methods, the proposed approach achieves a speedup of approximately 9.4x while maintaining high accuracy.

Acknowledgments

The author thanks Prof. Te-Sheng Lin (Department of Applied Mathematics, National Yang Ming Chiao Tung University, NYCU) for guidance on synthetic data generation and providing computational resources.

References

J. Curtin. TAP TONES AND WEIGHTS OF OLD ITALIAN VIOLIN TOPS. Violin Society of America 20(2), 162 (2006). [↩]
C. M. Hutchins. Catgut Acoustical Society, Inc. In SMAC 83: Proceedings of the Stockholm Music Acoustics Conference, July 28–August 1, 1983 42, 47. Kungl. Musikaliska akademien (1983). [↩]
C. R. Qi, H. Su, K. Mo, L. J. Guibas. PointNet: Deep learning on point sets for 3D classification and segmentation. arXiv preprint arXiv:1612.00593 (2016). [↩]
C. R. Qi, L. Yi, H. Su, L. J. Guibas. PointNet++: Deep hierarchical feature learning on point sets in a metric space. arXiv preprint arXiv:1706.02413 (2017). [↩]
G. Z. Jin. Leveraging graph neural networks and neural operator techniques for high-fidelity mesh-based physics simulations. AIP Advances 13(4), 046109 (2023). DOI: 10.1063/5.0167014 [↩]
L. Lu, P. Jin, G. E. Karniadakis. DeepONet: Learning non-linear operators for identifying differential equations based on the universal approximation theorem of operators. arXiv preprint arXiv:1910.03193 (2019). [↩]
J. He, S. Koric, D. Abueidda, A. Najafi, I. Jasiuk. Geom-DeepONet: A point-cloud-based deep operator network for field predictions on 3D parameterized geometries. arXiv preprint arXiv:2403.14788 (2024). [↩]
Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, A. Anandkumar. Fourier Neural Operator for Parametric Partial Differential Equations. arXiv preprint arXiv:2010.08895 (2020). [↩] [↩] [↩] [↩] [↩]
Z. Li, N. B. Kovachki, C. Choy, B. Li, J. Kossaifi, S. P. Otta, M. A. Nabian, M. Stadler, C. Hundt, K. Azizzadenesheli, A. Anandkumar. Geometry-Informed Neural Operator for Large-Scale 3D PDEs. (2023). arXiv preprint arXiv:2309.00583 (2023). [↩]
C. Zeng, Y. Zhang, J. Zhou, Y. Wang, Z. Wang, Y. Liu, L. Wu, D. Z. Huang. Point cloud neural operator for parametric PDEs on complex and variable geometries. arXiv preprint arXiv:2501.14475v2 (2025). [↩] [↩] [↩]
K. Yokota, T. Kurahashi, M. Abe. Physics-informed neural network for acoustic resonance analysis in a one-dimensional acoustic tube. The Journal of the Acoustical Society of America 156(1), 30–43 (2024). [↩]
D. Salvi, S. Gonzalez, F. Antonacci, A. Sarti. Parametric Optimization of Violin Top Plates using Machine Learning. arXiv preprint arXiv:2102.07133 (2021). [↩]
S. Gonzalez, D. Salvi, D. Baeza, F. Antonacci, A. Sarti. A Data-Driven Approach to Violin Making. arXiv preprint arXiv:2102.04254 (2021). [↩] [↩]
U. G. K. Wegst. Wood for sound. American Journal of Botany 93(10), 1439–1448. DOI: 10.3732/ajb.93.10.1439 (2006). [↩]