Boundary-Value PDEs Meet Higher-Order Differential Topology-aware GNNs
Authors: Yunfeng Liao
Yangxin Wu
Xiucheng Li
Venue: NeurIPS 2025 spotlight ()
Keywords: Partial differential equations
Neural operator
Higher-order graph neural network
Exterior differential calculus
Electromagneticism
Primary area: machine_learning_for_sciences

Abstract: Recent advances in graph neural network (GNN)-based neural operators have demonstrated significant progress in solving partial differential equations (PDEs) by effectively representing computational meshes. However, most existing approaches overlook the intrinsic physical and topological meaning of higher-order elements in the mesh, which are closely tied to differential forms. In this paper, we propose a higher-order GNN framework that incorporates higher-order interactions based on discrete and finite element exterior calculus. The time-independent boundary value problems (BVPs) in electromagnetism are instantiated to illustrate the proposed framework. It can be easily generalized to other PDEs that admit differential form formulations. Moreover, the novel physics-informed loss terms, integrated form estimators, and theoretical support are derived correspondingly. Experiments show that our proposed method outperforms the existing neural operators by large margins on BVPs in electromagnetism. Our code is available at https://github.com/Supradax/Higher-Order-Differential-Topology-aware-GNN.

Boundary-Valle PDEs Meet Higher-Order Differential Topology-aware GNNs
Yunfeng Liao Yangxin Wu Xiucheng Li
School of Computer Science and Technology, Harbin Institute of Technology (Shenzhen) lixiuchengGhit, edu.
Abstract
Recent advances in graph neural network (GNN)-based neural operators have demonstrated significant progress in solving partial differential equations (PDEs) by effectively representing computational meshes. However, most existing approaches overlook the intrinsic physical and topological meaning of higher-order elements in the mesh, which are closely tied differential forms. this paper, we propose a higher-order GNN framework that incorporates higher-order interactions based on discrete and finite element exterior calculus. The time-inne iemeiindeeendent boundary value problems (BVPs) in electromagnetism framework. It can be easily generalized to other PDEs that admit differential form formulations. Moreover, the novel physics-informed form estimators, and theoretical support are derived correspondingly. Experiments show that proposed method outperforms the existing neural operators by large margins on B VPs in electromagnetism. Our code is available at attpss//gi thub. com/ Supp/ Hi gher- Seteer-Diffeeentii Topollyy aware-GNN.
instantiated to illustrate the proposed
informed loss terms, integrated
1 Introduction
Solving partial differential equations (PDEs) accurately is fundamental in scientific computations. Traditional numerical solvers [1] and the emerging physics-informed neural networks [2] rely iterative computations. This is a significant bottleneck that precludes their application time- sensitive domains where slight inaccuracy is tolerable but speed matters, such gaming engines and interactive simulations. To address these limitations, neural operators [3-6] propose to directly learn the mapping between initial/boundary conditions and complete PDE solutions, climinating tine-consuming iterations while maintaining the capability to handle the PDEs whose exact analytical solutions are unattainable.
Convolutional neural network-based solvers are inherently constrained to regular, grid-like domains, whereas numerical solvers typically employ meshes to represent irregular solving regions approach naturally aligned with graph neural networks (GNNs). This compatibility has spurred growing interest in GNNs for time-dependent physics simulations, demonstrated through applications from fabric dynamics in wind [7] and granular particle systems [8] to neural mesh refinement schemes [9]. However, solving lime-independent PDES, especially boundary value problems (BVPs), presents greater challenges due to the absence of temporal guidance (initial data) and limited feature representation. Recent work [10] attempts apply GNNs BVPs in electromagnetism and has shown promising results, but is still limited to relatively simple cases.
Conventional GNN-based BVP solvers primarily utilize vertex adjacency in meshes while neglecting higher-order topological elements (edges, faces, cells), despite their fundamental physics interpre- tations from the perspective of differential forms [11]. While vector analysis has long dominated
Corresponding Author.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

3
D
Node
Eáge
Pece
Cell
Figure 1: In the proposed DEC-HOGNN, scalar and vector fields nodes, edges, faces and cells are encoded into k-simplex features, passing through HOGNN and decoded back as target fields.
naturally interpreted as differential forms on manifolds [12]. In electromagnetism, it reveals a hierarchical structure: potentials (f) manifest O-forms, field intensities (E, H) 1-forms, flux densities (D, B) as 2-forms and density distributions (p) as 3-forms [13]. This formalism, crystallized in the Maxwell's 's House representation of Maxwell's equations [14], treats traditional vector fields as mere proxies for underlying forms. Discrete exterior calculus [15] operationalizes this approach through De Rham mappings that represent k-forms as integrals over k-simplices [16]. Also, introducing geometric objects like differential forms allows us to naturally generalize this framework to PDEs on curved spaces [17] beyond merely incorporating topological structures [18, 19].
physical modeling, modern physics recognizes that many vector fields
In this paper, we propose higher-order GNN-based PDE solver framework by exploring the ideas in discrete exterior calculus (DEC) and finite element exterior calculus (FEEC) in principled manner. By encoding the integrals over k-simplices as k-simplex features within higher-order GNNs (HOGNNs), which explicitly model interactions between simplices of varying dimensions, we establish a principled framework for solving form-based PDEs while preserving the topological and physical structure inherent the problem domain. The resulting framework is dubbed DEC-HOGNN and illustrated in Figure 1. DEC-HOGNN enjoys better physical interpretation from the differential form perspective and can be extended higher-dimensie cases naturally. Our main contributions are summarized as follows: 1) We design differential topology and decodes PDE operators based on the principles of DEC and FEEC. 2) Various physies-informed loss terms are derived under DEC-HOGNN, including the boundary condition ones, which can enable solving the boundary-value PDEs more effectively. 3) The universal approximation property in solving Poisson problems (electrostatics and magnetostatics) is presented, and the performance excellence is demonstrated via empirical experiments.
HOGNN, which naturally encodes
Related Work
Higher-Order GNNs. HOGNN is an extended learning framework on generalized graphs, i.c.. hypergraphs [20]. A hypergraph G' allows hyperedge to contain more than two vertices [21]. HOGNN leverages the more abundant adjacencies on hypergraphs and mimics what GNN does on plain graphs via redefining various neighborhoods. If G has no further decorated structures, then one can define the k-node-tuple adjacency [22]; while boundary, co-boundary, upper, and lower adjacencies are available when G is simplical complex [18] cell complex [23]. Mechanisms in GNN are mostly based on special adjacency induced by edges and thus can be easily transplanted to hypergraphs. The graph-convolution [24], attention [25]. and generalized message passing [26] of HOGNN all fall into this category. HOGNN has been rell-studied in various regions, such as recommendation systems [27] and molecular predictions [28], where multi-body interactions are of significance, but beyond the expressive ability of plain graphs.
Neural Operators. Neural operator [3] learns function-tor coriion-tt-funnttit functiono mappings mainly using data- driven loss instead of PDE-based physies loss. Its original implementation is furnished with kernel convolution of (n22) complexity. which can be improved to O(n log n) in the spectral domain via discrete fast Fourier transform [3]. is then extended non-square-like regions [4] and spatial- spectral neural operators realized by wavelet transform [29-31]. This field later gradually shifts towards Transformer architectures with PDE-compatible attentions, from GNOT [5] employing boundary- cross-attention to Transolver [6] using attentions among domain slices.

---

