math.NAOct 12, 2026

A Method for Exact Enforcement of Boundary Conditions on General Curved Triangular Domains: Algorithm and Application to Scientific Machine Learning

Authors: Suchuan Dong, Yuchuan Zhang

Organizations: Center for Computational and Applied Mathematics Department of Mathematics Purdue University, USA

Read PDF
Ask AI
Cite

Copy citation

DOINot available for this paper

Abstract

We present a method for enforcing exactly Dirichlet, Neumann, and Robin conditions on general triangular domains with arbitrarily curved boundaries. The algorithmic constructions developed here complements those in a previous work for general quadrilateral domains. Our method combines (i) the exact mapping between a general curved triangular domain and the standard triangle, (ii) TFC constrained expressions, and (iii) Barnhill transfinite interpolations, to provide explicit closed-form trial ansatze that mathematically satisfy the imposed boundary conditions exactly. Corner compatibility induced by the boundary conditions is pivotal, and enforcing exactly the imposed boundary conditions hinges on the exact enforcement of the induced compatibility constraints. The corner compatibility associated with the acute vertices of the standard triangle and the coupling between gradient components when the hypotenuse corresponds to a Neumann/Robin boundary pose particular challenges. Two strategies are critical to resolving these issues, by reformulating conditional corner constraints into equivalent unconditional forms and by introducing a separate TFC-style auxiliary edge construction together with an associated linear solve on the hypotenuse. We analyze the corner compatibility constraints and present formulations for the general form of functions that exactly satisfy the imposed Dirichlet, Neumann, or Robin conditions on curved triangular domains. The proposed method has been implemented using extreme learning machine type randomized neural networks for solving PDE problems. We present ample numerical experiments with linear/nonlinear PDE problems on a variety of domains with complex boundary geometries. Numerical results show that the current method has enforced the Dirichlet, Neumann, and Robin conditions on curved domain boundaries to the machine precision.

Figures & tables

Figure 1: Sketch of (a) a general triangular domain Ω\Omega with curved boundaries, and (b) the standard triangular domain Ωst={(ξ,η) : −1⩽ξ,η⩽1, ξ+η⩽0}\Omega_{st}=\{(\xi,\eta)\ :\ -1\leqslant\xi,\eta\leqslant 1,\ \xi+\eta\leqslant 0\} .
 3 [g(ξ,η)−g(ξ,−1)]−[g(−1,η)−g(−1,−1)],\displaystyle\ \big[g(\xi,\eta)-g(\xi,-1)\big]-\big[g(-1,\eta)-g(-1,-1)\big], 
 H(ξ,−ξ)=\displaystyle H(\xi,-\xi)= [g(ξ,−ξ)−g(ξ,−1)]−[g(−1,−ξ)−g(−1,−1)],\displaystyle\ \big[g(\xi,-\xi)-g(\xi,-1)\big]-\big[g(-1,-\xi)-g(-1,-1)\big], 
 H(−η,η)=\displaystyle H(-\eta,\eta)= [g(−η,η)−g(−η,−1)]−[g(−1,η)−g(−1,−1)].\displaystyle\ \big[g(-\eta,\eta)-g(-\eta,-1)\big]-\big[g(-1,\eta)-g(-1,-1)\big]. 
Algorithm 1 Computing the NN solution (DBC on all boundaries)
Figure 2: Variable-coefficient Poisson equation: Domain geometries (top row) and distributions of the exact solution (bottom row) for five domains: (a,f) domain #1, (b,g) domain #2, (c,h) domain #3, (d,i) domain #4, (e,j) domain #5.
Figure 3: Variable-coefficient Poisson equation (DBC on all boundaries): Distributions of the NN solutions (top row) and their point-wise absolute errors (bottom row) on the five domains. Simulation parameters: (a,f) Rm=4.0R_{m}=4.0 , Q=70Q=70 , M=400M=400 , simple map; (b,g) Rm=3.5R_{m}=3.5 , Q=100Q=100 , M=500M=500 , simple map; (c,h) Rm=4.0R_{m}=4.0 , Q=100Q=100 , M=500M=500 , Winslow map ( pw=0p_{w}=0 , Qw=200Q_{w}=200 ); (d,i) Rm=4.5R_{m}=4.5 , Q=100Q=100 , M=700M=700 , Winslow map ( pw=2p_{w}=2 , Qw=200Q_{w}=200 ); (e,j) Rm=4.4R_{m}=4.4 , Q=120Q=120 , M=900M=900 , Winslow map ( pw=5p_{w}=5 , Qw=200Q_{w}=200 ).
 Domain #1Domain #2Domain #3Domain #4Domain #5
