arXiv CSopen access

SplineNet: An Isogeometric Deep Learning Method for Complex Shells

SplineNet: An Isogeometric Deep Learning Method for Complex Shells Shizhou Luo1 and Xiaodong Wei1* 1

Global College, Shanghai Jiao Tong University, Shanghai, 200240, China.

*Corresponding author(s). E-mail(s): [email protected]; Contributing authors: [email protected];

arXiv:2607.06026v1 [cs.LG] 7 Jul 2026

Abstract We present a novel isogeometric deep learning method, termed SplineNet, for the seamless design and analysis of shell structures with complex geometries. The proposed approach is built upon watertight spline representations, e.g., analysis-suitable unstructured T-splines, and features exact geometric descriptions of Computer-Aided Design (CAD) models in neural networks. Bézier extraction is used to build the network architecture, where Bernstein polynomials serve as the nonlinear activation functions. SplineNet can be applied in a data-free or data-driven way. In the data-free case, energybased formulations can be naturally incorporated as loss terms, which fulfill the need of ComputerAided Engineering (CAE) and can be accurately calculated. In particular, the Kirchhoff–Love (KL) model is adopted to solve for the mechanical behaviors of shell structures. This way, CAD and CAE can be tightly integrated in a deep neural network without the time-consuming model/data exchange process. In the data-driven case, SplineNet can be used as the trunk net of Deep Operator Networks (DeepONet) to provide interpretability. Given such a trained network and unseen input data, results can be immediately obtained without retraining the network or repeatedly performing the traditional workflow for analysis. In the end, a variety of numerical examples are studied to demonstrate the effectiveness of the proposed method, especially when real-world complex geometries are involved. Keywords: Deep Operator Networks, Interpretability, Isogeometric Analysis, Analysis-suitable Unstructured Splines, Bézier Extraction

1 Introduction

the finite element method (FEM) [2]. In particular, spline-based discretizations possess higherorder continuity, making them especially suitable for high-order partial differential equations (PDEs), such as the Kirchhoff–Love (KL) shell formulation [3]. Despite these advantages, traditional IGA methods become computationally expensive for large-scale problems and require repeated simulations in design optimization [4] and uncertainty quantification [5]. To overcome these drawbacks, leveraging deep learning to enhance the performance of traditional

Isogeometric Analysis (IGA) was proposed to bridge the gap between Computer-Aided Design (CAD) and Computer-Aided Engineering (CAE) by employing spline-based representations used in CAD as basis functions for numerical analysis [1]. By directly utilizing exact geometric descriptions, IGA enables higher geometric fidelity and improved numerical performance compared with

1

PDE solvers has attracted significant attention [6]. Based on learning objectives, these approaches can be broadly classified into two categories: solution learning and operator learning. In the solution learning paradigm, the neural network acts as a function approximator to solve a specific problem. Many methods have been studied accordingly based on both strong and energy/weak forms of PDEs, such as physicsinformed neural networks (PINNs) [7], the deep Ritz method (DRM) [8], the deep energy method (DEM) [9], and variational physics-informed neural networks (VPINNs) [10]. PINNs utilize strong form residuals to define loss functions. They provide a flexible, mesh-free approach for solving forward and inverse problems. However, they require high-order derivatives and lead to training difficulties [11, 12]. DRM, DEM, and VPINNs employ the energy/weak form of PDEs as loss functions, which can help reduce the order of derivatives [13] and improve computational efficiency [14]. Despite these improvements, in many implementations of these methods, essential boundary conditions are imposed through penalty terms that cause conflicting gradients during training [15, 16]. Another line of solution-learning methods is inspired by classical discretization techniques, such as hierarchical deep learning neural networks (HiDeNN) [17–20] and finite element neural network interpolation (FENNI) [21, 22]. These methods construct finite-element-like shape functions using neural network architectures. Such networks can enhance interpretability, accurately impose boundary conditions, and facilitate advanced strategies such as r -adaptivity [18, 21] and multigrid training [21]. However, these approaches remain instance-specific: any change in geometries, boundary conditions, or loadings requires retraining. In contrast, operator learning methods aim to learn mappings from input functions to output solution fields, allowing them to provide instant predictions for unseen inputs without retraining. For example, Fourier neural operators (FNO) [23] leverage the fast Fourier transform (FFT) to achieve efficient and resolution-independent operator learning. However, FNO typically requires functions to be represented on uniform grids, which limits its application to complex domains. Geo-FNO [24] addresses this issue by learning to transform an irregular physical domain into

a latent uniform grid. Graph kernel networks (GKNs) [25] approximate kernel integral operators on graph-based discretizations and can therefore handle irregular meshes. Deep Operator Networks (DeepONet) [26] employ a branch net to encode input functions and a trunk net to encode output locations, providing a flexible framework for learning nonlinear operators. Several variants exist to enhance the performance of DeepONet. For example, the residual U-Net is leveraged to encode complex geometries [27]. Point-DeepONet [28] integrates PointNet to learn geometric information from point clouds. Neural operators on Riemannian manifolds (NORM) [29] represent the input and output using the eigenfunctions of the Laplace–Beltrami operator defined on the underlying manifold, thereby providing a spectral basis for learning the input-output mapping. These methods extend operator learning beyond regular grids. However, there exists a fundamental gap for these methods to be applied in engineering design, because their geometric representations are recognized by sampled points in the computational domain rather than by CAD models. This inconsistency poses several limitations. First, the performance of neural nets may be sensitive to the choice of point locations, especially in problems where local features, such as stress concentration and cracks, are important [30–32]. Second, even with a carefully curated point set, analysis results still need to be mapped back to the original CAD model to facilitate design iterations, which can be error-prone, especially for large-scale problems. Third, visualization of analysis results requires additional care, including the choice of point locations for evaluating the design and the triangulation of such points to yield intuitive and high-quality visualizations, which often require domain knowledge. Therefore, incorporating the isogeometric concept, i.e., using the same geometric representation everywhere, into learning-based methods has significant potential to enhance conventional numerical solvers while seamlessly integrating CAD and CAE in neural nets. Recent studies have started to exploit this idea by integrating spline representations into solution learning methods. For instance, isogeometric neural networks (IGN) [33] train a neural network to predict control variables, analogous to 2

