Highlights:

• Geometry-Preserving ITD Framework: Enables large-scale simulation on voxel-native datasets by coupling mesh-free geometric representation with C-HiDeNN-TD solvers.

• Exact Dirichlet Enforcement: Strong boundary condition satisfaction via body-fitted functions eliminates the need for variational penalties or complex interface quadrature.

• Optimized Numerical Performance: Optimal \(O(h^{p+1})\) convergence is established through a priori error analysis

1 Introduction↩︎

Body-fitted finite-element methods (FEM) deliver high-order accuracy on smooth Computer-Aided Design (CAD) geometries, achieving the optimal \(\mathcal{O}(h^{p+1})\) convergence rate in the \(L^2\) norm for polynomials of degree \(p\)[1]. This accuracy, however, is predicated on the availability of a clean, watertight, conforming mesh. In modern computational-mechanics workflows, generating such a mesh has become the dominant cost: industrial CAD assemblies routinely contain slivers, gaps, and overlapping features, and the associated geometric clean-up and mesh generation consume the majority of the analysis cycle rather than the solution itself [2]. The pursuit of analysis frameworks that relax the conforming-mesh requirement is therefore a central theme of contemporary numerical analysis.

In several application areas the very premise of a clean parametric CAD file no longer holds. In biomedical engineering and geophysics, geometries such as trabecular bone microstructures, soft-tissue organs, and porous rock formations are acquired by Computed Tomography (CT) or Magnetic Resonance Imaging (MRI), which produce volumetric density fields rather than boundary representations (B-reps) [3]. In additive manufacturing—particularly multi-material printing and topology optimization—the design space is intrinsically voxelated, with spatially varying material properties and lattice infills that resist a clean B-rep description[4]. Converting such volumetric or voxel-native data into smooth NURBS surfaces is both computationally expensive and information-lossy. These trends motivate analysis frameworks that operate directly on volumetric representations, bypassing surface reconstruction altogether.

Immersed (or embedded) methods sidestep body-fitted meshing by discretizing the governing partial differential equations on a regular Cartesian background mesh, on which the physical domain \(\Omega\) is defined implicitly—most often through a level-set or signed-distance function. This idea underlies the finite cell method and its higher-order and image-based extensions [3], [5], [6], among other fictitious-domain techniques. On the immersed Dirichlet boundary \(\partial\Omega\), the constraint \(u = g\) can be imposed either weakly, through Nitsche-type penalization of the boundary residual [7]–[10], or strongly, through trial-function multiplication. The latter is exemplified by the \(\phi\)-FEM approach [11], [12], in which the trial function is constructed to vanish on \(\partial\Omega\) by multiplication with a level-set function, so that the boundary condition is satisfied exactly by construction rather than enforced variationally.

Classical implementations of these strategies share a common difficulty: they require numerical integration over geometrically irregular cut subdomains, and often along the immersed interface, typically via cut-cell or shifted-boundary quadrature [13]. This integration complicates assembly and, crucially, breaks the tensor-product structure of the underlying Cartesian mesh—the very structure that makes regular grids attractive for large-scale computation in the first place.

Preserving the tensor-product structure of the background grid is especially valuable because it is the prerequisite for reduced-order tensor methods. Fields defined on a regular Cartesian grid admit accurate low-rank Tucker and Canonical-Polyadic (CP) representations, and the trial space itself can be built in a separable format such as Proper Generalized Decomposition (PGD) [14], [15] or the more recent Convolution Hierarchical Deep-learning Neural Network with tensor decomposition (C-HiDeNN-TD) [16]–[18], which together turn \(d\)-dimensional problems into a sequence of one-dimensional ones and break the curse of dimensionality. Such separable solvers, however, require every operator in the weak form to retain a tensor-product form—precisely what cut-cell quadrature destroys.

A natural route to reconcile the two is to represent the geometry by a smooth implicit field [19], [20] and to enforce the Dirichlet condition by multiplying the trial space with a body-fitted mask that vanishes on \(\partial\Omega\) and itself admits a low-rank tensor representation. All integrals then remain volume integrals over the regular grid and reduce to tensor-product Gauss quadrature, with no cut-cell bookkeeping. Yet, despite this apparent compatibility, a systematic integration of strong, body-fitted immersed enforcement with separable tensor-decomposition solvers has not, to our knowledge, been reported. Several practical issues remain open, including the robust evaluation of the mask and its derivatives from voxelated signed-distance data and the control of conditioning in the fictitious (void) region.

In this paper we propose the Immersed Tensor Decomposition (ITD) framework, which couples a body-fitted function with the C-HiDeNN-TD reduced-order solver [16]–[18] to enable large-scale simulation on non-Cartesian domains using regular voxel meshes. The specific contributions are as follows.

• An exact Dirichlet immersed C-HiDeNN-TD formulation. The trial function is multiplied by a body-fitted function \(\Phi\) (BFF) based on the signed distance function (SDF) that vanishes on \(\partial\Omega\), yielding strong, exact enforcement of the Dirichlet condition without any variational penalty or interface quadrature. We derive the corresponding transformed weak form, in which a key algebraic cancellation reduces the geometric coefficients that arise naively (five in 2D, six in 3D) to only three Tucker-decomposable fields, \(\Phi\), \(\Phi^2\), and \(\Phi\,\nabla^2\Phi\). This reduction is what makes the formulation tractable in a separable tensor-decomposition setting.

• An a priori error analysis. We establish an a priori error estimate for the formulation (Theorem 1), characterizing its convergence behavior on immersed domains.

• Numerical study. The formulation is evaluated on canonical 2D circular and annular domains, a 3D spherical domain, and complex starfish- and gear-shaped domains, demonstrating optimal convergence and robustness on non-Cartesian geometries discretized by regular voxel meshes.

The remainder of the paper is organized as follows. Section 3 develops the exact Dirichlet immersed C-HiDeNN-TD, which constitutes the main contribution of this work. Section 5 presents numerical examples, Section 6 contains a discussion, and Section 7 outlines conclusions and future work. The proof of Theorem 1 is given in Appendix 8.

2 Geometric Representation using mesh-free body-fitted functions↩︎

Different geometric representations have been widely adopted across various computational domains based on their specific application requirements. In computer graphics and geometric modeling, researchers have extensively leveraged the signed distance function (SDF) as an efficient, continuous mathematical representation of shapes. By definition, an SDF calculates the shortest Euclidean distance from any spatial point \(\mathbf{x}\) to the domain boundary \(\partial \Omega\), with the sign designating whether the point resides inside (typically negative) or outside the domain. This implicit formulation is highly advantageous for parametric modeling in the latent space; a notable example is DeepSDF, which utilizes deep auto-decoders to learn continuous SDFs for high-fidelity 3D shape generation and completion [21].

Despite these useful mathematical characteristics, SDFs cannot be directly utilized as geometric descriptors for rigorous physics-based simulation, particularly in the field of Computer-Aided Engineering (CAE). In CAE, spatial fields such as material density or structural compliance possess distinct, conserved physical meanings within the domain \(\Omega\). Because the SDF strictly describes relative proximity to the boundary \(\partial \Omega\) without capturing interior field physics, it fundamentally lacks physical interpretation across the problem domain.