max uu -error (domain)2.027×10−72.027\times 10^{-7}1.989×10−61.989\times 10^{-6}3.433×10−73.433\times 10^{-7}2.185×10−52.185\times 10^{-5}2.008×10−52.008\times 10^{-5}
rms uu -error (domain)2.184×10−82.184\times 10^{-8}3.830×10−73.830\times 10^{-7}5.932×10−85.932\times 10^{-8}2.852×10−62.852\times 10^{-6}5.593×10−65.593\times 10^{-6}
max DBC-error ( AB‾\overline{AB} )2.776×10−162.776\times 10^{-16}0.00.02.220×10−162.220\times 10^{-16}2.220×10−162.220\times 10^{-16}6.661×10−166.661\times 10^{-16}
rms DBC-error ( AB‾\overline{AB} )1.687×10−161.687\times 10^{-16}0.00.07.697×10−177.697\times 10^{-17}3.445×10−173.445\times 10^{-17}1.750×10−161.750\times 10^{-16}
max DBC-error ( BC‾\overline{BC} )8.882×10−168.882\times 10^{-16}4.441×10−164.441\times 10^{-16}6.661×10−166.661\times 10^{-16}6.661×10−166.661\times 10^{-16}4.441×10−164.441\times 10^{-16}
rms DBC-error ( BC‾\overline{BC} )2.351×10−162.351\times 10^{-16}1.243×10−161.243\times 10^{-16}1.519×10−161.519\times 10^{-16}1.057×10−161.057\times 10^{-16}1.472×10−161.472\times 10^{-16}
Table 1: Variable-coefficient Poisson equation (DBC on all boundaries): Maximum and rms errors of the NN solution over the domain and the maximum and rms DBC errors on the three boundaries. Simulation parameters follow those of Figure 3 .
Figure 4: Variable-coefficient Poisson equation on domain #1 (NBCs and RBCs): Distributions of the NN solutions (top row) and their point-wise absolute errors (bottom row) under four boundary configurations. (a,e) Config #1: NBC on AC‾\overline{AC} and DBC on the other boundaries; (b,f) Config #2: NBC on AB‾\overline{AB} and AC‾\overline{AC} , and DBC on BC‾\overline{BC} ; (c,g) Config #3: RBC on AC‾\overline{AC} with αAC=1.0\alpha_{AC}=1.0 , and DBC on the other boundaries; (d,h) Config #4: RBC on AB‾\overline{AB} and AC‾\overline{AC} with αAB=0.5\alpha_{AB}=0.5 and αAC=1.0\alpha_{AC}=1.0 , and DBC on BC‾\overline{BC} . Simulation parameters: (a,e) Rm=3.3R_{m}=3.3 , (b,f) Rm=3.3R_{m}=3.3 , (c,g) Rm=3.4R_{m}=3.4 , (d,h) Rm=4.0R_{m}=4.0 . Q=100Q=100 and M=500M=500 , together with the simple domain map, are used in all four configurations.
 Config #1Config #2Config #3Config #4