nodal values in FEM, rather than pointwise solution values, thereby reducing the computational burden and keeping the learned solution in an isogeometric approximation space. Isogeometric convolution HiDeNN (C-IGA) [34] introduces convolution into HiDeNN to reproduce splines and support adaptive refinement. Later, it has been extended to the multi-patch case [35] by imposing compatibility conditions across patch interfaces. Alternatively, multi-patch isogeometric neural solvers [36] formulate neural ansatz functions on IGA reference domains, where interface conditions are enforced and an energy-based loss is used for training. Nonetheless, multi-patch representations are rather restrictive in modeling complex geometries, as they require a multi-block mesh structure that is generally difficult to achieve. Moreover, enforcing interface conditions can be tedious and challenging for high-order problems such as Kirchhoff–Love shells. To address these issues, we propose a novel isogeometric deep learning method, termed SplineNet, for the seamless design and analysis of thin shell structures with complex geometries; see Fig. 1 for an overview. Inspired by basis-function-aware networks, e.g., HiDeNN and C-IGA, SplineNet directly embeds analysis-suitable unstructured splines (ASUTS) [37], an advanced spline technique for modeling complex geometries in a watertight manner, into neural networks via Bézier extraction. The weights and biases of SplineNet are obtained from the extraction matrices, with Bernstein polynomials serving as the activation functions. The inputs are parametric coordinates, whereas the learnable parameters are control variables or control points, depending on the application of interest. In this way, the geometric model of SplineNet remains the same in both design and analysis. Due to the modeling capability of ASUTS, SplineNet can accommodate geometries from realworld applications that go beyond simple singleor multi-patch settings. SplineNet supports both solution learning and operator learning. In the former case, SplineNet acts as a physics-informed neural net to solve governing PDEs using an energy-based method for a single instance, and it is referred to as Physics-informed SplineNet (P-SNet ). In the latter case, the trunk net of DeepONet is replaced with SplineNet. As a result, the solution field is

represented as a linear combination of splines, thus enhancing the interpretability of operator networks. The variant in this case is called Operator SplineNet (O-SNet ). The main contributions of this work are summarized as follows: 1. We propose SplineNet to explicitly embed unstructured splines into neural nets, achieving watertight representations of complex geometries and seamless CAD/CAE integration in neural nets. 2. We employ SplineNet in both solution and operator learning to enhance IGA with learning capabilities and to provide interpretability for operator networks. 3. We apply SplineNet with Kirchhoff–Love shells to demonstrate its capability in handling realworld geometries and high-order PDEs. The remainder of the paper is organized as follows. Section 2 introduces ASUTS along with Bézier extraction. Section 3 presents the formulation of Kirchhoff–Love shells. Section 4 details P-SNet and O-SNet. Numerical examples are presented in Section 5. Finally, Section 6 draws conclusions and provides suggestions for future research directions.

2 Analysis-suitable unstructured T-splines Non-uniform rational B-splines (NURBS) [38] are the current industrial standard for representing CAD geometries. Complex geometries are commonly represented as collections of trimmed NURBS patches. Although trimmed NURBS provide great flexibility for geometric modeling, the trimming operation generally leads to gaps and overlaps between adjacent trimmed patches, hindering the model/data exchange between CAD and CAE. To overcome this fundamental issue, Tsplines [39] and related variants have been developed. Among them, analysis-suitable unstructured T-splines (ASUTS) [37] provide desirable properties for both design and analysis, such as watertightness for complex geometries, global smoothness, linear independence, partition of unity, optimal convergence, and adaptivity. We therefore adopt ASUTS as the foundation for the seamless integration of CAD and CAE in this

3

CAD/CAE integration in AI Unstructured Splines

+

Neural Network

SplineNet

Solution learning paradigm

SplineNet

+

Physics Loss

P-SNet

Operator learning paradigm Branch Net SplineNet (Trunk Net)

+

Data Loss

O-SNet

Fig. 1: Overview of SplineNet and its applications in both solution learning and operator learning. Unstructured splines are explicitly embedded into the neural network to construct the foundational SplineNet (top). It supports two distinct paradigms: a solution learning paradigm (P-SNet) driven by a physical loss (e.g., the total energy of Kirchhoff-Love shells) for solving PDEs, and an operator learning paradigm (O-SNet) that incorporates an additional branch net and a data loss to achieve instant predictions under varying input conditions. where Ni (ξ, η ) denotes the ASUTS basis function associated with the control point Pi , and n is the number of basis functions. Every ASUTS basis function is a bicubic piecewise smooth polynomial. They are defined elementwise and vary with element types. According to the influence of EPs, faces in an ASUTS control mesh are mainly classified into irregular, transition, and regular elements. An irregular element is a face that contains at least one EP. Transition elements are adjacent to irregular ones. The remaining faces are regular. ASUTS basis functions in regular elements are simply bicubic Bsplines. However, special treatments are required around EPs to ensure desired properties such as smoothness and convergence. More specifically, extra face-based points and associated splines are added in irregular elements to guarantee convergence, on top of which the D-patch method [40] is

work. In what follows, we briefly review the basics of ASUTS. Interested readers may refer to [37] for details. ASUTS are built from two components: an unstructured quadrilateral control mesh (possibly with T-junctions) and the corresponding spline basis functions. The control mesh consists of control points, edges, and faces (or elements). A key to ASUTS is the treatment of extraordinary points (EPs). An EP is an interior point of valence other than four, or a boundary point of valence other than three or two, where the valence of a point is the number of edges sharing the point. EPs are inevitable for representing complex geometries; see Fig. 2. A geometry is represented as



 n x(ξ, η ) X r (ξ, η ) = y (ξ, η ) = Ni (ξ, η )Pi , i=1 z (ξ, η )

(1)

4

(a)

(b)

Fig. 2: Car hull model (a) and its control mesh (b) with extraordinary points (red dots). adopted to build C 1 -continuity into splines. Subsequently, the truncation mechanism [41] is used in transition elements to make the newly added splines compatible with existing ones and also to recover partition of unity. ASUTS basis functions are represented through Bézier extraction for every element [42]. On element e, the local ASUTS basis functions N e (ξ, η ) are expressed as a linear combination of the 16 bicubic Bernstein polynomials B (ξ, η ):

N e (ξ, η ) = C e B (ξ, η ),

elements, which determine k . We observe that k varies in different elements. Moreover, we pick one spline function, NI (ξ, η ), near an EP to illustrate how it is represented elementwise; see the red dot in Fig. 5a, which is its corresponding control point. When restricted on element e, NI (ξ, η ) is expressed as a linear combination of the 16 Bernstein polynomials, 16 X e Nie (ξ, η ) = Cij Bj (ξ, η ), (5) j=1

(2)

where i is the local index of I on element e, and e the coefficients Cij are from C e . Their specific values corresponding to the three representative elements are shown in Fig. 5. Note that the original irregular face is locally split into 2 × 2 Bézier elements, and only the lower-left one is considered here for illustration. Such coefficients will be used to construct the proposed SplineNet. This elementwise representation is central to the proposed method. It converts the global ASUTS into a unified elemental data structure, where each element is described by Bernstein polynomials and a precomputed extraction operator. As will be shown in Section 4, this structure naturally leads to the construction of SplineNet.

where C e is the Bézier extraction operator and B (ξ, η ) is [B1 (ξ, η ), . . . , B16 (ξ, η )]T := [ b1 (ξ )b1 (η ), b2 (ξ )b1 (η ), . . . , b4 (ξ )b4 (η ) ]T . (3) Here bi (·) is defined by

 bi (s) =

 3 si−1 (1 − s)4−i , i−1

(4)

where i = 1, 2, 3, 4 and s ∈ [0, 1]. C e has dimension k × 16, where k is the number of ASUTS basis functions defined on element e. Note that k = 16 in regular elements, whereas k varies in irregular and transition elements. An example is illustrated in Fig. 3. For details, we consider three representative Bézier elements as shown in Fig. 4. The red, magenta, and green elements correspond to the irregular, transition, and regular Bézier elements, respectively. The splines associated with the highlighted control points have support on these three

