Neural-Schwarz Tiling for Geometry-Universal PDE Solving at Scale
Paolo Secchi1
Daniel S. Balint1
Marco Maurizi2
1
arXiv:2605.12343v1 [cs.LG] 12 May 2026
2
Imperial College London Italian Institute of Artificial Intelligence [email protected] [email protected]
Abstract Most learned PDE solvers follow a global-surrogate paradigm: a neural operator is trained to map full problem descriptions to full solution fields for a prescribed distribution of geometries, boundary conditions, and coefficients. This has enabled fast inference within fixed problem families, but limits reuse across new domains and makes large-scale deployment dependent on expensive problem-specific data generation. We introduce NEST (Neural-Schwarz Tiling), a local-to-global framework that shifts learning from full-domain solution operators to reusable local physical solvers. The central premise is that, although global PDE solutions depend on geometry, scale, and boundary conditions, the physical response on small neighborhoods can be learned locally and composed into global solutions through classical domain decomposition. NEST learns a neural operator on minimal voxel patches (3 × 3 × 3) with diverse local geometries and boundary/interface data. At inference time, an unseen voxelized domain is tiled into overlapping patches, the learned local solver is applied patchwise, and global consistency is enforced through iterative Schwarz coupling with partition-of-unity assembly. In this way, generalization is shifted from a monolithic neural model to the combination of local physics learning and algorithmic global assembly. We instantiate NEST on nonlinear static equilibrium in compressible neo-Hookean solids and evaluate it on large, geometrically complex 3D domains far outside the scale of the training patches. Our results show that local neural building blocks, coupled through Schwarz iteration, offer a reusable local-training path toward scalable learned PDE solvers that generalize across domain size, shape, and boundary-condition configurations.
1
Introduction
Machine learning has transformed domains in which data admit a universal representation. In language, large-scale models trained on token sequences generalize across a broad range of tasks, domains, and contexts [1]. This success rests on a common compositional substrate: sentences, documents, code, and instructions can all be represented as sequences of tokens, enabling a shared modeling paradigm across heterogeneous problems. No comparable foundation has emerged for physical simulation. Many important problems in science and engineering—fluid flow, heat transfer, elasticity, fracture, electromagnetism, and coupled multiphysics—are governed by partial differential equations (PDEs). These describe fields evolving over geometries, boundary conditions, material distributions, and forcing terms that vary dramatically from one problem to another. Unlike language, PDE problems are tied to their spatial domain, discretization, boundary representation, and physical coefficients. The central object of physical Preprint.
simulation is thus not a single task, but a family of operators whose structure changes with the geometry and physics of the problem. This lack of universality is a major bottleneck for scientific machine learning (SciML) [2–4]. Classical numerical solvers such as finite element modeling (FEM) are general and reliable, but their cost grows rapidly with problem size, resolution, and the number of repeated solves required by downstream workflows like inverse design, uncertainty quantification, digital twins, or design-space exploration [5, 6]. Learned surrogates can be fast after training, but are typically specialized: a model trained on one geometry, discretization, or coefficient distribution often fails to transfer outside that regime. Consequently, today’s PDE learning pipelines replace one expensive computation with another: generating large datasets and retraining specialized models for each problem family. Physicsinformed neural networks (PINNs) [7] sidestep data generation by embedding the PDE residual into the loss, but each new geometry, boundary condition, or coefficient field requires re-optimizing the network from scratch, so generalization across problems is essentially absent. The key challenge is therefore not to accelerate a fixed benchmark, but to learn reusable physical computation that generalizes across domains and scales with problem size [8–10]. A general learned solver should not need retraining whenever the geometry, boundary conditions, or domain size change; it should instead learn local physical laws composable into global solutions on previously unseen structures. This is the central question of this work: can we learn a small, reusable neural building block for PDE solution, and assemble it algorithmically to solve large, geometrically arbitrary physical systems? SciML [8] often formulates simulation as learning a parameterized surrogate operator Gθ ≈ G † : X → U, where G † maps functional inputs — such as boundary data or initial conditions — to solution fields. Neural operators [11] have emerged as a principled class of models for approximating infinite-dimensional mappings, typically of the form Gθ = Q ◦ σ (WL + KL ) ◦ · · · ◦ σ (W1 + K1 ) ◦ P, where P, Q are lifting and projection maps, Kℓ are learned integral kernel operators, Wℓ are local linear maps, and σ is a pointwise nonlinearity. They act on functions rather than fixed finitedimensional vectors and admit guarantees of discretization invariance and universal approximation [11]. However, these guarantees apply to operators between function spaces on fixed domains. In practice, a learned operator is trained for a prescribed geometry class, boundary representation, and coefficient distribution. While several architectures improve empirical transfer across meshes, resolutions, and parameter fields, the learning problem remains global: the model approximates an entire PDE solution operator over a given problem family. Changing the domain, enlarging the computational region, or modifying boundary-condition structure typically requires new data generation, retraining, or architecture-specific adaptation, limiting neural operators as reusable computational primitives. In this work, we take a different view. Rather than learning a global operator for each class of PDE problems, we learn a local physical solver once and reuse it compositionally. The core idea is to train a neural operator on small voxelized subdomains, exposing it to a rich distribution of boundary/interface conditions. At test time, arbitrary geometries are voxelized into overlapping local neighborhoods, and global consistency is enforced through Schwarz iteration [12, 13]. Generalization thus shifts from the neural model alone to the combination of a learned local solution operator and a classical domain-decomposition procedure. Our main contributions. We introduce NEST (Neural-Schwarz Tiling), a local-to-global framework for solving PDEs on previously unseen domains without retraining. NEST learns a reusable neural building block on small voxel neighborhoods and assembles these blocks through overlapping Schwarz iteration to obtain global solutions on large, geometrically complex structures. This construction separates local physics learning from global domain assembly: the neural operator approximates local solution updates, while the Schwarz procedure propagates information, enforces compatibility across overlapping subdomains, and enables scaling to domains much larger than those seen during training. Our contributions are threefold. First, we propose a neural-Schwarz formulation in which a single local operator, trained on small voxel patches with diverse boundary conditions, can be reused across arbitrary voxelized geometries. Second, we show that global PDE solutions can be recovered by 2
composing these learned local updates through Schwarz iteration, enabling inference on domains whose size, shape, and topology differ from the training setting. Third, we demonstrate that this paradigm provides a reusable local-training alternative to training specialized global surrogates for each new geometry class or problem family, opening a path toward learned PDE solvers that behave less like task-specific predictors and more like reusable computational building blocks.
2
Related Work
Neural operators for parametric PDEs. Operator learning has emerged as a central paradigm for scientific machine learning, aiming to approximate mappings between infinite-dimensional function spaces rather than finite-dimensional input-output pairs. Early architectures such as DeepONet [14] and the Fourier Neural Operator (FNO) [15] demonstrated that neural networks can learn solution operators for families of parametric PDEs, enabling fast inference after training. Subsequent work developed a broader theory of neural operators, including universal approximation and discretizationinvariance results [11], as well as physics-informed variants that incorporate PDE residuals or physical constraints during training [16]. These methods have significantly expanded the scope of learned PDE solvers, but they typically learn a global operator associated with a prescribed problem family. In particular, the training distribution usually fixes or strongly constrains the domain class, boundary representation, coefficient statistics, and input-output structure. As a result, changing the geometry or boundary-condition structure often requires new data, retraining, or architecture-specific adaptation. Learning PDEs on meshes and irregular geometries. A complementary line of work addresses the limitations of grid-based operator learning by representing physical systems on meshes, graphs, point clouds, or geometry-aware latent spaces. MeshGraphNets [17] introduced a graph-network simulator that performs message passing on adaptive mesh representations and demonstrated strong performance on deformable solids, fluids, and cloth. Graph-based and geometry-informed neural operators extend this idea to operator learning on irregular domains, using graph kernels, point-cloud representations, signed-distance functions, or learned mappings between irregular and regular domains [18]. More recently, transformer architectures have advanced PDE modeling; Transolver [19, 20] handles complex geometries via physics-aware attention mechanisms [21], whereas HAMLET [22] tackles parametric problems using graph attention. These approaches are important steps toward geometry-aware and parametric learned simulation. However, they still primarily follow a global learning paradigm: the model is trained to map from a full geometry and its associated physical inputs to a full solution field. Thus, generalization to geometries, topologies, domain sizes, or boundary-condition configurations far outside the training distribution remains difficult. Moreover, scaling to very large domains can require substantially larger models, more memory, or additional training data. In contrast, our goal is not to learn a single global operator over an entire geometry distribution, but to learn a reusable local solver that can be composed algorithmically into global solutions. Domain decomposition and Schwarz methods. Domain decomposition methods provide a classical route for solving PDEs on large or complex domains by decomposing the global problem into smaller subproblems. Among these, overlapping Schwarz methods iteratively solve local problems on overlapping subdomains and exchange interface information until a globally consistent solution is obtained [23]. These methods are attractive because they expose locality, parallelism, and modularity: global solution structure emerges from repeated local solves and communication across overlaps. Recent work has begun to combine domain decomposition with machine learning, for example by learning interface conditions, accelerating local solvers [24], or coupling neural operators with Schwarz-type iterations [25–27]. Such hybrid approaches suggest that classical numerical algorithms can provide useful structure for learned PDE solvers, especially when extrapolation beyond a fixed training domain is required. However, existing studies focus primarily on 2D domains, and extending these frameworks to 3D problems is not straightforward. Current methods often rely on learning non-minimal subdomains, which makes it difficult to tile complex 3D geometries [26]. Furthermore, the reliance on these specific subdomain configurations is computationally intensive and difficult to tailor to diverse problem sets such as PDE parameters. Positioning of NEST. NEST builds on this local-to-global perspective, but differs from existing neural operators and neural domain-decomposition approaches in its use of a small, reusable voxel3
level building block. Instead of training a model on full domains [28, 18, 29] or geometry-specific decompositions [24], NEST trains a local neural operator on small voxel neighborhoods with diverse boundary data and interface conditions. At inference time, arbitrary voxelized domains are tiled into overlapping neighborhoods, the same learned operator is applied repeatedly, and Schwarz iteration enforces global consistency. This separates the learning problem from the global geometry: the neural model learns local physical response, while the Schwarz procedure performs global assembly. As a result, NEST aims to provide a path toward learned PDE solvers that are reusable across domain size, shape, topology, and boundary-condition configurations, rather than specialized to a fixed benchmark or geometry class.
3
Preliminaries and Problem Statement
We consider a general boundary value problem on a bounded domain Ω ⊂ Rd with boundary ∂Ω: given a (possibly nonlinear) differential operator Lϕ parameterized by ϕ ∈ Φ (e.g. material constants, coefficients) and a prescribed Dirichlet boundary condition g : ∂Ω → Rm , find u : Ω → Rm such that Lϕ [u] = 0 in Ω, u = g on ∂Ω, (1) where the Dirichlet boundary condition is imposed on the entire boundary ∂Ω. Our goal is to learn a surrogate solution operator Sθ : (Ω, g) 7→ u (2) that generalizes across both g and the geometry Ω ∈ D, where D is a family of admissible domains rather than a fixed reference shape. Crucially, training Sθ must be cheap enough that the surrogate can be easily retrained for new operator parameters ϕ, enabling rapid tailoring to different PDE problems. To illustrate NEST, we instantiate (1) on quasilinear elliptic boundary value problems in divergence form, Lϕ [u] = − div σϕ (∇u) = 0 in Ω, u = g on ∂Ω, (3) where σϕ : Rm×d → Rm×d is a (generally nonlinear) constitutive map. Different choices of σϕ recover a wide range of stationary physical models such as heat conduction with temperaturedependent conductivity, nonlinear electrostatics, Darcy flow in heterogeneous media, and finite-strain solid mechanics. In this work, as a representative benchmark, we consider the nonlinear static equilibrium of a compressible neo-Hookean solid, with d = 3, m = 3, and ϕ = (λ, µ) the Lamé constants. Let u denote the displacement, F = I + ∇u the deformation gradient, and J = det F . The strain energy density ψϕ (F ) =
µ λ tr F ⊤ F − 3 − µ log J + (log J)2 2 2
(4)
yields the first Piola–Kirchhoff stress Pϕ (F ) = ∂ψϕ /∂F , and identifying σϕ (∇u) = Pϕ (I + ∇u) recovers the divergence-form structure of (3). In the absence of body forces and surface tractions, the strong form reduces to div Pϕ (I + ∇u) = 0 in Ω, u = g on ∂Ω, (5) which is the instance of (1) on which we benchmark NEST. Solutions are computed with the FEniCS finite-element library [30].
4
NEST — NEural-Schwarz Tiling
NEST is a two-stage framework for solving PDEs on large, geometrically complex, arbitrary domains by combining a learned local neural operator with a global Schwarz coupling strategy. Instead of training a neural operator directly on full-domain solutions for a prescribed family of geometries and boundary conditions, NEST learns a reusable local surrogate, Sθ , on a small canonical patch. The local operator maps patch-level problem data—including geometry and boundary/interface conditions—to the corresponding local solution field. At inference time, an arbitrary domain Ω ∈ D is decomposed into overlapping patches. The learned operator Sθ is applied independently on each patch, while neighboring predictions are coupled through an additive Schwarz iteration. Interface information is exchanged across overlapping regions, so that local solutions are progressively made consistent with one another and with the prescribed 4
global boundary conditions. In this way, NEST turns a neural operator trained only on local building blocks into a solver that can be applied to domains and boundary conditions that are much larger and more complex than those seen during training. The following subsections describe the local neural operator, the tiling of arbitrary domains into overlapping patches, and the Schwarz-based global coupling procedure. 4.1
Local Operator Learning
We restrict Sθ to domains carved out of a fixed canonical patch Ω0 = [0, 1]3 , which we discretize as a regular 3 × 3 × 3 grid of hexahedral cells {Ck }27 k=1 . This is the smallest cubic voxel patch containing a single fully interior cell, making it the minimal hexahedral patch that supports overlapping Schwarz with non-trivial interior updates. A training sample is a pair (Ω, g) obtained as follows. S Geometry sampling. Each training domain Ω = k∈A Ck ⊆ Ω0 is defined by an active set A ⊆ {1, . . . , 27} of solid cells. We draw A uniformly from the set of admissible configurations T = A ⊆ {1, . . . , 27} : A ̸= ∅, A ̸= {kc }, A is face-connected , (6) where kc is the central cell. Here, face-connected means that the graph (A, ∼), with k ∼ k ′ if and only if Ck and Ck′ share a full two-dimensional face, is connected. This guarantees that the solid is a single continuous component, while admitting arbitrary void configurations. Boundary data sampling. Given a sampled domain Ω, we prescribe the Dirichlet datum g on ∂Ω as a smooth, multiscale random field. Generalizing the procedure introduced in [27] to 3D, let x = (x1 , x2 , x3 ) ∈ Ω denote material coordinates and define the projection s(x) = (x1 + x2 + x3 )/3 together with its rescaled counterpart s̃(x) = (s(x) + δ)/ζ + s0 . The components of g, for i ∈ {1, 2, 3}, are then gi (x) = α
N X
h i 1 (i) (i) (i) (i) (i) a cos 2πn s̃(x) + b + c sin 2πn s̃(x) + d + e n n n n . (n + 1)κ n n=1 (i)
(i)
(i)
(i)
(7)
(i)
The amplitudes an , cn ∼ U[0.5, 1.0], phases bn , dn ∼ U[π/4, 5π/4], and biases en ∼ U[−0.25, 0.25] are sampled independently per component i and per mode n; N is the number of active modes, κ a fixed spectral-decay exponent, and α a global magnitude. The phase shift s0 ∼ U[0, 1], offset δ ∼ U[0, ζ], and zoom log ζ ∼ U [log ζmin , log ζmax ] are sampled once per domain. The pair (ζ, δ) selects a window of length 1/ζ from the field, so the patch sees boundary data consistent with what a 3 × 3 × 3 tile would receive inside a ζ-times-larger structure, therefore exposing the surrogate at training time to the multiscale gradients it must handle at Schwarz inference on generic domains. The field g is evaluated at the boundary nodes of the FEM discretization of Ω and assigned as Dirichlet data there, as illustrated in Figure 1. Architecture. We parameterize Sθ as a Graph Neural Operator (GNO) [31, 11] on the patch graph G = (V, E). The node set V consists of the corners of the 4 × 4 × 4 reference grid that are incident to at least one cell of Ω, and is partitioned into boundary nodes ∂V on ∂Ω and interior (active) nodes V \ ∂V. Edges E = {(i, j) : 0 < ∥xi − xj ∥ ≤ r} connect node pairs within a fixed radius r = 1.5 h, where h is the cell edge length. The Dirichlet datum g, prescribed on ∂V, is extended to every node by inverse-distance-squared interpolation [27], ∥xi − xj ∥−2 , −2 k∈∂V ∥xi − xk ∥
wij = P g̃(xi ) =
X
wij g(xj ) ∀i ∈ V \ ∂V,
(8)
g̃(xi ) = g(xi ) ∀i ∈ ∂V,
(9)
j∈∂V (0)
(0)
yielding the input feature hi = g̃(xi ) ∈ R3 at every node. A lifting MLP maps hi to a hidden width C, after which L kernel-integral layers update features by X 1 (ℓ+1) (ℓ) (ℓ) (ℓ) (ℓ) (ℓ) hi = hi + σ LNW (ℓ) hi + κ(ℓ) xj − xi , hi , hj ⊙ hj , (10) |N (i)| j∈N (i)
5
NEST prediction Training
Local Operator
Figure 1: Left: Assignment of displacement boundary conditions. A continuous displacement field (represented by blue streamlines) is evaluated at the boundary nodes of the voxelized domain to prescribe Dirichlet boundary conditions (red vectors). Right: A close-up highlights the local operator tiling geometries during inference.
where N (i) are the graph neighbours of i, κ(ℓ) is an MLP producing per-channel multiplicative weights from the relative position and the features of the two incident nodes, σ is the SiLU activation, (L) and LN denotes layer normalization. A node-wise projection MLP then maps hi to the prediction ŷ(xi ). We train two variants of Sθ that share this architecture but differ in their output target: a c 11 (xi ) ∈ displacement surrogate with ŷ(xi ) = û(xi ) ∈ R3 , and a gradient surrogate with ŷ(xi ) = ∇u R, the same logic can be potentially extended to the other gradient components. The displacement gradient can be then used to compute strain and stress fields. Predicting ∇u with a dedicated surrogate, rather than recovering it by differentiating û, is empirically less sensitive to small displacement errors. Training. We generate 15,000 samples and split them 90/10 into train and validation. Both surrogates are trained for 100 epochs with Adam [32] at learning rate 10−4 and the checkpoint with the lowest validation loss is retained. Inputs and outputs are channel-wise z-normalized using statistics computed on the training split. The gradient surrogate is trained with the MSE on its scalar output. The displacement surrogate is trained with the composite loss L = MSE(û, u) + MSE(∇û, ∇u). This term acts as a gradient-consistency regularizer, indeed, at Schwarz inference, the converged displacement prediction of each patch is fed into the gradient surrogate to produce the ∇u. 4.2
Schwarz Method
Having trained the local surrogate Sθ , we approximate the global solution on an arbitrary domain via an overlapping additive Schwarz method, employing Sθ as the local solver. We cover Ω with overlapping 3 × 3 × 3 patches {Ωp }P p=1 , each matching the canonical training geometry as illustrated int ext in Figure 1. We write the boundary of each patch as ∂Ωp = Γext p ∪ Γp . Here, Γp = ∂Ωp ∩ ∂Ω int inherits the macroscopic Dirichlet datum, and Γp is the interface with neighbouring patches. Let V denote the active nodes of Ω, and let Vp ⊆ V denote those of patch p. Finally, let {χp }P p=1 be a partition of unity over V such that each χp vanishes identically outside Vp . Starting from an initial iterate u(0) that matches g on ∂Ω, each Schwarz sweep n → n+1 proceeds, in parallel for every patch, by assembling local Dirichlet data and invoking the displacement surrogate, ( g gp(n) = u(n) Γint p
on Γext p , on Γint p , 6
û(n+1) = Sθ (Ωp , gp(n) ), p
(11)
followed by partition-of-unity assembly into a global field, u(n+1) =
P X
χp û(n+1) , p
u(n+1) ∂Ω = g.
(12)
p=1
Iteration stops once the relative change on interior active nodes, ∆(n+1) =
∥u(n+1) − u(n) ∥2, V\∂Ω , ∥u(n+1) ∥2, V\∂Ω
(13)
falls below a prescribed tolerance ε. Algorithm details can be found in Appendix D. Once the displacement iteration has converged to u⋆ , we recover ∇u on Ω by a single non-iterative pass of the gradient surrogate. For each patch p, we extract the converged displacement u⋆ ∂Ω on p the patch boundary and feed it as Dirichlet datum to a dedicated GNO trained, as described above, to predict one component of ∇u. Repeating the pass per component and assembling the patchwise predictions through the partition of unity {χp } yields the global gradient field.
5
Experiments
We design the experiments to evaluate three claims: (i) NEST trained on 3 × 3 × 3 patches generalizes to much larger, unseen, geometrically distinct domains without re-training, similarly to FEM; (ii) once trained, contrary to FEM, inference cost is largely independent of the complexity of the constitutive model; and (iii) local patch training provides a paradigm shift: from global training on full-domain solutions to reusable locally trained neural surrogates. All experiments use the displacement and gradient surrogates described in Paragraph 4.1. The local models are trained once and then frozen for all macro-scale evaluations. The test set contains four macro-scale geometries: two triply periodic minimal surface (TPMS) structures and two jet-engine brackets from the SimJEB dataset [33]. These geometries are not used during local training. Each geometry is voxelized at multiple resolutions n × n × n, where n denotes the number of voxels per axis. The results for n = 60 in Figure 2 highlight the generalization capabilities of NEST. Unless otherwise stated, NEST is run with the overlapping Schwarz iteration described in Section 4.2, using a convergence tolerance of ε = 10−5 for the relative iterate change in Eq. (13). An ablation isolating the role of Schwarz coupling is reported in Appendix B, while running time is detailed in Appendix C. Reference solutions are computed with the FEniCS finite-element solver. 5.1
Accuracy on unseen macro-scale geometries
Table 1 reports max-normalized mean-squared errors of NEST on unseen SimJEB bracket and TPMS geometries at increasing voxel resolutions, for displacement and for the 11-component of the displacement gradient; full error definitions are given in Appendix A. A single local model trained on 3 × 3 × 3 patches is used for all geometries and resolutions. No macro-scale data are used for training or fine-tuning. When multiple geometries belong to the same class, the reported values are averaged within that class. Table 1: Accuracy of NEST on unseen macro-scale geometries. A single local model trained on 3 × 3 × 3 patches is reused for all test geometries and resolutions. SimJEB brackets
TPMS
Resolution
emaxMSE u
emaxMSE ∇u11
emaxMSE u
emaxMSE ∇u11
30 45 60
8.147 × 10−3 6.670 × 10−3 6.929 × 10−3
3.102 × 10−2 2.566 × 10−2 2.096 × 10−2
5.574 × 10−3 2.539 × 10−3 1.715 × 10−3
3.376 × 10−3 5.155 × 10−3 5.997 × 10−3
NEST maintains stable displacement accuracy as the macro-scale resolution increases, despite being trained only on local voxel patches. The displacement error remains below 10−2 across all tested geometries and resolutions. Errors are higher on SimJEB brackets than on TPMS structures, which is 7
FEM ground truth ∇𝑢11
NEST prediction ∇𝑢11
FEM ground truth 𝑢
NEST prediction 𝑢
FEM ground truth ∇𝑢11
NEST prediction ∇𝑢11
FEM ground truth 𝑢
NEST prediction 𝑢
Figure 2: NEST vs. FEM on unseen macro-scale test geometries. The test set contains two SimJEB jet-engine brackets and two TPMS structures, voxelized at multiple resolutions, here shown at n = 60. NEST is trained only on random 3 × 3 × 3 local patches and is not fine-tuned on these geometries. expected because the bracket geometries contain sharper geometric features and less regular local topology. Nevertheless, the same local model remains applicable across both geometry classes, supporting the central premise that local neural building blocks can be composed through Schwarz iteration to solve geometrically distinct domains. 5.2
Comparison with a favorable global neural-operator baseline
We next compare NEST with a global neural-operator baseline in an intentionally favorable setting for the baseline. A standard global neural operator is trained to map full macro-scale geometries and boundary data directly to full solution fields. In practice, such models are known to be highly sensitive to the geometry distribution used for training: when trained on only a small number of macro-scale geometries and evaluated on substantially different unseen geometries, performance typically deteriorates dramatically. For this reason, we do not use the global baseline as a strict out-of-distribution geometry-generalization test. Instead, we ask a more conservative question: can NEST, trained only on local 3 × 3 × 3 patches and never on macro-scale solutions, match a global neural operator trained directly on the target macro-scale geometry class? We use Transolver [19] as a representative global neural-operator baseline. The comparison is deliberately favorable to Transolver: the model is trained on macro-scale TPMS and SimJEB geometries at resolution 603 , using full-domain FEM solutions. In contrast, NEST is trained only once on local voxel patches and is reused without fine-tuning on all macro-scale geometries and resolutions. Thus, while Transolver learns a geometry-specific full-domain surrogate, NEST learns a reusable local solver that is composed through Schwarz iteration at inference time. Transolver was trained for 1000 epochs with an initial learning rate of 10−3 , halved every 200 epochs. On the TPMS geometries at resolution 603 , Transolver reaches a displacement error of 1.740 × 10−3 and a gradient-component error of 9.392 × 10−3 for ∇u11 . On the same setting, NEST reaches a displacement error of 1.715 × 10−3 and a gradient-component error of 5.997 × 10−3 . Therefore, NEST matches the displacement accuracy of a global neural operator trained directly on macro-scale data, while producing a more accurate estimate of the displacement-gradient component. This comparison should be interpreted as an upper-bound-style baseline rather than as a standard out-of-distribution benchmark for global neural operators. Its purpose is to test whether local neural-Schwarz composition can be competitive with a global surrogate that has access to far more task-specific information during training. The result is notable because NEST does not require 8
macro-scale training data, does not observe the target geometries during training, and reuses the same local model across all tested geometry classes and resolutions. 5.3
Reusable local training versus geometry-specific global training
The previous comparison evaluates predictive accuracy, but it does not by itself capture the main data advantage of NEST. The key distinction is not only the cost of generating one training solve, but the type of information contained in the training set. Global neural operators are trained on full-domain examples tied to a particular macro-scale geometry distribution. Consequently, when the target geometry class, resolution, or boundary-condition structure changes substantially, new macro-scale FEM data may be required. In contrast, NEST is trained on local 3 × 3 × 3 patch problems. These local problems are not associated with a specific macro-scale geometry; they define reusable physical building blocks that can be composed through Schwarz iteration on different domains. To make this comparison concrete, we compare NEST and global neural-operator baselines under a similar data-generation budget. NEST is trained on 15,000 local patch solves, corresponding to 2.412 × 106 total degrees of freedom and 8.880 × 103 seconds of data-generation time. The global neural-operator baselines are trained on full macro-scale FEM solutions from TPMS and SimJEB geometries, corresponding to 2.575 × 106 total degrees of freedom and 9.062 × 103 seconds of datageneration time. Thus, the total amount of generated FEM data and wall-clock data-generation time are comparable across the two settings, but the resulting training sets have fundamentally different reuse properties. Table 2: Training-data comparison between local NEST training and global neural-operator training. The total generated degrees of freedom and data-generation wall time are comparable, but NEST uses local patch problems that can be reused across macro-scale geometries and resolutions, whereas the global baselines use full-domain examples tied to the target geometry families. Method
Training domains
# FEM solves
Tot DoFs
Data-gen. time [s]
NEST NOs
Local 3 × 3 × 3 patches Macro-scale geometries
15,000 13
2.412 × 106 2.575 × 106
8.880 × 103 9.062 × 103
For the global neural-operator baselines, we train Transolver and GNO on a macro-scale dataset consisting of 7 TPMS geometries, obtained from linear combinations of 8 primitive TPMS equations [34, 35], and 6 SimJEB bracket geometries. The dataset is split 90/10 into training and validation sets. For each architecture, two separate models are trained, one for the displacement field u and one for the gradient component ∇u11 . The models are trained for 500 epochs with an initial learning rate of 10−3 , halved every 100 epochs. Errors are evaluated on the test geometries shown in Figure 2 at resolution 603 . Table 3: Accuracy comparison between global neural-operator baselines and NEST at resolution 603 . The global baselines are trained on macro-scale data from the target geometry families, whereas NEST is trained only on local 3 × 3 × 3 patch problems and reused without macro-scale fine-tuning. Brackets
TPMS
Model
Parameters u model
emaxMSE ∥u∥
emaxMSE ∇u11
emaxMSE ∥u∥
emaxMSE ∇u11
Transolver GNO NEST
3.885 × 105 2.011 × 105 7.661 × 104
1.457E-02 8.601E-03 6.929E-03
3.932E-02 3.433E-02 2.096E-02
5.156E-02 2.449E-01 1.715E-03
1.286E-01 4.405E-02 5.997E-03
Under a comparable data-generation budget, NEST achieves the lowest errors across both geometry classes and both reported quantities. More importantly, the comparison highlights a qualitative difference in how the generated data are used. The global neural operators consume macro-scale solutions associated with specific geometry families, so their training data must cover the target distribution of full domains. NEST instead consumes local physical solves that are independent of any particular macro-scale geometry. The same trained local model is then reused across TPMS and SimJEB structures and across multiple resolutions through Schwarz assembly. Thus, the advantage 9
of NEST is best understood as reusable local training: comparable data-generation effort produces a composable local solver rather than a geometry-specific global surrogate.
6
Conclusion
We introduced NEST, a neural-Schwarz framework that learns local solvers on minimal voxel patches and composes them through overlapping Schwarz iteration to solve PDEs on larger, unseen domains. By separating local physics learning from global assembly, NEST avoids geometry-specific fulldomain training and enables reuse across sizes, shapes, and boundary conditions. On nonlinear 3D elasticity, the same locally trained model generalizes across geometries and surpasses state-of-the-art neural-operator baselines trained on macro-scale data, even under favorable evaluation for those baselines. Once trained, NEST also keeps inference largely independent of PDE nonlinearities. These results support local-to-global composition as a path toward reusable learned PDE solvers. Future work will improve Schwarz convergence and extend NEST to richer physics and boundary conditions.
7
Data and Code Availability
The datasets and source code supporting the findings of this work will be made publicly available upon acceptance of the manuscript.
8
Acknowledgements
This work was funded by the UKRI Engineering and Physical Sciences Research Council (EPSRC), and Tata Steel Research and Innovation Ltd (TSRIL), UK, through an Industrial Cooperative Awards in Science and Engineering (iCase) award, number 220109, to Paolo Secchi. Furthermore, the authors gratefully acknowledge Marco Benedetti for his valuable insights and constructive discussions.
References [1] Tom Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared D Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, et al. Language models are few-shot learners. Advances in neural information processing systems, 33:1877–1901, 2020. [2] Shashank Subramanian, Peter Harrington, Kurt Keutzer, Wahid Bhimji, Dmitriy Morozov, Michael W Mahoney, and Amir Gholami. Towards foundation models for scientific machine learning: Characterizing scaling and transfer behavior. Advances in Neural Information Processing Systems, 36:71242–71262, 2023. [3] Michael McCabe, Bruno Régaldo-Saint Blancard, Liam Parker, Ruben Ohana, Miles Cranmer, Alberto Bietti, Michael Eickenberg, Siavash Golkar, Geraud Krawezik, Francois Lanusse, et al. Multiple physics pretraining for spatiotemporal surrogate models. Advances in Neural Information Processing Systems, 37:119301–119335, 2024. [4] Alicia Tierz, Mikel M Iparraguirre, Icíar Alfaro, David González, Francisco Chinesta, and Elías Cueto. On the feasibility of foundational models for the simulation of physical phenomena. International Journal for Numerical Methods in Engineering, 126(6):e70027, 2025. [5] George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics-informed machine learning. Nature Reviews Physics, 3(6):422–440, 2021. [6] Marco Maurizi, Derek Xu, Yu-Tong Wang, Desheng Yao, David Hahn, Mourad Oudich, Anish Satpati, Mathieu Bauchy, Wei Wang, Yizhou Sun, et al. Designing metamaterials with programmable nonlinear responses and geometric constraints in graph space. Nature Machine Intelligence, 7(7):1023–1036, 2025. [7] Maziar Raissi, Paris Perdikaris, and George E Karniadakis. 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. 10
[8] Siying Ma, Mehrdad M Zadeh, Mauricio Soroco, Wuyang Chen, Jiguo Cao, and Vijay Ganesh. Learning data-efficient and generalizable neural operators via fundamental physics knowledge. arXiv preprint arXiv:2602.15184, 2026. [9] Mikhail Masliaev, Dmitry Gusarov, Ilya Markov, and Alexander Hvatov. Towards universal neural operators through multiphysics pretraining. arXiv preprint arXiv:2511.10829, 2025. [10] Qian-Ze Zhu, Paul Raccuglia, and Michael P Brenner. Generalizing pde emulation with equation-aware neural operators. arXiv preprint arXiv:2511.09729, 2025. [11] Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces with applications to PDEs. Journal of Machine Learning Research, 24(89):1–97, 2023. URL https://jmlr.org/papers/v24/21-1524.html. [12] Andrea Toselli and Olof Widlund. Domain decomposition methods-algorithms and theory, volume 34. Springer Science & Business Media, 2004. [13] Tarek Poonithara Abraham Mathew. Domain decomposition methods for the numerical solution of partial differential equations. Springer, 2008. [14] Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis. Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. Nature Machine Intelligence, 3(3):218–229, 2021. doi: 10.1038/s42256-021-00302-5. URL https://doi.org/10.1038/s42256-021-00302-5. [15] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. arXiv preprint arXiv:2010.08895, 2020. [16] Zongyi Li, Hongkai Zheng, Nikola Kovachki, David Jin, Haoxuan Chen, Burigede Liu, Kamyar Azizzadenesheli, and Anima Anandkumar. Physics-informed neural operator for learning partial differential equations. ACM/IMS Journal of Data Science, 1(3):1–27, 2024. doi: 10.1145/3648506. URL https://doi.org/10.1145/3648506. [17] Tobias Pfaff, Meire Fortunato, Alvaro Sanchez-Gonzalez, and Peter W. Battaglia. Learning mesh-based simulation with graph networks. In International Conference on Learning Representations, 2021. URL https://openreview.net/forum?id=roNqYL0_XP. [18] Zongyi Li, Nikola Kovachki, Christopher Choy, Boyi Li, Jean Kossaifi, Shourya Otta, Mohammad Amin Nabian, Maximilian Stadler, Christian Hundt, Kamyar Azizzadenesheli, and Animashree Anandkumar. Geometry-informed neural operator for largescale 3d PDEs. In Advances in Neural Information Processing Systems, volume 36, 2023. URL https://proceedings.neurips.cc/paper_files/paper/2023/hash/ 70518ea42831f02afc3a2828993935ad-Abstract-Conference.html. [19] Haixu Wu, Huakun Luo, Haowen Wang, Jianmin Wang, and Mingsheng Long. Transolver: A fast transformer solver for PDEs on general geometries. In Proceedings of the 41st International Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, pages 53681–53705. PMLR, 2024. URL https://proceedings.mlr.press/v235/wu24r. html. [20] Huakun Luo, Haixu Wu, Hang Zhou, Lanxiang Xing, Yichen Di, Jianmin Wang, and Mingsheng Long. Transolver++: An accurate neural solver for pdes on million-scale geometries. arXiv preprint arXiv:2502.02414, 2025. [21] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin. Attention is all you need. Advances in neural information processing systems, 30, 2017. [22] Andrey Bryutkin, Jiahao Huang, Zhongying Deng, Guang Yang, Carola-Bibiane Schönlieb, and Angelica I. Aviles-Rivero. HAMLET: Graph transformer neural operator for partial differential equations. In Proceedings of the 41st International Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, pages 4624–4641. PMLR, 2024. URL https://proceedings.mlr.press/v235/bryutkin24a.html. 11
[23] Victorita Dolean, Pierre Jolivet, and Frédéric Nataf. An Introduction to Domain Decomposition Methods: Algorithms, Theory, and Parallel Implementation, volume 144 of Other Titles in Applied Mathematics. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2015. doi: 10.1137/1.9781611974065. URL https://doi.org/10.1137/1.9781611974065. [24] Weihang Ouyang, Yeonjong Shin, Si-Wei Liu, and Lu Lu. Noem: efficient and scalable finite element method enabled by reusable neural operators. Nature Computational Science, 6(4): 417–429, 2026. [25] Ali Taghibakhshi, Nicolas Nytko, Tareq Uz Zaman, Scott MacLachlan, Luke Olson, and Matthew West. Learning interface conditions in domain decomposition solvers. In Advances in Neural Information Processing Systems, volume 35, 2022. URL https://proceedings.neurips.cc/paper_files/paper/2022/hash/ 2f8928efe957139e9c0efc98f173f4be-Abstract-Conference.html. [26] Jianing Huang, Kaixuan Zhang, Youjia Wu, and Ze Cheng. Operator learning with domain decomposition for geometry generalization in pde solving. arXiv preprint arXiv:2504.00510, 2025. [27] Rui Wu, Nikola Kovachki, and Burigede Liu. A learning-based domain decomposition method. Computer Methods in Applied Mechanics and Engineering, 453:118799, 2026. [28] Zongyi Li, Daniel Zhengyu Huang, Burigede Liu, and Anima Anandkumar. Fourier neural operator with learned deformations for pdes on general geometries. Journal of Machine Learning Research, 24(388):1–26, 2023. [29] Xinyi Li, Zongyi Li, Nikola Kovachki, and Anima Anandkumar. Geometric operator learning with optimal transport. arXiv preprint arXiv:2507.20065, 2025. [30] Anders Logg, Kent-Andre Mardal, and Garth Wells. Automated solution of differential equations by the finite element method: The FEniCS book, volume 84. Springer Science & Business Media, 2012. [31] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Graph kernel network for partial differential equations. arXiv preprint arXiv:2003.03485, 2020. [32] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014. [33] Eamon Whalen, Azariah Beyene, and Caitlin Mueller. Simjeb: simulated jet engine bracket dataset. In Computer Graphics Forum, volume 40, pages 9–17. Wiley Online Library, 2021. [34] Oraib Al-Ketan and Rashid K Abu Al-Rub. Mslattice: A free software for generating uniform and graded lattices based on triply periodic minimal surfaces. Material Design & Processing Communications, 3(6):e205, 2021. [35] Maxine Perroni-Scharf, Zachary Ferguson, Thomas Butruille, Carlos Portela, and Mina Konaković Luković. Data-efficient discovery of hyperelastic tpms metamaterials with extreme energy dissipation. In Proceedings of the Special Interest Group on Computer Graphics and Interactive Techniques Conference Conference Papers, pages 1–12, 2025.
12
Appendix A
Error Metric
The following error metrics have been employed: N 2 1 X ∥û (xi ) − u (xi )∥ maxMSE e∥u∥ = N i=1 max ∥u (x )∥2 j j
2 c 11 (xi ) − ∇u11 (xi ) N ∇u X 1 = emaxMSE ∇u11 2 N i=1 maxj ∇u11 (xj )
(14)
where N is the total number of points in the test dataset. Errors are normalized by the maximum value to avoid division by zero, as ∇u is close to zero at certain points within the domain.
B
Role of Schwarz coupling
To isolate the contribution of Schwarz iteration, we compare the full NEST solver against a one-shot local tiling baseline in Table 4. The one-shot baseline applies the trained local model independently on all patches and assembles the predictions with the same partition of unity, but does not iterate to enforce interface consistency. Table 4: Ablation of Schwarz coupling. The one-shot baseline uses the same trained local model and the same patch decomposition as NEST, but removes iterative interface coupling. The relative error on Gyroid TPMS with resolution 60 is evaluated. Method
Resolution
emaxMSE u
emaxMSE ∇u11
60 60
5.711 × 10−1 1.424 × 10−3
8.539 × 10−1 5.593 × 10−3
One-shot local tiling NEST with Schwarz
This ablation is important because it tests whether global accuracy comes from the local neural predictor alone or from the combination of local prediction and iterative Schwarz consistency.
C
Resolution scaling and runtime
We next evaluate the computational cost of NEST inference. Table 5 compares the wall-clock time of the proposed neural-Schwarz prediction pipeline against the FEniCS solver used to generate reference solutions. NEST inference is executed on a single RTX 6000 GPU, while FEM simulations are executed on a single CPU core. This comparison is therefore not intended as a hardware-normalized benchmark against optimized FEM solvers; rather, it quantifies the practical wall-clock speedup of the current NEST implementation relative to the FEM ground-truth generator used in our experiments. Table 5: Wall-clock runtime comparison on a Gyroid TPMS geometry. NEST is run on a single RTX 6000 GPU, while FEM is run on a single CPU core. Resolution
Iterations
FEM (s)
NEST total (s)
NEST / iter. (s)
Speedup
15 30 45 60
4 14 29 61
29.90 166.70 459.60 811.15
0.27 9.43 57.81 257.55
0.07 0.67 1.99 4.22
110.7× 17.7× 8.0× 3.2×
NEST provides substantial wall-clock speedups at all tested resolutions. The speedup decreases with increasing resolution because the number of Schwarz iterations grows from 4 iterations at resolution 153 to 61 iterations at resolution 603 . This behavior is consistent with one-level Schwarz methods, where global information propagates through local overlap communication. Thus, the current implementation should not be interpreted as demonstrating total O(N ) runtime. Rather, the per-iteration cost is expected to scale approximately linearly with the number of active patches/nodes, while the total runtime also depends on the number of Schwarz iterations required for global convergence. 13
To directly test this scaling claim, we decompose the total NEST runtime as TNEST ≈ Nit Tsweep ,
Tsweep = O(P α ),
(15)
where P is the number of active patches, Nit is the number of Schwarz iterations, and α is the per-sweep scaling exponent. The expected behavior Tsweep = O(P ) corresponds to α = 1. We estimate α empirically between consecutive resolutions as log Ti+1 /Ti , αi→i+1 = (16) log Pi+1 /Pi where Ti and Pi denote the per-iteration runtime and number of active patches at resolution i. Table 6 reports the number of active nodes, number of active patches, per-iteration runtime, and the resulting empirical exponent. Table 6: Empirical scaling of one NEST Schwarz sweep on a Gyroid TPMS geometry across resolutions. The empirical exponent α is computed pairwise between consecutive resolutions; values close to 1 indicate the predicted O(P ) per-sweep cost. Resolution 15 30 45 60
Active nodes 3
2.42 × 10 8.60 × 103 2.89 × 104 6.94 × 104
Active patches 3
1.11 × 10 1.37 × 104 4.02 × 104 8.93 × 104
NEST / iter. (s)
Empirical exponent α
0.07 0.67 1.99 4.22
– 0.90 1.01 0.94
The measured exponents lie in the narrow range 0.90–1.01, empirically supporting the predicted O(P ) per-sweep cost. Combined with Table 5, this confirms that the resolution-dependent decay in speedup is driven by the growth in Schwarz iterations Nit rather than by super-linear per-sweep work. These results identify global Schwarz convergence, rather than local neural inference, as the main scalability bottleneck of the current implementation. This motivates multilevel coarse corrections, learned relaxation strategies, or accelerated Schwarz variants as natural directions for further improving large-scale performance. A further expected advantage of NEST is that, once the local neural operator has been trained, inference does not require evaluating or linearizing the constitutive law. Therefore, the cost of a NEST patch evaluation is largely independent of whether the underlying PDE contains a simple or complex constitutive model. This contrasts with nonlinear FEM solvers, whose cost can increase substantially with more complex material laws due to quadrature, stress evaluation, tangent assembly, and nonlinear iterations.
D
NEST inference: pseudo-code and algorithmic details
This appendix complements Section 4.2 with a self-contained algorithmic description of the NEST inference pipeline. Algorithm 1, Phase 1 performs offline pre-processing of an arbitrary voxelized domain Ω, and Algorithm 1, Phase 2 performs the Schwarz iteration itself. We adopt the notation of Section 4: the macro voxel grid has cells {Ck } of edge length h with Nd cells along axis d ∈ {1, 2, 3}, active node set V ⊆ R3 , macro Dirichlet boundary ∂Ω ⊂ V, overlapping 3 × 3 × 3 patches {Ωp }P p=1 int ext with active nodes Vp , edges Ep , and boundary partition ∂Ωp = Γext ∪ Γ , Γ = ∂Ω ∩ ∂Ω. p p p p We write {χp }P for the partition of unity over V, S for the frozen local GNO of Section 4.1, θ p=1 g : ∂Ω → R3 for the macro Dirichlet datum, and ge for its inverse-distance−2 extension to interior nodes defined in Section 4.1. The per-axis partition-of-unity weight is the piecewise-linear ramp 0 ≤ ξ < 13 and s > 0, 3ξ, τ (ξ; s, N ) = 3(1 − ξ), 23 < ξ ≤ 1 and s < N − 3, (17) 1, otherwise, where ξ ∈ [0, 1] is the patch-reference coordinate along one axis, s the patch anchor, and N the number of macro cells on that axis. 14
Algorithmic clarifications. Partition-of-unity weights. χp is the tensor product of three per-axis ramps (17): linear from 0 to 1 on the first cell unless flush-left (sd = 0), constant 1 on the central cell, and linear from 1 to 0 on the last cell unless flush-right (sd = Nd − 3). Weights are floored at P 10−6 and renormalised at assembly so that p χp (v) ≡ 1 on every covered v ∈ V. Update style and damping. The sweep is additive and Jacobi-style: in iteration n + 1 all patches read the same u(n) , with no Gauss–Seidel sequential pass. The damping θ ∈ (0, 1] relaxes the swept field after partition-of-unity assembly; θ = 1 recovers undamped additive Schwarz and is the default in all reported experiments. Algorithm 1 NEST Framework: Pre-processing and Schwarz Inference Phase 1: Pre-processing of a voxelized domain Ω Require: macro voxel mask S ∈ {0, 1}N1 ×N2 ×N3 (largest face-connected component of the voxelizer output) ext int Ensure: patches {Ωp }P p=1 with Vp , Ep , Γp , Γp , χp 1: V ← {v : ∃ Ck ∈ adj(v) with Sk = 1} 2: ∂Ω ← {v ∈ V : ∃ Ck ∈ adj(v) with Sk = 0 or Ck ∈ / S} 3: for (s1 , s2 , s3 ) with 0 ≤ sd ≤ Nd − 3 for d ∈ {1, 2, 3} do 4: Ωp ← cells of S in the 33 block anchored at (s1 , s2 , s3 ) 5: Vp , ∂Ωp ← same 8-cell stencil restricted to Ωp 6: if Vp = ∅ or ∂Ωp ∩ Vp = ∅ then 7: continue 8: end if int ext 9: Γext p ← ∂Ωp ∩ ∂Ω, Γp ← ∂Ωp \ Γp
Ep ← {(i, j) ∈ Vp2 : 0 < ∥xi − xj ∥ ≤ 1.5 h} 11: χp (v) ← τ (ξ1 ; s1 , N1 ) τ (ξ2 ; s2 , N2 ) τ (ξ3 ; s3 , N3 ) for each v ∈ Vp with patch-reference (ξ1 , ξ2 , ξ3 ) ∈ [0, 1]3 12: end for 13: return {Ωp }P p=1 , V, ∂Ω
10:
Phase 2: NEST Schwarz inference (additive, Jacobi-style) Require: patches {Ωp }P p=1 from Phase 1; frozen Sθ ; macro datum g; tolerance ε; damping θ ∈ (0, 1]; max iterations Nmax Ensure: global displacement u⋆ : V → R3 1: u(0) ← g e on V 2: for n = 0, 1, . . . , Nmax − 1 do 3: a ← 0 ∈ R|V|×3 , w ← 0 ∈ R|V| 4: for p = 1, . . . , P do (n)
▷ u(0) ∂Ω = g
(n)
← u(n) ▷ Eq. (11) g (n+1) (n) 6: ûp ← Sθ Vp , Ep , gp (n+1) 7: a[Vp ] += χp ⊙ ûp , w[Vp ] += χp 8: end for 1 9: u(n+ 2 ) (v) ← a(v)/w(v) for v ∈ V with w(v) > 0 ▷ Eq. (12) 1 10: u(n+1) ← θ u(n+ 2 ) + (1 − θ) u(n) , u(n+1) ∂Ω ← g 11: Compute ∆(n+1) as in Eq. (13); if ∆(n+1) < ε then break 12: end for 13: u⋆ ← u(n+1) 14: Optional gradient pass: for each p, feed u⋆ ∂Ω to the dedicated ∇uij surrogate and assemble p per component through {χp }. 15: return u⋆ 5:
gp
Γext p
← g, gp
Γint p
15