max uu -error (domain)9.113×10−89.113\times 10^{-8}8.459×10−78.459\times 10^{-7}8.596×10−88.596\times 10^{-8}4.430×10−74.430\times 10^{-7}
rms uu -error (domain)1.462×10−81.462\times 10^{-8}1.873×10−71.873\times 10^{-7}1.954×10−81.954\times 10^{-8}2.128×10−72.128\times 10^{-7}
max bc-error ( AB‾\overline{AB} )8.882×10−168.882\times 10^{-16}4.441×10−164.441\times 10^{-16}8.882×10−168.882\times 10^{-16}3.747×10−163.747\times 10^{-16}
rms bc-error ( AB‾\overline{AB} )3.237×10−163.237\times 10^{-16}1.047×10−161.047\times 10^{-16}3.237×10−163.237\times 10^{-16}1.476×10−161.476\times 10^{-16}
max bc-error ( BC‾\overline{BC} )0.00.00.00.00.00.00.00.0
rms bc-error ( BC‾\overline{BC} )0.00.00.00.00.00.00.00.0
Table 2: Variable-coefficient Poisson equation on domain #1 (NBCs and RBCs): Maximum/rms NN domain errors and the maximum/rms boundary-condition errors (DBC, NBC, or RBC) on the three boundaries. Different configurations correspond to those in Figure 4 .
Figure 5: Variable-coefficient Poisson equation on domain #3 (NBCs and RBCs): Distributions of the NN solutions (top row) and their point-wise absolute errors (bottom row) under four BC configurations. (a,e) Config #1: NBC on AC‾\overline{AC} , and DBC on the other boundaries; (b,f) Config #2: NBC on AB‾\overline{AB} and AC‾\overline{AC} , and DBC on BC‾\overline{BC} ; (c,g) Config #3: RBC on AC‾\overline{AC} with αAC=1.0\alpha_{AC}=1.0 , and DBC on the other boundaries; (d,h) Config #4: RBC on AB‾\overline{AB} and AC‾\overline{AC} with αAB=0.5\alpha_{AB}=0.5 and αAC=1.0\alpha_{AC}=1.0 , and DBC on BC‾\overline{BC} . Simulation parameters: (a,e) Rm=3.9R_{m}=3.9 , (b,f) Rm=3.3R_{m}=3.3 , (c,g) Rm=3.9R_{m}=3.9 , (d,h) Rm=3.2R_{m}=3.2 . Q=100Q=100 and M=500M=500 , together with the Winslow map ( pw=0p_{w}=0 , Qw=50Q_{w}=50 ), are used in all four configurations.
 Config #1Config #2Config #3Config #4
max uu -error (domain)6.542×10−66.542\times 10^{-6}2.139×10−52.139\times 10^{-5}6.490×10−66.490\times 10^{-6}1.169×10−51.169\times 10^{-5}
rms uu -error (domain)1.286×10−61.286\times 10^{-6}1.317×10−61.317\times 10^{-6}1.196×10−61.196\times 10^{-6}4.804×10−64.804\times 10^{-6}
max bc-error ( AB‾\overline{AB} )4.441×10−164.441\times 10^{-16}4.441×10−164.441\times 10^{-16}4.441×10−164.441\times 10^{-16}8.882×10−168.882\times 10^{-16}
rms bc-error ( AB‾\overline{AB} )3.444×10−173.444\times 10^{-17}1.388×10−161.388\times 10^{-16}3.444×10−173.444\times 10^{-17}2.126×10−162.126\times 10^{-16}
max bc-error ( BC‾\overline{BC} )0.00.00.00.00.00.00.00.0
rms bc-error ( BC‾\overline{BC} )0.00.00.00.00.00.00.00.0
Table 3: Variable-coefficient Poisson equation on domain #3 (NBCs and RBCs): Maximum/rms NN domain errors and maximum/rms boundary-condition errors (DBC, NBC, or RBC) on the three boundaries. Different configurations correspond to those in Figure 5 .
Figure 6: Nonlinear Helmholtz equation: Domain geometries (top row) and distributions of the exact solution (bottom row) for five domains: (a,f) domain #1, (b,g) domain #2, (c,h) domain #3, (d,i) domain #4, (e,j) domain #5.
Figure 7: Nonlinear Helmholtz equation (DBC on all boundaries): Distributions of the NN solutions (top row) and their point-wise absolute errors (bottom row) on the five domains. Simulation parameters: (a,f) Rm=3.9R_{m}=3.9 , Q=100Q=100 , M=500M=500 , Winslow map ( pw=1p_{w}=1 , Qw=100Q_{w}=100 ); (b,g) Rm=4.0R_{m}=4.0 , Q=100Q=100 , M=500M=500 , Winslow map ( pw=1p_{w}=1 , Qw=100Q_{w}=100 ); (c,h) Rm=4.0R_{m}=4.0 , Q=80Q=80 , M=300M=300 , simple map; (d,i) Rm=4.8R_{m}=4.8 , Q=100Q=100 , M=500M=500 , simple map; (e,j) Rm=4.0R_{m}=4.0 , Q=100Q=100 , M=300M=300 , simple map.
 Domain #1Domain #2Domain #3Domain #4Domain #5