3 Kirchhoff–Love shell formulation The Kirchhoff–Love (KL) shell theory provides a classical model for thin-shell analysis by neglecting transverse shear deformation and requiring the shell director to remain normal to the middle surface during deformation [3]. The resulting 5

(a) Control mesh

(b) Bézier mesh for analysis

Fig. 3: Example of an ASUTS control mesh and its Bézier mesh. (a) The control mesh contains regular control points (black dots), extraordinary points (red dots), and face-based points (yellow dots), where irregular, transition, and regular elements are shaded dark gray, gray, and light gray, respectively. (b) The corresponding Bézier mesh through Bézier extraction; note that every irregular element in (a) yields four Bézier elements due to the 2 × 2 split, which is used by ASUTS to locally increase the degrees of freedom. weak formulation involves second-order derivatives of the displacement and therefore requires at least C 1 -continuous basis functions. This makes ASUTS particularly suitable for KL shell analysis. Given a thin shell X in the 3D Euclidean space, let its reference mid-surface be parameterized by

X (θ ) : Ω̂ ⊂ R2 → R3 ,

θ = (θ1 , θ2 ),

The membrane and bending strain measures are defined as

εαβ =

Aαβ = Aα · Aβ , Bαβ = Aα,β · a3 ,

(6)

(7)

δWint (u; δu) = δWext (δu) ∀δu,

α = 1, 2,

(11)

(12)

with

Z aα = x,α ,

aαβ = aα · aβ , bαβ = aα,β · a3 .

The weak form is obtained from the principle of virtual work: find u such that

where u is the mid-surface displacement. The reference and current covariant basis vectors are defined by

Aα = X,α ,

καβ = Bαβ − bαβ , (10)

where

and let its deformed configuration be

x(θ ) = X (θ ) + u(θ ),

1 (aαβ − Aαβ ) , 2

δWint =

(8)

δεαβ nαβ + δκαβ mαβ dA,

ZA δWext =

with the corresponding unit normals

(13)

δui fi dA. A

A1 × A 2 A3 = , ∥A1 × A2 ∥

a1 × a2 a3 = , ∥a1 × a2 ∥

Here nαβ and mαβ are the membrane stress resultants and bending moments, fi is the applied load,

(9)

6

(a)

(b) k = 31

(c) k = 17

(d) k = 16

Fig. 4: Supported splines on different types of elements. (a) Three types of elements under consideration: an irregular element (red), a transition element (magenta), and a regular element (green). (b–d) The splines associated with the highlighted control points have support on the element of interest, i.e., the irregular, transition, and regular elements, respectively.

7

0.027778 0.027778 0.024306 0.019965 0.055556 0.055556 0.048611 0.039931 0.111111

0.111111

0.097222 0.079861

0.222222 0.222222 0.194444 0.159722

(a)

(b)

(c)

(d)

Fig. 5: Extraction coefficients associated with a selected spline function. (a) The selected control point (marked in red). (b–d) The extraction coefficients associated with the selected spline in the irregular (dark gray), transition (gray), and regular (light gray) elements, respectively. Note that only the lowerleft Bézier sub-element is shown for the irregular element.

8

p and dA = |Aαβ |dθ1 dθ2 . The constitutive relations follow the St. Venant–Kirchhoff KL shell model. The present work employs the linear KL shell [43] as the model problem. Under the small-deformation assumption, the total potential energy is written as Π(u) = Em (u) + Eb (u) − Wext (u), where

Z 1 Em = εαβ nαβ dA, 2 A Z 1 καβ mαβ dA, Eb = 2 A Z Wext = fi ui dA.

parametric coordinates. It contains 16 neurons whose activation functions are the 16 bicubic Bernstein polynomials. All the weights are one, and the bias is zero. As a result, given the input parametric coordinate (ξ, η ), where 0 ≤ ξ, η ≤ 1, the first layer simply outputs the corresponding Bernstein polynomial values B (ξ, η ). This layer does not involve trainable parameters. The second hidden layer produces ASUTS basis functions as linear combinations of Bernstein polynomials. It contains k neurons, where k is the number of ASUTS basis functions defined on the element of interest (e.g., element e). The weights come from the extraction matrix C e , whereas the bias is zero. Recall that an ASUTS basis function Nie is expressed as a linear combination of BernP16 e stein polynomials Bj , Nie = j=1 Cij Bj , where e e Cij is an element of C and it is used to define the corresponding weight in this layer. ASUTS basis functions serve as the activation functions of this layer. As a result, given the input B (ξ, η ), this layer outputs ASUTS basis functions N e (ξ, η ); see Eq. 2. The final output combines ASUTS basis functions with control variables to predict the solution at (ξ, η ). The weights are the trainable control variables U e . U e is a k × d matrix, where d (=2 or 3) is the spatial dimension. Each column of U e is a vector of control variables along a certain direction. The bias is zero. No activation function is needed in this layer. The output is the predicted solution, and it is obtained by

(14)

(15)

A

This energy functional will be used in the next section to construct the physical training loss.

4 Isogeometric deep learning method In this section, we introduce the details of SplineNet, which serve as the building block of the proposed method. SplineNet can be used in both solution learning and operator learning. In the former case, the Physics-informed SplineNet (P-SNet) is presented with an energy-based loss to solve single-instance problems, whereas in the latter case, SplineNet serves as the trunk net of DeepONet to enable instant predictions in a CAD/CAE-integrated manner.

ue (ξ, η ) = (U e )T N e (ξ, η ) = (U e )T C e B (ξ, η ),

4.1 SplineNet SplineNet embeds spline basis functions into neural networks through Bézier extraction, thereby enabling a CAD/CAE-integrated representation. Since both the geometry and spline basis functions are evaluated for each Bézier element in IGA, the elementwise SplineNet is constructed. Its architecture is illustrated in Fig. 6. For each element, SplineNet is a fully connected network with two hidden layers. The input layer is the parametric coordinate (ξ, η ) and the output layer is the predicted solution ue (ξ, η ). The two hidden layers correspond to Bernstein polynomials and ASUTS basis functions, respectively. More specifically, the first hidden layer produces Bernstein polynomials evaluated at certain

(16)

which recovers the discretization of IGA. Remark 1 Note that, in this work, the weights of the second hidden layer are fixed as the Bézier extraction coefficients C e , and the bias is set to zero. With this choice, the second hidden layer exactly recovers the desired ASUTS basis functions. While these weights and biases can also be trained (which yields the r adaptivity), additional treatments are then required to preserve the C 1 -continuity of ASUTS and only marginal improvement can be achieved. Therefore, we keep the weights and the bias fixed.

9

Elementwise SplineNet

Inputs

Outputs

B1

B3

N2e

B4

N3e

N4e

B13

N5e

B14

…

B15

Nke

B16

16 neurons

Ue1 Ue2 Ue3 U4e U5e

ue (ξ, η)

…

η

N1e

…

ξ

B2

Uke coefficients of Bézier extraction operator Ce trainable control variables Ue

k neurons