n (M)
f (M)
(M) d (P(
d
d
0 (+(M) a (+(M)
B
n (+M) (P(
Figure 2: The primal manifold .M has 4 tetrahedrals. The i-cell or E M {0, 1, 3)) and its dual cell *0; are illustrated, associated with specific magnetic quantities. Identify the scalar potential as 0-form, vector potential A and magnetic intensity H as I-form and magnetic flux density B as 2-form, and the relations H V. B V A, V H 0, V B are encoded into exterior derivative d and its property d² 0. The electric case is similar.
GNN-based PDE Solvers. GNN-based neural operators have been studied in time-evolving mesh- based and particle-based physics simulations [7, 8]. The partieles appear in the form of point clouds, and adjacency is built on local neighborhoods, in which equivariance is introduced to enhance model performance, such as subequivariance [32] and IsoGCN [33], an equivariant data-driven neural differential operator. To align with Neumann boundary in PDEs, NIsoGCN [34] further introduces Neumann term into the differential kernel. This category requires data on fields and their differentials, but the latter is often intractable. It is found that message passing in GNNs can represent various numerical methods for time-dependent PDEs [35] and aligns with finite volume methods to achieve local mass conservation [36]. In addition, GNN-based approaches for BVPs [10] and inverse problems [37] are also explored.
3 Preliminary
De-Rham Mapping. In algebraic topology, a simplex chain complex C(X) on a topological space X is graded vector space of k-orderr simplices (x). decorated with the boundary operator a A k-cochain o² (X)) Ca (X) + I is real-valued function of σ E Cx(X). If there exists a diffeomorphism tp between and smooth manifold M, then for any k-form w, we obtain covariant functor mapping from Ca (X) to C* (X). namely, Re-Rham mapping [16]:
If w.
F F:C (*))
pplo.)
DEC and FEEC. Discrete Exterior Calculus (DEC) offers a comprehensive and differential topology- preserving toolkit to discretize operators on manifolds. In contrast to Graph Calculus [38], which treats a discretized manifold as a graph-aa at the cost of losing essential differential properties and thereby introducing inaccuracies -DCiaaaciis-eCCCCCCCC maintains differential topology properties by working with integrations. Specifically, it characterizes a k-form w E 12 (M) on a discrete manifold M through its integral over every k-simplex [16]. It yields high accuracy in applications where differential information matters, and is widely used in computational physics [39-42]. Another benefit of DEC is that it allows for processing manifolds from dual perspective. For instance, the Hodge star * (M) - -2"="(" (MM on M gives the important constitution relation between B and H in clectromagnetism. Even if whw is intractable the discrete case, we can still estimate its integral on the dual manifold +M (shown in Figure 2) without a priori the metric [15]. Finite Element Exterior Calculus (FEEC) generalizes the Finite Element Method (FEM) and a numerical implementation of Galerkin methods. FEM is in essence an interpolation method with nodal functions, a.k.a., Lagrangian element [1]. However, such node-wise interpolation is inherently incompatible with differential operators such as div and curl, whereas the adoption of higher-order finite elements in FEEC enables the exact representation these operators in the integral sense via straightforward linear combinations [14].
Whitney Element. Whitney elements are a type of finite elements used in FEEC and DEC. They provide way approximate differential forms discretized manifold that respects the geometry

---

dE
Convelution
Prinal Bdge Cobousdery Mijoomcy
Groph Convolation
Field Recovery Recovery-Deeading
D
Daal Edge Caboundary Mijomny
E
Field Enconding
Comvolution
Edgn Adpoeeyy
Feature Forwarding
Figure 3: The illustration of DEC-HOGNN for 2D-electrostaties BVPs.
and topology of the manifold. Let w and w denote scalar and vector fields, respectively. Given tetrahedralized mesh representation of manifold M with dim M 3, let [i, denote the edge containing Uj and so are [ k), [i, k, for face and cell, respectively. The canonical node elements tw; (v,(x) d.g} are functions M. in which 6,z is the Kronecker delta. Let (MM), W (²(M) and (MM) denote the function spaces spanned by node elements wi, edge elements W ((.d)) face elements Wi..... and cell elements ([..kk.) (scalar fields), respectively. Then the integral of W[sjj] along [E, .]. the flux of A]..j.A] through face [f, kk], and the volume integral of (k..kk.] within tetrahedral [T, J, k, 4] all equal to [14]. The element definitions are as follows:
W[oa = wi Vwy - W, Vw; W K. 2(w, Vw; X Vwg) + wy (Vw x Vw;) (W, Vw,))
(3)
[i Σ wi Vw, Vwk) Vwi Xxx(F.iA.1 vol(
(4)
in which wi, Wj, WR, wi E w (M). For ease of reference, the notation frequently used throughout the paper is summarized in Table Appendix A.
4 Methodology
Problem Setup and Motivation. In this work, we focus on PDEs that admit differential form formulation (many important PDEs in physies and engineering can be formulated using differential forms, such as Maxwell's, Navier-Stokes, Yang-Mills Equations, etc.). Given PDEs (often formulated in vector fields) on a discrete manifold M, we aim learn neural operator Ge that takes as inputs the observed scalar fields {s;} and vector fields {v;), and outputs the target vector field of interest:
u = Ge((s, 3m (on.:.(v... (Vi...., M).
(5)
In comparison to the vector field formulation, the differential-fome characterization of PDEs is coordinate-inde indeeende and explicitly reveals the geometric and topological aspects of the spaces. Motivated by this, we propose to transform the vector fields into the differential form formulation of PDEs, which enables us to develop the neural operators by combining the principles DEC and FEEC. The benefits are threefold: 1) they permit exploring the higher-order topological elements to facilitate PDE solving, 2) developing various chysics-consistent losses, especially for the boundary- value problems, and 3) differential operators can be interpreted as simple linear combinations in DEC and FEEC, aligning with the higher-order message passing framework. Once the PDEs are solved, we then translate the results into the vector field language for downstream analysis.
Method Overview. Since differential f-forms can be represented integrals on k-dimensional elements and further identified as higher-order element features, we introduce physics-informed higher-order interactions and aggregations into existing HOGNN frameworks in Section 4. Nev- ertheless, the initial input and expected output are often vectors instead of forms in practice. To circumvent this, we propose proper encoder-decoder to enable consistent transformations in Sec- tion 4.2, In nutshell, the input vectors encoded into forms (higher-order features), processed by

---

HOGNN, and then decoded back into vectors, as illustrated in Figure 1. The physics- loss and the universal approximation property of the proposed method are presented in Section 4.3.
As an illustrative example, we instantiate our framework by solving the classical BVPs in electro- magnetism, and the proposed method can be easily generalized other PDEs that can be formulated using differential forms. Recall that a typical Neumann electrostatic BVP is
V D p. E O in S2; E = Eo, D Do on 02;; D = E; E 92,
(6)
in which p. D, E are the charge density, displacement field, and electric field, respectively, whereas {22, is partition of the domain Ω. The permittivity tensor E; can vary in different 22 consisting of different materials. In this case, the learned operator Ge will take as input (fp). {Eo. Do} 92) and produce the complete fields (E, D). To this end, we rewrite Eq. 6 in its differential form formulation.
(7) in which d- denotes the exterior derivative and D, E are the corresponding k-forms. The input scalar field p and masked vector fields E, D are first encoded into integrated forms in the form of vertex, edge and face features. Three types of edge adjacencies are used in the higher-order GNN, which give important physics interpretations (as shown in Figure 3): the curl-free property dE is implicitly included in the primal edge co-boundary adjacency, while Gauss's s Law dD is involved in the dual edge co-boundary adjacency. These GNN layers can be stacked sequentially, and the aggregated features can be either decoded back as complete E, D or forwarded to other networks.
dD = dE = 0 in S2 E Eo D = Do on 0Ω; == E in S4g
4.1 Physics-Informer Higher-Order Interactions
In light of DEC, we identify the potential on node, the circulation along edge, the flux through face, and the mass within cell as the node, edge, face, and cell features, respectively. Given a discrete manifold M and its dual xM, since both of them are cell complexes, an element e onward have four types of neighborhood: boundary B(c)= == de and co-boundary C(c) == {d E B(d)), upper adjacency (+(c) := {d: 36, {c, d) C B and lower adjacency Ni(e) := fd 36, fc, d) c C())) Following the paradigm of [26], the element feature h' of e is updated by its original feature h. and the aggregated messages from four neighborhoods B(c), ((C)..
via:
h2 = he, m m. m No "."" m²i), CEM,
(8)
m
B A Aggregate (B( ha),
(9)
m
E
c
= Aggregatc (he, ha),
(10)
N₁
m
SggegatteueNNoo).eeeeee PPN (he, ha, hs),
(11)
N1
m
AggregatejeN,,e) dEBIC)) PAN (be, has hg).
(12)
Since forms in a PDE can be defined on both M and +M (e.g., E and D in Eq. 7), we also propose a higher-order MPNN on the dual manifold xM for dual forms by the fact:
B(*c) = (xd : d E ((c)), C(*c) = (*d : d E B(c)), (+c) = [xd d (c)),, Ni(xc) = (td: d e (G(c)).
Usually, not all neighborhoods and elements will be used in higher-order MPNN due to computational complexity and the absence of features on corresponding elements. For instance, 2D-electrostatics BVPs only involve E, D, p and thus only adjacencies about edges and faces are considered. In recent studies, there are various ways to implement the message passing on four different adjacencies, including generalized convolution, attention, and Transformer. We list several candidates of the adjacency layer backbones in Table and will study their impact on model efficacy in Section 5.2.
In Table 1, [,(k) D A 2D ITT is the higher-order Laplacian [21], defined by indicator matrix T,y == Xo,EN((a,) and diagonal degree matrix D; İ. is the Laplacian of an extended graph G = (C+ U Ck+1, [(c, d) : Ct. d E C(c))); he.d, is the concatenation of h ha, h& and he.d is likewise. In our implementation. the relative orientation sign(c, d) between two elements is further considered, via replacing , by sign (c, ddh..d-
Consistency with Conservation Law. It is beneficial for solving BVPs by introducing higher-order interactions, because the interactions have meaningful physics interpretations, e.g., the sum of edge
5

---

Table 1: Possible implementations for co-boundary/lower adjacency layers.
Layer
Co-boundary and lower adjacencies
h = (WN())) 1
Σ
( (A))))))))) + (ruttpattigghit
(6,8)EM(())
Convolution
h = =((o(c)]=
Σ
((10)))))))))
dEC(c)
h/ = o(Woho
(eaa W. a (d,8)eM(())
Attention
h/ = o(Wohe
softmaxet((o)Wxxh.ooo Wahhadd
dec(e)
h/g o(Wche
( ka)
(d.5)EMi(e)
Transformer
hi = =(Wohe
deC(e)
features around a face corresponds to the vorticity, while the sum of face features around a cell indicates the divergence that location. This enables us preserve structural conservation discretely by using integrated differential forms. Moreover, prior work also shows that adding conservation regularization can boost the performance of PDE solvers. For instance, a conservative GNN solver for 2D fluid dynamics [36] employs message passing based on face lower adjacency in the HOGNN framework and achieves local mass conservation (divergence-
through asymmetric aggregation. We can realize more conservation constraints with DEC, e.g., enforcing vorticity conservation (V x = 0) and divergence conservation (V B 0, V D = p) electromagnetism.
Consistency with Electromagnetic Constitution Law. As a side-product, we can estimate the electric and magnetic permeability E and μ of the medium while solving BVPs. The constitution law B = μH, D = zE can be equivalently interpreted by Hodge star: B = #' H E. It allows us estimate the μ. E along different directions via the definition DEC:
vol (xOx) vol (ox) S. w, WE)+(M).
(13)
4.2 Encoder-Decoder between Vector Fields and Forms
The encoder aims to transform vector fields into higher-order element features (integrated forms). In 3D-cases, 1-form Wi and 2-form w2 be realized by vector proxies u, and u22 with the aid of a unit tangent vector t on an edge and a unit normal vector n on a face a2, respectively; while 3-form W3 can be identified as scalar field ug : R R. and the integrated form is the usual volume integral. More formally.
Encoder.
[ [W1==
[1: : / u uddss
[ W2 = / U2 nds,
/ w3
wa u3dV
(14)
Jaz
Jr3
Ja3
Given the samples {x,) observed on o, (i (i=1) (1 2, 3), the integrals on the right-hand sides of Eq. 14 can be estimated yield the integrated forms numerically with either the Monte Carlo method or canonical quadrature rules on simplices (c.g., cubic quadrature on a triangle [43]). Similarly, the 2D case is given as follows and can be estimated similarly:
[.
w1
tds,
W2
u2dady.
(15)
In DEC, the differential forms in the primal manifold and dual manifold are related by the Hodge star x (e.g., E E Ω (M) and D E S22 (+M) Eq. 7). Since we only have observations primal nodes, it is simpler to estimate the integrated 2-forms on the primal manifold. However, DEC requires the integrals on the dual manifold, which poses a challenge for implementation. To sidestep this, we derive approximation of the integrated dual forms based on the integrated primal forms, presented Theorem 1. The details of the theorem can be found in Appendix Intuitively, the theorem offers us simpler way to calculate the integrated dual forms.

---

Theorem 1 (Non-Cycle Forms Estimation). For smooth k-form w that is not a cycle defined on bounded region, i.e. dw does not vanish identically in any k-simplex dge one can estimate the integrated dual forms [Jo... w On-k k. E +Cn+ (X)) by using the integrated primal forms (.a,, w: oh E [x)) up to accuracy ((ghtl) in which & == SUPoxECx(X) diam ok-
Decoder. The decoder intends to recover the vector fields from the derived integrated forms. The principle is that the integral of proper vector field proxy on a simplex should equal the integrated form onward. To find the proxy at v; for a I-form (e.g.. E' in Eq. 7) along edge [V;, Uy], note that Mes.sy]]the the final edge feature output by DEC-HOGNN indicates the circulation of the proxy and the circulation of Wpr,ey) in Eq. 2 along edge (U, 1y] is 1, and thus the Whitney electric field contributed by edge [o,, my] at v, is set to hi W Wp,.e,l- The contributions to v, from different edges can be aggregated by average pooling or attention-style weighted sum. Similarly, for the proxies of 2-forms (e.g.. D in Eq. 7) with final face feature [jp..,y.ou] and 3-forms (e.g., p in Eq. 7) with final cell feature Miu..,y.....n) the face and cell contribution are h [p W [p...y.mu) and Mps. (P respectively, with Wi [x) "[ver"s.Me.rr defined by Eq. 3 and Eq. 4.
Remark. The interpolation nature of the decoding approach offers a continuous field not only defined the nodes but also in the cells. Though other GNN operators can also estimate the inner fields by multi-linear interpolation based on values at nodes, such methods can destroy certain important physical properties compared with Whitney elements. For instance, the reflection law in clectrodynamics clucidates that the tangent component of E on the medium interface is continuous while the normal component is not. With prior knowledge of the interface positions, one can choose to average out contributions from each side respectively and obtain two different E's before and after reflection, which is consistent with physics on the interface, while usual interpolation based on node-wiee wise data cannot. Lastly, we present the theoretical analysis in Theorem 2, demonstrating that our encoding and decoding scheme preserves sufficient information for scalar and vector fields on sufficiently fine mesh
Theorem 2 (Proper Sncooer-Decodder Given a uniform partition {[z,, F24 111 on the
unit interval, i.e., 0= I1 T2... < T2 +1 1, the encoding operator En WW1.2 1)) + g²a is
defined as the 2 integrals on each [1,, 14+1] and the decoding operator D : R2² + C (00 1))
maps any 2 -dimensional feature to a function fxx) on [0, 1) such that | f (x) = h, Then
for any [20 1]) and E > 0, there exists an integer M > O such that for all N > M,
1l - Dx Ex(1))le> <e
(16)
By definition, we have EN DN id.
4.3 Informee Loss and Universal Approximation Property
Physics-Informn Loss. Purely data-driven methods are likely to go against physics. Hence, PINN introduces physics-constrainee loss regularization to guide models to learn beyond the data resolution [44]. In PINN, the PDE residual is usually estimated by point-wise sampling and cittoddifferentatiinn. This is incompatible with GNNs since vertices are discrete, making it hard to enforce physics-consistency constraints, e.g., constraining a field to be curl-free [10]. But in integrated forms, one can interpret differential version of PDEs into the integral version. For instance, magnetic flux density B is always divergence-free, inducing a DEC-version constraint:
/
L B (17)
dB + B 0
V
B-n Σ 02E0 43
B Dog B-
Σ
03E 72
In our proposed framework, J., B is indeed feature of face 02. Therefore, we can introduce physics- informed loss without sampling and differentiation, but by simply summing up corresponding higher-order features. More physics-informer loss terms are derived in Appendix B. covering the divergence, vorticity, and boundary conditions.
Universal Approxímation Property. Theorem shows the universal approximation ability of DEC-HOGNN in solving the Poisson problems. The proof is available in Appendix D.
Theorem 3 (Universal Approximation Property). Let (2) denote the Hilbert space on a bounded closed region Ω with C² -b-buddarry and P, D, Eg be the encoder- decoder projection with
7

---

Maroed actre R
Paiteed Dww.ccamenctoooocoo.. D30
Camplea Seurca JoX
Camplita Anerttriee Fead Ese
Conciraa Dediacament Dla
Figure 4: In the electrostatic field task, the model is required to output complete fields shown in the right-side two panels based on the left-side three panels.
Figure 5: Models are compared on various settings. From left to right: 2D electrostatics BVP on square 2D-mesh with two kinds of medium; 2D magnetostatics irregular holed mesh; 3D) electrostatics and magnetostatics on a 3D socket model and quarter-gear model, respectively.
resolution Given a simplex-partition (Xa) of with measure-zero intersections, for any f V VUU V.V.V= == ([o ([" (f +I flleo < M), there exists network in the form of D L.... L1 o with infinite neurons to solve the Dirichler Poission problem Q such that:
lim sup diam
11D LN... Li E(oE( - u|12> = 0,
(18)
in which Au(x) = f(x). Q2 (xx = g x an
As corollary, note that E Vu and the experiment setting below is equivalent to solving Neumann boundary condition in a linear medium, the universal approximation property also rings true by similar proof in Appendix D via single-layer potential method.
5
Experiment
5.1 Experiments and Benchmarks
Experiment Tasks. We assess the performance of the proposed model on 2D electric and magnetic boundary value problems. In the electrostatic case, the system is governed by potential yo:
( p(x), x E Ω2 (x(x) u(x), x E an,
which gives the electric intensity E and electric displacement D. In our settings, the domain Ω is partitioned into two regions 12 22 with measure-zen intersection. Q2 has isotropic permeability E while &t has non-isotropic linear one E2, shown in the blue and green-colored regions in Figure 4, and
D(x) = =xE(), x Ω, 1,2.
Let H2) { 11/112 + ||V <x) and H²( the vector-valued function space on Ω with ( components. The operator G: (22 x (an)) + H²((2)), E- Xxeen, D Xxeσn) + (E, D recovers the field E, D on Ω based on prior knowledge the boundary field data (E, D) in and the charge density distribution p on S2 as shown Figure 4.
the function space on Ω with L derivatives, i.e., H²
Dataset Generation. Both 2D and 3D meshes are adopted illustrated in Figure 5. The electrostatic data obtained by FEM-based electromagnetism PDE solvers. Ω is partitioned into many triangles X; with measure-7ero intersections. The charge density in each X, is randomly assigned following a uniform distribution. The PDEs are then solved in large enough vacuum region with Neumann boundary conditions. E, D boundary vertices and p all vertices are sampled as input while

---

Table 2: Performance of different neural operators in terms of MSE. The bold and underscored numbers indicate the best and recond-best in each column, respectively.
Model
2D Electrostatics
Magnetostatics
3D Electrostatics
3D Magnetostaties
DeepONet
1.866=0.025
1,438:0.021
0.842+0.09
0.483+0.000
MKGN
1.960=0.097
2.324:0.200
0.793:0.043
1.382+0.033
Galerkin- Typp
1.209=0,022 0,2
1.895+0. 144
1.176+0.055 030
0.120+0.003
GNOT
2.064==0.052
1.142+0.149
7.266+0.014
4.480=0.543
Transolver
2.751=0.009
7.43240.143
8.249+0.00 0099
7.089=0.460
GCN-based
1.36220.015
1.224:0.026
1.98220.89
0.18220,000
GAT-based
1.360=0.014
1.75540.040
1.02940.00 .49
1.22240.00
Graph UNet-based
1.917=0.024
1,283+0.00 004
0,809
0.375=0.000
GTTbbaeed
0,990=0.031
1,405+0.077
0,244+0.025
0.171=0.025
DEC-HOGNN (Ours)
0.623+0.058
0.87540.052
0.195+0.005
0.158+0.008
Table 3: Test loss and performance drop of different variants against the entire models.
Interactions
Lower adjacency
PD
D
None
PD
PD
PD None
None
Co-boundary
PD
PD
PD
PD
D
P
None
Test loss (mean square error)
0.652
0,676
0.678
0.719
0,954
0.959
1.244
2.077
Performance drop (%)
0.00
3.81
4.12
10.34
46.35
47.23
90.98
218.87
E, on all vertices are the target output. Scalar fields are normalized, and vector fields are shrunk with respect to the average norm to eliminate magnitude differences.
Note that the 2D magnetostatic field is also governed by potential A, (the vector potential A has only one non-vanishing component A₂ the 2D space identified as the zyy-llnne). Therefore, the data generation is similar. For more details, please refer to Appendix E.
5.2 Numerical Results
Baselines and Implementation Details. To showcase the necessity of introducing specifically devised solvers for BVPs and the efficacy the proposed method, we evaluate our model against the following general time-dependent PDE neural solvers, in which the input fields are masked accord-
ingly: DeepONet [45], MKGN [46], Galerkin-type Attention [47], GNOT [5], and Transolver [6]; in addition, various GNN-based solvers devised particularly for BVPs are also included (as in [10]): GCN [48], GAT [49], Graph U-Net [50], and Graph-Trranformert
BVP solvers [51].
Main Results. Table presents the MSE (mean square error) of different methods on 2D/3D clecrostaicc//mgnettsaaasss BVPs. The neural operators from DeepONet to Transolver, which are mostly designed for time-dependent PDEs, yield larger errors compared with GNN-based BVP solvers. because many time-evolving PDE solvers like DeepONet predict an increment based on previous observations, which is unfortunately intractable in static cases. The multi- design in MKGN and the heterogeneous cross-atte attention GNOT both suffer from the lack of features, as time-independent PDEs can be determined by the boundary conditions takes up a minority. These make them no better than simple GNN operators in [10]. The proposed DEC-HOGNN considers implicit higher-order interactions in the governing systems and thereby outperforms usual GNN operators.
boundary nodes, which
Ablation I: On Higher-Oree Interactions. Different types of interactions can contribute variably the performance of DEC-HOGNN. this experiment, the full model incorporates four types of interactions, namely, primal and dual lower/co-bounde ablated evaluate their impact. The resulting variants are denoted as Primal-Dual (PD), Primal-only (P), Dual-only (D), and None. Table 3 presents the results in 2D electrostatic BVPs, which shows that message passing on edge lower adjacency brings minor enhancement while co-boundary adjacency matters much more. This result coincides with the governing PDE in electrostatics, dE 0, = whose inducing loss terms are supposed be computed based on both primal and dual co-boundary adjacency (Appendix B).
adjacencies, which are systematically

---

Table 4: The improvements of DEC-HOGNNs equipped with different higher-order graph convolu- tions and interactions against the vanilla GNN on the 2D-maggetosstttics.
Backbone
GCN
GAT
Graph Transformer
Both co- boundary and lower adjacency
2,450=0.055
2.246=0.350
1.222#0.036
Lower adjacency only
3,463+0.069
3.491=0.081
1.629=0.068
Co-boundary only
2.159=0.044
2.161=0.030
1.386+00.04
Vanilla GNN Improvement (%)
1.876=0.038
2.747=0.541
1.64620.026
-15.07
21.33
25.75
2.00
Fourie/
Loss 1200 spochs)
Loss 1400 spochs)
1.75
Loss 1500 apochs)
Transelver
1.50
GNBT
Geo FNO
1.25
DeC Hn Ouss
0.75
Aspect
0.50
Equi-tingle Showness
Area Fanskion Raiio
0.25
0.40
16.21
1.00
150
200
250
300
Dropped Edge Amount
Figure 6: Left: Mesh quality and model performance of DEC-HOGNN variations as more edges are dropped; Right: Model comparison meshes with different topological characteristics.
Ablation I: On Different GNN Backbones. In this experiment, study the impact of different
variants of higher-order GNNs on the efficacy of DEC-HOGNN by selecting different backbones, including GCN, GAT, and Graph Transformer (GT), 2D electrostatic BVPs. Results in Table 4 manifest that higher-order interactions can enhance model performance (GAT and GT), which is consistent with the previous results. It further shows that adding more higher-order interactions does not necessarily imply more positive sides. For instance, the appearance of lower adjacency impedes GAT and GCN, while the co- ba-uounddry adjacency much more beneficial for all backbones. Also, not every convolution can fit the DEC-HOGNN framework well just like the counterexample GCN.
Ablation III: On Performance Degeneration due to Mesh Quality. We analyze the negative impact out of mesh degradation by randomly dropping a certain amount of edges hierarchically. To adopted with arrows implying degrading
measure the mesh quality quantitatively, three indicators directions since good quality usually comes with uniform and regular elements. Further details on these tailored meshes are covered in Appendix E. The left panel of Figure 6 reflects that dropping edges from a triangularized mesh is followed by mesh degeneration. Also, minor degeneration would affect the performance violently while major one leads to salient performance drop. Note that is also infeasible to adopt classical solvers using meshes with prominent quality issues. Thus these negative effects are tolerable.
Ablation IV:
On Different Topological Characteristics. As shown in Figure 12, eight 2D magne- totastics benchmarks used for the evaluation, which named after their different topological properties by (Connected Component Amount, Hole Amount). The right panel of Figure shows that the advantage of our approach persists as the underlying topology changes.
6 Conclusion
In this paper, we propose BVP solver via integrating higher-order topological interactions, which aligns with the discrete representation of differential forms, and the resulting model is aware of differential topology. Several novel chysics-informen are also developed. Both theoretical analysis and experimental results demonstrate the advantages of incorporating higher-order interactions via integrated differential forms. Similar to traditional mesh-based numerical solvers, the performance of DEC-HOGNN influenced by mesh quality, as Whitney elements may degrade on poorly shaped triangulated or tetrahedral meshes, potentially leading to convergence issues. Extending this approach more general time-dependent PDEs and developing methods to mitigate mesh quality dependence are left for our future work.
loss terms and integrated form estimators
10

---

Acknowledgements
This work is supported by the National Natural Science Foundation of China under Grant No. 62206074 and 62472125, Guangdong Basic and Applied Basic Research Founda- tion under Grant No. 2025A1515012922, Shenzhen College Stability Support Plan under Grant No. 1XWD2002081173223300 Shenzhen Science and Technology Program No.
1C11200241 and S 1200008.
References
[1] Gouri Dhatt, Emmanuel Lefrançois, and Gilbert Touzot. Finite element method. John Wiley & Sons, 2012.
[2] Maziar Raissi. Deep hidden physics models: Deep learning of nonlinear partial differential equations. Journal of Machine Learning Research (JMLR), 2018.
[3] Zongyi Li, Nikola Kovachki, Kamyar Arizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Animashree Anandkumar. Fourier neural operator for parametric partial differential equations. In International Conference on Learming Representations (ICLR), 2020.
[4] Zongyi Li, Daniel Zhengyu Huang, Burigede Liu, and Anima Anandkumar. Fourier neural operator with learned deformations for pdes general geometries. Journal of Machine Learning Research (JMLR), 2023.
[5] Zhongkai Hao, Zhengyi Wang, Hang Su, Chengyang Ying, Yinpeng Dong, Songming Liu, Ze Cheng, Jian Song, and Jun Zhu. GNOT: A general neural operator transformer for operator learning. In Proceedings of the 40th International Conference on Machine Learning (ICML), 2023.
[6] Haixu Wu, Huakun Luo, Haowen Wang, Jianmin Wang, and Mingsheng Long. Transolver: fast transformer solver for PDEs on general geometries. In Proceedings of the 41 International Conference on Machine Learning (ICML), 2024.
[7] Tobias Pfaff, Meire Fortunato, Alvaro Sanchez Gonzallz, and Peter Battaglia. Learning mesh- based simulation with graph networks. In International Conference Learning Representations (ICLR), 2021.
[8] Alvaro Sanchez- Boechez-Gozzzlez, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, and Peter Battaglia. Learning to simulate complex physics with graph networks. In Proceedings of the 37th International Conference on Machine Learming (ICML), 2020.
[9] Ferran Alet, Adarsh Keshav Jeewajee, Maria Bauza Villalonga, Alberto Rodriguez, Tomas Lozano-Perez, and Leslie Kaelbling. Graph element networks: Adaptive, structured computation and memory. In Proceedings of the 36th International Conference on Machine Learning (ICML), 2019.
[10] Winfried Lötzsch, Simon Ohler, and Johannes Otterbach. Learning the solution operator of boundary value problems using graph neural networks. In ICML 2022 2nd Al for Science Workshop, 2022.
[11] Barbara Burke Hubbard and John H. Hubbard. Vector Calculus, Linear Algebra and Differential Forms; A Unified Approach. Pearson, 2nd edition, 2001.
[12] Karl Warnick and Peter H. Russer. Differential forms and electromagnetic field theory. Progress In Electromagnetics Research (PIER), 2014.
[13] Alain Bossavit. On the geometry of electromagnetism. J. Japan Soc. Appl. Electromagn. & Mech, 1998.
[14] Alain Bossavit. Computational electromagnetism: variational formulaiion,,, complementarity, edge elements. Academic Press, 1998.
[15] Anil Nirmal Hirani. Discrete exterior calculus. California Institute of Technology, 2003.

---

[16] Mathieu Desbrun, Eva Kanso, and Yiying Tong. Discrete differential forms for computational modeling. In ACM SIGGRAPH 2006 Courses. Birkhäuser Basel, 2006.
[17] Yunfeng Liao, Jiawen Guan, and Xiucheng Li. Curvature-aware graph attention for PDEs on manifolds. In Forty-second International Conference on Machine Learning (ICML), 2025.
[18] Stefania Ebli, Michaël Defferrard, and Gard Spreemann. Simplicial neural networks. In 7DA & Beyond, 2020.
[19] Mathieu Alain, So Takao, Brooks Paige, and Marc Peter Deisenroth. Gaussian processes cellular complexes. In Proceedings of the s International Conference on Machine Learning (ICML), volume 235, pages 879-905. PMLR, 2024.
[20] Maciej Besta, Florian Scheidl, Lukas Gianinazzi, Grzegorz Kwasniewski, Shachar Klaiman, Jürgen Müller, and Torsten Hoefler. Demystifying higher-order graph neural networks. arXiv preprint arXiv:2406. 12841, 2024.
[21] Federico Battiston, Giulia Cencetti, Iacopo lacopini, Vito Latora, Maxime Lucas, Alice Patania, Jean-Gabriel Young, and Giovanni Petri. Networks beyond pairwise interactions: Structure and dynamics. Physics reports, 2020.
[22] Christopher Morris, Martin Ritzert, Matthias Fey, William L. Hamilton, Jan Eric Lenssen, Gaurav Rattan, and Martin Grohe. Weisfeiler and leman go neural: Higher-order graph neural networks. Proceedings of the AAAI Conference on Artificial Intelligence (AAAI), 2019.
[23] Cristian Bodnar, Fabrizio Frasca, Nina Otter, Yuguang Wang, Pietro Lio, Guido F Montufar, and Michael Bronstein. Weisfeiler and lehman go cellular: Cw networks. Advances in neural information processing systems (NeurIPS), 2021.
[24] Eric Bunch, Qian You, Glenn Fung, and Vikas Singh. Simplicial 2-complex convolutional neural networks. In TDA Beyond, 2020.
[25] Lorenzo Giusti, Claudio Battiloro, Lucia Testa, Paolo Di Lorenzo, Stefania Sardellitti, and Sergio Barbarossa. Cell attention networks. In 2023 International Joint Conference on Neural Networks (IJCNN, 2023.
[26] Cristian Bodnar, Fabrizio Frasca, Yuguang Wang, Nina Otter, Guido F Montufar, Pietro Lio, and Michael Bronstein. Weisfeiler and lehman topological: Message passing simplicial networks. International conference on machine learning (ICML). 2021.
[27] Shiwen Wu, Fei Sun, Wentao Zhang, Xu Xie, and Bin Cui. Graph neural networks in recom- mender systems: a survey, ACM Computing Surveys, 2022.
[28] Wei Ju, Zequn Liu, Yifang Qin, Bin Feng, Chen Wang, Zhihui Guo, Xiao Luo, and Ming Zhang. Few-shot molecular property prediction via hierarchically structured learning on relation graphs. Neural Networks, 2023.
[29] Gaurav Gupta, Xiongye Xiao, and Paul Bogdan. Multiwavelet-hasel operator learning for differential equations. In Advances in Neural Information Processing Systems (NeurlPS), 2021.
[30] Tapas Tripura and Souvik Chakraborty. Wavelet neural operator for solving parametric partial differential equations in computational mechanics problems. Computer Methods in Applied Mechanics and Engineering. 2023.
[31] Xiongye Xiao, Defu Cao, Ruochen Yang. Gaurav Gupta, Gengshuo Liu, Chenzhong Yin, Radu Balan, and Paul Bogdan. Coupled multiwavelet operator learning for coupled differential equations. In International Conference on Learning Representations (ICLR), 2023.
[32] Jiaqi Han, Wenbing Huang. Hengbo Ma, Jiachen Li, Joshua B. Tenenbaum, and Chuang Gan. Learning physical dynamics with subequivariant graph neural networks. In Advances in Neural Information Processing Systems (NeurIPS), 2022.
[33] Masanobu Horie, Naoki Morita, Toshiaki Hishinuma, Yu Ihara, and Naoto Mitsume. Isomet- ric transformation invariant and equivariant graph convolutional networks. In International Conference Learning Representations (ICLR), 2021.
12

---

[34] Masanobu Horie and Naoto Mitsume. solvers with mixed boundary conditions. Advances in Neural Information Processing Systems (NeurIPS), 2022.
neural networks: Graph neural pde
[35] Johannes Brandstetter, Daniel E. Worrall, and Max Welling. Message passing neural PDE solvers. In International Conference on Learning Representations (ICLR), 2022.
[36] Masanobu Horie and Naoto Mitsume. Graph neural PDE solvers with conservation and similarit-equiivaranne. In Forty-first International Conference on Machine Learning (ICML), 2024.
[37] Qingqing Zhao, David B. Lindell, and Gordon Wetzstein. Learning to solve PDE-constrained inverse problems with graph networks. In JCML 2022 2nd for Science Workshop, 2022.
[38] Matthias Hein, eves Audibert, and Ulrike von Luxburg. Graph laplacians and their convergence random neighborhood graphs. Journal of Machine Learning Research (JMLR), 2007.
[39] Sharif Elcott, Yiying Tong, Eva Kanso, Peter Schröder, and Mathieu Desbrun. Stable, circulation- preserving, simplicial fluids. ACM Transactions on Graphics (TOG), 2007.
[40] Pankaj Jagad, Abdullah Abukhwejah, Mamdouh Mohamed, and Ravi Samtaney. A primitive variable discrete exterior calculus discretization of incompressible navier-stokes equations over surface simplicial meshes. Physics of Fluids (PoF), 2021.
[41] Pankaj Jagad and Ravi Samtaney. Effects of rotation on vorticity dynamics sphere with discrete exterior calculus. Physics of Fluids (PoF), 2021.
[42] Patrick Mullen, Keenan Crane, Dmitry Pavlov, Yiying Tong, and Mathieu Desbrun. Energy- preserving integrators for fluid animation. ACM Transactions on Graphics (TOG), 2009.
[43] G. R. Cowper. Gaussian quadrature formulas for triangles. International Journal For Numerical Methods Engineering, 1973.
[44] Sifan Wang, Hanwen Wang, and Paris Perdikaris. Learning the solution operator of parametric partial differential equations with physics physics-informen informed deeponets. CoRR, 2021.
[45] 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. Nat. Mach. Intell., 2021.
[46] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Andrew Stuart, Kaushik Bhattacharya, and Anima Anandkumar. Multipole graph neural operator for parametric partial differential equations. In Advances in Neural Information Processing Systems (NeurIPS), 2020.
[47] Shuhao Cao. Choose a transformer: Fourier or galerkin. In Advances in Neural Information Processing Systems (NeurIPS), 2021.
[48] Xinzheng Xu, Xiaoyang Zhao, Meng Wei, and Zhongnian Li. A comprehensive review of graph convolutional networks: Approaches and applications. Electronic Research Archive (ERA), 2023.
[49] Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Liò, and Yoshua Bengio. Graph attention networks, In International Conference on Learning Representations (ICLR), 2018.
[50] Hongyang Gao and Shuiwang Ji. Graph u-nets. In Proceedings of the 36th International Conference on Machine Learning (ICML), 2019.
[51] Seongjun Yun, Minbyul Jeong, Rachyun Kim, Jaewoo Kang, and Hyunwoo J. Kim. Graph transformer networks. In Advances in Neural Information Processing Systems (NewrIPS), 2019.
[52] L.C. Evans. Partial Differential Equations. American Mathematical Society, 2nd edition, 2010.
13

---

[53] Loring W. Tu. Differential Geometry: Connections. Curvature, and Characteristic Classes. Springer, 1st edition, 2017.
[54] George Cybenko.
Approximation by superpositions of a sigmoidal function. Mathematics of Control, Signals and Systems (MCSS), 1989.
[55] Jonathan Richard Shewchuk. Delaunay refinement algorithms for triangular mesh generation. Computational Geometry, 22(1):21-74, 2002. ISSN 0925-7721. 16th ACM Symposium on Computational Geometry.
[56] Hang Si. Tetgen, a delaunay-based quality tetrahedral mesh generator. ACM Transactions on Mathematical Software (TOMS), 2015.
14

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
The abstract and introduction are consistent with our actual contributions.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: In the Conclusion section, we discuss the limitation the proposed framework, which can be affected by low mesh-quality.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g.. independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
the paper. The authors should use their best
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
15

---

Justification:
The complete statement and proof can be found in Appendix D.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: Details on how to generate datasets for the conducted experiments are men- tioned in Appendix E and experiment settings are also covered in Appendix F.
Guidelines:
The answer NA means that the paper does include experiments.
If the paper includes experiments, No answer this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm. (b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome deseribe the particular way they provide for reproducibility, In the case of closed-source models, it may be that access to the model is limited in some way (e.g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
16

---

reproducible, they Answer: [Yes]
Justification:
repository containing data and code is provided in the abstract.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (uttps://nips.cm cc/
cesssooooociill for more details.
While encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https: / /nips.
gttcc//ggee///eeeeeeeeeee Conee/Cmdduummmii ssi assionPoll cniionoooyy) for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: Details how generate datasets for the conducted experiments men- tioned in Appendix E and experiment settings are also covered in Appendix F.
Guidelines:
The answer NA means that the paper does include experiments.
The experimental setting should be presented the core the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: The main results are reported with average and standard deviation error, as shown in Table 2 and Table 4. We conduct these experiments with at least three random seeds under each setting.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results accompanied by emror bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
17

---

The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (e.g. negative enror rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: See Appendix F.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental rums as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported the paper (e.g., preliminary or failed experiments that didn't 't make it into the paper).
9. Code of ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics https: http///eurrrp./peeee//ee httppullicttiiicccsieeeeeeeeeeee
Answer: [Yes]
Justification:
We have read it and followed the NeurlPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make to preserve anonymity (e.g., if there is special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive socictal impacts and negative societal impacts of the work performed?
Answer: [No]
Justification: Our work mainly focus on topics over scientific computation. has no significant societial impacts.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
18

---

The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: Models and data are used for scientific computation. There is no risk of misuse.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g.. code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: Appendix E depicts data details. Code and models are either original or open sources.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
19

---

If this information is not available online, the authors are encouraged to reach out to the asset's creators.
13. New assets
Question: Are assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: The data on electromagnetic static fields will be uploaded to the repository. Their details can be found in Appendix E.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moddl part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: Our work focuses on solving partial differential equations arising math and physics and does not involve human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then much detail as possible should be included the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: Our work focuses on solving partial differential equations arising in math and physics and does not involve human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere the NeurIPS Code of Ethics and the guidelines for their institution.
20

---

For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods in this research? Note that the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology, scientific rigorousness, or originality the research, declaration is not required.
Answer: [No]
Justification:
The core of this paper does not involve LLM usage.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non non-standard components.
Please refer to our LLM policy (https / wwtpe//cooeoeeece/////a/ for what should or should not be described.
21

---

Table 5: Summary of symbols and their corresponding meanings.
Symbol
Meaning
2 (M)
The space spanned by all k-form on manifold M. The dual manifold of the primal manifold M.
*M
C(X)
The simplex chain complex of topology space X. The set of all k-simplices in C(X). Its element denoted by Ok. The set of all k-cochain on Cx(X). Vo E Cx(X), Cx(X) IR, oj + (ox)..
(X)
Ck (X)
*
The Hodge star acting on k-forms. ** (M) + Ω (MM In
DEC, maps k-forms the discrete primal manifold to the discrete dual manifold, i.e., ** 2 (M) - -0" (*M)
A n-simplex with vertices Vig 2 ... 1/ in order.
dh
The exterior derivative. dh (MM
12++1 (xM).
The boundary operator. 0t Cx(X) + Cx-1(X). The function space on region Q with I,2 derivatives.
(2))
H (2)
The vector-valued function space region Ω, each component having 2 derivative.
WM(M
The space spanned by k-th order Whitnyy element defined on M. Indicator function. XA if condition A is satisified and otherwise
Xa
XA O.
A
Notation
The notation used throughout the paper is summarized in Table 5.
B Training Loss Function Derivation
The loss below is under the assumption that the dual manifold is constructed based on cireumcenters, and fields are recovered by interpolation. For brevity, we only focus on discrete manifolds in E
Div-free and Curl-free. V B = 0 + dB = 0, E (+MM. By Stokes theorem,
dB
B= Σ
B=
Σ sign (o2)h
(19)
Das
in which the partial-order relation a d' implies that is face of And thus the physics-informed 2 loss can be:
Ldiv-froe == Σ
Σ
sign 03)
(20)
Then for V x = 0+ dl 0, EM likewise, have:
Ceurt-free == Σ
Σ
sign(o1,
(21)
OGE*M 01 202
0, Ts E M, 02 xM by 03 E *M,0g M. And for brevity, will only cover the case +M since that on M is almost the same. For non-divergence- free scenarios like . = p, then the 2 loss becomes:
As for forms D, E such that -
0, one can obtain the formula by modifying
Laiv = Σ
plos) vol(os) - Σ sign (22 F3 0z)ho
(22)
MEERM
Boundary Condition: n n-B =B=0. Consider tetrahedral 03 [A, B, DJ with positive-oriented faces [A, B, C), (C, D, A), [A, D, B), [B, D, C), as shown in Figure Let hi denote the altitude to
the opposite face [B, D, C] in d3, WA the nodal function at vertex A, and similarly for the others.
22

---

C
nn
B
hn
D
A
Figure 7: Tetrahedral [A, D,B,,,CK C,
AB indicates the unit vector pointing from to B. Let H -hahphe, H'' -hghghchp Hhp and V = vol(o3). Without loss of generalityWLLOG), set [A, CJC to be the unique boundary face in Jg with outwards normal np.
can be easily shown that - A VWA = np [14] and the interpolation fields on [A, B, C1 only come from itself and its upper- adjacent neighbors. We first check whether its neighbors contribute normal components. WLOG, observe [C, D, A] and we have:
W[C.D.A] nD = WA (Vwc x VwD) up + we (WD x AWA) np + wD (WWA X Vwc) np =
(23)
The first and second terms vanish since - hD VwD = no while the last term vanishes since wp = on face [A, B, C). Therefore, the only normal component contributor is [A, B, CJ itself and have:
f
[A,B.C] W[A,B,C] ny dS 2
Σ
Waww Vwc) nn ds
(24)
|A,B,O
eye
[Os 0 ligii
cos
Bs we
[A.B.C]
WA
loDA WB
A cos 0 cos cos
dS
(25)
2
Cos OA +
Vo,,[A, C1)
(26)
linli
lehl
haho
2vol([A, B, CI)
Σ eye
hA 0A
(27)
3H
which 0 is the intersection angle of DĀ and np, сye indicates the cyclic summation with respect to the tuple (A, B, C). Thus in the sense of average approximation, the loss boundary condition n B 0 is
vol(a2)
C.
Σ
h2
Σ eye
ha Co
(28)
Ho
GGE0AM
Boundary Condition: H O. WLOG, assume only [A, BJ is on the boundary in as and it is obvious that the fields on [A, B] comes from edges {o| 3og, a2, of < <o' < 0y A 01 < < 0s)- First, observe that none in this set, except [A, B] itself, has tangential component contribution. WLOG,
Wac. AB (WA Vwc - we VwA) AB =
(29)
The first term vanishes since Vwg. -B 0 while the second vanishes since we on edge [A, B). And thus the circulation on edge [A, B]
is:
/
W(A AB
NOABug- Wp VWA) .AB
(30)
AA
vol
BJ)
(VwB AB + VWa BÅ)
(31)
vol B))
cos cos
(32)
h hc
23

---

A
A
D
B
B D C
E
C
(a) Estimate the integrated dual 1-form on [A, DJ via integrated primal 1-forms on ö[A, B, C).
(b) Estimate the integrated dual 1- form on [A, E, D) via integrated primal I-forms on 0[A, B, C, D).
in which cos AB = ng AB. Whence the corresponding loss is:
LaxH = 1 Σ h2 vol(o1)
COS 0g
(33)
h
BEdo
Boundary Condition B = 0. This condition is less common than B 0 in electromagnetism. involves more terms since all faces of By have contributions. The term from [A, B, C) is:
[A, = f [AB..C
C
cos HD + WA hARB
cos
cos Oc HC
WA
HÂ + WA
(34)
habi
lohh
([,, B, CT)
h cos 0 HA
(35)
3H
eyc
The terms from [C, D, B), D, D1 C are respectively:
[.B.C.I...... = hDll
V
(h bc
(36)
V
[A.BB.C,A..D.III =
hDH (BBCB haCh)
(37)
V
I...B.C.BB.D.C-- hpll (heAC hBAB).
(38)
Thus, the loss is:
CaxB =
Σ
sign 02 (en(or)hooto Io +
Σ
sign(o)) I To
(39)
G1E0-M
C Estimates among Integrated Forms
For triangle region [A, B, C] whose diameter is dominated by and a smooth function f. then direct expansion gives: Vx [A, B,C.
JABO ((... = S(A,8,c
+ ( = S(nn (40)
And likewise,
f B.O f(x)dz = [CI (
(41)
Let ABB = vol([A, B)) and assume that |BD/ > =), then
/
/
w(x) vo, DJ) + -(E3) [AB.DD]
dw
[A.B] +
w +
(42)
[B.D]
(D.
24

---

J/c +
1 S.2.)
w(x)) vo, + 0(5²)
J.....D du
+
(43)
[G,D]
J.D.A]
Note that w(x)= = w(x) + O(e) and we have:
H.B.DD dw dw B,)[
+ BD
+ +(())
(44)
[ dw Vol DI) |DC+
DC
and
J.A.D.D du JAA. u (BD)/ /
S.D.A + O
(45)
Ju..c.m) dw Jaq w (ICD)/ BCI) Jnn.q w [[D.A] + +
Combining Eq. 44 and Eq. 45 and rearranging the terms gives:
JAPI" w
BD] BC|
f
|CD| [BC| [A.B]
w+
(46)
|A,D]
[A,C]
Likewise, one can mimic the process above and obtain the estimation of dual 2-forms by primal 2-forms:
J A.E.D.I w
BE BC| J.C.D
|CE]
w
w O(E))
(47)
[C
A.B.DI
This can be generalized to high- digh-dimensional cases. For brevity, let I(or Sa w and
[Uig .. 22 w). == [Pior .... 22 W, e) and we have:
( + vol( [[) Sr, dw
Son w
(48)
vol( u|] Sr. dw Sors w
the right-hand side (RHS) equals:
ΣØ= n 1=0 (-1))//(uuu (1,, (-1)m-1 [1 2 w)) + (-1)₂ 10/0//1/1 ( (1ii.. 21 w/)
(-10 w) + (-1)²- n- u() + (-1) (( Pi) ul)
w))
RHS
(49)
> -= (- voll 10 (,, w) (-1)n- ((Ui w|) (-1)"
0(s"-1)
Σ voll (0.0
10 15(1,,) + (-1)n- (( 21 uI) + (-1)a
2, w)) + (6)
(50)
in which (u, w) =1 w)), 1 (w, v) ( =[(0 w. 2l). This gives:
vol((u; w|)
(-1)a-1 ((10) w) + (-1) (1/11122))) uI)
vol([w,)) vl) (-1)n- 1 (( u)) + 1) (1i w]) + 0
(51)
I(Weenn-g w|) I( [M/mm-- u|)
0(5a~1)
(52)
T( T(Vijnn-2 w))
combine with the volume decomposition identity vol((u, w)) vol((w, v/) vol(uu, v)) and have:
I I(Viin- w|) vol([w, e)) (3)/((en I((0imm-a u/) + vol [u, w|) vol([u, ul) I( I((utm-29 b)) +0)))=0000 O(E"-1).
(53)
vol([u, e))
C.1 Estimate Integrated Dual Forms from Integrated Primal Forms
We only propose ways to estimate dual I-forms in IR? and R and 2-forms from primal I-forms and 2-forms since usual integration quadrature can be leveraged directly. Apart from that, we further show how to estimate dual I-forms from primal 1-forms in R". Higher-order estimates and higher dimensions are left for future work.
O is a non-overlappling vertex with vertices A, B,C s.t. O span A, B, C). Specifically, 0 the dual node and the circumcenter of [A, B, C), then ([O,+[B, CI) (T(JO, B) - ((O,Cl))/2 Therefore, reduced to show the relation between T((O, X]) and T(XXy]) where X, Y {A, B, C).
25

---

A
B D
C
Figure 9: Estimate D2' (+M) by O2 (M) in
vo ([A, DJ) and I[AA, DJ) be estimated by I[AA, B)) and T[[A, C)). And appears done. However, it is possible to avoid introducing new vertex D, which case it is easier to extend to higher-dimensional cases. For brevity, let S == abs(vol (B, S = vol([A, B, CT) and SB, Sc is likewise defined.
Note that ([AA))
vol((A.0)
vol(IC, D))
vol([A, C, D)) vol([B,]) [, vol( Se + SBB S SA
vol(
SB
SB
(54)
And thus:
O[[A,0)) = vol([A,0)
|vol[[C, DJ) I([A, BJ) vol([B, D)) I([A,C)
(ol DJ)
voo[[B, CJ)
vol([B, CT)
S- SA
Sn I[[A, BJ) + Sc I[A,, CT)
S
S SA
S-
SB
I[A, B)) + Sc
I([A,T)
(55)
8
Combine with I([O+([B (T([O, Bl) + I[OO, C1))/2, we have:
SA I(AA, Bl) + SB Sc T(BB, S
SA ((C( AI)
[[[O,{[B,C]]
(56)
S
namely,
100, +[B, C()
SA -SB
SB Sc
SSA
I[[A, B)) IiB. C I[CC,
=
SB -Sc
- SA Sc
+O(=?)
(57)
28
SA
SB
Sc
However, this linear transform is not invertible since its determinant zero. This implies that one cannot recover the integrated primal forms simply by this transform formula. We will later turn this topic in Appendix C.2.
Fortunately, this formula can be easily extended to dimension with the volume decomposition identity. For on = [to. ..., Un] and u E span Tn we have:
Σ(-1) vol((T(o-x) ul) vol(o.)
(58)
in which - on {vx) and is the permutation of the rest vertices. This be easily verified by the identity Eq. 59, in which there are only n non-zero terms within 2n terms in total.
vol(a,,) = det(v1 (vi vo: 1/ uo....., On - vo)
det((e) - u) + (u - oo).....
(Un - u) + (u - vo))
= Σ (-1) vol([T(a
(59)
Then we claim:
Theorem 4, For simplex (u,, U. Wi,..., Wr)). vo int
vol((o,o_wl)
[(u,o])= Σ
(Iu, w]) + O(s3).
(60)
vol(o)
WEO-N
26

---

Proof. WLOG, consider a simplex o = (u, Wi... Wn), then for & span(w;), then Elp
[Wi 2) n [u, E [m]m]] pl.
o( vol((u, o))
[(uu, pl)
(61)
vol((u, pl
vol([u, o[)
/vo[(pp, 9l)
vol([u, q)) I( [[u,, o]) + vol(lu, pl)
I([u, 9])
(62)
vol((u, pl)
vol(lu, gl)
vol(uu, o[) vol((u,pl)
vol((p. gl)
((u, v/) + vol(ee.p)
vol(lo ( 4)) (v, w|
[(u, w|)
vol((u,g]) q])
vol((u,g]) wEr-v
(63)
Whence we only need to check the coefficients Eq. 60. We will use property of vol. If span span and vol vol r, then for any simplex 5 whose vertex is not in span we have vol([o, 6]) - vol([7,6]).
The first term is:
vol((u, o|) vol( g))
vol([u, o, q)) vol [u. p. g))
(64)
vol( (u, pl) vol( g))
vol((u, z)) vol (u, v, g))
vol([u, o, q. {w;}-x])
(65)
vol([u, v, {wo]-x])
vol((u, o, (w,)l) vol((o,
(66)
vol((u, (,<e))
vol(o)
The second term is: vol([u, ol) vll( vol((o-o-u 9))
vol((u, ol) vol((v,p. a- vol((o 4)
(67)
vol((u, pl) vol([u, d)) vol((o-u-w w|)
vol((u, pl) vol((e, q, o-w-ll) vol((u, ol) vol((p, o-w))
w|)
(68)
vol((u. (I) vol vol((o, a-w]) vol
(69)
This completes the proof.
By Eq. 60, we obtain all o(),, E a. But what about T[ o))? Note that [x o) C
[o, and thus apply Eq. 60 once again in (o, au), denoted by Tu Consequently, have:
vol([[Ta) o)) [[[e, o)) + o(E)
(o0, * = =- Σ
(70)
vol
vol(((Tu)-,,)) (olllo-w,he I[[u, w)) + O(=2)
(71)
vol
vola
VεT WEO-
Note that for simplex a (x (m u, U, w,.... Wr). there exists and only exists one intersection between (u, o) and (v, wim) for all o E into and also only one q intersected at (v, p) and (wI, ..., Wm).
Therefore, combine this result with Eq. 53, repeat the proof above and we have:
Corollary 5. For simplex (I1, ...Tm U, U, WT, ... Wn). Vo into,
vol((o, [([21mm u, w)) + (
[[ o))
(72)
vol(a)
wg (
The estimation of dual 2-forms by primal 2-forms, special case, is:
vol([B, Fl)
[[[A,C, DJ) + vol((C, F?))
I[A,E, E, DJ) vol([B, C1)
I([A, B, D))
(73)
vol([B, CT)
(IA, vol([F., D)) vol([D, H))
I[AA, H, DJ) = vol((FE, D)) vol (A,, H)) vol([,HH) vol([A, Ol) ([A, OJ
(74)
27

---

B
A
D
B
A
E
(a) Estimate primal forms by neighbors.
(b) Estimate primal forms by sub-division. Figure 10: Illustrations for integrated primal form estimation.
I[AA, DJ) vol([A, OJ) vol( H))
El)
ET) vol( I ([A, B, DJ)
vo I[[A, DJ)
vol([A, H)) vol( D])
(75)
vol([A,, 0J) [A, H)) VB + Ve
Ve B,)
VB
[[A,, C, D))
vol(1O, H]) vol (l),O)
Ve+ VB
Vc + VB
(76)
((A,,, + VB Ve Ve V
V ([A,C,D)
(77)
C.2 Estimate Integrated Primal Forms from Integrated Dual Forms
One may ask why we need to estimate the integrated primal forms by dual forms. After all, the dual manifold not tetrahedralized and thus causes much more complexity estimation. stems from the intrinsic requirement in our previous proof. The approximation So dw/ Jr dus as vol / vol r, fine if w not exact, i... di + 0. Any exact forms cannot be properly estimated by the aforementioned quadratures.
But, fortunately, the converse estimation remedy this limitation. Since See d'w/ Ser d'w se vol #oT vol *7 fails if d' O. In smooth cases, forms on M are also forms *.M. If this also holds in the discrete cases, then by Hodge decomposition, the approximation valid for any form up harmonic form.
Two methods can be applied: sub-division and estimating from neighbors. The latter has a default: it fails if this dual edge is the boundary, i.c.. it has no neighbors. For brevity, only estimating primal 1-forms from dual 1-forms TR2 is discussed, which, however, is not challenging to extend to higher-dimensional spaces and forms.
As shown in Figure 10a, given integrated dual forms, integrations of [A, [B, D), [A,
[C, D), [B, C] are all tractable. Then again by Stoke's theorem:
SAB.D du [a. w (B.D) w ST..A vol (AA JIA.C.D) du Jin. la w I...A w
vol([A, B, DJ)
(78)
in which JJA.DY such a quadrilateral AB DC does not exist. Figure 10b offers another approach. By subdividing the primal triangles, note that [A, D) is contained in both [A, D, B] and [A, D, E] where the integrated forms are known on every edge except [A, D). Therefore, replace C in Eq. 78 by E. and we obtain the estimate.
be estimated. But this approach fails when [4, D] is on the boundary since
28

---

D Universal Approximation on Poission Equation
In this section, we show that for any properly exposed Dirichlet boundary value problems, there exists a proposed Higher-order GNN network that approximates the continuous operator. Our proof follows the famous result on the universal approximation property of neural networks for operator.
In our higher-order architecture, we encode the field vectors into differential forms and decode them via Whitney interpolation. Formally,
NN D LN... 01 3
(79)
in which L; E, D are the i-th layer, encoder layer and decoder layer, respectively. We first show that this encoding decoding method has a pseudo-inverse property and an asymptotic identity property in proper function space. First, we consider simple but heuristic example, real-valued function op (T) on a closed interval C TR. WI.OG, := [0,1].
Proof of Theorem 2.
Proof. The existence of Dx is trivial since for all E R, there exists E R such that fii1 x)dr + ((ztzz44)/2 a(x - (z; + zz+1)/2) c. Since there are only finite non-smooth points, i.e., the interval endpoints and the midpoints, one can use a smooth mollifier there without changing the integral.
[u (()+ a(
Recall the canonical mollifier 7 := /x////// [52] and ffe := -7=* ff. Since y E (I) is continuous, *== - &p uniformly on the compact set I in the sense of L Thus, we only need to show that f + + 2 Since (e: & E (T), we have M1 == supre/ (x). oo and Mo == supzer max (f(xx),((x)))
Direct expansion gives:
Ri11
(x141 s)?
Ps (2)d
(s(()) + ()(())(z s)dx (s)2 N (((Ce)
(80)
2
2z
Since fc= + ip uniformly, there exists E small enough such that VEn > 0,
50
+ (=(2)2 + +C(((s) (ri41 s)²
xz,)²
((1) + ((G))
2M,
2
which gives:
-s) x(t) s $2 1 [(((sa)
(rit1 s)²
+
(I++1 - +t²
+
Eo
(81)
2
2
2M
2 M22 NN
Eo
(82)
2.M
< 6 ED
(83)
2M1
Thus NN > log, / such that Vs, E [0,, (s) ((() Whence with 2o/2Mo
i (v(x) - (xx) = / (1((=)=))(d])) (==))lp((a) ((x))(() |dx + / [Vec(o) ((x))l(((a) ((x)
(84)
< 2M + 8o ) + 2Mo 6n
(8.5)
2M,
2Mo
350
(86)
Corollary 6. For uniformly bounded space By (00, 1]) = {i E L)) ||/| < M).
Vip (0, 1)). exists M such that for all N M.
|| - DN Ex(r)ll <..
The proof can be analogously extended to higher dimension with C2 boundary by hypercube- partition.
29

---

Corollary 7. Given a [2, 1)4) n (u: ||u||x < M), EN encodes the line integral along every edge of Xa while D preserves them. Then there exists ò such that for supa diamXa ö,
of [0. 1] [Xa]with measure-zen intersections, for v e
u Do Es(uillea <.
Proof. Consider hypercube partition where the edges align with the coordinate axes. Therefore u t can be viewed a function on Ii -axis, which is reduced to Theorem and we have:
u u- DN o Ex( S
Σ I - DN o Ex())) tull2 < dME ==1
(87)
Specifically, in IR' encoding u onto faces in the form of flux through the face is also proper since n, is equivalent to encoding that to the edges. It because faces of codimension in essence.
Define the projection operator PN == DN o Ex- For Dy. the degree of freedom of Dx is 2N and the interpolation is usually continuous with the encoded integral features in -e-see. And the integral operator yn is continuous since EN is linear and bounded, as the integral domain is compact. Therefore, Px is continuous with proper interpolation.
Now consider the continuous functional family C(K) on compact set K. Let V cC (K) BMr in which B BKK) == {f K +
< M). By continuity, PxV is still compact. Consider the space N U PNV and W == = W V VUU P.V.
R
Lemma 8 Arzela-Ascoll X is Banach space and K cc X. A subset V C() has compact closure if and only V is uniformly bounded and equicontinuous.
Theorem 9. W is closed, uniformly bounded and equicontinous. By Arzela-Ascoli's 's lemma, W is compact.
Proof. Consider a function sequence {.n]. If [fr} has subsequence in V, then it converges to in V since V is compact. If [Sn} has subsequence (i.,) such that sup, &n < M, then it is again completed as the finite union preserves compactness. Therefore, we only need to consider the case (P.,, in which fa E V, in increasing to infinity. Note that { must have subsequence converging to in V and thus:
IPrr, - 0ll2 I - + I1f - 2ll2 < 8
(88)
Thus lim, Jn 9 a.c. W is bounded due to |f < M and Poincaré inequality. W is cquicontinuous if the interpolation Lipschitz. For instance, the linear interpolation has a Lipschitz constant 2/C..
( f(x)dx
f E.+1
$
f(x)dxx
12
(89)
< C2-N
fx(x)ddx
(90)
C2
(91)
Whence functions in W are all Lipschitz and have a common upper bound C/CM and are therefore equicontinuous. By Arzela-Ascoli Lemma, W is compact.
Consider continuous convolution operator
[[g](r) = 1 f(x y-g
(92)
and suppose that G V
W, then have:
[[و[ NN
g]1l = ||وو|| P.oll
(93)
< |[][g]
||(1 - I..
(94)
S sup
E(o))l
(95)
New
30

---

Note that W is compact and so is GW == [G[g] W). therefore
lim sup || - P. Pb)][[]] =0 R=000 NeW
Note that in the linear interpolation scheme, [|D(h, h2)ll2 |ل h2/l1. Thus, the only thing to show is that o glol - o E[g] can be bounded with a large enough N.
Theorem 10 (Green Function [52]). The Green functionn of Poisson problem in R is:
G(x,y) = (n 2)
||x - y x y
n 223
(96)
In
2x
in which on is the surface area of a n-ddmmnsiooall unit sphere and the solution for Dirichlet problem
Au(x) = f(x), x 12; = g(x) = 0, E an
(97)
given by
/ G(x G(x-y)f(yddy / O,, y)gly)dy / G(x - y) (98) 12
u(x)
Jan
Lemma 11 (Extension Theorem [52]). For bounded closed with C -boundary, there exists an extension operator E (.. + WR*,(R"), such that E(u) in Ω and E(u) is compactly
supported and continuous:
Eu (W S c..allulllxx=0000))
(99)
For bounded closed Ω C R", it must be paracompact and thus has smooth partition of unity {p.] [53]. For square-like uniform mesh, one can construct finite open cover {O,} } such that 45 > 0, there exists a finite open cover suPaæt diam(O, O Og) < & and O₂ contains the partitioned hypercube a. For those hypercubes not intersected with the boundary, the corresponding open covers can also have no intersections.
Recall that the extension theorem can be shown by the reflection method, in which can locally extend ulo,ro = peu, an c U O due to the boundary compactness. It obvious that
p(supp(Eu)) 2 p(supp(u)).
(100)
in which y is the Lebesgue measure. And also note that
u(supp(Ew) s 1 supp(u) U supp
Σ R. [peul]
(101)
< (suppp(u) +
(102)
==(spp(()) +
(supp(aa))
(103)
u(supp(u))
(supp(a))
(104)
S (supp(()) + 60
(105)
in which R, is the reflection operator concerning the open set O,, and the last holds because ass the partition gets finer, the total volume of boundary hypercubes tends to be zero. Therefore, f defined on Ω can be extended to R such that Jq G[y) y)dy fg- E[f(x y)dy + s. For brevity, we omit and identify Elf] as in the following discussion. Consider component of
31

---

Eo G[g), i.e., the integral on I, and we have:
J. (ik--)f(y)dduue J. G(y)f(x- -
(106)
Σ I1I,
-1,,,as)/((x-- ((y)((x-
y)dydz
(107)
= Σ J.J, -G))//x--)))-->VVxxx-- - -y--yyy-y -yy)dddr
(108)
E1, 11x - ys)de J, Gly)dy - J.J, C(y)VS(X- ( (i) -y
(109)
T1m f(x)dx G(y)dyy J,J, [(x)V ()). ys yo)dydr
(110)
E +
(111)
Note that Jj fy)dy are components of Elg] and thus the neural network L is required to learn N² coefficients in the form of S1 f(z y )dr, continuous function of æ. But here the integral is on the translated J, i.e., 7 + YJ the sense of Minkowski sum. Since f extended to R. such that = 0 a.e. on R 52 we can omit those + y. not in the mesh. Therefore, encoding f(y N integrals is still enough to estimate Ti,r,J-
For T2,,.,, , in sub-regions with fxx is not constant 0, as N + co, we must have [f(x) > 8 due to continuity. Then we have
12, S J, |G()) -1V (0((x-x55)) -c)lll- Iyy yolldy s M diam(J)
| G(r)dy
(112)
Thus as N + xo, dJa + 0 = T2,1, oLS, G(y)dy) while Tı,,, = OUJ G(y)dy). We conclude that 72,,,J = o(Ti,r.)) and therefore, we only need to focus on the main part. The first term can be interpreted (N*)² coefficients. This can be approximated by the canonical universal approximation theorem [54]. This seemingly holds for special meshes, but in fact, one can adopt a bounded Jacobian deform into a more generalized mesh.
matrix multiplication where receives
B as the input and outputs
Note that the operator norm of G[ == Jn G(x y)f(y)d dy is finite long Ω is closed and bounded. Then by Proper Encoder-Decode theorem, we reach Theorem 3. Following an analogous proof, can see that it also holds for non-homogeneous Dirichlet conditions. Thus, we have the following corollary:
Corollary 12. For oon-homogeneous Dirichlet problems, Au(x) = fxx), x E Q u(x) = g(x). x e 82, there also exists a network universally approximates Take gradients on both sides, there exists a network that universally approximates Vu, i.e., the electric intensity E.
E
Dataset Details
Table 6 lists out details of several benchmarks. 2D meshes are obtained by adding random points into a 2D convex hull and then triangulated via the open Python library Triangle [55]. 3D meshes come from example meshes provided by free C++ tetrahedraliization library TerGen [56].
2D cases. The 2D elelectrosstatics case involves triangle pieces constituting square, including two media: silicon and a non non-isotropic linear material. Each piece is assigned random charge density excitation, and the region boundary condition is set as equal voltage. The 2D-magnetostatics case includes mu-metal (a material with a non-linear B-H curve) triangle pieces, composing an irregular mesh with two holes. Each triangle randomly assigned a current density. The region boundary condition is balloon, implying infinitely large region.
32

---

Table 6: Various meshes adopted in different BVPs. ELE and MAG are for electrostatics and magnetostatics, respectively. NA means not applicable.
Mesh
Square
Holes
Socket
Gear
Vertex amount
105
180
2,407
2,453
Edge amount
277
487
12,865
12.833
Face amount
173
306
19.032
18,843
Cell amount
NA
NA
8,580
8,462
Train set size
120
120
240
240
BVP type
2D-ELE
2D-MAG
3D-ELE
3D-MAG
Figure 11: Illustrations degraded meshes. from top to bottom.
The edge amount decreases from left right and then
Dropped Edge Amount
0
25
50
100
150
200
250
300
Vertex amount
180
180
180
180
177
168
155
114
Edge amount
487
462
437
387
334
275
216
139
Face amount
306
281
256
206
156
106
60
24
AR
1.286
1.301
1.318
1.365
1.401
1.458
1.595
1.929
EAS ATR
0.595
0.594
0.607
0.670
0.783
0.811
0.975
0.816
0.197
0.251
0.316
0.404
0.487
0.584
0.657
0.860
Table 7: Details about the generated degraded meshes.
3D cases. The 3D-electrostatics case is an AlzOg-based socket, in which four SiO₂ cubic blocks are embedded with uniformly distributed random charge density. The solving region boundary condition is set uniform electric potential to simulate enclosed empty cell in metal. The 3D-magnttsstatics case a cobalt quarter-geal with three copper circuits with uniformly distributed random currents. The region boundary condition is that H is tangential to the boundary and the integral on the boundary vanishes.
All coordinates are normalized to [o, 1] d = 2, 2, scalar fields are standardized, and vector fields are divided by average norm eliminate scaling differences during pre-processing.
Degenerated Meshes via Dropping Edges. By randomly dropping a certain amount of edges and further deletion fix the topology, these degraded meshes illustrated in Figure 11 and depicted in Table 6 are obtained. Usually, mesh is considered having good quality if its elements are uniform and regular. We briefly introduce these three indicators here:
Aspect Ratio (AR). There are many versions of AR. We define it as the ratio between the maximum edge length over the minimum.
33

---

Figure 12:
Illustrations meshes with different topological characteristics. It is clear that those colored identically have the same connect compone
Equi-Angle Skewness (EAS). An N-polygon 's ideal angle is 0 = x(N 2)/N. For an actual angle 9; its offset := |8 / max 0). And its EAS max; 6,.
Area Transition Ratio (ATR). It depicts the uniformity of element size. Let Ai, Aj be the area of two neighboring elements. Then its ATR is |A; Ajl/ max Aj}.
F Implementation Details
The hardware configuration consists of an Ubuntu 20.04 LTS operating system platform powered by four NVIDIA RTX A6000 Graphics Processing Units with 48GB memory. We use the PyTrrch framework version 1.13 alongside Python 3.12 interpreter. Model optimization involved selection of systematically chosen from the discrete value set (i x 10-3 : E [1,5). E {1, 2, 3, 4,5} to achieve optimal performance metrics. Datasets partitioned into training, validation, and test sets with the ratio 8:1:1. Each model in the main result is trained until the validation loss converges within 1,000 epochs.
the initial learning parameter 7 for the Adam optimizer, which
34


Reviews:
[{"review_id": "OKxp7EgjQR", "reviewer": "Reviewer_5x2r", "summary": "solve boundary value PDE problems. The authors incorporate discrete and finite element exterior calculus into GNNs to deal with geometrical relationships between differential forms (e.g., charge density and electric field). They also derived the physics-informed loss to measure the consistency of the output against the considered physical law. In addition, the authors provided the universal approximation theorem for Poisson problems, justifying the use of the model for electromagnetic phenomena. The experimental results suggest that the proposed method has a high capacity to predict electromagnetic phenomena for 2D and 3D settings.", "questions": "* In the experiments, how much was the error of PDE constraints, e.g., divergence-free condition?\n* What would be necessary to extend the method for other type of PDEs, e.g., Navier–Stokes equations?", "limitations": "yes", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n* The writing is well-construcgted and easy to follow. The authors reveal the connection of the present method to existing ones from the viewpoint of DEC (discrete exterior calculus), clarifying the position of the research.\n* The derivation of physics-informed loss for GNN is useful. Thanks to the DEC formulation, a certain physical constraint turned out to be the sum of integration, which is easy to compute.\n* The method is demonstrated to have high accuracy compared to the considered baseline models, showing the expressibility of the model.\n\nWeakness:\n* The relation to existing exterior-calculus-based or simplicial-complex-based methods is not clearly stated. For instance, differences and advantages against research on cellular complexes [Alain+ ICML 2024] and simplicial complexes [Ebli+ NeurIPS Workshop 2020] could be stated, and possibly, compared experimentally.\n* There is no evaluation of computation time. Since the classical solvers have high accuracy (with possibly a lot of computation resources), the evaluation of the accuracy alone is not enough. Therefore, the reviewer recommends evaluating the speed-accuracy tradeoff with changing spatiotemporal resolution (at least for classical solvers).\n\nMinor points:\n* Simplex complex (l.63 p. 2) could be called simplicial complex instead.\n* Since x is not used in the PDE, $x \\in \\Omega$ in Equation 6 is confusing. Either using x in the PDE or writing $\\mathrm{in} \\ \\Omega$ would be preferable.", "was_revised": true, "final_justification": "Although most of the concerns are addressed by the authors, the evaluation of the speed-accuracy tradeoff compared to classical solvers, which is the most essential one, is yet to be addressed. Therefore, I keep the original score.", "rebuttal": "We are grateful for your thoughtful review and constructive suggestions. All your comments have been carefully addressed in the point-by-point responses below. \n\n> **Q1**. For instance, differences and advantages against research on cellular complexes [Alain+ ICML 2024] and simplicial complexes [Ebli+ NeurIPS Workshop 2020] could be stated, and possibly, compared experimentally.\n\nThanks for pointing out these important works on incorporating higher-order topology. A crucial observation is that one can well define the neighborhood of a $k$-cell as a set of $(k-1),k$ or $(k+1)$-cells, respectively. Since lots of canonical results of graphs are based on a certain type of adjacency, one can generalize the usual graph Laplacian to higher-order Laplacians, based on which *Ebli et al.* proposed a higher-order GCN while *Alain et al.* derived a higher-order graph Gaussian process. **The message passing mechanism in our paper actually covers Ebli et al.'s work**. These methods take advantage of multi-body interactions in different scales in a topological structure since a $k$-cell can be interpreted as the result of the interactions among $k$ individuals (nodes).  \n\nOur approach is a step beyond merely leveraging topological inductive bias. Indeed, by connecting higher-order features with differential forms, it **paves the way for further leveraging geometry inductive bias** since forms are more than commonplace in differential geometry. By encoding input data as forms and building their relations properly, our method is able to turn differential operators into higher-order feature interactions. In all, our approach is **an attempt to open the door towards the geometry world**. \n\n> **Q2**. Since the classical solvers have high accuracy (with possibly a lot of computation resources), the evaluation of the accuracy alone is not enough. Therefore, the reviewer recommends evaluating the speed-accuracy tradeoff with changing spatiotemporal resolution (at least for classical solvers).\n\nThanks for your suggestion. The neural operators often outperform classical solvers by **several orders of magnitude** in terms of speed and thus the time cost comparison is less talked about. As stated in FNO[1] (Sec. 1), \"*On a 256×256 grid, the Fourier neural operator has an inference time of only **0.005s** compared to the **2.2s** of the pseudo-spectral method used to solve Navier-Stokes.*\"\n\n**In our setting, the Ansys 3D electrostatics solver takes 14 seconds in total  and 427MB memory for adaptive meshing to reach 0.1% energy loss, let alone time-consuming work on preliminary mesh checking.** Conversely, most prevailing neural operators are on the milisecond level. Furthermore, our experiments mainly concentrate on the effectiveness of differential-form neural representations, and thus only time-independent PDEs are considered, i.e., it **does not involve spatiotemporal resolution**. \n\nBut still, to validate our approach on time efficiency, we compare our methods with neural baselines. We kindly invite you to refer to the response to **Reviewer RTKj 's Q1** where a comprehensive comparison on computation cost against neural operators is presented. \n\n> **Q3**. In the experiments, how much was the error of PDE constraints, e.g., divergence-free condition?\n\nLet us further clarify some technical details on this topic. \n\nIn the magnetostatics experiments, we optimize the divergence-free loss $\\mathcal L_{\\text{{div}}}$ weighted by face area (in 2D cases) or cell volume (in 3D cases), and the data loss $\\mathcal L_{\\text{{data}}}$ simultaneously. Note that divergence-free loss **should be viewed as a requirement instead of a regularization term**. Otherwise, it may overfit the data at the vertices and the complete field recovered via interpolation can thereby go against physics laws. \n\nA further experiment is conducted: we set the total loss $\\mathcal L:=\\mathcal L_{\\text{{data}}}+10^2\\mathcal L_{\\text{{div}}}$ and then observe the training curves with and without this extra divergence-free condition. \n\nWe define the ratio $R:=10^2\\mathcal L_{\\text{{div}}}/\\mathcal L$. The higher $R$ is, the more likely that the model overfits the data at the vertices and breaks the physics law. And also, we study whether this extra term would have a negative impact on fitting given data points by observing the validation data loss $\\mathcal L_{\\text{{data}}}$. The table below shows two validation curves within 500 epochs. The result shows that **merely fitting the data points is likely to violate the divergence-free condition**. \n\n| Train with $\\mathcal L_{\\text{{div}}}$ | Item                            | 0     | 100    | 200    | 300    | 400    | 500    |\n| -------------------------------------- | ------------------------------- | ----- | ------ | ------ | ------ | ------ | ------ |\n|                                        | $R$                             | 1.05% | 51.76% | 74.89% | 85.88% | 90.09% | 92.47% |\n|                                        | $\\mathcal L_{\\text{{data}}}$    | 9.44  | 4.39   | 2.97   | 2.26   | 1.81   | 1.54   |\n|                                        | $10^2\\mathcal L_{\\text{{div}}}$ | 0.10  | 4.71   | 8.86   | 13.74  | 16.45  | 18.91  |\n| √                                      | $R$                             | 0.73% | 1.56%  | 4.00%  | 5.23%  | 6.67%  | 8.98%  |\n| √                                      | $\\mathcal L_{\\text{{data}}}$    | 9.47  | 5.06   | 3.60   | 2.90   | 2.52   | 2.23   |\n| √                                      | $10^2\\mathcal L_{\\text{{div}}}$ | 0.07  | 0.08   | 0.15   | 0.16   | 0.18   | 0.22   |\n\nThus the conclusion is: **without this divergence-free loss term, the model is prone to overfit the data at sampled points, hence, greatly violating the divergence-free condition.** Though sometimes fitting slower at given data points, some properties of the field are preserved, like *divergence-free*, coinciding with the insights of preserving symmetries in *Lagrangian Network* and *Hamiltonian Network*. \n\n> **Q4**.  What would be necessary to extend the method for other type of PDEs, e.g., Navier–Stokes equations?\n\nAs stated, our current method can cover various linear PDEs including common operators like $\\text{grad, curl, div}$. While some geometric neural representations are needed to achieve non-linear PDEs. Below, we briefly demonstrate some possible technical approaches for solving Navier-Stokes Equation (NSE) and leave them for future work. \n- The viscosity $\\mu\\Delta \\mathbf{u}$ term in NSE can be decomposed into $\\mu \\text{grad } (\\text{div }\\mathbf{u})$ and $\\mu\\text{curl }(\\text{curl } \\mathbf{u})$. In 2D cases, the integration of $\\mathbf{u}\\cdot \\mathbf{n}$ on an edge $e$ ($\\mathbf{n}$ is the normal of a primal edge and also the tangent of its dual counterpart $\\star e$) can be interpreted as the flux through $e$ and thus a 1-form on $\\star e$. This again represents $\\text{div }\\mathbf{u}$ properly via the *edge coboundary* in HOGNN, while $\\text{grad}$ can be implemented on the dual graph based on graph calculus. The other term $\\mu\\text{curl }(\\text{curl } \\mathbf{u})$ can be likewise handled.\n- Directly handling the non-linear term $\\mathbf{u}\\cdot \\nabla \\mathbf{u}$ may be difficult. However, this term is associated with Lie derivative $\\mathcal L_{X}$, which is more profound in differential geometry and allows generalization to non-Euclidean spaces. There already exist some preliminary implementations of Lie derivative and vector fields onward. Nevertheless, how to represent them effectively for neural networks remains open and is worthwhile to explore in the future. \n\n> **Q5**. On some minor issues. \n> (i) Simplex complex (l.63 p. 2) could be called simplicial complex instead. \n> (ii) Since x is not used in the PDE, $x \\in \\Omega$ in Equation 6 is confusing. Either using x in the PDE or writing in $\\Omega$ would be preferable.\n\nThanks for your correction. We have modified both in our manuscript. For (ii), it is indeed a bit confusing though $\\mathbf{x}$ is implicitly used in the fields."}, {"review_id": "NJPzBSqFno", "reviewer": "Reviewer_1rUG", "summary": "This paper proposes a neural operator framework for solving boundary value problems (BVPs) governed by partial differential equations (PDEs). The key innovation lies in leveraging Discrete Exterior Calculus (DEC) and Finite Element Exterior Calculus (FEEC) to explicitly encode higher-order topological structures (nodes, edges, faces, cells) as k-simplex features in graph neural networks (HOGNNs). This allows for principled modeling of scalar and vector fields as differential forms, enabling a topologically-aware neural solver.\n\nThe approach includes physics-informed encoder-decoder mappings between vector fields and forms, integrated physical loss functions that correspond to conservation laws (e.g., Gauss’s law), and theoretical guarantees of universal approximation for certain BVP classes. There are limited experimental evaluations on 2D and 3D electrostatics and magnetostatics.", "questions": "- The formulation of the neural operator in Eq. (5) is too abstract and hard to interpret until the reader reaches the concrete instantiation in Eq. (7). The authors could either move the motivating example earlier or introduce a simplified case (e.g., scalar Poisson) first.\n\n- Figures 1 and 2 are critical but lack self-contained descriptions. For a technical audience unfamiliar with DEC or the specific PDE being solved, these need extensive legends explaining what each symbol and arrow denotes.\n\n- More focussed experiments: Pick 2-3 shapes in 2D and 3D (exhibiting different topological characteristics), and comprehensively evaluate performance against other state of art. Additionally, while the paper briefly notes that performance is affected by mesh shape and quality, this is an important limitation for deployment in practical engineering settings. Additional empirical analysis (e.g., mesh convergence study) would strengthen the impact.", "limitations": "Yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 4, "originality": 4, "strengths_and_weaknesses": "*Strengths*\n\n- The methodology is technically sound and grounded in well-established mathematical frameworks (DEC and FEEC), with rigorous theoretical support (e.g., encoder-decoder preservation, universal approximation).\n- This work is an important step toward building structure-preserving, physically grounded neural solvers for BVPs, providing a new approach. The proposal to use higher-order GNNs aligned with exterior calculus opens a new direction for learning-based PDE solvers in domains like electromagnetism, fluid dynamics, and elasticity.\n- The synthesis of DEC/FEEC with HOGNNs for BVPs appears novel and nontrivial. While GNNs and physics-informed learning are well-trodden, the topological interpretation via differential forms and the dual-manifold message passing are particularly unique contributions (for solving PDEs).\n- The physics-informed loss design is very elegant and potentially efficient! I love how the work leverages integrated features (which appear naturally in DEC) directly, which avoids the need for expensive quadrature.\n\n*Weakness*\n- This is a very dense paper, with some parts written with excessive jargon (particularly early on), making it hard for a broader audience to parse. Key ideas (like Eq. 5) are underexplained at first reading. \n- Figures (e.g., Figure 1 and 2) are central to the understanding but require much more extensive captions and annotations to be fully self-contained.\n- It is unclear why the 3D magnetostatic results are not competitive (relative to other baselines) even though the model is designed to generalize naturally to 3D.\n-  The comparison misses modern BVP-specific baselines such as FNO, CNO, SCOT, etc., which are more appropriate than DeepONet or Transolver (designed for time-evolving systems).\n- The experimental results are very perfunctory, and should be strengthened. \n- While the stated focus is on time-dependant PDEs (in the abstract), all analysis and results are for the time-**independant** PDE. This is rather confusing.", "was_revised": true, "final_justification": "Thank you to the authors for actively engaging during this process. \nI also appreciate the additional results that the authors added. \nI am happy to increase my score to 5", "rebuttal": "We are sincerely grateful for the time and effort you have dedicated to reviewing our manuscript. Below, we address each of your comments in detail. \n\n> **Q1**. While the stated focus is on **time-dependent PDEs (in the abstract)**, all analysis and results are for the time-independent PDE. This is rather confusing.\n\nWe actually describe **in the abstract** that the **time-independent** boundary value problems (BVPs) in electromagnetism are instantiated to illustrate the proposed framework, not time-dependent.\n\nThough this framework can be easily generalized to time-evolving PDEs, we evaluate our method on time-independent cases, aiming to verify its capability of capturing local behaviors by introducing higher-order elements.\n\n> **Q2**. More focused experiments: Pick 2-3 shapes in 2D and 3D (exhibiting different topological characteristics), and comprehensively evaluate performance against other state of art. \n\nDue to character limit, only the results on 2D magnetotastics are provided. 8 relative benchmarks (from A1 to C2) are used for evaluation, each having different topological properties on connected component amount and hole amount:\n\n|Benchmark Id|A1|A2|A3|B1|B2|B3|C1|C2|\n|-|-|-|-|-|-|-|-|-|\n|Connected Component Amount|1|1|1|2|2|2|4|4|\n|Hole Amount|2|4|8|2|4|8|2|4|\n\nPerformance on these benchmarks are compared with other baselines. The table below shows the test loss of each model. It turns out that the advantage of our approach persists as the underlying topology changes. \n\n| Model/Benchmark Id | A1                  | A2                  | A3                  | B1                  | B2                  | B3                  | C1                  | C2                  |\n| ------------------ | ------------------- | ------------------- | ------------------- | ------------------- | ------------------- | ------------------- | ------------------- | ------------------- |\n|FourierType|3.969|3.670|3.532|4.963|4.639|5.561|3.457|4.211|\n|GalerkinType|2.457|2.883|3.207|**1.907**|2.511|2.912|**2.093**|3.275|\n|MKGN|9.177|5.208|6.769|8.428|5.813|4.689|4.259|6.722|\n|Transolver|2.733|$\\underline{2.479}$|$\\underline{2.167}$|2.301|$\\underline{2.273}$|$\\underline{2.112}$|2.234|**1.901**|\n|Geo-FNO|9.423|9.448|9.453|9.448|9.450|9.451|9.448|9.451|\n|GNOT|$\\underline{2.406}$|3.277|3.009|2.799|3.006|2.588|3.373|2.213|\n|DEC-HOGNN(Ours)|**1.573**|**1.487**|**1.498**|$\\underline{2.088}$|**1.901**|**1.727**|$\\underline{2.181}$|$\\underline{2.176}$|\n\n**Rmk**. Baselines achieving top 1 and top 2 performance are marked via bolding and underscoring, respectively. \n\n> **Q3**. Additionally, while the paper briefly notes that performance is affected by mesh shape and quality, this is an important limitation for deployment in practical engineering settings. Additional empirical analysis (e.g., mesh convergence study) would strengthen the impact.\n\nThanks for your suggestion. We then analyze the negative impact out of mesh degradation by randomly dropping a certain amount of edges on a mesh hierarchically. Several extra edges are further removed to fix the topology. Due to character limit, only partial data can be demonstrated. \n\nTo measure the mesh quality quantitatively, three indicators are adopted where the arrows imply the direction in which the mesh degrades. It is reflected that dropping edges from a triangularized mesh usually comes with mesh degeneration. \n\nThen, we train models on the meshes suffering from different levels of degeneration. The average mesh quality and validation loss within 500 epochs are reported. It is found that minor degeneration would not affect the model performance violently while a major one leads to salient performance drop. Note that it is also infeasible to adopt classical solvers using these meshes with prominent quality issues. Thus these negative effects are tolerable. \n\n| Edge Drop | #E | #V| #F | AR$↑$ | EAS $↑$ | ATR$↑$ | 0     | 100  | 200  | 300  | 400  | 500  |\n|-|-|-|-|-|-|-|-|-|-|-|-|-|\n|0|487|180|306|1.286|0.595|0.197|9.75|4.60|2.95|2.14|1.85|1.44|\n|25|462|180|281|1.301|0.594|0.251|9.32|4.24|2.89|2.18|1.65|1.42|\n|50|437|180|256|1.318|0.607|0.316|9.32|4.44|3.26|2.42|2.01|1.51|\n|100|387|180|206|1.365|0.670|0.404|9.33|4.58|3.52|2.66|2.32|2.04|\n|150|334|177|156|1.401|0.783|0.487|9.41|4.90|3.76|2.99|2.51|2.11|\n|200|275|168|106|1.458|0.811|0.584|9.59|5.22|3.89|3.35|3.04|2.70|\n|300|139|114|24|1.929|0.816|0.860|10.67|6.46|5.78|5.27|4.93|4.82|\n\n**Rmk.**  #E, #V, #F are the amount of edges, vertices and faces, respectively. Usually, a mesh is considered having good quality if its elements are uniform and regular. We briefly introduce these three indicators:\n- **Aspect Ratio (AR)**. There are many versions of AR. We define it as the ratio between the maximum edge length over the minimum.\n- **Equi-Angle Skewness (EAS)**. An N-polygon's ideal angle is $\\theta:=\\pi (N-2)/N$. For an actual angle $\\theta_{i}$, its offset is $\\delta_{i}:= |\\theta_{i}-\\theta|/\\max\\{\\theta,\\pi-\\theta\\}$. And its EAS is $\\max_{i}\\delta_{i}$. \n- **Area Transition Ratio (ATR)**. It depicts the uniformity of element size. Let $A_i,A_j$ be the area of two neighboring elements, then its ATR is $|A_i-A_j|/\\max\\{A_i,A_j\\}$. \n\n> **Q4.** It is unclear why the 3D magnetostatic results are not competitive (relative to other baselines) even though the model is designed to generalize naturally to 3D.\n\nThis seemingly abnormal result mainly comes with different complexities of the datasets. Take Poisson equation ($\\Delta u(\\mathbf{x})=f(\\mathbf{x}),\\mathbf{x}\\in\\Omega;u(\\mathbf{x})=g(\\mathbf{x}),\\mathbf{x}\\in \\partial\\Omega$) as an example, this BVP is dominated by both source term $f(\\mathbf{x})$ and boundary condition $g(\\mathbf{x})$. As stated in Appendix G, the behavior of the source term $f(\\mathbf{x})$ in 2D cases is more complicated in the sense that every triangle is endowed with random charges. While in 3D cases, only several embedded cubes with charge or circuits are involved. And as we know, such higher-order structure only captures local features since these higher-order features are originated by differential operators, which only reflect local behaviors. Consequently, the more the BVP is dominated by local restriction $f(\\mathbf{x})$ instead of global restriction $g(\\mathbf{x})$, the better performance our proposed method would achieve. In this perspective, the result that 2D cases appear better than the 3D counterparts, is consistent with both the data and the architecture. \n\nBesides, our method focuses more on how to capture local behaviors governed by $f(\\mathbf{x})$ rather than the global by $g(\\mathbf{x})$, due to the nature of integrated form representation. We will leave how to further combine global behaviors reflected by boundary conditions with integrated form representations for future work. \n\n> **Q5.** The comparison misses modern BVP-specific baselines such as FNO, CNO, SCOT, etc., which are more appropriate than DeepONet or Transolver (designed for time-evolving systems).\n\nThanks for your suggestion. However, these baselines **cannot directly fit our experiment settings due to their requirement of a regular mesh**. Both FNO and CNO requires the spatial domain or its parametric space to be a regular grid (e.g., *Elliptic mesh for the airfoil problem*(p.55 fig.23)[1], Appendix 3[2]). This fact is also supported by *\"state-of-the-art neural operators such as FNO, CNO, and scOT which are tailored for Cartesian grids\"* [3] (p.8). \n\nThus these baselines are improper to handle irregular meshes as shown in fig.5, which are more common in practical use. Also, it is worthwhile mentioning that involved baselines like *GNOT* also outperform earlier baselines like interpolated FNO, Geo-FNO in regular-mesh settings (Sec.4 Tab.1)[4]. \n\nStill, we can understand your concerns on insufficient experiments. Therefore, we further evaluate the performance against *Geo-FNO* and *FourierType Transformer Operator*, which can handle arbitrary meshes. We will clarify why baselines like *FNO, CNO, and scOT* are not adopted in our manuscript and add extra comparisons in the updated version.\n\n| Model| 2D Eletrostatics | 2D Magnetostatics |\n|-|-|-|\n|Geo-FNO|2.017|3.289|\n|FourierType|1.553|1.986|\n|GNOT|2.064|1.142|\n|DEC-HOGNN(Ours)|**0.623**|**0.875**|\n\nRef.\n\n[1] (CNO) Convolutional Neural Operators for robust and accurate learning of PDEs. arxiv.\n\n[2] (FNO) Fourier Neural Operator for Parametric PDE. \n\n[3] RIGNO: A Graph-based framework for robust and accurate operator learning for PDEs on arbitrary domains. arxiv. \n\n[4] GNOT:A General Neural Operator Transformer for Operator Learning. \n\n> **Q6.** *Suggestions on writing*. \n> (i) This is a very dense paper, with some parts written with excessive jargon (particularly early on), making it hard for a broader audience to parse. Key ideas (like Eq. 5) are underexplained at first reading.\n> (ii) Figures (e.g., Figure 1 and 2) are central to the understanding but require much more extensive captions and annotations to be fully self-contained.\n> (iii) The formulation of the neural operator in Eq. (5) is too abstract and hard to interpret until the reader reaches the concrete instantiation in Eq. (7). The authors could either move the motivating example earlier or introduce a simplified case (e.g., scalar Poisson) first.\n\nThanks for your advice. For (i) and (ii), we have added more detailed explanations on these terms, equations and figures, making it more friendly to broader readers. For (iii), we have exchanged the abstract presentation of the neural operator and the concrete example properly."}, {"review_id": "NsDFQ55iR0", "reviewer": "Reviewer_dQYA", "summary": "This paper proposed a novel GNN neural operator (NO) that is topology-aware by leveraging discrete exterior calculus (DEC). The advantage of this novel NO includes: 1. preserving the topological structure of the manifold which is also the problem domain of PDE; 2. It can turn some equations involving divergence and curl into boundary conditions by Stokes' theorem, such that the physical loss from equation no longer needs taking derivatives, e.g., divergence and curl, like PINN does.\n\nIn order to leverage the DEC, the authors designed a GNN that handles 4 types of neighborhood (see eq (8)-(12)). In experiments of 2D and 3D problems, the proposed method outperforms all baselines and showing strong advantage.", "questions": "N.A.", "limitations": "N.A.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 3, "originality": 4, "strengths_and_weaknesses": "**Strength**\n\nThis paper is highly innovative and origin by introducing DEC to learning operators with GNN. As mentioned in summary, there are two advantages: 1. Preserving topology structure of the manifold; 2. Simplifying the equation constraints to boundary constraints (see eq (19)-(22)). In prior works, topology structure is either ignored (e.g., transformer-based NO and FNO) or partially encoded (with graph adjacency of GNN). In contrast, DEC utilizes 4 types of neighborhood, which is a complete methodology. Ablation study (Table 2) shows that the complete setting outperforms incomplete ones. The second advantage naturally enables the proposed method to incorporate physics-informed loss without taking derivatives (for some equations, e.g., Maxwell, Navier-Stokes, etc.).\n\nThe paper is well written and elegantly presented. Also, for self-completeness, the paper provides universal approximation analysis.\n\n**Weakness**\n\n1. In Table 3, I think the bold numbers are the best of each column. Please clarify accordingly to avoid confusion.\n2. The current data is generated with finite element method, which may not be compatible with DEC (as the discussion in line 104-105). If data generation by DEC is not so convenient, this may limit the application of the method to some extent. \n3. Typo? On line 476, $H\\in\\Omega^1(\\star M)$ -> $ H\\in\\Omega^1(M)$ to align with Fig. 2.", "was_revised": true, "final_justification": "After reading other reviews and rebuttals, I see that all reviews are positive and concerns are addressed. Therefore, my score remains.", "rebuttal": "We sincerely appreciate your constructive feedback and meticulous evaluation of our work. Below, we provide responses to each point raised.\n\n> **Q1**. The current data is generated with finite element method, which may not be compatible with DEC (as the discussion in line 104-105). If data generation by DEC is not so convenient, this may limit the application of the method to some extent.\n\nThanks for your comments. Reviewer RTKj has similar concerns. We can fully understand your concerns on incompatibility between data and model design. However, FEM-based solvers are significantly more mature and widely adopted than their DEC and FEEC counterparts. Theoretically, all of these solvers can **converge to the same solution**, provided that the PDE is well-posed. As such, FEM, FEEC, and DEC are all valid choices for data generation when no numerical issues arise. This is why data generation using DEC is not strictly necessary.\n\n> **Q2**. *On minor issues*.\n> (i) Typo? On line 476,  $H\\in \\Omega^1(\\star M)\\to H\\in\\Omega^1(M)$  to align with Fig. 2.\n> (ii) In Table 3, I think the bold numbers are the best of each column. Please clarify accordingly to avoid confusion.\n\nThanks for your careful reading and corrections. For (i),  $H$ should be indeed defined on $M$ since $B$ is on $\\star M$; for (ii), we will clarify its meaning in the caption."}, {"review_id": "glsk3H31qh", "reviewer": "Reviewer_RTKj", "summary": "This paper introduces a novel method called DEC-HOGNN in the category of neural operators. It makes use of higher-order graph neural networks that can handle interactions between simplices of different orders, not just same orders. It also designs an encoder-decoder structure for differential forms of different orders involved to make predictions. By writing down fields as differential forms, the physical losses can also be easily implemented. This approach is tested on boundary value problems in electrostatics and magnetostatics for both 2D and 3D.", "questions": "1. How are the boundary conditions $\\nabla \\phi = u(\\bf{x}), \\bf{x} \\in \\partial \\Omega$ implemented as an input for this neural operator?\n2. Related to point 1 in weaknesses, why are the data still generated with traditional FEM, instead of FEEC? If data are generated with FEEC, then no encoding is required and all differential forms are naturally solved out.", "limitations": "Yes, the limitations are discussed in this paper.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 4, "originality": 4, "strengths_and_weaknesses": "Strengths:\n1. This is an interesting attempt to implement physical losses elegantly and effectively in neural operators, which an inspiring and meaningful first exploration for AI for PDE community. This method does not require differentiation for physical losses.\n2. The mathematical motivations of this work are demonstrated in an easy to follow manner.\n3. There are clear proofs for important theorems mentioned in this paper, including one universal approximation property for this proposed neural operator.\n\nWeaknesses:\n1. Though a large part of this paper is devoted to FEEC, the data generation is still achieved with traditional FEM methods, but not FEEC. This point makes this work not in a \"pure\" FEEC manner.\n2. The datasets explored in the experiments of this work do not cover a wide enough range of problems. Only electrostatics and magnetostatics are covered, which does not provide enough evidence for the superiority of this approach.\n3. The inference speed is not discussed and compared with other methods.", "was_revised": true, "final_justification": "All of my concerns have been resolved from this rebuttal period. This paper is of good quality, so I do not have further concerns. Therefore, I recommend acceptance for this paper.", "rebuttal": "We sincerely appreciate your insightful feedback and valuable assessment of our work. Below we address each comment in detail.\n\n> **Q1**. The inference speed is not discussed and compared with other methods.\n\nThanks for your comment. We have updated the manuscript to compare our model with baselines on computation cost. Below we provide various relative indices on the 2D electrostatics benchmark and 3D electrostatics benchmark respectively. \n\nSince our experiments mainly focused on the effectiveness of the model, there is still much room to reduce the overheads of the model. Also, higher-order elements come with extra expenditures. The adjacency amount soars up when lots of tetrahedrals appear in a finer tetrahedralization, among which some might be unnecessary. An interesting direction for future research is how to fully leverage higher-order topological adjacencies without being hindered by the proliferation of small tetrahedrals. Nevertheless, our results show that the inference efficiency of the proposed model is acceptable and remains significantly faster than classical solvers, which typically require several seconds per inference.\n\nTab 1: Comparison on the 2D electrostatics benchmark. \n\n| Model            | FLOP(M) | Memory(MB) | Inference Time(ms) |\n| ---------------- | ------- | ---------- | ------------------ |\n| DeepONet         | 0.10    | 2.74       | 0.33               |\n| MKGN             | 0.22    | 19.80      | 3.73               |\n| Galerkin-Type    | 2.96    | 18.20      | 1.67               |\n| GNOT             | 1.99    | 25.30      | 1.55               |\n| Transolver       | 1.34    | 19.05      | 1.84               |\n| GAT-based        | 0.10    | 7.69       | 1.48               |\n| Graph UNet-based | 0.14    | 3.36       | 4.69               |\n| GT-based         | 10.35   | 53.96      | 1.54               |\n| DEC-HOGNN(Ours)  | 2.85    | 55.59      | 5.00               |\n\nTab 2: Comparison on the 3D electrostatics benchmark. \n\n| Model            | FLOP(M) | Memory(MB) | Inference Time(ms) |\n| ---------------- | ------- | ---------- | ------------------ |\n| DeepONet         | 0.10    | 93.41      | 0.28               |\n| MKGN             | 0.53    | 957.25     | 5.50               |\n| Galerkin-Type    | 0.95    | 1122.03    | 5.29               |\n| GNOT             | 1.99    | 281.99     | 2.13               |\n| Transolver       | 1.35    | 235.51     | 2.46               |\n| GAT-based        | 0.11    | 264.42     | 2.08               |\n| Graph UNet-based | 0.78    | 257.22     | 6.76               |\n| GT-based         | 10.35   | 367.42     | 3.41               |\n| DEC-HOGNN(Ours)  | 5.37    | 9380.33    | 94.09              |\n\n> **Q2**.  How are the boundary conditions $\\nabla \\phi=u(x)$ implemented as an input for this neural operator?\n\nUsually a scalar potential $\\phi$ is unavailable in many physics scenarios. In electromagnetism, $\\nabla \\phi$ may be represented by electric intensity $\\mathbf{E}$ or magnetic intensity $\\mathbf{H}$, which can be easily measured by sensors in practice. Therefore, in our approach, we take these vector fields as both input and final output. In fact, we introduce $\\phi$ only to show the theorems.\n\nBut still, given a scalar field $\\phi$, we can leverage the gradient operator defined in graph calculus [1] or simply learn one like NIsoGCN [2] to yield a vector field $\\nabla \\phi$. It can then be interpreted as a vector proxy of a 1-form and thereby represented as integrated forms on boundary edges. \n\nRef.\n\n[1] Calculus on Graphs. Joel Friedman et al. \n\n[2] Physics-embedded neural networks: Graph neural pde solvers with mixed boundary conditions. \n\n> **Q3**. Though a large part of this paper is devoted to FEEC, the data generation is still achieved with traditional FEM methods, but not FEEC. This point makes this work not in a \"pure\" FEEC manner. Related to point 1 in weaknesses, why are the data still generated with traditional FEM, instead of FEEC? If data are generated with FEEC, then no encoding is required and all differential forms are naturally solved out.\n\nReviewer dQYA has similar concerns. We can fully understand your concerns on the compatibility between traditional FEM and its FEEC and DEC counterparts. For data generation, any of these methods is acceptable, as they yield the same solution for a well-posed PDE. Therefore, the method used for data generation does not compromise the purity or integrity of this work.\n\nRegarding your second point, if input data is presented in a FEEC manner, then we can indeed cast away the encoding process. However, such an assumption can be too ideal to fit more applications. Because it is easy to deploy sensors to sample vector intensity like $\\mathbf{B},\\mathbf{E}$ while it is hard to measure a differential form down there. Hence, the encoding process is still necessary for a broader applicability. \n\n> **Q4**. The datasets explored in the experiments of this work do not cover a wide enough range of problems. Only electrostatics and magnetostatics are covered, which does not provide enough evidence for the superiority of this approach.\n\nWe use Maxwell's equations as a case study, as they encompass a broad range of boundary value problems (BVPs) and exhibit a well-structured formulation within the framework of exterior calculus. Specifically, they involve differential operators such as $\\text{div}$ and $\\text{curl}$, which also appear in BVPs like *Darcy Flow*, as well as in time-dependent PDEs such as the *Heat and Wave equations*. Upon closer examination, these equations share structural similarities with Maxwell's equations, suggesting they follow the same underlying patterns. Consequently, experiments on electrostatics and magnetostatics may provide evidence for the framework’s broader applicability to a variety of linear PDEs. We will include them in the next version of our paper.\n\nWe also mention that this framework can also potentially be applied to non-linear PDEs like *Navier-Stokes Equations (NSE)*. This is mainly because $\\text{div},\\text{curl}$ can be encoded as higher-order element features. We briefly sketch a possible approach and kindly invite you to the response to **Reviewer 5x2r's Q4**."}]

Author rebuttal:

Meta-review: The paper presents a GNN-framework for solving partial differential equations that incorporates higher-order interactions based on discrete and finite element exterior calculus. The method outperforms existing neural operators and  provies theoretical guarantees of universal approximation for certain voundary value problem classes. Reviewers noted that the methodology is technically sound, highly innovative, and appreciated the rigorous theoretical support. They also commented that the work opens up a new direction for learning-based PDE solvers, that the physics-informed loss is elegant, and they appreciated the improved accuracy compared to baselines. On the other hand, reviewers raised concerns about clarifying the relationship to existing methods based on exterior calculus or simplicial complexes, the clarity of the writing, the relatively limited number of experimental results, and the omission of certain baselines (e.g., FNO, CNO, SCOT).

The authors responded to the reviewer concerns with a rebuttal, which provided additional experimental evaluation, clarified why certain baselines were not included (e.g., due to their requirement for a regular grid), added an assessment of computational costs, and replied to a number of other technical questions. The reviewers agreed that all or most of their concerns were addressed by the rebuttal, and there was consensus to accept the paper. 

I recommend the paper be accepted, and I also note that it stands out from other papers because of the thoroughness of the work, the technical novelty of the approach, and the potential for broad impact in the area of neural operators and physics-informed neural networks.
Decision: Accept (spotlight)