max uu -error (domain)1.238×10−61.238\times 10^{-6}2.242×10−82.242\times 10^{-8}2.521×10−72.521\times 10^{-7}1.887×10−51.887\times 10^{-5}4.867×10−84.867\times 10^{-8}
rms uu -error (domain)2.965×10−72.965\times 10^{-7}4.869×10−94.869\times 10^{-9}6.013×10−86.013\times 10^{-8}3.586×10−63.586\times 10^{-6}9.449×10−99.449\times 10^{-9}
max DBC-error ( AB‾\overline{AB} )4.441×10−164.441\times 10^{-16}4.441×10−164.441\times 10^{-16}2.220×10−162.220\times 10^{-16}4.441×10−164.441\times 10^{-16}2.220×10−162.220\times 10^{-16}
rms DBC-error ( AB‾\overline{AB} )1.340×10−161.340\times 10^{-16}1.592×10−161.592\times 10^{-16}6.438×10−176.438\times 10^{-17}1.689×10−161.689\times 10^{-16}9.541×10−179.541\times 10^{-17}
max DBC-error ( BC‾\overline{BC} )5.130×10−165.130\times 10^{-16}6.245×10−166.245\times 10^{-16}3.192×10−163.192\times 10^{-16}5.551×10−165.551\times 10^{-16}4.441×10−164.441\times 10^{-16}
rms DBC-error ( BC‾\overline{BC} )1.372×10−161.372\times 10^{-16}1.277×10−161.277\times 10^{-16}1.010×10−161.010\times 10^{-16}1.420×10−161.420\times 10^{-16}1.252×10−161.252\times 10^{-16}
Table 4: Nonlinear Helmholtz equation (DBC on all boundaries): Maximum/rms domain errors of the NN solution and the maximum/rms DBC errors on the three boundaries. Simulation parameters follow those of Figure 7 .
Figure 8: Nonlinear Helmholtz equation (NBCs/RBCs on boundaries): Distributions of the NN solutions (top row) and their point-wise absolute errors (bottom row) with domains #1 to #4 under four settings. (a,e) Case #1: NBC on AC‾\overline{AC} and DBC on the other boundaries with domain #1; (b,f) Case #2: NBC on AB‾\overline{AB} and AC‾\overline{AC} , and DBC on BC‾\overline{BC} , with domain #2; (c,g) Case #3: RBC on AC‾\overline{AC} with αAC=1.0\alpha_{AC}=1.0 , and DBC on the other boundaries, with domain #3; (d,h) Case #4: RBC on AB‾\overline{AB} and AC‾\overline{AC} with αAB=0.5\alpha_{AB}=0.5 and αAC=1.0\alpha_{AC}=1.0 , and DBC on BC‾\overline{BC} , with domain #4; Simulation parameters: (a,e) Rm=3.5R_{m}=3.5 , Q=100Q=100 , M=500M=500 , Winslow map ( pw=1p_{w}=1 , Qw=100Q_{w}=100 ); (b,f) Rm=3.7R_{m}=3.7 , Q=100Q=100 , M=500M=500 , Winslow map ( pw=1p_{w}=1 , Qw=100Q_{w}=100 ); (c,g) Rm=3.7R_{m}=3.7 , Q=100Q=100 , M=500M=500 , simple map; (d,h) Rm=4.2R_{m}=4.2 , Q=100Q=100 , M=500M=500 , simple map.
 Case #1Case #2Case #3Case #4