Fig. 6: Architecture of the elementwise SplineNet. The input parametric coordinate (ξ, η ) is first mapped to the Bernstein polynomial layer, and then transformed into ASUTS basis functions through the Bézier extraction matrix C e . The final layer combines the ASUTS basis functions with the elemental control variables to predict the displacement field on element e.

4.2 Physics-informed SplineNet

operation using the elementwise index list from ASUTS. Each sub-network predicts the displacement ue (ξ, η ) restricted to its Bézier element, which is then used to compute the total potential energy of the KL shell. The derivatives of the geometry and displacement fields are also required to evaluate the energy. These derivatives can be obtained in two ways. One option is to use automatic differentiation in PyTorch with respect to the parametric coordinates (ξ, η ). Alternatively, we explicitly construct derivative layers for the spline basis functions. The first hidden layer is replaced by the derivatives of the 16 Bernstein polynomials instead of the Bernstein polynomials themselves, which gives B,α (ξ, η ), where α ∈ {ξ, η}. The extraction matrix is then applied in the same way as in SplineNet:

Physics-informed SplineNet (P-SNet) applies SplineNet to solve a single-instance problem without precomputed training data. As shown in Fig. 7, the overall SplineNet for the entire domain is obtained by stacking elementwise SplineNets for all the Bézier elements. In other words, each elementwise SplineNet serves as a sub-network in the overall architecture, where there are no connections among different elements, and thus sparsity is achieved to facilitate training. Note that different sub-networks generally share common control variables. In other words, a single control variable appears in several elements/sub-networks. To guarantee their uniqueness, similar to the conventional FEM/IGA, a local-to-global index map is introduced to build the correspondence between the local degrees of freedom U e to their global ones U . In practice, U e is directly gathered from U by, for example, the PyTorch [44] indexing

e N,α (ξ, η ) = C e B,α (ξ, η ).

10

(17)

Inputs

SplineNet

Outputs

Loss functions

element 1

…

…

u1 (ξ, η) …

…

ξ

element e

…

…

Boundary Loss

…

…

η

Π (u; U) +

ue (ξ, η)

element n

…

…

un (ξ, η) Minimize Update U Fig. 7: Schematic of P-SNet.

Therefore, the obtained by

displacement

derivatives

are

boundary condition:

(18)

L(U ) = Π(u; U )+

M

e ue,α (ξ, η ) = (U e )T N,α (ξ, η ).

1 X 2 ∥û(ξb , ηb ) − u(ξb , ηb )∥ , M b=1 (19) where the energy Π(u; U ) is defined in Eq. 14, M is the number of predefined boundary points, and û(ξb , ηb ) and u(ξb , ηb ) are given and predicted displacements of the chosen boundary points, respectively. Usually, these boundary points are equally spaced in the parametric domain. In this work, the prescribed boundary displacements are imposed directly on the corresponding boundary control variables. To accurately evaluate Π, the input parametric coordinates are chosen to be the Gaussian points during training. This way, PSNet features a spline representation and enforces the shell governing physics through the energy formulation.

Higher-order derivatives are computed in the same manner by replacing the first hidden layer with the corresponding higher-order derivatives of the Bernstein polynomials. The same procedure is used for the geometric mapping by replacing the displacement control variables U e with the control point coordinates. Thus, the geometric derivatives aα and aα,β in Eq. 11 are obtained with the displacement derivatives. The Jacobian of the geometry mapping is further constructed from these geometric derivatives. Moreover, the derivatives of ue with respect to the physical coordinates (x, y, z ) can be obtained by the chain rule if needed. During training, the global control variables U are updated by minimizing the loss term that combines the total potential energy and the Dirichlet

4.3 Operator Spline Net Operator SplineNet (O-SNet) extends SplineNet to learn mappings from input functions (e.g.,

11

loadings) to shell structural responses. The architecture is developed based on DeepONet [26], which can be used to provide instant predictions for new input without retraining. A standard DeepONet contains two networks: a branch net and a trunk net. The branch net is used to encode the input function, while the trunk net is used to encode the query location where the output function is evaluated. The outputs of them are then coupled to approximate the target operator. Different kinds of neural architectures can be used in both networks, such as multi-layer perceptrons (MLP) and convolutional neural networks. In practice, the input function is represented by discrete values at a group of sensor points. For function f , these values can be T written as f = [f (x1 ), f (x2 ), . . . , f (xm )] , where xi is a sensor point and m is the number of sensor points. The branch net takes f as input and outputs a latent vector. Meanwhile, the trunk net takes the query locations y as input and outputs another latent vector. The two latent vectors are coupled through an inner product to predict the output function value at y , i.e., G(f )(y ), where G is the mapping/operator from the input function to the output function. A schematic of DeepONet is shown in Fig. 8.

Branch Net

…

f1 f2 f3

Ultimately, the trunk net of O-SNet provides the ASUTS basis functions N (ξ, η ). Doing so can not only explicitly embed high-fidelity geometric representations in a learning-based method, but also enhance IGA with learning capabilities. Instead of letting a neural net learn the basis of the output function, as is done in DeepONet, O-SNet introduces a prior to fulfill the same task. As a result, O-SNet provides predictions in the same discretized form as the conventional IGA, and thus enhances the interpretability of DeepONet. Moreover, the performance of O-SNet is independent of query point locations. Indeed, the trunk net of OSNet remains the same no matter where the query points are chosen. The branch net of O-SNet is simply MLP. It remains mostly the same as DeepONet except the fact that it encodes the input function f into control variables U (θ), where θ denotes the trainable parameters of the branch net. When combining with the output of the trunk net, O-SNet yields the desired prediction, u(ξ, η ; θ) = U(θ)T N (ξ, η ). This way, O-SNet preserves the overall structure of DeepONet while keeping the output field in the same spline space for CAD/CAE integration. The operator learning process in O-SNet is carried out in a data-driven manner. A set of labeled (i) data {f (i) , û(i) }N denotes i=1 is given, where f a discretized input function evaluated at sensor points, and û(i) is the corresponding displacement field obtained from an FEM/IGA solver. The predicted displacement field corresponding to f (i) is denoted by u(i) . The loss function consists of two terms: the mean squared error (MSE) of the displacement at sampling points (ξ, η ) and the boundary conditions,

…

fm

Discretized input function

G(f)(y) Output function G(f) at location y Trunk Net

…

y

Location point

L(θ) = Ldata + Lbc

Fig. 8: Architecture of DeepONet. The branch net encodes the discretized input function evaluated at sensor points, while the trunk net encodes the query location. The two latent vectors are coupled to approximate the target operator and predict the output function value at the query location.

=

N

K

N

M

2 1 X X (i) û (ξj , ηj ) − u(i) (ξj , ηj ; θ) N K i=1 j=1

2 1 X X (i) û (ξb , ηb ) − u(i) (ξb , ηb ; θ) , N M i=1 b=1 (20) where K denotes the number of sampling points for the output, and M denotes the number of boundary points. In this work, boundary conditions are enforced through control variables.

+

In O-SNet, this structure is adapted to the isogeometric setting mainly through the trunk net, as shown in Fig. 9. The trunk net of O-SNet features SplineNet rather than other options (e.g., MLP). 12