To bridge this gap, we propose a new geometric representation termed body-fitted functions (BFF), whose support is precisely the computational domain \(\Omega\), defined by \[\Phi(\mathbf{x}) = \left\{\begin{array}{cc} =0, & \text{exterior} \\ =0, & \text{on } \partial\Omega \\ >0, & \text{interior} \end{array}\right.\] with a constraint of non-vanishing interior boundary gradient \[\lim_{\mathbf{x}\in \partial \Omega, \eta \to 0^+} |\nabla \Phi(\mathbf{x}-\eta \boldsymbol{n}) | \geq a >0,\] where \(\boldsymbol{n}\) is the unit outward normal vector, and \(a\) is a positive lower bound constant. The BFF thus acts as a smooth indicator of the physical domain: it is positive within the phyical domain, identically zero outside and on the boundary, and possesses a non-vanishing normal gradient at the boundary.

Specifically, we develop a formulation in which Dirichlet boundary conditions are enforced strongly by multiplying the trial function \(v^h\) with a body-fitted mask \(\Phi(\mathbf{x})\) that strictly vanishes on \(\partial \Omega\) there regardless of the discrete representation of \(v^h\). Historically, meshfree methods and physics-informed machine learning have relied on approximate distance functions to satisfy kinematic admissibility a priori [22]. Our approach guarantees exact geometric enforcement, entirely eliminating the boundary-residual penalty terms, and the associated hyperparameter tuning, which are required by traditional Nitsche-type methods. This geometric enforcement is achieved at the price of an additional pre-computation effort to assemble the geometric tensors that arise in the bilinear form.

2.1 SDF-based construction of body-fitted functions↩︎

Importantly, the definition of the BFF does not prescribe a specific functional form: any sufficiently smooth function meeting these conditions qualifies as a body-fitted function. Different applications may call for different constructions, for instance, a polynomial level-set, a tanh-based phase field, or a data-driven neural implicit.

While the general definition admits multiple constructions, a particularly convenient route exploits the SDF, which is often already available from the geometry. Here, we adopt the following SDF-based exponential mask as the BFF: \[\Phi_\beta(\mathbf{x}) \;=\; \begin{cases} 1 - \mathrm{e}^{\,\mathrm{SDF}(\mathbf{x})/\beta}, & \mathbf{x} \in \Omega,\\ 0, & \mathbf{x} \notin \Omega, \end{cases} \label{eq:Phi95def}\tag{1}\] which decays exponentially toward the boundary at a rate set by a scaling parameter \(\beta\). By operating on the standard convention where \(\mathrm{SDF}(\mathbf{x}) \le 0\) inside the domain, this formulation guarantees that \(\Phi(\mathbf{x})\) smoothly approaches zero at the boundary while avoiding the gradient singularities often introduced by absolute distance fields. The parameter \(\beta\) controls how sharply \(\Phi\) transitions from 0 to approximately 1. The proposed body-fitted geometry is illustrated in Figure 2.

2.2 Tucker decomposition of body-fitted functions↩︎

Building upon the proposed body-fitted geometric representation, we next seek to solve the governing partial differential equations (PDEs) through a reduced-order model (ROM) constructed in a purely data-free manner. Traditional ROM techniques, such as Proper Orthogonal Decomposition (POD), typically require an expensive offline phase to generate high-fidelity snapshot data to form a reduced basis. To circumvent this computational bottleneck, we leverage Tucker decomposition to formulate a low-rank representation of the BFF \(\Phi(\mathbf{x})\) using C-HiDeNN Tucker Decomposition as shown in Eq. 2 . \[\Phi(x,y,z)\approx\; \sum_{p=1}^{P}\sum_{q=1}^{Q}\sum_{k=1}^{K} \mathcal{G}_{pqk}\,\Phi_x^{(p)}(x)\,\Phi_y^{(q)}(y)\,\Phi_z^{(k)}(z), \label{eq:tucker95Phi}\tag{2}\] where \(P,Q,K\) are the number of Tucker modes; \(\Phi_{x_i}^{(p)}(x_i)\) refers to the Tucker factor in each mode \(p\) and dimension \(x_i\); \(\mathcal{G}_{pqk}\) is the core tensor that has a reduced dimension of \(P\times Q \times K\). Each factor is interpolated using C-HiDeNN shape function \(\boldsymbol{\widetilde{N}}(x_i)\) [18].

\[u^{h}(x_i)=\sum_{k=1} ^ {nnode}\widetilde{N}_{k}(x_i;a_{k},s_{k},p_{k}, \mathcal{A}_k)u_{k} \label{full95chidenn}\tag{3}\]

where \(u^{h}(x_i)\) is the approximated function; \(nnode\) is the total number of nodes and \(k\) is the nodal index; \(\widetilde{N}_{k}(x_i)\) is the C-HiDeNN shape function at the \(k\)-th node in dimension \(x_i\); \(a_k,s_k,p_k, \mathcal{A}_k\) are the node-wise hyperparameters in the C-HiDeNN shape function, which controls the smoothness of the interpolation. It should be noted that these hyperparameters can vary across different nodal patches since C-HiDeNN can optimize these hyperparameters like machine learning parameters, rendering an adaptable functional space without altering the number of global nodes or hidden layers. Detailed formulation of C-HiDeNN shape function can be found in [16], [23].

As shown in Fig 3, by expressing the high-dimensional discretized space as a compact core tensor multiplied by corresponding factor matrices along each spatial mode, Tucker decomposition systematically reduces the geometric complexity. Operating directly within this low-rank tensor format allows us to evaluate and solve the governing equations in the low-rank representation. Most importantly, projecting the formulation into this latent space drastically reduces the total number of unknowns required to solve the system, yielding significant computational savings while preserving the rigorous enforcement of boundary conditions and the underlying physics.

3 Exact Dirichlet immersed method↩︎

With the body-fitted geometric representation established, we now construct the exact Dirichlet immersed formulation. The trial solution is built by multiplying the unknown field with \(\Phi\), thereby satisfying the Dirichlet condition by construction. The weak form is then derived, and a key algebraic cancellation drastically reduces the number of geometric coefficients.

We consider the following Poisson equation defined in an arbitrary domain \(\Omega\): \[\nabla^2u(\mathbf{x}) + b(\mathbf{x}) = 0 \quad \text{in } \Omega, \label{eq:poisson95strong}\tag{4}\] which is subjected to the Dirichlet boundary condition: \[u(\mathbf{x}) =\; g(\mathbf{x}) \quad \text{on } \partial\Omega, \label{eq:dirichlet95bc}\tag{5}\] where \(u(\mathbf{x})\) is the solution; \(b(\mathbf{x})\) is the forcing term; \(g(\mathbf{x})\) is the boundary condition.

3.1 Strong Dirichlet enforcement using body-fitted function↩︎

We introduce a Cartesian background region \(\mathcal{D}\) that covers the computational domain \(\Omega\), as illustrated in Fig. 4. A structured Cartesian mesh is laid over \(\mathcal{D}\), independent of the geometry of \(\Omega\). The mesh nodes do not necessarily lie on \(\partial\Omega\), and the elements are in general incompatible with the geometry. This incompatibility is resolved by the BFF \(\Phi\), which encodes \(\Omega\) exactly as defined in Section 2. Using a regular Cartesian grid offers two principal advantages: (1) it eliminates mesh generation at essentially zero cost, since the grid is constructed by tensor-product subdivision of \(\mathcal{D}\); and (2) the tensor-product structure is naturally compatible with the Tucker-decomposition-based reduced-order formulation developed in Section 4, where the separation of variables along each Cartesian direction enables significant computational savings.

The trial solution defined on the Cartesian background mesh is constructed as \[u^h(\mathbf{x}) \;=\; \Phi(\mathbf{x})\,v^h(\mathbf{x}) + g^h(\mathbf{x}), \label{eq:strong95split}\tag{6}\] where \(g^h\) is a sufficiently smooth extension of the Dirichlet boundary condition into the background domain \(\mathcal{D}\), and \(\Phi\) is a BFF satisfying \(\Phi(\mathbf{x}) > 0\) in \(\Omega\), \(\Phi = 0\) on \(\partial \Omega\), and \(\Phi \equiv 0\) in the void \(\mathcal{D}\setminus\overline{\Omega}\). Equation 6 therefore enforces \(u^h = g^h = g\) exactly on \(\partial \Omega\), regardless of the discrete representation of \(v^h\), even though the background Cartesian mesh may not conform to the geometry \(\Omega\).

3.1.1 Weak form↩︎

Substituting Eq. 6 into Eq. 4 , the auxiliary variable \(v\) satisfies \[\nabla^2(\Phi v^h) \;=\; -\bigl(\nabla^2 g^h + b\bigr) \quad \text{in } \Omega, \label{eq:strong95v95pde}\tag{7}\] where \(\nabla^2 g^h + b\) become a new source term.

Using \(\Phi w^h\) as the test function (so that \(\Phi w^h = 0\) on \(\partial \Omega\) and the boundary contributions on \(\partial \Omega\) vanish), multiplying Eq. 7 by \(\Phi w^h\) and integrating over \(\mathcal{D}\) yields, after integration by parts, \[\int_{\mathcal{D}} \nabla(\Phi w^h)\cdot \nabla(\Phi v^h) \,\mathrm{d}\mathcal{D} \;=\; \int_{\mathcal{D}} \Phi\, w^h\,\bigl(\nabla^2 g^h + b\bigr) \,\mathrm{d}\mathcal{D}. \label{eq:weak95strong95pre}\tag{8}\]

Expanding the bilinear form, \[\nabla(\Phi w^h)\cdot \nabla(\Phi v^h) \;=\; \Phi^2\,(\nabla w^h \cdot \nabla v^h) + \Phi\, w^h\,(\nabla v^h \cdot \nabla \Phi) + \Phi\, v^h\,(\nabla w^h \cdot \nabla \Phi) + |\nabla \Phi|^2\,w^h\, v^h. \label{eq:phiw95phiv95expanded}\tag{9}\] The two cross terms combine via \(\Phi w^h(\nabla v^h \cdot \nabla \Phi) + \Phi v^h(\nabla w^h \cdot \nabla \Phi) = \tfrac{1}{2}\,\nabla(\Phi^2)\cdot \nabla(w^h v^h)\), and a further integration by parts, together with cancellation with the \(|\nabla \Phi|^2 w^h v^h\) term, yields \[\int_{\mathcal{D}} \Phi^2\,(\nabla w^h \cdot \nabla v^h) \,\mathrm{d}\mathcal{D} \;-\; \int_{\mathcal{D}} (\Phi\,\nabla^2 \Phi)\, w^h\, v^h \,\mathrm{d}\mathcal{D} \;=\; \int_{\mathcal{D}} \Phi\, w^h\,\bigl(\nabla^2 g^h + b\bigr) \,\mathrm{d}\mathcal{D}. \label{eq:weak95strong95unstabilized}\tag{10}\]

The transformed weak form 10 is algebraically equivalent to the standard weak form obtained by substituting 6 into the Poisson equation; both describe exactly the same continuum problem. The difference lies in the structure of the bilinear form. In a naive expansion of the bilinear form 9 , five (2D) or six (3D) distinct geometric coefficient fields appear: \(\Phi^2\), \(\Phi\,\nabla\Phi\) (two or three vector components), \(|\nabla\Phi|^2\), and \(\Phi\) on the right-hand side, each requiring an independent low-rank Tucker approximation in the reduced-order setting. After the algebraic cancellation, the bilinear form involves only three geometric fields, and \(v^h\) couples to the geometry through only two bilinear integrals: a weighted stiffness term (\(\Phi^2 \nabla w^h \cdot \nabla v^h\)) and a weighted mass-like term (\(\Phi\,\nabla^2\Phi\, w^h v^h\)). The right-hand side involves a single geometric coefficient \(\Phi\). Consequently, in the subsequent tensor-decomposition ROM (Section 4), only three coefficient fields must be Tucker-approximated, directly reducing the offline pre-computation and the per-mode assembly cost by nearly a half.

3.1.2 Stabilization in the void↩︎

In the void \(\mathcal{D}\setminus\overline{\Omega}\), both \(\Phi\) and \(\Phi^2\) vanish, so the bilinear form in Eq. 10 is degenerate and yields a singular stiffness matrix. To stabilize the discrete system, we add an isotropic stiffness in \(\mathcal{D}\setminus\overline{\Omega}\): \[\int_{\mathcal{D}} \bigl(\Phi^2 + \gamma(\mathbf{x})\bigr)\,(\nabla w^h \cdot \nabla v^h) \,\mathrm{d}\mathcal{D} \;-\; \int_{\mathcal{D}} (\Phi\,\nabla^2 \Phi)\, w^h\, v^h \,\mathrm{d}\mathcal{D} \;=\; \int_{\mathcal{D}} \Phi\, w^h\,\bigl(\nabla^2 g^h + b\bigr) \,\mathrm{d}\mathcal{D}, \label{eq:weak95strong}\tag{11}\] where \(\gamma(\mathbf{x}) \ll 1\) is a small stabilization function.

Compared with existing immersed methods that enforce Dirichlet conditions through boundary terms, the stabilization strategy in Eq. 11 is structurally simpler. Nitsche-type methods [7]–[10], [24], in both their sharp-interface and diffuse-interface forms, require a consistency term, an adjoint consistency term, and a penalty term. Each of these involves the boundary measure and each demands a separate Tucker compression in the reduced-order setting. The \(\phi\)-FEM approach [11], [12], while sharing the idea of trial-function multiplication with the present work, still relies on a boundary stabilization to control the normal-flux mismatch across the immersed interface, introducing additional surface-localized integrals. In contrast, the proposed formulation stabilizes the system entirely through the volumetric stiffness shift \(\gamma(\mathbf{x})\) in the void. Since the Dirichlet condition is satisfied strongly by construction, no boundary residual, no penalty parameter, and no interface quadrature are needed. This has two practical consequences for tensor-decomposition reduced-order models: (1) only volume fields are Tucker-compressed, rather than multiple boundary-localized fields; (2) the global domain decomposition in TD proceeds uniformly across the entire background \(\mathcal{D}\) without partitioning at the immersed interface, and the 1D sub-problems inherit only smooth volumetric coefficient fields.

3.2 Error estimate↩︎

The exact Dirichlet formulation 11 admits the following a priori error analysis. A sketch of the proof is referred to  8. We emphasize that Eq. 1 is merely one concrete realization of the BFF concept. The theoretical framework developed in this Subsection is independent of such construction, requiring only the defining properties and sufficient regularity of \(\Phi\).

To facilitate the analysis we recall two extension operators. The Stein extension [25], [26] is a total extension operator \(E: H^s(\Omega) \to H^s(\mathbb{R}^d)\) that extends Sobolev functions from a Lipschitz domain \(\Omega\) to the whole space while preserving the \(H^s\) regularity. For the error estimate we require a Stein extension \(\tilde{v} \in H^{p+1}(\mathcal{D})\) that coincides with \(v\) on \(\Omega\) and possesses the same smoothness as the exact solution across the background domain. The harmonic extension \(\hat{v}\) of \(v\) into the void \(\mathcal{D}\setminus\overline{\Omega}\) is defined as the solution of \(\Delta \hat{v} = 0\) in \(\mathcal{D}\setminus\overline{\Omega}\) with \(\hat{v} = v\) on \(\partial\Omega\), and \(\hat{v} = v\) in \(\Omega\), whose weak form is precisely the stabilization term with the weighting coefficient \(\gamma(\mathbf{x})\). While the harmonic extension exactly maintains the Laplacian of the interior solution in the void that differs from the original equation, it is generally only \(C^0\) across \(\partial\Omega\), whereas the Stein extension preserves the full \(H^{p+1}\) regularity at the expense of not satisfying any PDE in the void. The mismatch \(\|\hat{v} - \tilde{v}\|_{H^1(\mathcal{D}\setminus\overline{\Omega})}\) appearing in the error bound quantifies this regularity gap.

In Eq. (?? ), the first term represents the C-HiDeNN interpolation error. The second and third terms both arise from the stabilization: the second term accounts for the deviation from the original governing equation inside the computational domain \(\Omega\), while the third term captures the smoothness issue across the boundary \(\partial \Omega\) caused by the harmonic extension outside \(\Omega\).

Corollary 1 shows that setting \(\gamma = 0\) inside \(\Omega\) eliminates the second term of ?? entirely, leaving only the C-HiDeNN interpolation error (first term) and the extension mismatch from the void stabilization (third term). This is the most natural choice for the exact Dirichlet formulation: since the BFF already enforces the boundary condition strongly, there is no need for interior stabilization to restore coercivity, and the error is dictated solely by the approximation properties of the trial space and the regularity gap across \(\partial\Omega\).

Corollary 2 provides a practical guideline for choosing the void stabilization. When \(\gamma\) scales as \(\mathcal{O}(h^{2p})\), the third error term, which arises from the smoothness issue of the harmonic extension, converges at the same rate \(\mathcal{O}(h^p)\) as the interpolation term. Balancing the two sources of error in this way recovers the overall \(h\)-convergence rate of the underlying C-HiDeNN approximation and prevents the void stabilization from becoming the accuracy bottleneck. The requirement that the fictitious-domain stabilization parameter be mesh-dependent to preserve high-order convergence is not specific to the present formulation; it is a general constraint across immersed and fictitious-domain methods. The \(\phi\)-FEM framework [11], [12], though enforcing Dirichlet conditions strongly through trial-function multiplication, requires a boundary stabilization whose magnitude must be adapted to the mesh size to control the normal-flux mismatch across the immersed interface. Similarly, in CutFEM [9], [10], the ghost-penalty form carries an intrinsic \(h\)-dependent scaling as well. In the finite cell method [3], [5], [6], the void stiffness parameter \(\alpha\) is chosen sufficiently small (between \(10^{-5}\) and \(10^{-10}\)) for high-order basis functions. In practice, setting \(\gamma\) to a small constant (\(10^{-10}\) in our numerical experiments, Section 5) satisfies this scaling for all mesh sizes of interest, since \(\gamma^{1/2} \ll h^p\) for moderate \(h\).

3.3 Regularity of the body-fitted function↩︎

Theorem 1 assumes the BFF \(\Phi\) to be sufficiently smooth on \(\mathcal{D}\), so that the geometric coefficients \(\Phi^2\) and \(\Phi\,\nabla^2\Phi\) appearing in the weak form 11 are well-defined and the interpolation estimates used in the proof are valid. In practice, however, the signed distance function from which \(\Phi\) is constructed is typically only \(C^0\) continuous, as the SDF is continuous everywhere but its gradient is discontinuous across the medial axis (the skeleton of \(\Omega\)), where the nearest point on \(\partial\Omega\) jumps between distinct boundary segments [27].

The exponential map defined in Eq. 1 serves as a regularization of the original SDF as well. To see this, compute the gradient of \(\Phi_\beta\): \[\nabla \Phi_\beta \;=\; -\frac{1}{\beta}\,\mathrm{e}^{\mathrm{SDF}/\beta}\,\nabla\mathrm{SDF} \;=\; \frac{1}{\beta}\,(1 - \Phi_\beta)\,\nabla\mathrm{SDF}, \label{eq:grad95Phi95smooth}\tag{12}\] where we have used \(|\nabla\mathrm{SDF}| = 1\) almost everywhere. At any point where \(|\mathrm{SDF}| \gg \beta\), the factor \(1 - \Phi_\beta = \mathrm{e}^{\mathrm{SDF}/\beta}\) is exponentially small, so \(\nabla\Phi_\beta \approx 0\) regardless of the direction of \(\nabla\mathrm{SDF}\). Consequently, if the parameter \(\beta\) is chosen small enough that the medial-axis locus lies deep in the interior where \(|\mathrm{SDF}| \gg \beta\), then \(\Phi_\beta \approx 1\) and \(\nabla\Phi_\beta \approx 0\) at those points, and the \(C^0\) irregularity of the SDF is exponentially attenuated in \(\Phi_\beta\) and its derivatives. A practical guideline is to set \(\beta \ll d_{\mathrm{med}}\), where \(d_{\mathrm{med}}\) is the minimum distance from \(\partial\Omega\) to the medial axis, so that the kinks of the SDF are effectively invisible to the body-fitted formulation. This regularization property is a key practical advantage of the exponential construction over, e.g., a linear ramp \(\Phi = \max(-\mathrm{SDF}, 0)\), for which \(\nabla\Phi\) inherits the full derivative discontinuity of the SDF.

4 Immersed C-HiDeNN-TD using exact Dirichlet immersed form↩︎

For large-scale problems, the computational cost can be reduced by model order reduction. The C-HiDeNN-Tensor Decomposition (TD) method, that leverages a special form of tensor decomposition called canonical polyadic decomposition, can be applied directly to the Cartesian background \(\mathcal{D}\).

The discrete trial and test functions follow the same TD structure as in the C-HiDeNN-TD framework [16], [18]: \[v^h(x,y,z) \;=\; \sum_{m=1}^{M} v_x^{(m)}(x)\,v_y^{(m)}(y)\,v_z^{(m)}(z), \qquad w^h \;=\; \delta v^h, \label{eq:cp95v}\tag{13}\] with the boundary lift \(g^h\) represented as a rank-\(M_g\) TD approximation \[g^h(x,y,z) \;=\; \sum_{m=1}^{M_g} g_x^{(m)}(x)\,g_y^{(m)}(y)\,g_z^{(m)}(z), \label{eq:cp95g}\tag{14}\] where each factor \(v_\xi^{(m)}, g_\xi^{(m)}\) (\(\xi \in \{x,y,z\}\)) is represented using 1D C-HiDeNN shape functions.

To improve the efficiency of multi-dimensional integrals in the weak form, the key coefficients \(\Phi\), \(\Phi^2\), and \(\Phi\,\nabla^2\Phi\) are also approximated by Tucker forms: \[\begin{align} && \Phi \approx \hat{\Phi}^{(1)} = \sum_{p=1}^{P_1} \sum_{q=1}^{Q_1} \sum_{k=1}^{K_1} \mathcal{H}^{(1)}_{pqk} \, \phi_{x}^{(p)}(x) \, \phi_{y}^{(q)}(y) \, \phi_{z}^{(k)}(z), \\ && \Phi^2 \approx \hat{\Phi}^{(2)} = \sum_{p=1}^{P_2} \sum_{q=1}^{Q_2} \sum_{k=1}^{K_2} \mathcal{H}^{(2)}_{pqk} \, \varphi_{x}^{(p)}(x) \, \varphi_{y}^{(q)}(y) \, \varphi_{z}^{(k)}(z), \\ && \Phi\,\nabla^2\Phi \approx \hat{\Phi}^{(3)} = \sum_{p=1}^{P_3} \sum_{q=1}^{Q_3} \sum_{k=1}^{K_3} \mathcal{H}^{(3)}_{pqk} \, \theta_{x}^{(p)}(x) \, \theta_{y}^{(q)}(y) \, \theta_{z}^{(k)}(z). \label{tucker95phi} \end{align}\tag{15}\] It is important to note that \(\hat{\Phi}^{(1)}\), \(\hat{\Phi}^{(2)}\) and \(\hat{\Phi}^{(3)}\) are decomposed independently using separate Tucker factorizations with distinct factor matrices and core tensors. Consequently, the algebraic relations \(\hat{\Phi}^{(2)} = (\hat{\Phi}^{(1)})^2\) and \(\hat{\Phi}^{(3)} = \hat{\Phi}^{(1)}\,\nabla^2\hat{\Phi}^{(1)}\) that hold for the exact geometric fields are generally not preserved in the Tucker-compressed representation. This deviation introduces an additional geometric consistency error. Enforcing these transformation relations by deriving \(\hat{\Phi}^{(2)}\) and \(\hat{\Phi}^{(3)}\) directly from \(\hat{\Phi}^{(1)}\) would eliminate this source of error but may increase the tensor rank required to represent the derived fields.

Substituting 13 into 11 and contracting against the Tucker cores 15 produces three contributions to the assembled 1D operator \(\mathbf{K}^{[\xi]}\) in each PGD direction \(\xi \in \{x,y,z\}\): a stiffness term from \(\Phi^2 + \gamma\), a mass-like term from \(\Phi\,\nabla^2 \Phi\), and a coupling term contributing to the right-hand side from \(\Phi\) and \(g^h\).

To circumvent the curse of dimensionality, the trial space is enriched mode by mode using PGD. At enrichment step \(M\), the solution is split into a known part (the first \(M-1\) modes) and an unknown rank-1 increment, \[v^h(x,y,z) \;=\; \underbrace{\sum_{m=1}^{M-1} v_x^{(m)}(x)\,v_y^{(m)}(y)\,v_z^{(m)}(z)}_{v^{<}(\mathbf{x})} \;+\; \bar{v}_x(x)\,\bar{v}_y(y)\,\bar{v}_z(z), \label{eq:pgd95split95exact}\tag{16}\] and the three 1D factors \(\bar{v}_x, \bar{v}_y, \bar{v}_z\) are computed by an alternating fixed-point iteration. Restricting the variation to the \(x\) component, \(\delta v^h = \delta\bar{v}_x(x)\,\bar{v}_y(y)\,\bar{v}_z(z)\), and substituting into 11 yields a 1D linear system \[\mathbf{K}^{[x]} \mathbf{v}_x \;=\; \mathbf{q}^{[x]}, \label{eq:pgd951d95system95exact}\tag{17}\] in which the operator \(\mathbf{K}^{[x]}\) and right-hand side \(\mathbf{q}^{[x]}\) are assembled by contracting the Tucker cores 15 against scalars obtained from 1D quadrature in the fixed directions \(y\) and \(z\). The procedure is repeated cyclically over \(x \to y \to z\) until the relative change in the rank-1 mode falls below a prescribed tolerance; the new mode is then appended to \(v^{<}\) and the next enrichment step is initiated.

4.1 Error estimate for immersed TD↩︎

In the immersed TD formulation, the BFF \(\Phi\) and its derived fields \(\Phi^2\) and \(\Phi\,\nabla^2\Phi\) are approximated respectively by independent Tucker decompositions \(\hat{\Phi}^{(1)}\), \(\hat{\Phi}^{(2)}\) and \(\hat{\Phi}^{(3)}\) as defined in Eq. 15 . This geometric compression introduces additional deviations whose impact is quantified by the following error estimate. Define the Tucker perturbation fields \[\label{eq:delta95def} \delta_1 = \hat{\Phi}^{(1)} - \Phi,\qquad \delta_2 = \hat{\Phi}^{(2)} - \Phi^2,\qquad \delta_3 = \hat{\Phi}^{(3)} - \Phi\,\nabla^2\Phi.\tag{18}\]

Comparing Theorem 2 with Theorem 1, the perturbation \(\delta_2\) enters the estimate in exactly the same manner as the void stabilization parameter \(\gamma\) (second and third terms). This reveals a fundamental structural analogy: the Tucker compression error in \(\hat{\Phi}^{(2)}\) acts as an effective additional stabilization whose \(L^\infty\)-norm combines additively with \(\gamma\). Consequently, the accuracy of the Tucker-compressed formulation is controlled by the total stabilization magnitude \(\gamma + \delta_2\), and the Tucker rank must be chosen such that \(\|\delta_2\|_{L^\infty}\) remains commensurate with \(\gamma\) to avoid degrading the convergence rate. The \(\delta_3\) perturbation, arising from the independent compression of \(\Phi\,\nabla^2\Phi\), contributes a mass-like error term (fourth term) that scales with \(\|\delta_3\|_{L^\infty}\). The \(\delta_1\) perturbation propagates directly through the right-hand side (fifth term). A complete proof is given in Appendix 9.

5 Numerical examples↩︎

In immersed methods, the diffuse interface (the narrow band where \(0 < \Phi < 1\)) requires higher-order quadrature to resolve the sharp gradients of the BFF and its derivatives. Applying uniformly dense Gauss quadrature over all elements is unnecessarily expensive: elements deep inside \(\Omega\) (where \(\Phi \approx 1\) and its derivatives are smooth) and elements in the void (where \(\Phi = 0\)) do not benefit from the extra integration points.

We adopt a two-level adaptive Gauss quadrature strategy. A pre-classification step evaluates the nodal values of \(\Phi\) on each element and labels elements intersecting the transition zone (those with at least one node satisfying \(\varepsilon_{\text{ext}} < \Phi < \varepsilon_{\text{int}}\)) as boundary elements; all others are non-boundary elements. The implementation follows Algorithm 5.

In the assembly of the stiffness matrix and force vector, elements in \(\mathcal{E}_B\) are integrated with \(n_g^{\text{bnd}}\) Gauss points per direction (\(20 \times 20\) in the present computations), while elements in \(\mathcal{E}_N\) use a coarser \(n_g^{\text{int}}\) rule (\(6 \times 6\)). The error computation employs a slightly enriched rule (\(n_g^{\text{bnd}}+2\), \(n_g^{\text{int}}+2\)) to avoid quadrature-dominated error estimates. For the multiply-connected void domain, the boundary element set \(\mathcal{E}_B\) typically accounts for approximately \(10\)–\(15\%\) of the total elements; the adaptive strategy therefore reduces the quadrature cost by roughly a factor of \(5\)–\(7\) compared with a uniform \(20 \times 20\) rule, with no measurable loss of accuracy. Table 1 summarizes the quadrature configuration.

5.1 2D circular / annular computational domain↩︎

The exact Dirichlet immersed formulation is evaluated on two 2D Poisson benchmark problems using the BFF \(\Phi_\beta\) defined in Eq. 1 , embedded in the square background domain \(\mathcal{D} = [-1,1] \times [-1,1]\). The void stabilization is set to \(\gamma(\mathbf{x})=10^{-10}\) in \(\mathcal{D}\), which is sufficiently small to be considered negligible within the range of mesh resolutions of interest and does not affect the observed convergence rates. The C-HiDeNN radial basis expansion parameter is \(r_a=50\) throughout.

5.1.1 2D circular computational domain↩︎

The first benchmark is the Poisson equation on a circular domain \(\Omega = \{ (x,y) \mid \sqrt{x^2+y^2} < R\}\) of radius \(R=0.7\), subject to homogeneous Dirichlet boundary conditions. The signed distance function for this geometry admits the closed form \[\mathrm{SDF}(x,y) = \sqrt{x^2 + y^2} - R, \label{eq:sdf95circle}\tag{19}\] with the sign convention illustrated in Fig. 1 (negative inside \(\Omega\), positive outside). The BFF \(\Phi_\beta\) is then constructed from Eq. 1 . Figure 6 shows the SDF and the resulting BFF \(\Phi_\beta\) for the circular domain.

The exact solution is chosen as \[u^{\text{Ext}}(x,y) = R^2 - x^2 - y^2, \label{eq:circle95exact}\tag{20}\] which vanishes on \(\partial\Omega\) and yields a constant source term \(b(x,y)=4\).

Figure 7 shows the convergence of the relative \(L^2\) and energy errors for \(p = 1, 2, 3\) with \(\beta = 0.01\), under uniform mesh refinement (\(n_x = 8, 16, 32, 64, 128, 256\)). The estimated convergence rates are summarized in Table 2.

For \(p = 1\) and \(p = 2\), the estimated rates are close to their respective theoretical optima. For \(p = 3\), the rates are \(\mathcal{O}(h^{3.4})\) in \(L^2\)-norm and \(\mathcal{O}(h^{2.5})\) in energy-norm, below the nominal \(\mathcal{O}(h^{4})\) and \(\mathcal{O}(h^{3})\). This suboptimality is attributed to the combined effects of the reduced regularity of \(\Phi_\beta\) in the transition layer, the finite accuracy of the adaptive boundary-integration quadrature, and the fact that, on the finer meshes, the error is sufficiently small that these previously subdominant error sources become visible.

We isolate the role of the void stabilization parameter \(\gamma\) by fixing the mesh (\(128 \times 128\)) and sweeping \(\gamma\) over nine decades, from \(10^{-2}\) to \(10^{-12}\). All other parameters are held constant (\(p=3\), \(s=3\), adaptive Gauss quadrature with \(20\times20\) boundary points and \(6\times6\) interior points).

Figure 8 plots the relative \(L^2\) and energy errors as functions of \(\gamma\) on a log-log scale. Two regimes are evident: (1) Stabilization-dominated regime (\(\gamma \gtrsim 10^{-7}\)). For large \(\gamma\), the void stabilization distorts the solution throughout the physical domain \(\Omega\). The energy error grows as \(\mathcal{O}(\gamma^{0.5})\), consistent with the \(\sqrt{\gamma}\) scaling of the third term in Corollary 2, where a linear regression over points 2-5 (corresponding to \(\gamma = 10^{-3}\) through \(10^{-6}\)) yields a fitted slope of approximately \(0.5\) for the energy norm. The \(L^2\) error exhibits a fitted slope of approximately \(1\) over the same interval. (2) Discretization-dominated regime (\(\gamma \lesssim 10^{-7}\)). As \(\gamma\) is further reduced, both the \(L^2\) and energy errors saturate to a plateau. In this regime the \(\gamma\)-induced error has fallen below the intrinsic discretization error of the \(128\times128\) mesh, so further reduction of \(\gamma\) yields no improvement in accuracy. This shows that choosing \(\gamma = 10^{-10}\) is appropriate in the above convergence study.

5.1.2 2D annular domain with an eccentric hole↩︎

The second benchmark introduces a non-convex, multiply-connected geometry: a circular outer boundary of radius \(R_1=0.7\) containing an eccentric circular hole of radius \(R_2=0.2\) centered at \((x_0, y_0)=(0.15, 0)\). For a multiply-connected domain defined by the set difference of an outer and an inner region, the signed distance function combines the individual distance fields via the maximum operator: \[\mathrm{SDF}(x,y) = \max\!\Bigl(\sqrt{x^2 + y^2} - R_1,\; R_2 - \sqrt{(x - x_0)^2 + y^2}\Bigr), \label{eq:sdf95void}\tag{21}\] which is negative when both \(\sqrt{x^2+y^2} < R_1\) (inside the outer circle) and \(\sqrt{(x-x_0)^2+y^2} > R_2\) (outside the inner hole) are satisfied. See Fig. 1 for the sign convention. The BFF follows from Eq. 1 . Figure 9 shows the SDF and the resulting BFF \(\Phi_\beta\) for the annular domain.

The exact solution, \[u^{\text{Ext}}(x,y) = -(R_1^2 - x^2 - y^2)\,\bigl(R_2^2 - (x - x_0)^2 - y^2\bigr), \label{eq:void95exact}\tag{22}\] is constructed to vanish on both the outer boundary and the inner hole, satisfying the homogeneous Dirichlet condition. The corresponding source term \(b(x,y)=-\Delta u^{\text{Ext}}\) is obtained analytically.

The key challenge in this example is the accurate resolution of the narrow ligament between the outer boundary and the eccentric hole, which tests the geometric fidelity of the body-fitted description under aggressive mesh refinement. Figure 10 shows the convergence of the relative \(L^2\) and energy errors for polynomial orders \(p = 1, 2, 3\), all using the adaptive Gauss quadrature of Table 1. The mesh size is defined as \(h = 2/n_x\), where \(n_x\) is the number of elements in the \(x\)-direction.

For \(p = 1\), the \(L^2\) error converges approximately as \(\mathcal{O}(h^{2.2})\) and the energy error as \(\mathcal{O}(h^{1.2})\), matching the expected optimal rates \(\mathcal{O}(h^{p+1})\) and \(\mathcal{O}(h^{p})\) respectively. For \(p = 2\), the \(L^2\) error decays as \(\mathcal{O}(h^{2.8})\) and the energy error as \(\mathcal{O}(h^{1.7})\), again in agreement with the theory. For \(p = 3\), the \(L^2\) rate is approximately \(\mathcal{O}(h^{3.4})\) and the energy rate is \(\mathcal{O}(h^{2.5})\); the slight suboptimality relative to the nominal \(\mathcal{O}(h^4)\) and \(\mathcal{O}(h^3)\) rates is attributed to the narrow-ligament geometry, where the element size must be sufficiently small to resolve the curvature of both the outer and inner boundaries simultaneously; this condition is met only on the finest meshes. The convergence rates are summarized in Table 3.

5.2 2D complex geometries↩︎

5.2.1 2D starfish computational domain with a central hole↩︎

The next benchmark considers the Poisson equation on a more complex geometry: a 7-petaled starfish domain with a central hole. The domain is embedded within a square background box \(\mathcal{D} = [0,2] \times [0,2]\) centered at \(\mathbf{c}=(1,1)\). Using local polar coordinates \(\rho = \sqrt{(x-1)^2+(y-1)^2}\) and \(\theta = \operatorname{atan2}(y-1, x-1)\), the computational domain is defined as \[\Omega = \Big\{ (x,y) \mid r_{\text{hole}} \le \rho \le R(\theta) \Big\},\] where the outer boundary is given by \(R(\theta) = r_{\text{base}} + \alpha \cos(7\theta)\). The geometric parameters are set to a mean petal radius \(r_{\text{base}} = 0.70\), petal amplitude \(\alpha = 0.24\), and a central hole radius \(r_{\text{hole}} = 0.16\). Consequently, the domain boundary \(\partial\Omega\) consists of two distinct components: the outer wavy curve and the inner circular void.

Homogeneous Dirichlet boundary conditions (\(u=0\)) are prescribed on all of \(\partial\Omega\), with constant body force \(b = 1.0\) as in Eq. 4 .

The Dirichlet conditions are enforced geometrically through the BFF \(\Phi\), utilizing the ansatz \(u^h = \Phi v^h\) (homogeneous case) with \(v^h\) pinned to \(0\) on the boundary of the background box. The BFF is defined via Eq. 1 with transition parameter \(\beta = 0.002\). The immersed weak form follows Eq. 11 : \[A_\gamma(w^h, v^h) = \int_{\mathcal{D}} (\Phi^2 + \gamma) \nabla w^h \cdot \nabla v^h - \Phi \nabla^2\Phi \, w^h v^h \,\mathrm{d}\mathcal{D},\] with the corresponding load functional \(L(w^h) = \int_{\mathcal{D}} \Phi w^h b \,\mathrm{d}\mathcal{D}\). The void stabilization parameter is set to \(\gamma = 10^{-10}\).

Unlike the circular benchmark, this problem does not possess a closed-form exact solution. Instead, the numerical results are compared against a high-fidelity reference solution obtained using a body-fitted \(P_1\) finite element solve on an identical triangulated starfish domain via FEniCSx. Figure 11 visualizes this comparison, displaying the reference field, the immersed C-HiDeNN solution, and the point-wise absolute error.

5.2.2 2D gear computational domain with a central hole↩︎

To further evaluate the framework’s capability in resolving sharp geometric features, the third benchmark considers the Poisson equation on a 20-tooth gear domain with a central hole. The domain is embedded within the same square background box \(\mathcal{D} = [0,2] \times [0,2]\) centered at \(\mathbf{c}=(1,1)\). Using the local polar coordinates \(\rho\) and \(\theta\) defined previously, the computational domain is defined as \[\Omega = \Big\{ \text{inside the 20-tooth gear polygon} \Big\} \setminus \Big\{ \rho < r_{\text{hole}} \Big\}.\]

The outer boundary is defined as a piecewise-linear polygon constructed from \(80\) vertices located at angles \(\theta_i = 2\pi i/80\). The corresponding radii cycle with a period of \(4\) in the sequence \([r_{\text{in}}, r_{\text{out}}, r_{\text{out}}, r_{\text{in}}]\). This generates square-like teeth, where each tooth rises from the root radius to a flat tip at the tooth radius and back via straight segments. The geometric parameters are set to a tooth-tip radius \(r_{\text{out}} = 0.90\), a tooth-root radius \(r_{\text{in}} = 0.70\), and a central hole radius \(r_{\text{hole}} = 0.50\). Thus, the domain boundary \(\partial\Omega\) consists of the outer 20-tooth gear polygon and the inner circular void.

In direct contrast to the smooth analytic curve of the starfish domain, the gear’s outer edge introduces multiple sharp, re-entrant corners at the \(80\) polygon vertices. These discontinuities present a significantly more stringent geometric challenge for resolving the BFF \(\Phi\).

The physical problem remains identical to the starfish benchmark (constant body force \(b = 1.0\), homogeneous Dirichlet boundary conditions), with the BFF, immersed weak form, and stabilization parameters (\(\gamma = 10^{-10}\), \(\beta = 0.002\)) retained from the previous examples.

Similar to the starfish benchmark, the numerical results are compared against a high-fidelity reference solution obtained via a body-fitted \(P_1\) finite element solve on an identical triangulated gear domain using FEniCSx. Figure 12 displays the reference solution, the immersed C-HiDeNN solution, and the point-wise absolute error, illustrating the framework’s performance around the sharp polygonal features.

5.2.3 3D sphere with a concentric spherical void↩︎

The final benchmark extends the immersed body-fitted formulation to three dimensions on a multiply-connected solid: a spherical shell formed by an outer sphere of radius \(R_1=0.7\) with a concentric spherical void of radius \(R_2=0.3\), both centered at \(\mathbf{x}_c=(1,1,1)\) inside the background box \([0,2]^3\). Writing \(r=\lVert\mathbf{x}-\mathbf{x}_c\rVert\), the signed distance function again combines the outer and inner fields through the maximum operator, \[\mathrm{SDF}(\mathbf{x}) = \max\!\bigl(r - R_1,\; R_2 - r\bigr), \label{eq:sdf95void3d}\tag{23}\] which is negative precisely in the shell \(R_2 < r < R_1\). The BFF \(\Phi_\beta\) follows from Eq. 1 . The exact solution, \[u^{\text{Ext}}(\mathbf{x}) = \bigl(R_1^2 - r^2\bigr)\bigl(r^2 - R_2^2\bigr), \label{eq:void3d95exact}\tag{24}\] vanishes on both the outer sphere and the inner void, enforcing the homogeneous Dirichlet condition, and the source term \(b=-\Delta u^{\text{Ext}}\) is obtained analytically.

Figure 13 shows the relative \(L^2\) and energy errors for \(s=p=1\) and \(s=p=2\), with mesh size \(h = 2/n_x\). The measured rates over the finest meshes, summarized in Table 4, recover the expected optimal orders \(\mathcal{O}(h^{p+1})\) in \(L^2\) and \(\mathcal{O}(h^{p})\) in energy, confirming that the body-fitted description retains full accuracy in three dimensions.

6 Discussions↩︎

The Immersed Tensor Decomposition (ITD) framework successfully integrates non-conforming geometries into reduced-order models by intertwining a continuous body-fitted representation with a separable C-HiDeNN-TD solver. This approach introduces several key advantages and practical considerations:

Strong Boundary Enforcement without Interface Quadrature: Unlike classical immersed methods that rely on weak Nitsche-type enforcement, ITD satisfies Dirichlet conditions strongly via the BFF \(\Phi\). This eliminates penalty parameter tuning and complex interface quadratures, preserving the purely Cartesian tensor-product structure required for efficient reduced-order modeling.

Algebraic Complexity Reduction: A direct substitution of the body-fitted trial function into the weak form would generate six distinct geometric coefficient fields in 3D (five in 2D). Our exact algebraic cancellation reduces this to just three scalar fields (\(\Phi\), \(\Phi^2\), and \(\Phi\,\nabla^2\Phi\)). This reduction cuts the tensor compression overhead by nearly a half. Balancing Hyperparameters (\(\beta\) and \(\gamma\)): Practical implementation requires balancing the exponential transition scaling \(\beta\)—which must smooth the SDF without demanding excessively dense quadrature—against the void stabilization \(\gamma\). As established numerically and theoretically, \(\gamma\) must scale as \(\mathcal{O}(h^{2p})\) to effectively stabilize the void without polluting the solution in the physical domain.

7 Conclusions and Future Work↩︎

This paper presented the Immersed Tensor Decomposition (ITD) framework, which successfully integrates mesh-free geometric representations with the separable C-HiDeNN-TD reduced-order solver. By utilizing a BFF (\(\Phi\)) derived from a signed-distance field, the formulation achieves strong, exact enforcement of Dirichlet boundary conditions on regular Cartesian background grids. This approach circumvents the need for conforming body-fitted meshes, variational penalty parameters, and irregular interface quadratures.A central mathematical contribution of this work is the algebraically transformed weak form, which collapses the geometric complexity from six initial geometric coefficients in 3D (five in 2D) into just three volume fields (\(\Phi\), \(\Phi^2\), and \(\Phi\,\nabla^2\Phi\)). This reduction is critical for maintaining the efficiency of the low-rank Tucker decomposition and avoiding the curse of dimensionality. We established a rigorous a priori error estimate demonstrating that the ITD formulation preserves the optimal \(\mathcal{O}(h^{p+1})\) convergence rate of the underlying high-order C-HiDeNN interpolation, provided the void stabilization parameter and the Tucker compression perturbations are properly scaled. Numerical evaluations on canonical 2D domains and 3D geometries confirmed these theoretical bounds, demonstrating robust sub-domain resolution and optimal energy norm convergence. Future research will focus on extending the ITD framework in the following directions:

Vector-field and non-linear PDEs: While the current framework was derived and validated for the scalar Poisson equation, extending the formulation to coupled vector-field problems, such as linear elasticity, is necessary for broader computational mechanics and computer-aided engineering applications.

Direct ingestion of voxel-native data: We intend to deploy the ITD framework directly on raw volumetric datasets, such as medical imaging (CT/MRI) and additive manufacturing lattice structures. In these scenarios, the generation of signed-distance fields and subsequent tensor compression can completely bypass traditional, information-lossy CAD surface reconstruction workflows.

Adaptive rank selection: The derived error estimates indicate that Tucker perturbation fields act as an effective additional stabilization. Developing an adaptive rank-selection algorithm to dynamically balance the tensor compression error against the intrinsic discretization error will further optimize offline pre-computation costs.

8 Proof of Theorem 1↩︎

For simplicity, we consider homogeneous Dirichlet boundary condition (\(g^h = 0\)); the non‑homogeneous case follows by a standard lifting argument. The stabilized weak form can be written as \[\label{eq:general95weak} A_\gamma(W^h, u^h) = L(W^h, b),\tag{25}\] where the bilinear form \(A_\gamma(\cdot,\cdot)\) and the linear functional \(L(\cdot,b)\) in Eq. (11 ) are defined by \[\begin{align} A_\gamma(W^h, u^h) &=& \int_{\mathcal{D}} \nabla W^h \cdot \nabla u^h \,\mathrm{d}\mathcal{D} + \int_{\mathcal{D}} \gamma(\mathbf{x}) \nabla (W^h/\Phi) \cdot \nabla (u^h/\Phi) \,\mathrm{d}\mathcal{D}, \\ & = & \int_{\mathcal{D}} \bigl(\Phi^2 + \gamma(\mathbf{x})\bigr)\,(\nabla (W^h/\Phi) \cdot \nabla (u^h/\Phi)) \,\mathrm{d}\mathcal{D} \;-\; \int_{\mathcal{D}} (\nabla^2 \Phi / \Phi)\, W^h\, u^h \,\mathrm{d}\mathcal{D}, \\ L(W^h, b) &=& \int_{\mathcal{D}} W^h\,b \,\mathrm{d}\mathcal{D}. \end{align}\] The test and trial spaces coincide: \[\mathcal{V}^h = \mathcal{S}^h = \bigl\{ W^h = \Phi w^h :\; w^h = \mathcal{I}^h w,\; w \in H^{p+1}(\mathcal{D}),\; w|_{\partial\mathcal{D}} = 0 \bigr\},\] where \(\mathcal{I}^h\) denotes the C‑HiDeNN interpolation operator. The energy norm is \(\|w\|_\gamma = \sqrt{A_\gamma(w,w)}\).

Proof of Theorem 1. Let the exact solution be extended to the whole background domain by \(u = \Phi \hat{v}\), where \(\hat{v}\) is the harmonic extension. Since \(\hat{v}\) may have limited smoothness across \(\partial\Omega\), we introduce a Stein extension \(\tilde{v} \in H^{p+1}(\mathcal{D})\) with \(\tilde{v}|_\Omega = v|_\Omega\). Using ?? , we estimate the two terms separately.

8.0.0.1 Approximation error

\[\begin{align} 2\inf_{z^h \in \mathcal{S}^h} \|u - z^h\|_\gamma &\le 2 \|u - \Phi\tilde{v}\|_\gamma +2 \inf_{z^h} \|\Phi\tilde{v} - z^h\|_\gamma \\ &= 2 \Bigl( \int_{\mathcal{D}\setminus\Omega} \gamma\, \bigl|\nabla(\hat{v}-\tilde{v})\bigr|^2 \,\mathrm{d}\mathcal{D} \Bigr)^{1/2} + C_1 h^p\|\tilde{v}\|_{H^{p+1}(\mathcal{D})} \\ &\le C_3 \sqrt{\sup_{\mathbf{x}\in\mathcal{D}\setminus\Omega} \gamma(\mathbf{x})}\,\|\hat{v} - \tilde{v}\|_{H^1(\mathcal{D}\setminus\overline{\Omega})} + C_1 h^{p}\, \|\tilde{v}\|_{H^{p+1}(\mathcal{D})}, \end{align}\] where the first term reflects the mismatch between the harmonic and smooth extensions in the void, and the second term follows from the approximation properties of C‑HiDeNN.

8.0.0.2 Consistency error

The exact solution satisfies the unstabilized weak form \(A_0(W^h, u) = L(W^h, b)\). Consequently, \[A_\gamma(\eta^h, u) - L(\eta^h, b) = \int_{\mathcal{D}} \gamma\, \nabla(\eta^h/\Phi) \cdot \nabla(\hat{v}) \,\mathrm{d}\mathcal{D}.\] As \(\int_{\mathcal{D}\setminus\Omega} \nabla w \cdot \nabla(\hat{v}) \,\mathrm{d}\mathcal{D}\) for any weighting function \(w\in H^1(\mathcal{D}\setminus\Omega)\) and \(w|_{\partial \Omega} = 0\), this residual reduces to \[A_\gamma(\eta^h, u) - L(\eta^h, b) = \int_{\Omega} \gamma\, \nabla(\eta^h/\Phi) \cdot \nabla(\hat{v}) \,\mathrm{d}\Omega.\] By the Cauchy–Schwarz inequality and the definition of the \(\gamma\)‑norm, \[\begin{align} \bigl| A_\gamma(\eta^h, u) - L(\eta^h, b) \bigr| &\le \sup_{\mathbf{x}\in\Omega \gamma(\mathbf{x})} \bigl\| \nabla(\eta^h/\Phi) \bigr\|_{L^2(\Omega)} \, \bigl\| \nabla(\hat{v}) \bigr\|_{L^2(\Omega)}. \end{align}\] Using the inequality \[\|\eta^h\|^2_\gamma \geq A_0 (\eta^h, \eta^h)=\|\nabla \eta^h\|^2_{L^2(\Omega)},\] we obtain \[\sup_{\eta^h \in \mathcal{V}^h\setminus\{0\}} \frac{|A_\gamma(\eta^h, u) - L(\eta^h, b)|}{\|\eta^h\|_\gamma} \le C_2 \sup_{\mathbf{x}\in\Omega} \gamma(\mathbf{x}) \sup_{\eta^h \in \mathcal{V}^h\setminus\{0\}} \dfrac{\|\nabla \eta^h/\Phi\|_{L^2(\Omega)}}{\|\nabla \eta^h\|_{L^2(\Omega)}} \cdot \|u\|_{H^1(\Omega)}.\]

Combining the two bounds in ?? gives \[\|u - u^h\|_\gamma \le C_1 h^{p}\, \|\tilde{v}\|_{H^{p+1}(\mathcal{D})} + C_2 \sup_{\mathbf{x}\in\Omega} \gamma(\mathbf{x}) \sup_{\eta^h \in \mathcal{V}^h\setminus\{0\}} \dfrac{\|\nabla \eta^h/\Phi\|_{L^2(\Omega)}}{\|\nabla \eta^h\|_{L^2(\Omega)}} \cdot \|u\|_{H^1(\Omega)} + C_3 \sqrt{\sup_{\mathbf{x}\in\mathcal{D}\setminus\Omega} \gamma(\mathbf{x})}\, \|\hat{v} - \tilde{v}\|_{H^1(\mathcal{D}\setminus\overline{\Omega})},\] which is precisely the statement of Theorem 1.

9 Proof of Theorem 2↩︎

We follow the structure of the proof of Theorem 1 (Appendix 8), treating the Tucker perturbation fields \(\delta_1, \delta_2, \delta_3\) as additional stabilization-like terms. For simplicity, we consider homogeneous Dirichlet boundary condition (\(g^h = 0\)); the non-homogeneous case follows by a standard lifting argument.

Substituting the Tucker-compressed geometric coefficients \(\hat{\Phi}^{(1)}, \hat{\Phi}^{(2)}, \hat{\Phi}^{(3)}\) into the exact Dirichlet weak form 11 yields the Tucker-perturbed bilinear form and linear functional: \[\begin{align} \hat{A}_\gamma(w, v) &=& \int_{\mathcal{D}} \bigl(\hat{\Phi}^{(2)} + \gamma(\mathbf{x})\bigr)\,(\nabla w \cdot \nabla v) \,\mathrm{d}\mathcal{D} \;-\; \int_{\mathcal{D}} \hat{\Phi}^{(3)}\, w\, v \,\mathrm{d}\mathcal{D}, \tag{26} \\ \hat{L}(w, b) &=& \int_{\mathcal{D}} \hat{\Phi}^{(1)}\, w\, b \,\mathrm{d}\mathcal{D}. \tag{27} \end{align}\] Using the perturbation fields \(\delta_1 = \hat{\Phi}^{(1)} - \Phi\), \(\delta_2 = \hat{\Phi}^{(2)} - \Phi^2\), \(\delta_3 = \hat{\Phi}^{(3)} - \Phi\,\nabla^2\Phi\), the Tucker forms relate to the exact forms 25 as \[\begin{align} \hat{A}_\gamma(w, v) &=& A_\gamma(w, v) + \int_{\mathcal{D}} \delta_2\,(\nabla w \cdot \nabla v) \,\mathrm{d}\mathcal{D} \;-\; \int_{\mathcal{D}} \delta_3\, w\, v \,\mathrm{d}\mathcal{D}, \tag{28} \\ \hat{L}(w, b) &=& L(w, b) + \int_{\mathcal{D}} \delta_1\, w\, b \,\mathrm{d}\mathcal{D}. \tag{29} \end{align}\] The discrete Tucker solution \(u^h\) satisfies \(\hat{A}_\gamma(\eta^h, u^h) = \hat{L}(\eta^h, b)\) for all \(\eta^h \in \mathcal{V}^h\), while the exact solution \(u\) satisfies \(A_0(\eta, u) = L(\eta, b)\) for all \(\eta \in \mathcal{V}^h\).

Define the Tucker energy norm \(\|w\|_{\hat{\gamma}} = \sqrt{\hat{A}_\gamma(w,w)}\). From 28 , \[\bigl| \|w\|_{\hat{\gamma}}^2 - \|w\|_\gamma^2 \bigr| \;\leq\; \|\delta_2\|_{L^\infty(\mathcal{D})} \|\nabla w\|^2_{L^2(\mathcal{D})} \;+\; \|\delta_3\|_{L^\infty(\mathcal{D})} \|w\|^2_{L^2(\mathcal{D})}. \label{eq:norm95diff}\tag{30}\] Since \(\Phi \in [0,1]\) and \(w = \Phi\,\tilde{w}\) for some \(\tilde{w} \in H^1(\mathcal{D})\), the \(L^2\) term is controlled by the gradient term via the Poincaré–Friedrichs inequality. Assuming the Tucker compression is sufficiently accurate that \(\|\delta_2\|_{L^\infty}\) and \(\|\delta_3\|_{L^\infty}\) are small, the two norms are equivalent: there exist constants \(c_0, c_1 > 0\) independent of \(h\) and of the Tucker rank such that \[c_0 \|w\|_\gamma \;\leq\; \|w\|_{\hat{\gamma}} \;\leq\; c_1 \|w\|_\gamma, \qquad \forall\, w \in \mathcal{V}^h. \label{eq:norm95equiv}\tag{31}\] In the estimates below, we absorb \(c_0^{-1}\) and \(c_1\) into the generic constants \(C_2, \dots, C_5\).

Applying Strang’s first lemma (Lemma 1) to the Tucker form \(\hat{A}_\gamma\): \[\label{eq:strang95tucker} \|u - u^h\|_{\hat{\gamma}} \;\leq\; 2 \inf_{z^h \in \mathcal{S}^h} \|u - z^h\|_{\hat{\gamma}} \;+\; \sup_{\eta^h \in \mathcal{V}^h\setminus\{0\}} \frac{\bigl| \hat{A}_\gamma(\eta^h, u) - \hat{L}(\eta^h, b) \bigr|}{\|\eta^h\|_{\hat{\gamma}}}.\tag{32}\] We estimate the approximation error and the consistency error separately.

Let \(u = \Phi \hat{v}\) be the exact solution extended harmonically to \(\mathcal{D}\), and let \(\tilde{v} \in H^{p+1}(\mathcal{D})\) be the Stein extension. Decompose the approximation error: \[\begin{align} 2\inf_{z^h \in \mathcal{S}^h} \|u - z^h\|_{\hat{\gamma}} &\le 2 \|u - \Phi\tilde{v}\|_{\hat{\gamma}} + 2 \inf_{z^h} \|\Phi\tilde{v} - z^h\|_{\hat{\gamma}}. \end{align}\] For the first term, using the definition of \(\hat{A}_\gamma\) and noting that \(u - \Phi\tilde{v} = \Phi(\hat{v} - \tilde{v})\) vanishes on \(\Omega\): \[\begin{align} \|u - \Phi\tilde{v}\|_{\hat{\gamma}}^2 &= \int_{\mathcal{D}\setminus\Omega} \bigl(\Phi^2 + \delta_2 + \gamma\bigr) \bigl|\nabla(\hat{v} - \tilde{v})\bigr|^2 \,\mathrm{d}\mathcal{D} \;-\; \int_{\mathcal{D}\setminus\Omega} (\Phi\,\nabla^2\Phi + \delta_3)\, |\hat{v} - \tilde{v}|^2 \,\mathrm{d}\mathcal{D} \\ &\le \bigl(\sup_{\mathbf{x}\in\mathcal{D}\setminus\Omega} \gamma(\mathbf{x}) + \|\delta_2\|_{L^\infty(\mathcal{D}\setminus\Omega)}\bigr) \|\hat{v} - \tilde{v}\|^2_{H^1(\mathcal{D}\setminus\overline{\Omega})} \;+\; \|\delta_3\|_{L^\infty(\mathcal{D})} \|\hat{v} - \tilde{v}\|^2_{L^2(\mathcal{D})} \\ &\le \bigl(\sup_{\mathbf{x}\in\mathcal{D}\setminus\Omega} \gamma(\mathbf{x}) + \|\delta_2\|_{L^\infty(\mathcal{D}\setminus\Omega)} + \|\delta_3\|_{L^\infty(\mathcal{D})}\bigr) \|\hat{v} - \tilde{v}\|^2_{H^1(\mathcal{D}\setminus\overline{\Omega})}. \end{align}\]

For the second term, the C-HiDeNN interpolation estimate in the \(\hat{\gamma}\)-norm is identical to that in Theorem 1 up to the perturbation of the norm, which is controlled by 31 . Hence, \[2\inf_{z^h} \|\Phi\tilde{v} - z^h\|_{\hat{\gamma}} \;\leq\; C_1 h^{p}\, \|\tilde{v}\|_{H^{p+1}(\mathcal{D})}.\]

Decompose the consistency residual using 28 –29 : \[\begin{align} \hat{A}_\gamma(\eta^h, u) - \hat{L}(\eta^h, b) &= \bigl[A_\gamma(\eta^h, u) - L(\eta^h, b)\bigr] \tag{33} \\ &\quad + \int_{\mathcal{D}} \delta_2\, \nabla w^h \cdot \nabla \hat{v} \,\mathrm{d}\mathcal{D} \tag{34} \\ &\quad - \int_{\mathcal{D}} \delta_3\, w^h\, \hat{v} \,\mathrm{d}\mathcal{D} \tag{35} \\ &\quad - \int_{\mathcal{D}} \delta_1\, w^h\, \tilde{b} \,\mathrm{d}\mathcal{D}, \tag{36} \end{align}\] where \(w^h = \eta^h / \Phi\), \(\hat{v} = u / \Phi\), and \(\tilde{b} = \nabla^2 g^h + b\) (with \(b\) the body force).

9.0.0.1 Term 33 : \(\gamma\)-stabilization

From the proof of Theorem 1 ( 8), the exact stabilized consistency error satisfies \[\bigl| A_\gamma(\eta^h, u) - L(\eta^h, b) \bigr| \;\leq\; \sup_{\mathbf{x}\in\Omega} \gamma(\mathbf{x})\, \|\nabla w^h\|_{L^2(\Omega)}\, \|\nabla \hat{v}\|_{L^2(\Omega)}.\] Dividing by \(\|\eta^h\|_{\hat{\gamma}}\) and using the norm equivalence 31 together with the inequality \(\|\eta^h\|_\gamma \ge \|\nabla \eta^h\|_{L^2(\Omega)}\) proved in Lemma 1 yields the bound with \(\sup_{\Omega} \gamma(\mathbf{x})\) as in the second term of Eq.  ?? .

9.0.0.2 Term 34 : \(\delta_2\) stiffness perturbation

The perturbation \(\delta_2\) enters as a coefficient of the stiffness integral, exactly as \(\gamma\) does. Splitting the integral over \(\Omega\) and \(\mathcal{D}\setminus\Omega\): \[\int_{\mathcal{D}} \delta_2\, \nabla w^h \cdot \nabla \hat{v} \,\mathrm{d}\mathcal{D} \;=\; \int_{\Omega} \delta_2\, \nabla w^h \cdot \nabla \hat{v} \,\mathrm{d}\Omega \;+\; \int_{\mathcal{D}\setminus\Omega} \delta_2\, \nabla w^h \cdot \nabla \hat{v} \,\mathrm{d}(\mathcal{D}\setminus\Omega).\] In the interior \(\Omega\), the Cauchy–Schwarz inequality gives \[\biggl| \int_{\Omega} \delta_2\, \nabla w^h \cdot \nabla \hat{v} \,\mathrm{d}\Omega \biggr| \;\leq\; \|\delta_2\|_{L^\infty(\Omega)}\, \|\nabla w^h\|_{L^2(\Omega)}\, \|\nabla \hat{v}\|_{L^2(\Omega)}.\] In the void \(\mathcal{D}\setminus\Omega\), the harmonic extension satisfies \(\nabla^2 \hat{v} = 0\). Integrating by parts and using \(w^h|_{\partial\Omega} = 0\) (since \(w^h = \eta^h/\Phi\) and \(\eta^h|_{\partial\mathcal{D}} = 0\)), the void integral vanishes under the same argument as the \(\gamma\) term in Theorem 1: \[\int_{\mathcal{D}\setminus\Omega} \delta_2\, \nabla w^h \cdot \nabla \hat{v} \,\mathrm{d}(\mathcal{D}\setminus\Omega) = 0.\] Thus, the \(\delta_2\) contribution combines additively with \(\gamma\) in the interior bound, yielding the term \(\|\delta_2\|_{L^\infty(\Omega)}\) inside the second line of ?? .

9.0.0.3 Term 35 : \(\delta_3\) mass perturbation

By the Cauchy–Schwarz inequality, \[\begin{align} \biggl| \int_{\mathcal{D}} \delta_3\, w^h\, \hat{v} \,\mathrm{d}\mathcal{D} \biggr| &\;\leq\; \|\delta_3\|_{L^\infty(\mathcal{D})}\, \|w^h\|_{L^2(\mathcal{D})}\, \|\hat{v}\|_{L^2(\mathcal{D})} \\ &\;\leq\; \|\delta_3\|_{L^\infty(\mathcal{D})}\, \|\eta^h/\Phi\|_{L^2(\mathcal{D})}\, \|\hat{v}\|_{L^2(\mathcal{D})}. \end{align}\] Dividing by \(\|\eta^h\|_{\hat{\gamma}}\), taking the supremum over \(\eta^h\), we obtain the fourth term of ?? .

9.0.0.4 Term 36 : \(\delta_1\) right-hand side perturbation.

\[\begin{align} \biggl| \int_{\mathcal{D}} \delta_1\, w^h\, \tilde{b} \,\mathrm{d}\mathcal{D} \biggr| &\;\leq\; \|\delta_1\|_{L^\infty(\mathcal{D})}\, \|w^h\|_{L^2(\mathcal{D})}\, \|\tilde{b}\|_{L^2(\mathcal{D})} \\ &\;\leq\; \|\delta_1\|_{L^\infty(\mathcal{D})}\, \|\eta^h/\Phi\|_{L^2(\mathcal{D})}\, \|\tilde{b}\|_{L^2(\mathcal{D})}. \end{align}\]

Combining the approximation error and the four consistency error contributions into the Strang decomposition 32 and converting from \(\|\cdot\|_{\hat{\gamma}}\) to \(\|\cdot\|_\gamma\) via 31 , we obtain the error estimate given in Eq.  ?? . \(\square\)

References↩︎