max uu -error (domain)1.039×10−51.039\times 10^{-5}1.524×10−71.524\times 10^{-7}2.472×10−92.472\times 10^{-9}2.517×10−42.517\times 10^{-4}
rms uu -error (domain)1.069×10−61.069\times 10^{-6}2.615×10−82.615\times 10^{-8}5.068×10−105.068\times 10^{-10}7.118×10−57.118\times 10^{-5}
max bc-error ( AB‾\overline{AB} )8.882×10−168.882\times 10^{-16}8.882×10−168.882\times 10^{-16}2.220×10−162.220\times 10^{-16}1.554×10−151.554\times 10^{-15}
rms bc-error ( AB‾\overline{AB} )3.361×10−163.361\times 10^{-16}2.688×10−162.688\times 10^{-16}8.345×10−178.345\times 10^{-17}3.589×10−163.589\times 10^{-16}
max bc-error ( BC‾\overline{BC} )0.00.00.00.00.00.00.00.0
rms bc-error ( BC‾\overline{BC} )0.00.00.00.00.00.00.00.0
Table 5: Nonlinear Helmholtz equation (NBCs/RBCs on boundaries): Maximum/rms domain errors of the NN solution and the maximum/rms boundary-condition errors (DBC, NBC, or RBC) on the three boundaries. Different cases correspond to those in Figure 8 .
Figure 9: Burgers equation (time-evolving spatial domains): Domain geometries (top row) and distributions of the exact solution (bottom row) in the space-time plane: (a,e) Config #1, DBCs on AB‾\overline{AB} and AC‾\overline{AC} , no boundary condition on BC‾\overline{BC} ; (b,f) Config #2, DBCs on AB‾\overline{AB} , BC‾\overline{BC} and AC‾\overline{AC} ; (c,g) Config #3, DBCs on AB‾\overline{AB} , BC‾\overline{BC} and AC‾\overline{AC} ; (d,h) Config #4, NBC on AB‾\overline{AB} , DBC on AC‾\overline{AC} , and no boundary condition on BC‾\overline{BC} .
Figure 10: Burgers equation: Distributions of the NN solutions (top row) and their point-wise absolute errors (bottom row) in the space-time plane for the four configurations. Simulation parameters: (a,e): Config #1, Rm=4.0R_{m}=4.0 , Q=100Q=100 , M=400M=400 , Winslow map ( pw=2p_{w}=2 , Qw=200Q_{w}=200 ); (b,f): Config #2, Rm=3.5R_{m}=3.5 , Q=100Q=100 , M=400M=400 , Winslow map ( pw=1p_{w}=1 , Qw=100Q_{w}=100 ); (c,g) Config #3, Rm=4.2R_{m}=4.2 , Q=100Q=100 , M=400M=400 , Winslow map ( pw=0p_{w}=0 , Qw=100Q_{w}=100 ); (d,h) Config #4, Rm=2.9R_{m}=2.9 , Q=100Q=100 , M=500M=500 , Winslow map ( pw=0p_{w}=0 , Qw=100Q_{w}=100 ); Different configurations correspond to those as defined in Figure 9 .
 Config #1Config #2Config #3Config #4
max uu -error (domain)6.205×10−76.205\times 10^{-7}3.975×10−53.975\times 10^{-5}4.210×10−74.210\times 10^{-7}3.006×10−53.006\times 10^{-5}
rms uu -error (domain)1.898×10−81.898\times 10^{-8}9.017×10−69.017\times 10^{-6}1.049×10−71.049\times 10^{-7}2.100×10−62.100\times 10^{-6}
max bc-error ( AB‾\overline{AB} )2.220×10−162.220\times 10^{-16}2.220×10−162.220\times 10^{-16}4.441×10−164.441\times 10^{-16}1.776×10−151.776\times 10^{-15}
rms bc-error ( AB‾\overline{AB} )6.101×10−176.101\times 10^{-17}6.809×10−176.809\times 10^{-17}9.322×10−179.322\times 10^{-17}3.823×10−163.823\times 10^{-16}
max bc-error ( BC‾\overline{BC} )–3.331×10−163.331\times 10^{-16}4.441×10−164.441\times 10^{-16}–
rms bc-error ( BC‾\overline{BC} )–8.534×10−178.534\times 10^{-17}1.046×10−161.046\times 10^{-16}–
Table 6: Burgers equation: Maximum/rms domain errors of the NN solution and the maximum/rms boundary-condition errors (DBC, NBC) for the four configurations. Different configurations correspond to those in Figures 9 and 10 . In Configs #1 and #4, note that no boundary condition is imposed on BC‾\overline{BC} .
More in this topic

PDE Solving

PDE solving develops numerical and learned methods to approximate solutions of partial differential equations across scientific applications.

See all 249 papers

Explore similar work

CardsList
Mar 23, 2026

A Novel Method for Enforcing Exactly Dirichlet, Neumann and Robin Conditions on Curved Domain Boundaries for Physics Informed Machine Learning