f(x1)

U1

…

… Branch net

Control variables

element 1

index mapping

Input functions

index mapping

Un

f(xm)

Displacement field

…

…

ξ η

η

u(ξ, η)

Minimize

×

Loss functions

element n

× ξ Parametric coordinates

SplineNet as the trunk net

Fig. 9: Architecture of O-SNet. The branch net maps discretized loading functions to global control variables U , while SplineNet serves as the trunk net to provide the ASUTS basis functions. Their combination yields the displacement field in the same spline space as IGA, enabling operator learning in a CAD/CAE-integrated representation.

5 Numerical examples

are considered where EPs are intentionally introduced to test ASUTS in the context of SplineNet; see Fig. 10. To keep the presentation concise, Fig. 11 reports the displacement field comparison on the finest mesh, whereas the free-edge midpoint displacements and errors for all meshes are summarized in Table 1. For all the three meshes, the relative errors of the free-edge midpoint displacement remain below 1%. In particular, on the finest mesh, P-SNet achieves a relative error of 0.679%, and the predicted displacement field closely matches the IGA reference solution. These results demonstrate that P-SNet can accurately predict the solution of the Scordelis–Lo roof with unstructured ASUTS meshes, validating its effectiveness as an energy-based solver for singleinstance shell problems.

In this section, numerical examples are presented to assess the proposed method. The first example verifies P-SNet as a data-free energy-based solver using the classical Scordelis–Lo roof benchmark. The remaining examples evaluate O-SNet for operator learning on three shell structures with increasing geometric complexity: a roof with extraordinary points, a B-pillar, and a plane nose. These tests examine the accuracy and generalization capability of the proposed method.

5.1 Solution learning with P-SNet To verify the performance of P-SNet, we start with the Scordelis–Lo roof problem, which has been widely used as a benchmark for shells. The roof is subjected to gravity and supported at the two curved edges. Three consecutively refined meshes

13

(a) Geometry and boundary conditions

(b) Mesh 0

(c) Mesh 1

(d) Mesh 2

Fig. 10: Scordelis–Lo roof and its unstructured Bézier meshes. (a) Problem settings: radius R = 25 m, length L = 50 m, thickness t = 0.25 m, uniform gravity load g = 90 N/m2 , Young’s modulus E = 432 MPa, and Poisson’s ratio ν = 0.0. (b) Initial mesh with 48 elements and 124 control points. (c) Refined mesh with 120 elements and 432 control points. (d) Refined mesh with 408 elements and 1624 control points.

4.5e-2 0.0

2.2e-3

0.0 -3.0e-1

-1.3e-3 (a) IGA result

(b) Prediction by P-SNet

(c) Error

Fig. 11: Displacement (z component) comparison for the Scordelis–Lo roof on Mesh 2.

14

Table 1: Free-edge midpoint displacement comparison. umid (IGA) umid (P-SNet) |∆umid | -0.25716 -0.25723 0.00007 -0.29781 -0.29780 0.00001 -0.30038 -0.30242 0.00204

Mesh Mesh 0 Mesh 1 Mesh 2

5.2.2 Roof

Table 2: Parameters of KL shells. Parameters E ν t

Meaning Young’s modulus Poisson’s ratio thickness

We begin with a simple roof problem. The geometric parameters are identical to those of the Scordelis–Lo roof, except that a different thickness is used and all four edges are fully clamped, as shown in Fig. 12. The roof is subjected to distributed surface loadings in the z -direction, generated using GRF with an amplitude factor of 1.5 × 108 . The roof is represented by ASUTS, where two EPs are intentionally introduced. A total of 1250 loading functions are generated, with 1000 used for training and the remaining 250 for testing. The branch net is implemented as an MLP with two hidden layers. We compare three network widths with 100, 300, and 500 neurons per hidden layer. All models are trained for 3000 epochs using the Adam optimizer [46] with the staged learning rate schedule shown in Fig. 13. The mean relative L2 error of the z -displacement is evaluated on both training and test sets at Gaussian quadrature points (GPs) and element corners (ECs). Only the GP values are used for training, whereas the EC errors are reported to assess the predicted control variables. Fig. 13 summarizes the training process and the best mean errors obtained with different network widths. As shown in Fig. 13a, the training and test errors decrease consistently as the number of epochs increases. The small gap between the training and test curves suggests that no obvious overfitting/underfitting occurs. Fig. 13b further compares different network widths. Increasing the width generally improves the prediction accuracy at both GPs and ECs, and the model with 500 neurons per hidden layer gives the lowest errors. Therefore, this model is selected for the following comparisons. The EC errors are slightly larger than the GP errors because the EC values are not directly used in the loss function. Table 3 reports the average L2 errors. To further illustrate the performance, three representative cases are selected from the test

Value 2 × 1011 0.3 0.02

5.2 Operator learning with O-SNet In this section, we investigate operator learning using O-SNet. Unlike P-SNet, which solves one fixed boundary value problem at a time, O-SNet learns the mapping from loading functions to structural responses. Several shell structures with EPs are studied to evaluate the performance of O-SNet.

5.2.1 Data generation and parameter settings In this work, the loading functions are generated by the mean-zero Gaussian random field (GRF):

f (x) ∼ γG 0, k (x, x′ ) ,

Relative error (%) 0.027 0.003 0.679

(21)

  ′ 2 ∥ where k (x, x′ ) = exp − ∥x−x is the Gaus2l2 sian kernel, l denotes the length scale controlling the smoothness of loading functions, and γ is the amplitude factor that amplifies the magnitude of loading functions. In this work, the length scale is set to l = 0.2, and the amplitude factor varies in different problems. The loading functions are sampled at the 16 Gaussian quadrature points of each Bézier element, which are also used as the inputs to SplineNet. For the linear Kirchhoff–Love shell tests (roof and B-pillar), the parameters are listed in Table 2. The structural displacement data are obtained from the in-house IGA solver [45] based on the ASUTS. The O-SNet model is implemented in PyTorch [44] and trained on an NVIDIA RTX 6000 Ada Generation GPU.

15

Relative L2 error

ux = uy = uz = 0 ux = uy = uz = 0

ux = uy = uz = 0

ux = uy = uz = 0

Epochs

(a)

(a)

Relative L2 error (%)

(b)

Network widths

Fig. 12: The roof problem for operator learning. (a) Problem setting with four fully clamped edges, subjected to z -direction loading over the entire surface, and (b) the corresponding Bézier mesh consisting of 120 elements and 432 control points.

(b)

Fig. 13: Training curve and network-width comparison for the roof problem. (a) Mean relative L2 error versus epochs, where solid and dashed lines denote test and training errors, respectively (with learning rates 10−3 , 10−4 , and 10−5 , each used for 1000 epochs). (b) Minimum relative L2 errors for different network widths, evaluated at Gaussian quadrature points (GPs) and element corners (ECs).

dataset corresponding to the first (Case 1), middle (Case 2), and last sample (Case 3). The predicted z -displacement fields are compared with the IGA reference solutions in Fig. 14. The relative L2 errors at GPs and ECs are 1.27% and 1.46% for Case 1, 1.72% and 2.20% for Case 2, and 1.90% and 2.95% for Case 3, respectively. The predictions closely match the reference solutions. The larger EC errors are a result of the fact that EC values are not directly enforced during training.

Table 3: Average relative L2 errors of the z direction displacement in the roof problem with 500 neurons. Error type GPs ECs

16

Training set 1.19% 2.21%

Test set 1.66% 2.56%

Case 1 5.1e-2

1.1e-3

0.0 0.0

-6.1e-3

-7.1e-2

Case 2 5.1e-2

2.5e-4 0.0

0.0

-7.1e-2

-8.0e-4

Case 3 3.5e-2

6.8e-4

0.0

0.0 -4.1e-2

-9.3e-4

(a) IGA result

(b) Prediction by O-SNet

(c) Error

Fig. 14: Comparison of z -displacement fields for three representative test cases of the roof problem.

17

Relative L2 error

5.2.3 B-pillar We next investigate O-SNet on a complex engineering structure, the B-pillar represented by ASUTS. As shown in Fig. 15, both ends of the structure are fully clamped. Distributed surface loadings are applied in the y -direction. The corresponding Bézier mesh is shown in Fig. 16. A total of 1250 loading functions are generated by GRF using an amplitude factor of 1.5 × 107 , with 1000 cases for training and 250 cases for testing. The branch net is an MLP with four residual blocks [47]. We also compare three network widths with 250, 500, and 750 neurons per hidden layer. Fig. 17 summarizes the training curves and the influence of network width. The mean relative L2 error of the y -displacement is evaluated on both training and test sets at GPs and ECs. The 750neuron network is selected, and its average errors are reported in Table 4. Representative cases (the first, middle, and last sample in the test dataset) are shown in Fig. 18. The relative L2 errors at GPs and ECs are

Epochs (a)

Relative L2 error (%)

ux = uy = uz = 0

Network widths (b)

Fig. 17: Training curve and the network width comparison for the B-pillar problem. (a) Mean relative L2 error versus epochs, where solid and dashed lines denote test and training errors, respectively (with learning rates 10−3 , 10−4 , and 10−5 , each used for 1000 epochs). (b) Minimum relative L2 errors for different network widths, evaluated at GPs and ECs.

ux = uy = uz = 0

Fig. 15: B-pillar model with two ends fully clamped, subjected to y -direction loading over the entire surface.

Table 4: Average relative errors of the y -direction displacement in the B-pillar problem with 750 neurons. Error type GPs ECs

Fig. 16: Bézier mesh of the B-pillar, consisting of 2881 elements and 3507 control points.

18

Training set 1.68% 2.38%

Test set 1.90% 2.51%

Case 1 2.4e-4

1.7e-2

1.0e-2 0.0

-1.1e-4

0.0

Case 2 1.6e-4

4.9e-3

0.0

0.0 -1.9e-4

-1.7e-3

Case 3

1.3e-4

0.0

0.0 -1.0e-2

-3.6e-4

-2.8e-2

(a) IGA result

(b) Prediction by O-SNet

(c) Error

Fig. 18: Comparison of y -displacement fields for three representative test cases of the B-pillar problem.

19

Relative L2 error

0.70% and 1.08% for Case 1, 2.55% and 4.05% for Case 2, and 0.76% and 0.93% for Case 3, respectively. These results show that O-SNet maintains accurate predictions for a large-scale engineering shell structure with complex geometry.

5.2.4 Geometrically nonlinear case: plane nose In the last case, we study O-SNet on a realworld engineering structure, the plane nose. This example is governed by the geometrically nonlinear Kirchhoff–Love shell formulation [3, 48], which introduces nonlinear load–displacement behavior and thus poses a more challenging operator learning task. As shown in Fig. 19, the circular edge on the left side, corresponding to the yz -plane cut, is fully clamped. Distributed surface loadings are applied over the entire surface in the x-direction. The corresponding Bézier mesh is shown in Fig. 20. Similar to the roof and B-pillar cases, a total of 1250 loading functions are generated using GRF

Epochs (a)

Relative L2 error (%)

ux = uy = uz = 0

Network widths (b)

Fig. 21: Training curve and the network width comparison for the geometrically nonlinear plane nose problem. (a) Mean relative L2 error versus epochs, where solid and dashed lines denote test and training errors, respectively (with learning rates 10−3 , 10−4 , and 10−5 , each used for 1000 epochs). (b) Minimum relative L2 errors for different network widths, evaluated at GPs and ECs.

Fig. 19: Plane nose model with the left circular edge fully clamped, subjected to x-direction loading over the entire surface.

with an amplitude factor of 5 × 108 . The dataset is split into 1000 training cases and 250 test cases. The displacement responses are obtained from the isogeometric solver tIGAr [49]. Unlike the linear Kirchhoff–Love shell cases, the plane nose uses E = 2.06 × 1011 , ν = 0.3, and t = 0.03. The branch net adopts an MLP with four residual blocks. We compare three network widths with 250, 500, and 750 neurons per hidden layer.

Fig. 20: Bézier mesh of the plane nose, consisting of 3011 elements and 3141 control points.

20

Fig. 21 summarizes the training curves and the influence of the network width for the nonlinear case. The errors are evaluated at both GPs and ECs. The average errors of the best case are reported in Table 5.

Table 5: Average relative errors of the x-direction displacement in the plane nose problem with 750 neurons. Error type Training set Test set GPs 1.44% 3.23% ECs 2.37% 3.93%

Case 1 0.0

2.5e-4 0.0

-1.0e-2

-4.5e-4

-1.9e-2

Case 2 1.6e-2

2.5e-4

1.0e-2 0.0

0.0

-2.1e-4

Case 3 0.0

2.5e-4

0.0 -1.0e-2

-1.8e-2

-3.9e-4

(a) IGA result

(b) Prediction by O-SNet

(c) Error

Fig. 22: Comparison of x-displacement fields for three representative test cases of the plane nose problem.

21

Representative cases (the first, middle, and last sample in the test dataset) are presented in Fig. 22. The relative L2 errors at GPs and ECs are 1.94% and 2.00% for Case 1, 0.70% and 1.11% for Case 2, and 1.55% and 1.64% for Case 3, respectively. These results demonstrate that O-SNet can learn the load-response operator for geometrically nonlinear shell problems with complex geometries.

CAD and FEA. John Wiley & Sons, Chichester, UK (2009). https://doi.org/10.1002/ 9780470749081 [2] Reddy, J.N.: An Introduction to the Finite Element Method, 2nd edn. McGraw-Hill, New York (1993) [3] Kiendl, J., Bletzinger, K.-U., Linhard, J., Wüchner, R.: Isogeometric shell analysis with Kirchhoff–Love elements. Computer methods in applied mechanics and engineering 198(4952), 3902–3914 (2009) https://doi.org/10.1 016/j.cma.2009.08.013