We present a systematic method for exactly enforcing Dirichlet, Neumann, and Robin type conditions on general quadrilateral domains with arbitrary curved boundaries. Our method is built upon exact mappings between general quadrilateral domains and the standard domain, and employs a combination of TFC (theory of functional connections) constrained expressions and transfinite interpolations. When Neumann or Robin boundaries are present, especially when two Neumann (or Robin) boundaries meet at a vertex, it is critical to enforce exactly the induced compatibility constraints at the intersection, in order to enforce exactly the imposed conditions on the joining boundaries. We analyze in detail and present constructions for handling the imposed boundary conditions and the induced compatibility constraints for two types of situations: (i) when Neumann (or Robin) boundary only intersects with Dirichlet boundaries, and (ii) when two Neumann (or Robin) boundaries intersect with each other. We describe a four-step procedure to systematically formulate the general form of functions that exactly satisfy the imposed Dirichlet, Neumann, or Robin conditions on general quadrilateral domains. The method developed herein has been implemented together with the extreme learning machine (ELM) technique we have developed recently for scientific machine learning. Ample numerical experiments are presented with several linear/nonlinear stationary/dynamic problems on a variety of two-dimensional domains with complex boundary geometries. Simulation results demonstrate that the proposed method has enforced the Dirichlet, Neumann, and Robin conditions on curved domain boundaries exactly, with the numerical boundary-condition errors at the machine accuracy.

Suchuan Dong, Yuchuan Zhang

Center for Computational and Applied Mathematics Department of Mathematics Purdue University, USA

PDE SolvingPhysics-Informed ML
AskPDF
Sep 8, 2026

Selective boundary condition reduction via learned error gating

Parametric PDEs can admit different boundary conditions with different accuracy and computational cost. We introduce a framework for learning when one reduced boundary condition can replace another: paired solutions train a neural network to estimate the resulting domain and boundary errors, and the simpler condition is used only when both predicted errors meet prescribed tolerances. We focus on singular limits in applications, in which a stiff Robin or nonlinear boundary law is replaced by its limiting Dirichlet form. We evaluate the method on a galvanic corrosion problem and other nonlinear stationary and evolution problems.

Daniel Fernández, Dominik Penk, Dominik Riedelbauch

Chair for Dynamics, Control, Machine Learning, and Numerics (Alexander von Humboldt Professorship), Department of Mathematics, Friedrich-Alexander-Universität Erlangen-Nürnberg, Cauerstraße, 91058 Erlangen, Germany. · Schaeffler Technologies AG & Co. KG, Industriestraße 1–3, Herzogenaurach, Germany.

PDE Surrogate ModelingMultifidelity ModelingPartial Differential Equations
AskPDF
Aug 27, 2026

Enforcing Dirichlet Boundary Conditions in Operator Learning

Operator learning in scientific machine learning is concerned with approximation of maps between infinite-dimensional function spaces; such maps frequently arise as the solution operators of partial differential equations (PDEs). Neural operators have demonstrated broad empirical success at approximating such maps from data. However, most existing neural operator architectures enforce boundary conditions indirectly through training from data even though the boundary condition is often known exactly. Furthermore, existing modifications and approaches that do enforce boundary conditions explicitly suffer from impractical restrictions, including boundary smoothness, uniform grids, and separable, box-like domains. In this work, we propose an architecture which, independently of training, satisfies homogeneous Dirichlet boundary conditions, whilst simultaneously retaining the expressivity of existing kernel-integral neural operator architectures. This is achieved by enforcing the property that the output of each layer is contained in the span of a subset of the homogeneous Dirichlet eigenfunctions of the Laplacian on the output domain. The method requires only that the output domain be bounded with Lipschitz boundary and places no restriction on the choice of discretization, making it applicable to arbitrary mesh data and general geometries. We prove universal approximation for the resulting architecture; furthermore the approach we adopt in the analysis proves universality for a broad class of kernel-integral neural operators thereby uniting existing theory for a variety of operator learning methods. We validate the proposed method on maps defined by the coefficient to solution map in 2D PDEs: Darcy flow on a square domain and the Helmholtz equation on a circular domain. Comparisons are made with alternative methods.

Andrew M. Stuart, Margaret Trautner

Caltech · ETH Zürich Department of Mathematics

PDE Operator LearningNeural OperatorsPDE SolvingNeural Network Approximation Theory
AskPDF