6 Conclusion and future work In this work, we propose a novel isogeometric deep learning method, SplineNet, which unifies scientific deep learning with CAD/CAE integration. Bézier extraction provides a unified mechanism for embedding ASUTS into neural networks. SplineNet has been applied to both solution learning (P-SNet) and operator learning (O-SNet). We first verified P-SNet on the Scordelis– Lo roof benchmark by employing the KL shell energy to solve a single-instance problem. This physics-informed version highlights the feasibility of using a spline-based neural representation for shell analysis. Building on this foundation, we further introduced O-SNet, where the displacement fields are instantly predicted given unseen loading functions for complex geometries. In contrast to sampling-based operator learning methods, our method operates on control variables. In the future, a particularly promising direction is to explore the proposed method in gradientbased shape optimization, where the position of control points can be treated as design variables within a CAD/CAE-integrated loop. Extending the approach to more general nonlinear settings (e.g., advanced constitutive models, contact, and buckling) and improving computational efficiency will further enhance its applicability to practical engineering design.

[4] Nagy, A.P., Abdalla, M.M., Gürdal, Z.: Isogeometric sizing and shape optimisation of beam structures. Computer Methods in Applied Mechanics and Engineering 199(1720), 1216–1230 (2010) https://doi.org/10.1 016/j.cma.2009.12.010 [5] Beck, J., Tamellini, L., Tempone, R.: IGAbased multi-index stochastic collocation for random PDEs on arbitrary domains. Computer Methods in Applied Mechanics and Engineering 351, 330–350 (2019) https://do i.org/10.1016/j.cma.2019.03.042 [6] Wang, Y., Bai, J., Lin, Z., Wang, Q., Anitescu, C., Sun, J., Eshaghi, M.S., Gu, Y., Feng, X.-Q., Zhuang, X., et al.: Artificial intelligence for partial differential equations in computational mechanics: a review. Applied Mechanics Reviews, 1–81 (2026) https://doi.org/10.1115/1.4071710 [7] Raissi, M., Perdikaris, P., Karniadakis, G.E.: Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics 378, 686–707 (2019) ht t p s : //doi.org/10.1016/j.jcp.2018.10.045

Acknowledgments. S. Luo and X. Wei are partially supported by National Natural Science Foundation of China (No. 12494550/12494555 and No. 12571408).

[8] E, W., Yu, B.: The deep Ritz method: a deep learning-based numerical algorithm for solving variational problems. Communications in Mathematics and Statistics 6(1), 1–12 (2018) https://doi.org/10.1007/s40304-018-0127-z

References [1] Cottrell, J.A., Hughes, T.J., Bazilevs, Y.: Isogeometric Analysis: Toward Integration of

[9] Samaniego, E., Anitescu, C., Goswami, S.,

22

[16] Berrone, S., Canuto, C., Pintore, M., Sukumar, N.: Enforcing Dirichlet boundary conditions in physics-informed neural networks and variational physics-informed neural networks. Heliyon 9(8), 18820 (2023) h t t p s : //doi.org/10.1016/j.heliyon.2023.e18820

Nguyen-Thanh, V.M., Guo, H., Hamdia, K., Zhuang, X., Rabczuk, T.: An energy approach to the solution of partial differential equations in computational mechanics via machine learning: concepts, implementation and applications. Computer Methods in Applied Mechanics and Engineering 362, 112790 (2020) https://doi.org/10.1016/j.cm a.2019.112790

[17] Saha, S., Gan, Z., Cheng, L., Gao, J., Kafka, O.L., Xie, X., Li, H., Tajdari, M., Kim, H.A., Liu, W.K.: Hierarchical deep learning neural network (HiDeNN): an artificial intelligence (AI) framework for computational science and engineering. Computer Methods in Applied Mechanics and Engineering 373, 113452 (2021) https://doi.org/10.1016/j.cm a.2020.113452

[10] Kharazmi, E., Zhang, Z., Karniadakis, G.E.: hp-VPINNs: variational physics-informed neural networks with domain decomposition. Computer Methods in Applied Mechanics and Engineering 374, 113547 (2021) https://doi.org/10.1016/j.cma.2020.113547

[18] Zhang, L., Cheng, L., Li, H., Gao, J., Yu, C., Domel, R., Yang, Y., Tang, S., Liu, W.K.: Hierarchical deep-learning neural networks: finite elements and beyond. Computational Mechanics 67(1), 207–230 (2021) https://do i.org/10.1007/s00466-020-01928-9

[11] Basir, S.: Investigating and mitigating failure modes in physics-informed neural networks (PINNs). Communications in Computational Physics 33(5), 1240–1269 (2023) https://do i.org/10.4208/cicp.OA-2022-0239 [12] Chickering, K.R.: A quasilinear algorithm for computing higher-order derivatives of deep feed-forward neural networks. arXiv preprint arXiv:2412.09752 (2024) https://doi.org/10 .48550/arXiv.2412.09752 arXiv:2412.09752

[19] Liu, Y., Park, C., Lu, Y., Mojumder, S., Liu, W.K., Qian, D.: HiDeNN-FEM: a seamless machine learning approach to nonlinear finite element analysis. Computational mechanics 72(1), 173–194 (2023) https://doi.org/10.1 007/s00466-023-02293-z

[13] Fuhg, J.N., Bouklas, N.: The mixed deep energy method for resolving concentration features in finite strain hyperelasticity. Journal of Computational Physics 451, 110839 (2022) https://doi.org/10.1016/j.jcp.2021.1 10839

[20] Lu, Y., Li, H., Zhang, L., Park, C., Mojumder, S., Knapik, S., Sang, Z., Tang, S., Apley, D.W., Wagner, G.J., et al.: Convolution hierarchical deep-learning neural networks (C-HiDeNN): finite elements, isogeometric analysis, tensor decomposition, and beyond. Computational Mechanics 72(2), 333–362 (2023) https://doi.org/10.1007/s0 0466-023-02336-5

[14] Bai, J., Lin, Z., Wang, Y., Wen, J., Liu, Y., Rabczuk, T., Gu, Y., Feng, X.-Q.: Energybased physics-informed neural network for frictionless contact problems under large deformation. Computer Methods in Applied Mechanics and Engineering 437, 117787 (2025) https://doi.org/10.1016/j.cma.2025.1 17787

[21] Škardová, K., Daby-Seesaram, A., Genet, M.: Finite element neural network interpolation: Part I—interpretable and adaptive discretization for solving PDEs. Computational Mechanics, 1–21 (2025) https://doi.or g/10.1007/s00466-025-02677-3

[15] Wang, S., Teng, Y., Perdikaris, P.: Understanding and mitigating gradient flow pathologies in physics-informed neural networks. SIAM Journal on Scientific Computing 43(5), 3055–3081 (2021) https://doi.org/10.1137/20m1318043

[22] Daby-Seesaram, A., Škardová, K., Genet, M.: Finite element neural network interpolation: Part II—hybridisation with the proper generalised decomposition for non-linear surrogate 23

360/nso/20240001

modelling. Computational Mechanics, 1–26 (2025) https://doi.org/10.1007/s00466-025-0 2676-4

[30] Wu, C., Zhu, M., Tan, Q., Kartha, Y., Lu, L.: A comprehensive study of non-adaptive and residual-based adaptive sampling for physicsinformed neural networks. Computer Methods in Applied Mechanics and Engineering 403, 115671 (2023) https://doi.org/10.1016/ j.cma.2022.115671

[23] Li, Z., Kovachki, N., Azizzadenesheli, K., Liu, B., Bhattacharya, K., Stuart, A., Anandkumar, A.: Fourier neural operator for parametric partial differential equations. In: International Conference on Learning Representations (2021). https://doi.org/10.48550/arXiv .2010.08895

[31] Gu, Y., Zhang, C., Zhang, P., Golub, M.V., Yu, B.: Enriched physics-informed neural networks for 2D in-plane crack analysis: theory and MATLAB code. International Journal of Solids and Structures 276, 112321 (2023) ht tps://doi.org/10.1016/j.ijsolstr.2023.112321

[24] Li, Z., Huang, D.Z., Liu, B., Anandkumar, A.: Fourier neural operator with learned deformations for PDEs on general geometries. Journal of Machine Learning Research 24(388), 1–26 (2023)

[32] Visser, C., Heinlein, A., Giovanardi, B.: PACMANN: point adaptive collocation method for artificial neural networks. Computer Methods in Applied Mechanics and Engineering 452, 118723 (2026) https://doi.org/10.1 016/j.cma.2025.118723

[25] Li, Z., Kovachki, N., Azizzadenesheli, K., Liu, B., Bhattacharya, K., Stuart, A., Anandkumar, A.: Neural operator: graph kernel network for partial differential equations. arXiv preprint arXiv:2003.03485 (2020) ht t p s : / / doi.org/10.48550/arXiv.2003.03485 arXiv:2003.03485

[33] Gasick, J., Qian, X.: Isogeometric neural networks: a new deep learning approach for solving parameterized partial differential equations. Computer Methods in Applied Mechanics and Engineering 405, 115839 (2023) https://doi.org/10.1016/j.cma.2022.1 15839

[26] Lu, L., Jin, P., Pang, G., Zhang, Z., Karniadakis, G.E.: Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. Nature machine intelligence 3(3), 218–229 (2021) ht tps://doi.org/10.1038/s42256-021-00302-5

[34] Zhang, L., Park, C., Lu, Y., Li, H., Mojumder, S., Saha, S., Guo, J., Li, Y., Abbott, T., Wagner, G.J., et al.: Isogeometric convolution hierarchical deep-learning neural network: isogeometric analysis with versatile adaptivity. Computer Methods in Applied Mechanics and Engineering 417, 116356 (2023) https://doi.org/10.1016/j.cm a.2023.116356

[27] He, J., Koric, S., Kushwaha, S., Park, J., Abueidda, D., Jasiuk, I.: Novel DeepONet architecture to predict stresses in elastoplastic structures with variable complex geometries and loads. Computer Methods in Applied Mechanics and Engineering 415, 116277 (2023) https://doi.org/10.1016/j.cm a.2023.116277 [28] Park, J., Kang, N.: Point-DeepONet: predicting nonlinear fields on non-parametric geometries under variable load conditions. Neural Networks 198, 108560 (2026) https: //doi.org/10.1016/j.neunet.2026.108560

[35] Zhang, L., Park, C., Hughes, T.J., Liu, W.K.: Multi-patch isogeometric convolution hierarchical deep-learning neural network. Computer Methods in Applied Mechanics and Engineering 434, 117582 (2025) https://doi. org/10.1016/j.cma.2024.117582

[29] Chen, G., Liu, X., Meng, Q., Chen, L., Liu, C., Li, Y.: Learning neural operators on Riemannian manifolds. National Science Open 3(6), 20240001 (2024) https://doi.org/10.1

[36] Tresckow, M., Ion, I.G., Loukrezis, D.: Multipatch neural solver with isogeometric mappings for partial differential equations on 24

a.2023.116523

computer-aided design domains. Engineering with Computers 42(4), 116 (2026) h t t p s : //doi.org/10.1007/s00366-026-02351-z

[44] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., et al.: PyTorch: an imperative style, high-performance deep learning library. In: Advances in Neural Information Processing Systems, vol. 32, pp. 8024– 8035 (2019)

[37] Wei, X., Li, X., Qian, K., Hughes, T.J., Zhang, Y.J., Casquero, H.: Analysis-suitable unstructured T-splines: multiple extraordinary points per face. Computer Methods in Applied Mechanics and Engineering 391, 114494 (2022) https://doi.org/10.1016/j.cm a.2021.114494

[45] Sheng, J., Wei, X.: Isogeometric topology optimization of thin-walled structures with complex design domains. Computer Methods in Applied Mechanics and Engineering 444, 118114 (2025) https://doi.org/10.1016/j.cm a.2025.118114

[38] Rogers, D.F.: An Introduction to NURBS: with Historical Perspective. Morgan Kaufmann, San Francisco, CA, USA (2001) [39] Sederberg, T.W., Zheng, J., Bakenov, A., Nasri, A.: T-splines and T-NURCCs. ACM transactions on graphics (TOG) 22(3), 477– 484 (2003) https://doi.org/10.1145/1201775. 882295

[46] Kingma, D.P., Ba, J.: Adam: a method for stochastic optimization. In: International Conference on Learning Representations (2015). https://doi.org/10.48550/arXiv .1412.6980 . https://arxiv.org/abs/1412.6980

[40] Reif, U.: A refineable space of smooth spline surfaces of arbitrary topological genus. Journal of Approximation Theory 90(2), 174–199 (1997) https://doi.org/10.1006/jath.1996.30 79

[47] He, K., Zhang, X., Ren, S., Sun, J.: Deep residual learning for image recognition. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 770–778 (2016). https://doi.org/10.1109/CV PR.2016.90

[41] Wei, X., Zhang, Y.J., Toshniwal, D., Speleers, H., Li, X., Manni, C., Evans, J.A., Hughes, T.J.: Blended B-spline construction on unstructured quadrilateral and hexahedral meshes with optimal convergence rates in isogeometric analysis. Computer Methods in Applied Mechanics and Engineering 341, 609–639 (2018) https://doi.org/10.1016/j.cm a.2018.07.013

[48] Kiendl, J., Hsu, M.-C., Wu, M.C., Reali, A.: Isogeometric Kirchhoff–Love shell formulations for general hyperelastic materials. Computer Methods in Applied Mechanics and Engineering 291, 280–303 (2015) https://do i.org/10.1016/j.cma.2015.03.010 [49] Kamensky, D., Bazilevs, Y.: tIGAr: automating isogeometric analysis with FEniCS. Computer Methods in Applied Mechanics and Engineering 344, 477–498 (2019) https://do i.org/10.1016/j.cma.2018.10.002

[42] Scott, M.A., Borden, M.J., Verhoosel, C.V., Sederberg, T.W., Hughes, T.J.: Isogeometric finite element data structures based on Bézier extraction of T-splines. International Journal for Numerical Methods in Engineering 88(2), 126–156 (2011) https://doi.org/10.1002/nm e.3167 [43] Casquero, H., Mathews, K.D.: Overcoming membrane locking in quadratic NURBSbased discretizations of linear Kirchhoff–Love shells: CAS elements. Computer Methods in Applied Mechanics and Engineering 417, 116523 (2023) https://doi.org/10.1016/j.cm 25