arXiv CSopen access

A Non-intrusive Approach for the Imposition of Strong Dirichlet Boundary Conditions in Unfitted Boundary Meshes

· arxiv_cs

arXiv CS · Papers · License: Open Access

software-architecture software-engineering testing

software engineering, software architecture, testing

A Non-intrusive Approach for the Imposition of Strong Dirichlet Boundary Conditions in Unfitted Boundary Meshes Juan Ignacio Camarottia , Ricky Aristioa , Riccardo Rossib,c , Rubén Zorrillab,c , Roland Wüchnera a Chair of Structural Analysis, Technical University of Munich, Arcisstr. 21, 80333 München, Germany b Departament d’Enginyeria Civil i Ambiental, Universitat Politècnica de Catalunya, Barcelona, 08034, Spain

arXiv:2609.15398v1 [cs.CE] 14 Sep 2026

c International Center for Numerical Methods in Engineering (CIMNE), Barcelona, 08034, Spain

Abstract The enforcement of essential boundary conditions is a fundamental challenge in unfitted boundary methods. This paper presents a non-intrusive, black-box strategy for imposing such conditions in unfitted meshes. The approach is intended for situations where the user does not have access to the solver’s source code or its mathematical formulation, which is often the case when using commercial software. The proposed algorithm allows solvers originally designed for body-fitted meshes to be used in unfitted cases, provided that four conditions are satisfied: (i) the solver must support user customization by means of scripting, (ii) allow the imposition of Dirichlet boundary conditions at the node level through scripting, (iii) permit the deactivation of elements outside the physical domain, and (iv) provide access to the solution gradient within active elements. The last condition can also be satisfied by externally reconstructing the gradient from nodal values and connectivity information, provided the element formulation is known, making it optional in practice. These requirements are very fair demands and are satisfied by the vast majority of productionready, possibly commercial, codes. In the current work, we show the application of this non-intrusive algorithm in the context of the Finite Element Method (FEM) and Isogeometric Analysis (IGA) discretizations, demonstrating optimal L2 -norm error convergence. This is demonstrated using the Kratos Multiphysics code (release v10.1) from the user API, simply leveraging the capabilities mentioned above. Keywords: Unfitted boundary methods, Black-box solver, Strong dirichlet boundary conditions, Finite element method (FEM), Isogeometric Analysis (IGA), Trimmed domain

1. Introduction Unfitted or embedded boundary methods have gained significant importance in computational mechanics, offering an efficient alternative to body-fitted (Fig. 1a) approaches for handling complex geometries as well as for arbitrarily large boundary motions. These methods embed the physical boundary within a structured background mesh (Fig. 1b), removing the need for conformal meshing. This characteristic is particularly advantageous in scenarios involving evolving geometries such as crack propagation [1, 2, 3], shape optimization [4, 5, 6], fluid-structure interaction [7, 8, 9], additive manufacturing [10], and biomedical applications [11, 12, 13, 14, 15, 16, 17], where re-meshing introduces substantial computational costs. Despite their geometric flexibility, unfitted methods present significant challenges related to the integration over cut elements, the enforcement of boundary conditions along immersed interfaces, and the stability of the formulation in the presence of small-cut cells [18]. Among these, the imposition of essential (Dirichlet) boundary conditions is particularly problematic, as the lack of alignment between mesh nodes and the physical boundary prevents direct enforcement. Variational approaches, such as Nitsche’s method [19, 20], Lagrange multipliers [21, 22], and penalty formulations [23, 24], impose these conditions weakly but often require careful tuning of penalty parameters and ∗ Corresponding author

Email addresses: [email protected] (Juan Ignacio Camarotti), [email protected] (Ricky Aristio), [email protected] (Riccardo Rossi), [email protected] (Rubén Zorrilla), [email protected] (Roland Wüchner) Preprint submitted to Engineering with Computers

September 15, 2026

may lead to ill-conditioned systems. Alternatively, some methods aim for strong enforcement by modifying basis functions so they become interpolatory along the boundary, as in WEB-splines [25, 26] and i-splines [27]. While these techniques avoid modifying the weak form, they introduce significant algorithmic and analytical complexities due to the need for constructing specialized basis functions. Notably, none of the methods reviewed in [18] support a black-box strategy for Dirichlet boundary condition enforcement, an increasingly important limitation in scenarios where users lack access to the solver internals. Several prominent embedded boundaries methodologies have been developed, including the Finite Cell Method (FCM) [28, 29, 30, 31, 32, 33], the extended Finite Element Method (XFEM) [34], the Cut Finite Element Method (CutFEM) [7, 35, 36, 37], Immersed Boundary Method (IBM)[16, 38], Immersogeometric Analysis (IMGA) [14], Isogeometric B-Rep Analysis (IBRA)[39, 40, 41, 42, 43], the Shifted Boundary Method (SBM) [44, 45, 46, 47, 48], each addressing specific challenges in accurate integration schemes, enforcement of boundary conditions and maintaining numerical stability. Each of the existing methodologies approaches the challenges mentioned above differently, but a common numerical difficulty across most immersed methods is the so-called small cut cell problem, which occurs when only a small fraction of an element lies inside the physical domain. This situation can lead to severe ill-conditioning of the system matrix, and different methods address it through different stabilization strategies. The Finite Cell Method (FCM) [28, 29] embeds the physical domain into a larger, structured computational mesh and distinguishes physical and fictitious regions using an indicator function. It enables high-order analysis on non-boundary-fitted meshes and typically enforces boundary conditions weakly (e.g., via Nitsche’s method). Accurate integration near the boundary is achieved through advanced quadrature schemes, but the method can suffer from ill-conditioning in cells with small physical volume fractions. The extended Finite Element Method (XFEM) enriches the standard finite element space with additional functions to capture discontinuities introduced by immersed boundaries, such as jumps or singularities. While this allows for accurate representation of complex features without mesh conformity, the method introduces extra degrees of freedom (DOFs), increasing computational cost. In both FCM and XFEM, the small cut cell problem can be effectively alleviated using the Eigenvalue Stabilisation Technique pioneered by Löhnert [49], which regularizes ill-conditioned element stiffness matrices via eigenmode filtering [50, 51]. The family of Immersed Boundary Methods (IBM), originally developed for fluid dynamics [38, 16], represents boundaries within a fixed Cartesian grid and applies boundary conditions via additional body forces in the governing equations. While this method simplifies mesh generation and allows for large deformations, it often suffers from mass conservation issues and difficulties in accurately enforcing boundary conditions at the interface. IBRA, on the other hand, employs the CAD model’s boundary representation (B-Rep) along with isogeometric basis functions to approximate solution fields, reducing pre-processing time but suffering from computational challenges associated with integrating trimmed elements and conditioning issues due to small cut cells. The Shifted Boundary Method (SBM) addresses boundary enforcement by shifting the physical boundary onto a surrogate location while modifying the imposed boundary conditions, effectively bypassing small cut-cell instabilities but introducing other challenges, such as the non-trivial modification of boundary conditions and a geometrically non-obvious treatment of the surrogate boundary [18]. CutFEM classifies elements as either fully inside the physical domain or intersected by the boundary. For cut elements, it employs sub-cell integration to accurately compute contributions over the physical region. Boundary conditions are enforced variationally, often using Nitsche’s method, and the method incorporates ghost penalty stabilization [52] to address issues related to small cut elements and ensure numerical stability. Despite this rich landscape of methods, explicit (strong-like) imposition of Dirichlet boundary conditions in unfitted meshes, without modifying the weak form or the basis functions used for analysis, has received limited attention, particularly in black-box settings where the solver’s internals are inaccessible. Motivated by this gap, we propose a physics-agnostic, non-intrusive iterative approach that strongly enforces Dirichlet boundary conditions in unfitted meshes without requiring access to the solver’s formulation or source code. This is especially relevant in modern computational practice, where black-box solvers such as commercial software packages are widely used and cannot be altered by the end user. The proposed method iteratively refines the BC imposition by minimizing the L2 -norm error between the numerical and prescribed boundary values. The stabilization mechanism ensures smooth normal gradient transitions within intersected and fully active elements. This is achieved through interpolation techniques such as Radial Basis Functions (RBF) and Moving Least Squares (MLS), which provide a robust approximation of gradient fields. Our approach presents several potential benefits over existing unfitted methods. Its physics-agnostic and non2

intrusive nature ensures broad applicability across various engineering fields, seamlessly integrating with existing black-box body-fitted solvers without requiring fundamental modifications. Furthermore, it eliminates the reliance on penalty parameters, thereby removing the need for empirical tuning, and enhances system conditioning, mitigating numerical instabilities associated with small cut cells. Nevertheless, some challenges persist, including potential less accuracy in the solution field due to the gradient approximation (especially in high-order geometries), increased computational cost due to the iterative enforcement process, and a dependency on interpolation schemes, which must be carefully chosen to achieve an optimal balance between accuracy and stability. The remainder of the paper is structured as follows: In Section 2, we introduce the governing equations and boundary conditions of the Poisson problem, which serves as the PDE for validating the method, and review its weak formulation. In Section 3, we present the concept of enforcing strong Dirichlet boundary conditions as an L2 -error minimization problem. Furthermore, we show that applying this approach to unfitted meshes results in a singular minimization problem. In Section 4, we show that this singularity issue can be mitigated by incorporating an additional stabilization term into the original functional, which is based on the solution normal gradient within the trimmed elements. In addition, we outline the resulting solution algorithm and describe the classification of elements within this framework. Different interpolation techniques for approximating the normal gradient within trimmed elements, such as Radial Basis Functions (RBF) and Moving Least Squares (MLS), are proposed in Section 5. In Section 6, we assess the validity of the proposed approach by numerically solving the Poisson problem using FEM and IGA discretizations. In addition, we demonstrate how the proposed method can be integrated into an existing black-box solver by implementing it in the Kratos Multiphysics [53, 54] framework (release v10.1) without modifying the underlying solver algorithms. This practical implementation confirms that the method can be deployed non-intrusively on top of production-ready, body-fitted solvers. Finally, in Section 7, we provide concluding remarks and discuss possible future applications of this method.

ΓD ΓD

Ω

y x

(a) Body-fitted triangular mesh

(b) Unfitted cartesian mesh

Figure 1: Comparison between body-fitted and unfitted meshes for a circular domain. (a) The body-fitted case, where the mesh conforms to the circular boundary using triangular elements. (b) The unfitted case, where the circular domain is embedded in a structured background mesh, allowing for simpler meshing but requiring special techniques for boundary condition enforcement.

2. Governing equations of Poisson’s problem In this paper, we consider the numerical solution of the Poisson equation. Let Ω ⊂ Rn be a bounded computational domain, where n = 2, 3 represents the spatial dimension. The domain boundary, ∂Ω = Γ, consists of two disjoint subsets: the Dirichlet boundary ΓD and the Neumann boundary ΓN , satisfying ΓD ∪ ΓN = Γ and ΓD ∩ ΓN = ∅. The

strong form of the Poisson equation is given by: −∇2 ϕ = σ

in Ω,

(1a)

ϕ = ϕD ∂ϕ ∇ϕ · n = = tn ∂n

on ΓD ,

(1b)

on ΓN .

(1c)

Here, ϕ is the primary variable, ϕD represents the prescribed value on the Dirichlet boundary ΓD , tN denotes the imposed normal gradient on the Neumann boundary ΓN , and σ acts as the forcing or source term. When the forcing term vanishes, i.e., σ = 0 throughout the entire domain Ω, the equation simplifies to the Laplace equation: ∇2 ϕ = 0

in Ω.

(2)

Finally, in Eq. 1c, n denotes the outward-pointing normal vector to the boundary Γ. 2.1. Weak form of the Poisson’s problem The Sobolev space H 1 (Ω) is a fundamental function space in the weak formulation of partial differential equations. To define it, we first introduce the space of square-integrable functions, denoted by L2 (Ω), which is defined as ( ) Z 2 2 L (Ω) := v : Ω → R | |v| dΩ < ∞ (3) Ω

This means that a function v belongs to L (Ω) if its square is integrable over the domain Ω. Building on this, the Sobolev space H 1 (Ω) extends L2 (Ω) by also requiring the first-order weak derivatives of v to be in L2 (Ω). Formally, it is defined as n o H 1 (Ω) := v ∈ L2 (Ω) | ∇v ∈ L2 (Ω) (4) 2

where ∇v represents the weak gradient of v. This ensures that functions in H 1 (Ω) not only belong to L2 (Ω) but also have weak derivatives that are square-integrable, providing a certain degree of smoothness. We introduce a special subspace of H 1 (Ω), known as H01 (Ω). This space is defined as n o H01 (Ω) := v ∈ H 1 (Ω) | v = 0 on ΓD (5) This definition implies that H01 (Ω) consists of functions in H 1 (Ω) that vanish on the Dirichlet boundary ΓD . It serves as the natural function space for test functions in the weak formulation when Dirichlet boundary conditions are imposed. Introducing the notation Z Z (α, β)Ω =

αβ dΩ

Ω

⟨α, β⟩Γ =

αβ dΓ

(6)

and applying Green’s theorem allows to write the weak formulation for the Poisson problem as: find ϕ ∈ H 1 (Ω) such that for all w ∈ H01 (Ω), (∇ϕ, ∇w)Ω = (σ, w)Ω + ⟨tN , w⟩ΓN (7) 3. Imposition of strong-like Dirichlet BCs as an L2 -norm error minimization problem For this section, we refer to the Poisson problem introduced in Section 2, defined on a domain Ω ⊂ Rn with boundary ∂Ω = ΓD ∪ ΓN , where ΓD and ΓN denote the Dirichlet and Neumann boundaries, respectively. The Dirichlet boundary condition ϕ = ϕD on ΓD can be interpreted as the minimization of the L2 -norm error between the field ϕ and the imposed value ϕD . The corresponding functional representing the error is defined as Z 1 (ϕ(x) − ϕD (x))2 dΓ ψ(ϕ) = (8) 2 Γ 4

This equation defines a least-squares fit problem and constitutes a projection of the Dirichlet boundary condition onto the discretized function space. It serves as the foundation for the auxiliary algebraic system used to impose the boundary conditions in our method. In a FEM-like approach, the solution is approximated by ϕh (x) =

n X

N j (x)ϕ j = N(x)ϕ

(9)

j=1

where N j (x) are the finite element shape functions, ϕ j are the corresponding nodal values and n represents the number of nodes (or control points in IGA). Substituting the approximation 9 into the error functional 8 yields 1 ψ(ϕ) = 2

Z ΓD

2  nΓ D  X  N j (x)ϕ j − ϕD (x) dΓ.

(10)

j=1

Remark 1. In Eq. 9, the index n refers to the total number of nodes (or control points) in the discretized domain. However, only those nodes whose associated shape functions satisfy N j (x) , 0 on the Dirichlet boundary ΓD actually contribute to the integral in Eq. 8. As a result, the effective linear system involves only a subset of these nodes. We denote this reduced number by nΓD . To determine the degrees of freedom ϕi , we take the derivative of ψ(ϕ) with respect to each ϕi and set it equal to zero:  nΓ Z X   D ∂ψ  (11) = N j (x)ϕ j − ϕD (x) Ni (x) dΓ = 0 i = 1, . . . , nΓD . ∂ϕi ΓD

j=1

This set of equations leads to the linear system nΓ D Z X j=1

ΓD

! Z Ni (x)N j (x) dΓ ϕ j =

ΓD

ϕD (x)Ni (x) dΓ,

i = 1, . . . , nΓD .

(12)

In matrix form, the linear system can be written as MΓD ϕ = b, where the mass matrix MΓD and the right-hand side vector b have entries defined by Z MΓ D i j = Ni (x)N j (x) dΓ, ΓD Z bi = ϕD (x)Ni (x) dΓ. ΓD

(13)

(14) (15)

This system is solved independently of the weak form of the original PDE and is used solely to compute boundary values for subsequent enforcement. In the following part we will demonstrate that, for unfitted discretizations, the resulting linear system is ill-posed, leading to infinite solutions. 3.1. Singularity problem for unfitted discretizations In this section, we analyze a simplified example involving two different FEM discretizations of a 1D bar structure using linear elements (Fig. 2): a body-fitted discretization (Fig. 3) and an unfitted discretization (Fig. 4). Our objective is to impose a Dirichlet boundary condition at the left end of the bar, denoted by ΓD . In the body-fitted discretization, ΓD coincides with the left-end of the first element, whereas in the unfitted case, it is located at the midpoint of the first element. 5

N1 (ξ) = 1−ξ 2

N2 (ξ) = 1+ξ 2

ξ −1

Figure 2: 2-node linear element in natural coordinates space [−1, 1] with corresponding shape functions.

u(x) = αx,

where α ∈ R

uΓD = 0 e1

Ω 8

Node Figure 3: Body-fitted discretization of a bar structure subjected to an external force F at its free end. The bar is fixed at the left support and undergoes a displacement u(x) due to the applied load. The discretization is achieved using a structured mesh, where boundary nodal points (marked in red) conform to the geometry of the bar. The horizontal axis represents the spatial coordinate x, while the vertical axis shows the displacement field u(x) for illustrative purposes. u

u(x) = αx + β,

uΓ D = 0 e3

7 Ω

where α, β ∈ R;

Node Figure 4: Unfitted discretization of a bar structure subjected to an external force F at its free end. The bar is fixed at the left support and undergoes a displacement u(x) due to the applied load. The discretization is achieved using a structured mesh, where boundary nodal points (marked in red) do not conform to the geometry of the bar. The horizontal axis represents the spatial coordinate x, while the vertical axis shows the displacement field u(x) for illustrative purposes.

The governing equation for the longitudinal displacement u(x) is given by ! d du EA = −p(x), in Ω, dx dx

(16)

subject to the boundary conditions u = uD

on ΓD ,

(17)

du EA = F on ΓN . (18) dx In Eq. 16, EA denotes the axial stiffness of the bar, p(x) represents a longitudinally distributed load, and uD corresponds to the imposed displacement. For this particular case, we assume p(x) = 0 and uD = 0. We now analyze the resulting linear system used for the imposition of the Dirichlet boundary condition (Eq. 13), considering both body-fitted and unfitted discretizations. Both in the body-fitted and unfitted cases, the system is given by      N1 (x = 0)N1 (x = 0) N1 (x = 0)N2 (x = 0) u1  uD (x = 0)N1 (x = 0) (19)     =   N2 (x = 0)N1 (x = 0) N2 (x = 0)N2 (x = 0) u2 uD (x = 0)N2 (x = 0)

Substituting the shape function values and the imposed displacement at the ΓD boundary for the body-fitted case, we obtain      1 0 u1  0    =   ⇒ u1 = 0, u2 ∈ R  (20) 0 0 0 u2 Thus, in the body-fitted case, the Dirichlet boundary condition is imposed strongly on the corresponding degree of freedom, enforcing u1 = 0. The remaining unknown, u2 , is not constrained by the boundary condition but will be determined by solving the system of equations arising from the weak form of the governing equation. Conversely, in the unfitted case, substituting the shape function values and the imposed displacement at the ΓD boundary leads to the following linear system  1 2    1 2     2  u1  0 2     (21) =      ⇒ Infinite solutions. 1 2 1 2  u  0 2 2 2 Since this system of equations is underdetermined, it admits infinitely many solutions. In other words, attempting to project the Dirichlet boundary condition onto the discretization space of an unfitted boundary leads to an ill-posed problem. 3.2. Physical interpretation of the singularity problem for unfitted discretizations Having presented the body-fitted and unfitted bar structure examples, we now turn our attention to the physical interpretation of this singularity problem observed in unfitted discretizations. Focusing on the first element of the bar structure, we can observe from Fig. 5 that there are infinite possible assignments for u1 and u2 that still satisfy the condition uΓD = 0. As observed in Fig. 5, the infinite solutions differ only in the gradient of the solution field within the trimmed element. Moreover, it is crucial to note that many of these solutions, such as those with a negative gradient inside the trimmed element, do not correspond with the physical behaviour of the problem. This raises a fundamental question whether this issue could be resolved by imposing an additional condition on the gradient within the first element (the trimmed element). This concept serves as the foundation for the methodology proposed in this work. Different du dx

u 13

u24

uΓD = 0

u 12

u21 e1

x Ω continues →

u 11

u22

u 14

u23 Node

Figure 5: Visualization of an unfitted discretization for a bar structure, illustrating how different displacement gradients ( du dx ) emerge from various combinations of nodal displacements in the unfitted boundary that satisfy the Dirichlet boundary condition uΓD = 0. The dashed lines represent possible displacement solutions u(x), and the domain Ω continues to the right, as indicated by the gray dashed line. Nodes are marked in red, and the structured mesh elements are annotated as e1 and e2 .

4. Stabilization of the L2 -norm error functional for unfitted meshes In this section, we present a novel approach to stabilize the L2 -norm error functional for unfitted meshes, enabling a strong-like, external imposition of Dirichlet boundary conditions without altering the underlying weak form of the 7

partial differential equation. This approach results in a non-intrusive iterative method suitable for body-fitted and unfitted discretizations, effectively addressing the challenges associated with trimmed elements. Before presenting the method, we clarify the concept of strong-like imposition of Dirichlet boundary conditions as used in this work. By strong-like, we refer to the explicit modification of both the left-hand side (LHS) and right-hand side (RHS) of the linear system arising from the weak form of the physical problem to enforce the boundary condition. The nodal values used for this enforcement are computed by solving an auxiliary algebraic system that is independent of the weak form governing the physics. We introduce a novel functional ψ(ϕ) to define strong-like Dirichlet boundary conditions in unfitted meshes. This functional consists of two key components and is given by Z Z 2 2 1 1 ϕk+1 (x) − ϕD (x) dΓ + ∇(ϕk+1 ) · n − e ∇(ϕk ) · n dΩ (22) ψ(ϕk+1 ) = 2 ΓD 2 Ωtrim In Eq. 22, ϕk+1 represents the solution at the Dirichlet boundary at iteration k + 1, Ωtrim denotes the intersected elements (or knot spans in the case of Isogeometric Analysis, ∇(ϕk+1 ) is the gradient inside an intersected element at iteration k+1, and e ∇(ϕk ) is an approximation of the gradient inside the intersected element, obtained using information computed at iteration k. The second term in Eq. 22 acts as a regularization mechanism that addresses the underdetermined nature of the projection problem in unfitted discretizations. This stabilization term introduces an additional constraint to the projection that guides the gradient within the trimmed region to be compatible with the gradient field in adjacent, fully-resolved elements. The term involves the normal component of the gradient, where e ∇(ϕk ) is an approximation computed at iteration k using information from neighboring non-intersected elements. Notably, the term is integrated over the entire trimmed element (not just the active portion) and has shown to produce accurate and stable results in practice. Remark 2. To approximate the gradient inside intersected elements, e ∇(ϕ), the gradient in non-intersected elements must first be determined. However, since the gradient in these elements depends on the nodal values, unknown prior to imposing boundary conditions, the procedure is inherently iterative. Notably, if the exact gradient at the boundary were known, the algorithm would converge in a single iteration. Remark 3. It is important to mention that at the first iteration (k = 0), the normal gradient within the intersected elements is initialized as ∇(ϕ0 ) · n = 0 (23) Now, we would like to introduce the new discrete form of the problem based on the proposed functional in Eq. 22. Taking the derivative of Eq. 22 with respect to each ϕk+1 and setting it to zero yields i     Z X Z n  X   n ∂ψ k+1 k  Ni (x) dΓ +  (∇N j (x) · n) ϕk+1 − e     (∇Ni (x) · n) dΩ = 0, = N (x)ϕ − ϕ (x) ∇(ϕ ) · n  j D j j    ∂ϕk+1 (24) ΓD j=1 Ωtrim j=1 i i = 1, . . . , n. In matrix form, the algebraic system can be expressed as a residual equation Aϕk+1 − f k = 0, where the matrix A and the right-hand-side vector f are defined as Z Z Ai j = Ni (x)N j (x) dΓ + (∇Ni (x) · n)(∇N j (x) · n) dΩ = (Ai j )γ + (Ai j )Ωtrim , fik =

ΓD

Ωtrim

ΓD

ϕD (x)Ni (x) dΓ +

Ωtrim

(e ∇(ϕk ) · n)(∇Ni (x) · n) dΩ = ( fi )γ + ( fik )Ωtrim

(25)

(26) (27)

This formulation extends the classical discrete system in Eq. 13 by incorporating a regularization (or augmentation) term within the intersected elements, thereby ensuring a well-posed problem. Notably, each individual contribution to the left-hand side is inherently ill-posed; only their combined effect yields a stable and solvable linear system. 8

Remark 4. It is important to emphasize that Equations 26 and 27 are not incorporated into the weak form of the physical problem. Instead, Equation 25 defines an independent algebraic problem that is solved iteratively to determine the boundary conditions to be imposed on the nodes of the cut elements. This procedure remains decoupled from the weak formulation given in Equation 7. Remark 5. In Equations 26 and 27, only the contribution ( fik )Ωtrim , associated with the intersected (or trimmed) elements needs to be updated at each fixed-point iteration. The boundary-related terms are computed only once at the beginning of the time step, as they remain constant throughout the iterative process.

Remark 6. Although the focus of this work is on the strong-like imposition of Dirichlet boundary conditions in unfitted meshes, the proposed framework can accommodate Neumann conditions in a fully black-box fashion. Given the weak form contribution Z tN w dΓ, (28) ΓN

the associated force vector is computed externally, using the geometry of the Neumann boundary, and then passed to the solver as an additional forcing term in its global right-hand side vector. The stiffness matrix assembled by the solver is left completely unchanged. For each cut element K, a boundary quadrature is constructed on ΓN ∩ K, the outward unit normal is computed from the exact geometry, and the contribution Z fN,i = tN Ni (x) dΓ (29) ΓN

is evaluated externally, where tN is the prescribed Neumann traction, Ni (x) is the i-th basis function of the host solver, and fN,i is the corresponding entry in the global load vector. The vector fN = { fN,i } is then added to the solver’s existing load vector via its public API. At this stage, it would be highly beneficial for the reader to re-examine a simple example based on the unfitted bar structure shown in Fig. 4 to better understand the iterative algorithm. This illustrative example is presented in Fig. 6, where we attempt to enforce a zero-displacement (u = 0) boundary condition at the midpoint of the first element. In the first iteration (green dashed line), the normal gradient inside the trimmed element is initially set to zero, leading the solution of Eq. 25 to yield u1 = u2 = 0. In the second iteration (cyan dashed line), information from the previous step regarding the normal gradients within the domain is incorporated, enabling a correction of the boundary condition values by approximating the normal gradient inside the trimmed element. The algorithm continues iterating until a specified absolute or relative tolerance is met. u(x) = αx + β, y

where α, β ∈ R

on 2

ti Itera

1 tion

Itera

F 1 e1

8 Ω

Figure 6: Illustration of the iterative correction process for enforcing a zero-displacement (u = 0) boundary condition at the midpoint of the first element in an unfitted bar structure. The bar, fixed at the left end, is subjected to an external force F at its free end. The structured mesh discretization is represented by nodal points (marked in red).

To further clarify the discussion, we now outline the general solution algorithm (Alg. 1) used to non-intrusively impose Dirichlet boundary conditions on both body-fitted and unfitted discretizations. The proposed method follows an iterative procedure that is applied at each time step for transient problems. First, an auxiliary problem is solved to 9

determine the nodal values on the Dirichlet boundary. These values are then imposed explicitly as boundary conditions in the main equations system. The physical problem is subsequently solved for the remaining degrees of freedom, the relative error between iterations is evaluated, and the process is repeated until convergence is achieved. Algorithm 1: Algorithm for the strong imposition of Dirichlet BCs in unfitted meshes Data: kmax : Max. iterations tol : Residual tolerance T : Simulation end time ∆t : Time step increment M : Background mesh EMB : Embedded geometry // Classify the elements (or knot spans) in the background mesh as active, inactive or intersected (note that this is a preprocessing step) ClassifyElements(M, EMB)→ A, I ; // A: subset of active elements, I: subset of intersected elements while t ≤ T do // Assemble the system matrices for the physical problem, considering only the contribution from active elements AssembleSystemMatrices(A) // Non-linear solution strategy loop while k ≤ kmax do // Solve the auxiliary algebraic problem to compute nodal boundary values in the Dirichlet boundary ΓD SolveAuxiliaryProblem(A, I) ; // Assemble and solve Eq.25 // Enforce the computed boundary values explicitly by altering the system matrix and right-hand side ModifySystemWithComputedBoundaryValues(I); // Solve for the remaining DOFs of the original problem with updated BCs Solve(A); // Compute the relative error between iterations erel ComputeRelativeErrorBetweenIterations(A, I); // Check convergence if erel ≤ tol then break; else k += 1 end end // Advance in time t += ∆t; end

The most critical preprocessing step in the algorithm is the elements classification, which determines their role in the computation and assembly of the global system. Elements are categorized into three types: active elements, intersected elements, and inactive elements. • Active elements: These elements are fully assembled into the global system of equations, contributing to both the left-hand side (LHS) and right-hand side (RHS) of the system. • Intersected elements: These elements participate in the auxiliary algebraic problem required to compute the coefficients ϕΓDi , but they are not included in the assembly of the global system. This exclusion arises because the nodal values of the active nodes in these elements are explicitly imposed. 10

• Inactive elements: These elements do not contribute to the global system and are excluded from the assembly process. Figure (7) illustrates an example of element classification applied to a plate with a hole.

ΓD

Ω y x Figure 7: Classification of elements within the Cartesian grid: Active (inside), Intersected (crossing the circle), and Inactive (outside).

5. Normal gradient approximation inside trimmed elements The key aspect of the proposed iterative algorithm for the strong-like imposition of Dirichlet BCs is the gradient approximation within the intersected elements. In this section, we present the fundamental principles of this approximation. As previously mentioned, the gradient within intersected elements is approximated using the computed gradients from neighboring non-intersected elements. In finite element discretizations, interpolation points are typically selected at the element centers of the N closest unperturbed elements (Def.1). Conversely, in isogeometric discretizations, they are chosen as the N nearest integration points within unperturbed knot spans. The number of interpolation points, N, is specified by the user and can be adapted based on the discretization density or problem setup. Definition 1. An unperturbed element or knot span is an element in which no degrees of freedom are influenced or modified by the auxiliary algebraic problem 25. Fig. 8 illustrates the concept of an unperturbed element within a Cartesian finite element background mesh. In this case, the physical domain consists of a square plate with a hole in the center. Through extensive numerical experiments, it has been observed that for the algorithm to achieve convergence, the gradient must be approximated using information exclusively from these unperturbed elements. These tests confirmed that including data from perturbed elements introduces inaccuracies that hinder the iterative process, reinforcing the necessity of relying solely on unperturbed regions for stable and accurate gradient reconstruction. Remark 7. Gradient reconstruction inside intersected elements should be performed exclusively using information from unperturbed elements, as contributions from perturbed elements may hinder convergence

ΓD

Ω y x Figure 8: Illustration of the unperturbed element concept: The integration point where the gradient is approximated is shown in purple, while the unperturbed elements used for the gradient approximation are highlighted in blue.

5.1. Radial basis function interpolation with polynomial extension Radial Basis Function (RBF) interpolation is a well-established technique for scattered data approximation and field transfer, widely used in computational mechanics. In immersed and embedded settings, it has been employed, for example, to transfer solution fields in the Finite Cell Method (FCM) after remeshing [55]. In the present work, RBF interpolation is adopted to approximate the gradient within intersected elements, using gradient information from the surrounding unperturbed active elements. Given N data points xi , i = 1, . . . , N, with known gradient values ∇ϕi , the RBF interpolation approximates the gradient field ∇ϕ(x) as Np N X X wi β(∥x − xi ∥) + λk qk (x), (30) ∇ϕ(x) = i=1

k=1

PN where β(r) is a chosen radial basis function, wi are the RBF weights, and p(x) = i=1p λ j q j (x) is an optional polynomial term used to improve approximation accuracy and ensure the reproduction of polynomial fields up to a given degree. Two conditions are imposed to determine the unknowns: 1. Interpolation condition: At each known point x j , N X

wi β(∥x j − xi ∥) +

i=1

Np X

λk qk (xj ) = ∇ϕ(x j ).

(31)

k=1

2. Polynomial orthogonality condition: The RBF component must be orthogonal to the chosen polynomial space, N X wi q(xi ) = 0 for all basis functions q(x) of p(x), (32) i=1

which ensures that the polynomial term is exactly reproduced, removes rank deficiencies in the system, and guarantees uniqueness of the solution. Combining these two conditions yields the augmented linear system " #" # " # B P w g = , PT 0 λ 0 12

(33)

where B is the interpolation matrix with entries B ji = β(∥x j − xi ∥), P is the polynomial matrix with rows p(x j ), w = (w1 , . . . , wN )T contains the RBF weights, λ are the polynomial coefficients, and g = (∇ϕ(x1 ), . . . , ∇ϕ(xN ))T contains the known gradients. The choice of β(r) influences accuracy and stability. Common options include: 2

• Gaussian: β(r) = e−(ϵr) , • Multiquadric (MQ): β(r) =

√

r2 + ϵ 2 ,

√ • Inverse Multiquadric (IMQ): β(r) = 1/ r2 + ϵ 2 , • Thin-Plate Spline (TPS): β(r) = r2 log r, • Cubic: β(r) = r3 . It is important to note that the shape parameter ϵ governs the flatness or sharpness of the radial basis function, directly influencing both the interpolation accuracy and the numerical stability of the solution. However, selecting an appropriate value for ϵ typically requires careful fine-tuning, which we aim to avoid in the present method. This motivates the investigation of alternative interpolation techniques that require no parameter tuning, such as the Moving Least Squares (MLS) method, which will be discussed later in this work. For an in-depth discussion of RBF interpolation and its polynomial extension, the reader is referred to [56] and [57]. 5.2. MLS interpolation The Moving Least Squares (MLS) method is a well-known numerical technique used for function approximation, surface reconstruction, and meshless methods in computational mechanics. Unlike traditional least squares approximation, MLS provides a localized and adaptive approach to fitting a function to scattered data points. The primary goal of MLS is to construct a smooth function approximation P(x) that closely follows a given set of data points while being flexible enough to adapt to local variations. This is achieved by minimizing a weighted least squares error that emphasizes nearby points more than distant ones. Given N a set of scattered data points xi , i = 1, . . . , N, with known gradient values ∇ϕi , the MLS method constructs a local polynomial approximation at each evaluation point x. The polynomial is defined as P(x) =

m X

a j (x)β j (x) = a(x) · β(x)

(34)

j=0

where β j (x) are basis functions (typically polynomials), and a j (x) are the unknown coefficients. The coefficients a j (x) are obtained by minimizing the weighted least squares error J(a) =

1X W(x − xi ) (P(xi ) − ∇ϕi )2 2 i

(35)

W(x − xi ) is a weight or kernel function introduced to ensure that points closer to x have greater influence. A common choice is the Gaussian weight function ! |x − xi |2 W(x − xi ) = exp − (36) h2 where h is a smoothing parameter controlling the influence range of each point. This leads to a system of linear equations, which is solved to determine the optimal values of a j (x). The system is given by M(x)a(x) = H(x) (37) being M(x) = BT (x)W(x)B(x) 13

(38)

and

  P  N W ∇ϕ  i i     i=1  N   P T H(x) = B (x)W(x)∇ϕ =  xi Wi ∇ϕi   i=1   P  N y W ∇ϕ  i

(39)

i=1

Here, M(x) is the weighted moment matrix and H(x) is defined as BT (x)W(x)∇Φ, being B(x) the matrix of basis functions and ∇Φ a vector containing the function values to be interpolated. The term W(x) is a diagonal weighting matrix, with diagonal entries Wii = W(x − xi ). Once the coefficients are determined, the final approximation function is given by: ∇ϕ(x) ≈ P(x)

(40)

which smoothly adapts to local data variations. For a comprehensive mathematical review of MLS interpolation, the reader is referred to [58] and [59]. 6. Application of the method to the numerical solution of the Poisson problem for FEM and IGA discretizations In this section, we assess the proposed algorithm for the Poisson problem 1a across various scenarios, considering different discretization methods (low- and high-order FEM basis functions, as well as B-Spline basis), different element geometries (triangular and quadrilateral elements), and different gradient approximation techniques within trimmed elements. The error convergence is evaluated using the following manufactured solution for the Poisson problem: ϕ(x, y) = sin(πx) cos(πy) For the given manufactured solution, the source term σ(x, y) is defined as: σ(x, y) = −2π2 sin(πx) cos(πy) The domain consists of a [1 × 1] square with a circular hole centered at C = (0.5, 0.5) and a radius of r = 0.25 (Fig. 9).We impose Dirichlet boundary conditions across the entire boundary, ΓD = Γ, by explicitly enforcing the manufactured solution. When the boundary is embedded (trimmed), the proposed algorithm strongly enforces these conditions. For body-fitted boundaries, where the geometry aligns with the underlying discretization, the same enforcement approach is used. The error is measured using the L2 -norm of the difference between the analytical solution and the numerical approximation: sZ (ϕ − ϕh )2 dΩ ∥ϕ − ϕh ∥L2 (Ω) = Ω

The proposed strategy was implemented using the Kratos Multiphysics API, based on the release version v10.1 of the open-source framework Kratos Multiphysics [53, 54]. This implementation demonstrates that it is possible to adapt existing body-fitted solvers to operate in unfitted mesh scenarios, provided that the solver satisfies four key conditions: (i) it supports user scripting, (ii) it allows Dirichlet boundary conditions to be imposed at the node level through scripting, (iii) it permits deactivation of elements outside the physical domain, and (iv) it provides access to the solution gradient within active elements. In Kratos, the structure and flow of a simulation are typically managed by a class called AnalysisStage. This class governs the full simulation lifecycle, including model initialization, time stepping, and solution loop execution. The proposed algorithm was implemented by extending AnalysisStage, overriding only the necessary member functions of this class to inject the Dirichlet enforcement step at the appropriate stage of the loop. To illustrate this, we present in Listing 1 a minimal but complete implementation of this algorithm using the Kratos API. The class shown integrates seamlessly with Kratos’ solver infrastructure and provides a reusable mechanism for solving unfitted problems without modifying the internals of the solver. 14

1 2 3 4 5

# --- Import Kratos core and the convection - diffusion solver wrapper --import Kr a to s Mu l t i p h y s i c s from K r at os M ul t i p h y s i c s . C o n v e c t i o n D i f f u s i o n A p p l i c a t i o n import p y t h o n _ s o l v e r s _ w r a p p e r _ c o n v e c t i o n _ d i f f u s i o n as so lver_wr apper from K r at os M ul t i p h y s i c s . anal ysis_st age import AnalysisStage from scipy . interpolate import Rbf # Used to intra / extrapolate gradients from the physical domain

6 7 8 9

# --- Main analysis class inheriting from Kratos ’ AnalysisStage inf rastruc ture --# This class wraps the entire solution process and integrates the proposed strategy . class E x t e n d e d G r a d i e n t M e t h o d C o n v e c t i o n D i f f u s i o n A n a l y s i s ( AnalysisStage ) :

10 11 12 13

def RunSolut io n Lo op ( self ) : # Retrieve the primary unknown ( e . g . temperature or concentration ) from solver settings unknown _ v a r i a b l e = self . _ G e t U n k n o w n V a r i a b l e ()

14 15 16

# Initialize the solution field from the previous time step ( used for gradient extrapolation ) self . _ I n i t i a l i z e O l d S o l u t i o n F i e l d ( u nk n o w n _v a r i a b l e )

17 18 19 20 21

# Loop over time steps while self . K e e p A d v a n c i n g S o l u t i o n L o o p () : self . _AdvanceTime () self . I n i t i a l i z e S o l u t i o n S t e p ()

22 23 24 25 26

# Begin fixed - point iterations for applying Dirichlet BCs strongly on an unfitted mesh self . i t e r a ti o n _ n u m b e r = 0 while self . _NotConverged () : self . _GetSolver () . Predict () # Kratos predictor step

27 28 29

# Proposed method : apply Dirichlet BCs using auxiliary system self . A p p l y D i r i c h l e t B o u n d a r y C o n d i t i o n s ()

30 31 32

# Solve the physical problem after applying boundary conditions self . _GetSolver () . S o l v e S o l u t i o n S t e p ()

33 34 35

# Check error between iterations to assess convergence self . _ U p d a t e C o n v e r g e n c e S t a t u s ()

36 37 38 39 40

# Update stored field and finalize current step self . _ U p d a t e O l d S o l u t i o n F i e l d ( u nk n o w n _ v a r i a b le ) self . F i n a l i z e S o l u t i o n S t e p () self . O u t p u t S o l u t i o n S t e p ()

41 42 43 44 45 46 47

def A p p l y D i r i c h l e t B o u n d a r y C o n d i t i o n s ( self ) : # Extract relevant submodel parts from Kratos ’ data structures self . model_part = self . _GetSolver () . G e t C o m p u t i n g M o d e l P a r t () self . i n t e r s e c t e d _ e l e m e n t s _ s u b _ m o d e l _ p a r t = self . model_part . Ge tS u bM od el P ar t ( " i n t e r s e c t e d _ e l e m e n t s ") self . a c t i v e _ e l e m e n t s _ s u b _ m o d e l _ p a r t = self . model_part . Ge tS ub M od el Pa r t ( " a ct i ve _e le m en ts " ) self . e m b e d d e d _ b o d y _ b o u n d a r y _ m o d e l _ p a r t = self . model_part . GetModel () . GetModelPart ( " embedded_body_boundary ")

48 49 50 51 52

# Assemble the LHS matrix and the RHS vector from the boundary only once per time step if self . i te r a t i o n _ n u m b e r == 0: self . LHS = self . CalculateLHS () # Compute the A matrix ( Includes both boundary and intersected - elements terms ) self . R H S B o u n d a r y C o n t r i b u t i o n = self . C a l c u l a t e R H S B o u n d a r y C o n t r i b u t i o n () # Compute f_gamma

53 54 55 56

# Compute the RHS contribution from the intersected elements ( updated at each iteration ) self . R H S T r i m m e d E l e m e n t s C o n t r i b u t i o n = self . C a l c u l a t e R H S T r i m m e d E l e m e n t s C o n t r i b u t i o n () # Compute ( f ^ k ) _trim self . RHS = self . R H S B o u n d a r y C o n t r i b u t i o n + self . R H S T r i m m e d E l e m e n t s C o n t r i b u t i o n

57 58 59

# Solve the auxiliary linear system to compute nodal Dirichlet values solution_dir = self . _ S o l v e L i n e a r S y s t e m ( self . LHS , self . RHS )

60 61 62 63 64

# Apply the computed Dirichlet values to the nodes belonging to the intersected elements # ( This typically calls node . Fix ( VARIABLE ) and node . S e t S o l u t i o n S t e p V a l u e ( VARIABLE , value ) # for each node in the intersected elements sub - model part ) self . _ A p p l y B o u n d a r y V a l u e s ( solution_dir )

65 66 67 68 69 70

# Assemble the global LHS matrix from both boundary and intersected element contributions def CalculateLHS ( self ) : LHS_bc = self . C a l c u l a t e L H S B o u n d a r y C o n t r i b u t i o n () LHS_trimmed = self . C a l c u l a t e L H S T r i m m e d E l e m e n t s C o n t r i b u t i o n () return self . _AssembleLHS ( LHS_bc , LHS_trimmed )

72 73 74

# Compute the RHS vector corresponding to boundary integration ( fixed per time step ) def C a l c u l a t e R H S B o u n d a r y C o n t r i b u t i o n ( self ) : return self . _ I n t e g r a t e B o u n d a r y B C s ()

75 76 77 78 79

# Use gradient extrapolation ( via RBF or MLS ) to compute trimmed - element RHS at each iteration def C a l c u l a t e R H S T r i m m e d E l e m e n t s C o n t r i b u t i o n ( self ) : g r a d i e n t _ i n t e r p o l a t o r = self . _ B u i l d R B F G r a d i e n t I n t e r p o l a t o r () return self . _ I n t e g r a t e T r i m m e d E l e m e n t C o n t r i b u t i o n s ( g r a d i e n t _ i n t e r p o l a t o r )

80 81 82 83 84

# --- Main script execution --# Loads parameters , constructs the analysis class , and runs the simulation if __name__ == " __main__ " : p r o j e c t _ p a r a m e t e r s _ f i l e _ n a m e = " P r o j e c t P a r a m e t e r s . json "

85 86 87 88 89

# Read solver and problem configuration from file . # This JSON file contains geometry and mesh definition , material properties , boundary conditions , solver settings and output configuration with open ( proj ect_par ameters _file_na me , ’r ’) as pa rameter _file : parameters = K r a t o s M u l t i p h y s i c s . Parameters ( parame ter_file . read () )

90 91 92 93 94 95 96 97

# Create model and run the simulation model = K r a t o s M u l t i p h y s i c s . Model () simulation = E x t e n d e d G r a d i e n t M e t h o d C o n v e c t i o n D i f f u s i o n A n a l y s i s ( model , parameters ) simulation . Initialize () # The elements cl assific ation and deactivation ( element . Set ( ACTIVE , False ) ) # is performed here simulation . R u nS ol ut i on Lo o p () simulation . Finalize ()

Listing 1: Core implementation of the proposed algorithm in Kratos.

Ω

ΓD r = 0.25

L = 1.0

ΓD

y x L = 1.0 Figure 9: Illustration of the computational domain, consisting of a [1×1] square with an embedded circular hole of radius r = 0.25, centered at C = (0.5, 0.5).

6.1. FEM discretizations In this section, we present numerical results for the application of the proposed method to the Poisson problem using different Finite Element Method (FEM) discretizations. The primary objective is to evaluate the effectiveness of the strong Dirichlet boundary condition enforcement technique for unfitted meshes and compare its performance across different element types and polynomial orders. We consider the manufactured solution already defined above for the Poisson problem described in Eq. 1a, which allows us to evaluate the L2 -norm error for different discretizations. As described in Section 6, the computational domain consists of a square region with an embedded circular hole, with Dirichlet BCs applied to the entire boundary. 16

The numerical experiments explore the performance of the method under different FEM discretizations, specifically considering linear triangular elements, as well as quadratic quadratic elements. For each case, we assess the convergence behaviour and accuracy of the numerical solution. Additionally, we analyze the impact of gradient approximation techniques (Radial Basis Functions (RBF) and Moving Least Squares (MLS)) on the accuracy of the strong Dirichlet BCs enforcement. 6.1.1. FEM Discretization with linear triangular elements This section examines the accuracy and convergence behaviour of the proposed strong Dirichlet BCs enforcement approach in unfitted triangular meshes. Figs.(10) and (11) present the L2 -norm error convergence for different gradient approximation strategies. The analysis considers two interpolation methods: Radial Basis Function (RBF) interpolation and Moving Least Squares (MLS) approximation. The RBF method is evaluated with both a linear and a multiquadric basis, while the MLS approach is tested with quadratic and cubic polynomial orders. Additionally, the performance of enforcing the exact gradient is compared against the MLS-based approximations (Fig. 11). RBF (linear) Gradient Interpolation RBF (multiquadric) Gradient Interpolation MLS (quadratic) Gradient Interpolation 10−2

kφ − φh kL2 (Ω)

MLS (cubic) Gradient Interpolation

2 1

10−3

2 × 10−2

3 × 10−2

4 × 10−2

6 × 10−2

10−1

Figure 10: Comparison of L2 -norm error convergence for linear triangular elements under different gradient approximation techniques. The error is computed for the manufactured solution ϕ(x, y) = sin(πx) cos(πy) using an unfitted boundary mesh. The study evaluates four gradient interpolation methods: RBF with linear basis (blue triangles), RBF with multiquadric basis (green triangles), MLS with quadratic basis (red squares), and MLS with cubic basis (cyan circles). The dashed black line represents the theoretical second-order convergence rate (O(h2 )).

All methods exhibit the expected second-order convergence rate (O(h2 )) in the L2 -norm error, confirming theoretical predictions for first-order FEM discretizations. However, multiquadric RBF interpolation consistently produces higher errors, emphasizing the impact of basis function selection on accuracy. In contrast, MLS cubic interpolation closely matches exact gradient imposition (Fig. 11), demonstrating its effectiveness in gradient approximation. The shown results confirm that MLS interpolation, both the quadratic and cubic variant, represent a robust and accurate approach for the gradient approximation in this algorithm. Finally, Fig. 12 provides a visualization of the numerical solution and the corresponding error distribution for the Poisson problem. The results correspond to an average element size of h = 0.0317, using MLS cubic interpolation for gradient approximation inside trimmed elements. 6.1.2. FEM Discretization with quadratic quadrilateral elements The performance of the strong-like Dirichlet BCs enforcement is analyzed for unfitted quadratic quadrilateral discretizations. The discretization employs 9-node quadratic quadrilateral elements. 17

10−2

kφ − φh kL2 (Ω)

Exact Gradient Imposition MLS (cubic) Gradient Interpolation

10−3

2 1

2 × 10−2

3 × 10−2

4 × 10−2

6 × 10−2

10−1

Figure 11: Comparison of L2 -norm error convergence for linear triangular elements discretization. The study compares the L2 -norm error using MLS with cubic basis functions gradient approximation (cyan pentagons) and and the L2 -norm error imposing the exact gradient in the trimmed elements (blue triangles). The dashed black line represents the theoretical second-order convergence rate (O(h2 )).

1.0

1.000

1.0 0.0006

0.796 0.8

0.592

0.0004

0.388

−0.020 0.4

−0.224

0.0002 0.6 0.0000

0.184

Numerical Solution

0.6

0.4 −0.0002

−0.428 0.2

−0.0004

0.2

−0.632

Error field φh (x, y) − φ(x, y)

0.8

−0.837

−0.0006 0.0

0.0 0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

(a) Numerical solution

(b) Error distribution

Figure 12: Comparison of the numerical solution and the corresponding error distribution for the plate with a hole problem using linear triangular elements discretization. (a) Computed numerical solution. (b) Error distribution, illustrating the deviation of the numerical solution from the exact solution.

Figs. (13) and (14) present the L2 -norm error convergence for different gradient approximation strategies. The study evaluates again MLS interpolation, RBF-based methods with linear and multiquadric basis functions, and direct exact gradient enforcement. The convergence behaviour differs among the tested methods. MLS cubic interpolation exhibits optimal third-order convergence (O(h3 )), consistent with theoretical expectations for quadratic FEM discretizations. In contrast, the RBF-based approaches show a slower error decay, indicating slightly suboptimal convergence. This suggests that while MLS effectively reconstructs the gradient field, the accuracy of RBF interpolation is more sensitive to the choice of basis function. These results highlight the superior performance of MLS, which 18

stems from its local minimization-based formulation. RBF (linear) Gradient Interpolation RBF (multiquadric) Gradient Interpolation 10−2

MLS (quadratic) Gradient Interpolation MLS (cubic) Gradient Interpolation

3 1

kφ − φh kL2 (Ω)

10−3

10−4

10−5

2 × 10−2

10−2

3 × 10−2

4 × 10−2

6 × 10−2

Figure 13: Comparison of L2 -norm error convergence for 9-node quadratic quadrilateral elements under different gradient approximation techniques. The error is computed for the manufactured solution ϕ(x, y) = sin(πx) cos(πy) using an unfitted boundary mesh. The study evaluates four gradient interpolation methods: RBF with linear basis (blue triangles), RBF with multiquadric basis (green triangles), MLS with quadratic basis (red squares), and MLS with cubic basis (cyan circles). The dashed black line represents the theoretical third-order convergence rate (O(h2 )).

10−2

Exact Gradient Imposition MLS (cubic) Gradient Interpolation

kφ − φh kL2 (Ω)

10−3

10−4

10−5

10−6 2 × 10−2

3 × 10−2

4 × 10−2

6 × 10−2

Figure 14: Comparison of L2 -norm error convergence for 9-node quadratic quadrilateral elements discretization. The study compares the L2 -norm error using MLS with cubic basis functions gradient approximation (cyan pentagons) and and the L2 -norm error imposing the exact gradient in the trimmed elements (blue triangles). The dashed black line represents the theoretical third-order convergence rate (O(h3 )).

The comparison between MLS cubic interpolation and exact gradient imposition (Fig. 14) reveals a greater discrepancy than in the previously analyzed cases. While MLS maintains the expected third-order convergence, its 19

absolute error remains higher than that of direct gradient enforcement. This suggests that for quadratic quadrilateral elements, the MLS approach introduces a larger numerical error, possibly due to increased sensitivity in gradient reconstruction when higher-order approximations are used. Finally, Fig. 15 provides a visualization of the numerical solution and error distribution for the Poisson problem, computed with an average element size of h = 0.014.

1.000

×10−6

1.0

0.796

0.75 0.8

0.6

0.184 −0.020

0.4

−0.224 −0.429

0.2

0.50 0.25

0.6 0.00

0.388

Numerical Solution

0.592

−0.25

0.4

−0.50 0.2

−0.633

−0.75

−0.837

−1.00

0.0 0.2

0.4

0.6

Error field φh (x, y) − φ(x, y)

0.8

0.0

0.2

0.4

0.6

0.8

1.0

(a) Numerical solution

(b) Error distribution

Figure 15: Comparison of the numerical solution and the corresponding error distribution for the plate with a hole problem using 9-node quadratic quadrilateral elements discretization. (a) Computed numerical solution. (b) Error distribution, illustrating the deviation of the numerical solution from the exact solution.

6.2. IGA discretizations In this section, we present the L2 -norm error convergence study for B-Spline basis functions, which form the foundation of Isogeometric Analysis (IGA) [60]. Rather than revisiting the mathematical formulation of B-Splines, typically defined via the Cox–de Boor recursion formula [61, 62], we refer the reader to the standard reference The NURBS Book by Piegl and Tiller [63], which offers a thorough and accessible introduction to the theory and practical aspects of B-Splines and NURBS. 6.2.1. Quadratic IGA discretizations This section studies the accuracy and convergence properties of the proposed strong Dirichlet BCs enforcement technique in unfitted (trimmed) quadratic B-Spline discretizations. Figs. (16) and (17) present the L2 -norm error convergence for different gradient approximation strategies. In the section discussing FEM discretizations, we demonstrated that the accuracy gain from using a cubic basis in MLS, as opposed to a quadratic basis, is negligible, while the computational cost of cubic basis functions is significantly higher. Given this observation, in this section, we reassess MLS interpolation using linear and quadratic basis, RBF-based methods with linear basis, and direct exact gradient enforcement. The convergence behavior varies among the tested methods. Both MLS quadratic and MLS linear interpolations exhibit the expected third-order convergence rate (O(h3 )), in agreement with theoretical predictions for quadratic IGA discretizations. In contrast, the RBF-based approach shows a slower error decay, indicating suboptimal convergence. This suggests that while MLS effectively reconstructs the gradient field, the reduced accuracy of RBF interpolation is likely due to the global nature of the interpolation, which can introduce approximation errors and reduced stability, particularly in regions with non-uniform point distributions. It is worth to mention that in this case, among the available RBF methods, only the linear basis (ϕ(r) = r) demonstrated proper accuracy in gradient reconstruction inside trimmed elements, highlighting its relative robustness in this specific setting. 20

RBF (linear) Gradient Interpolation 10−1

MLS (linear) Gradient Interpolation MLS (quadratic) Gradient Interpolation

10−2

kφ − φh kL2 (Ω)

10−3

10−4

10−5 10−2

2 × 10−2

3 × 10−2

4 × 10−2

6 × 10−2

Figure 16: Comparison of L2 -norm error convergence for a quadratic IGA discretization under different gradient approximation techniques. The error is computed for the manufactured solution ϕ(x, y) = sin(πx) cos(πy) using an unfitted boundary mesh. The study evaluates three gradient interpolation methods: RBF with linear basis (blue triangles), MLS with linear basis (red squares), and MLS with quadratic basis (cyan circles). The dashed black line represents the theoretical third-order convergence rate (O(h3 )).

Exact Gradient Imposition MLS (quadratic) Gradient Interpolation

10−2

10−3 3

kφ − φh kL2 (Ω)

10−4

10−5

10−6

10−2

2 × 10−2

3 × 10−2

4 × 10−2

6 × 10−2

Figure 17: Comparison of L2 -norm error convergence for a quadratic B-Spline discretization. The study compares the L2 -norm error using MLS with quadratic basis functions gradient approximation (cyan pentagons) and and the L2 -norm error imposing the exact gradient in the trimmed elements (blue triangles). The dashed black line represents the theoretical third-order convergence rate (O(h3 )).

Analyzing the results of gradient interpolation using quadratic MLS and exact gradient enforcement (Fig. 17) reveals some difference in accuracy. While MLS achieves the expected third-order convergence, its error remains larger than that of enforcing the exact gradient. This suggests that for quadratic IGA elements, the gradient approximation 21

process in MLS introduces greater inaccuracies compared to lower-order discretizations. It is important to highlight that the interpolation points used for gradient reconstruction belong to the so-called unperturbed elements (Def. 1), where no degrees of freedom are modified by this method. In the case of IGA with a quadratic discretization (p = 2), where a basis function spans over p + 1 = 3 knot spans, this implies that the gradient information is extracted from interpolation points located at a distance of approximately O(3h), with h being the knot span size. Conversely, in the FEM case, where shape function support is strictly local within each element, the gradient information used for interpolation is, in the worst-case scenario, at a distance of only O(2h) from the evaluation point. This fundamental difference in gradient information clearly contributes to the observed accuracy discrepancies between MLS-based and exact gradient enforcement methods. An important aspect to analyze is how the L2 -norm error varies with the iterative solver tolerance ϵ and how the number of iterations needed for convergence depends on the mesh size h for different target accuracies. The corresponding results are presented in Fig. 18 for the embedded circle case. From these plots, it can be observed that reducing the iterative solver tolerance ϵ leads to a decrease in the L2 -norm error, as expected. However, this reduction in tolerance also results in a higher number of iterations N required for convergence. As in most engineering problems, an optimal balance must be struck between computational cost and the required accuracy for a given application. In this particular case, a tolerance of ϵ = 0.001 provides satisfactory results. It should also be noted that only a few iterations are typically needed to reach convergence, highlighting the efficiency of the proposed algorithm. Additionally, the first iteration is always performed under the assumption of a zero normal gradient within each trimmed element (∇ϕ · n = 0), so there is no associated computational cost for the gradient approximation inside the trimmed elements during this step. Tolerance = 0.01 Tolerance = 0.001 Tolerance = 0.0001

Tolerance = 0.01 Tolerance = 0.001 Tolerance = 0.0001

Number of iterations for convergence N

−2

10−3

kφ − φh kL2 (Ω)

10−4

10−5

2 × 10−2

3 × 10−2

4 × 10−2

6 × 10−2

0.01

0.02

0.03

0.04

0.05

0.06

(a) L2 -norm error convergence plot for different algorithm tolerances

(b) Number of iterations N for convergence

Figure 18: Analysis of solver accuracy and convergence behavior for the embedded circle geometry using quadratic IGA discretization. (a) Convergence of the L2 -norm error for different iterative solver tolerances. (b) Required number of iterations N to achieve convergence for varying solver tolerances.

Finally, Fig. 19 provides a visualization of the numerical solution and error distribution for the Poisson problem, computed with an element size of h = 0.021. 6.2.2. Cubic IGA discretizations In this section, we present an evaluation of our developed method for strongly imposing Dirichlet BCs, focusing on its performance in terms of accuracy and convergence when applied to unfitted (trimmed) cubic B-Spline 22

1.000

1.0

0.796

0.00010 0.8

0.592

−0.020 0.4

−0.224 −0.429

0.2

0.6 0.00000

0.184

0.00005

Numerical Solution

0.388 0.6

0.4 −0.00005 0.2

−0.633

Error field φh (x, y) − φ(x, y)

0.8

−0.00010

−0.837 0.0 0.2

0.4

0.6

0.8

0.0

0.2

0.4

0.6

0.8

1.0

(a) Numerical solution

(b) Error distribution

Figure 19: Comparison of the numerical solution and the corresponding error distribution for the plate with a hole problem using quadratic IGA discretization. (a) Computed numerical solution. (b) Error distribution, illustrating the deviation of the numerical solution from the exact solution.

discretizations. Figs. (20) and (21) present the L2 -norm error convergence for different gradient approximation strategies. Among the methods evaluated, convergence characteristics differ noticeably. Both the quadratic MLS and linear MLS interpolation techniques achieve a fourth-order convergence rate (O(h4 )), aligning with theoretical expectations for cubic IGA discretizations. In contrast, the RBF-based method demonstrates a slower reduction in error, suggesting that its convergence performance is suboptimal. Once again, the quadratic MLS method delivers the best accuracy, underscoring its robustness as the most reliable approach for gradient approximation. RBF (linear) Gradient Interpolation MLS (linear) Gradient Interpolation 10−1

MLS (quadratic) Gradient Interpolation 4 1

kφ − φh kL2 (Ω)

10−2

10−3

10−4

10−5 2 × 10−2

3 × 10−2

4 × 10−2

Figure 20: Comparison of L2 -norm error convergence for a cubic IGA discretization under different gradient approximation techniques. The error is computed for the manufactured solution ϕ(x, y) = sin(πx) cos(πy) using an unfitted boundary mesh. The study evaluates three gradient interpolation methods: RBF with linear basis (blue triangles), MLS with linear basis (red squares), and MLS with quadratic basis (cyan circles). The dashed black line represents the theoretical fourth-order convergence rate (O(h4 )).

10−2

Exact Gradient Imposition MLS (quadratic) Gradient Interpolation 1

10−3

kφ − φh kL2 (Ω)

10−4

10−5

10−6

10−7

10−8

2 × 10−2

3 × 10−2

4 × 10−2

Figure 21: Comparison of L2 -norm error convergence for a cubic B-Spline discretization. The study compares the L2 -norm error using MLS with quadratic basis functions gradient approximation (cyan pentagons) and and the L2 -norm error imposing the exact gradient in the trimmed elements (blue triangles). The dashed black line represents the theoretical third-order convergence rate (O(h4 )).

An analysis of the gradient interpolation results using quadratic MLS versus exact gradient enforcement (see Fig. 21) indicates a noticeable accuracy difference. Although the MLS method achieves the anticipated 4th-order convergence (O(h4 )), its error level remains higher than that obtained with exact gradient enforcement. This is once again attributed to the fact that the gradient information for interpolation is sourced from the so-called unperturbed elements, which, for cubic discretizations, are located approximately 4h away from the interpolation point. An analysis of the gradient interpolation results using quadratic MLS versus exact gradient enforcement (Fig. 21) indicates a noticeable accuracy gap. Although the MLS method achieves the anticipated O(h4 ) convergence rate, its absolute error remains consistently higher than that obtained with exact gradient enforcement. This is attributed to the fact that the gradient information for interpolation is sourced exclusively from the so-called unperturbed elements, which, for cubic discretizations, are located approximately 4h away from the interpolation point (see 1). This distance amplifies the approximation error in trimmed regions. At present, no specific countermeasure has been incorporated to mitigate this loss of accuracy for high-order B-Splines discretizations; addressing this will be part of future work. We note, however, that the impact is significantly reduced for lower-order elements (including quadratic), for which our method achieves accuracy comparable to established techniques. As a final illustration, Fig. 22 presents the numerical solution and error distribution for the Poisson problem using quadratic MLS gradient approximation, alongside the error distribution plot for exact gradient enforcement, computed with an element size of h = 0.019. 6.3. Comparison of the proposed strategy, IBRA and SBM In this section, we compare our proposed method with the Isogeometric B-Rep Analysis (IBRA) methodology [39, 40, 41] and the Shifted Boundary Method (SBM) [44, 45, 46], evaluating both the accuracy in terms of the L2 -norm error and the conditioning of the resulting linear system. To this end, we consider the geometry shown in Fig. 23a and prescribe the following manufactured solution of the Poisson problem as a Dirichlet boundary condition on the entire boundary, Γ = ΓD : ϕ(x, y) = sin(x) sinh(y) where the corresponding source term is identically zero throughout the domain, i.e., σ(x, y) = 0. 24

(41)

1.000

1.0

0.796

0.0002 0.8

0.592

−0.020 0.4

−0.224 −0.429

0.2

0.6 0.0000

0.184

0.0001

Numerical Solution

0.388 0.6

0.4 −0.0001 0.2

−0.633

Error field φh (x, y) − φ(x, y)

0.8

−0.0002

−0.837 0.0 0.2

0.4

0.6

0.8

0.0

0.2

0.4

0.6

0.8

1.0

(a) Numerical solution

(b) Error distribution quadratic MLS as the gradient approximator

×10−8

1.0

0.8

2 0.6

0 0.4 −2 0.2

−4

Error field φh (x, y) − φ(x, y)

−6

0.0 0.0

0.2

0.4

0.6

0.8

1.0

x (c) Error distribution for exact gradient imposition

Figure 22: Comparison of the numerical solution and the corresponding error distribution for the plate with a hole problem using cubic IGA discretization. (a) Computed numerical solution. (b) Error distribution, illustrating the deviation of the numerical solution from the exact solution. c) Error distribution, for exact gradient imposition

The primary distinction between the methods lies in how each handles the imposition of boundary conditions. Our approach employs a non-intrusive framework that enforces Dirichlet boundary conditions strongly by solving an auxiliary algebraic problem, independent from the physics-based weak form. This strategy not only simplifies implementation but also avoids the need for integrating over trimmed knot spans, as the stabilization term is computed across the entire intersected element. In contrast, the Isogeometric B-Rep Analysis (IBRA) imposes Dirichlet boundary conditions weakly using techniques such as penalty methods, Lagrange multipliers, or Nitsche’s method. IBRA maintains optimality by triangulating and integrating the active portions of cut knot spans, as illustrated in Fig. 24. This often leads to a high concentration of integration points near trimmed boundaries, where accurate numerical integration is achieved through fine triangulation. To address the computational burden and potential over-integration in these regions, reduced integration techniques using fewer Gauss points have been proposed [42]. However, the presence of very small trimmed elements can introduce the so-called “small-cut cell instability”, which may significantly degrade the conditioning of the resulting linear systems. In the following examples, IBRA employs a penalty-free Nitsche method for boundary 25

2.0 3.315

1.8

ΓD y

L = 2.0

ΓD

y x

1.6

2.947

1.4

2.578

1.2

2.210

1.0

1.842

0.8

1.473

0.6

1.105

0.4

0.737

0.2

0.368

0.0

Numerical Solution

Ω

0.000 0.0

0.2

0.4

0.6

L = 2.0

0.8

1.0

1.2

1.4

1.6

1.8

2.0

(a) Computational domain

(b) Solution field

Figure 23: Computational setup and solution field for the problem defined by the field ϕ(x, y) = sin(x) sinh(y). (a) The computational domain consists of a square domain Ω = [0, 2] × [0, 2] with an embedded quadrilateral hole where Dirichlet BCs are imposed. (b) Solution field ϕ(x, y). 2.086

1.825

1.564

1.304

1.043

Integration points

0.782

0.521

0.261

0.000 0.000

0.261

0.521

0.782

1.043 x

1.304

1.564

1.825

2.086

Figure 24: Distribution of IBRA integration points for quadratic B-splines, with 9 integration points per knot span. The plot is shown in the parametric space (ξ, η), which in this case coincides with the physical domain (x, y). A high concentration of integration points is observed near the trimmed boundary, where triangulation is applied to enforce accurate integration in cut elements.

condition enforcement, as introduced in [64] and further validated in [46]. The Shifted Boundary Method (SBM), on the other hand, avoids direct integration over trimmed elements by shifting the boundary conditions from the true boundary to a surrogate boundary located entirely within the computational domain. This enables the use of standard integration schemes, sidestepping the geometric complexities of cut elements. Dirichlet conditions are weakly enforced on the surrogate boundary using techniques like Nitsche’s method. In Figs. (25a) and (25b), we compare the three approaches for the embedded square computational domain in terms of the L2 -norm error. For degree p = 1, the three methods yield very similar results in terms of the L2 -norm error. However, for degree p = 2, IBRA and SBM exhibit better accuracy than the proposed approach, likely due to inaccuracies introduced by the gradient approximation within the trimmed elements. A similar trend is observed for degree p = 3. Interestingly, as shown in Fig. 25b, when the exact gradient is imposed in the trimmed elements, the accuracy of the proposed method remains slightly lower than that of IBRA and SBM, but the difference in accuracy is significantly reduced. This suggests that the method remains optimal in terms of the L2 -norm error, with the loss 26

IBRA Proposed method - Cubic B-Splines (p = 3) Linear MLS gradient approximation Proposed method - Cubic B-Splines (p = 3) Exact gradient imposition SBM

IBRA Method (p = 1 and p = 2) 10−1

10−1

Proposed method - Linear B-Splines (p = 1) Quadratic MLS gradient approximation Proposed method - Quadratic B-Splines (p = 2) Quadratic MLS gradient approximation SBM (p = 1 and p = 2)

10−2 4

10−3

2 1

kφ − φh kL2 (Ω)

10−3

10−5

10−4

10−7

10−5

10−6 10−9

10−7 10−2

10−1

10−2

(a) Comparison of the L2 -norm error between IBRA, SBM and the proposed method

2 × 10−2

3 × 10−2

4 × 10−2

6 × 10−2

(b) Comparison of the L2 -norm error between IBRA, SBM and the proposed method

Figure 25: Error convergence comparison between IBRA, SBM and the proposed approach. (a) Convergence results for IBRA, SBM and the proposed method using linear (p = 1) and quadratic (p = 2) B-Spline discretizations. In the proposed method, the gradient is approximated using a quadratic MLS scheme. (b) Convergence results for IBRA, SBM and the proposed method using cubic B-Spline discretizations (p = 3). The proposed method is evaluated using both a linear MLS gradient approximation and exact gradient imposition.

of accuracy primarily originating from the gradient reconstruction in the trimmed elements. This issue becomes increasingly complex for higher polynomial degrees. The condition number of the system matrix κ(A) is crucial in FEM-like simulations, as it affects numerical stability and iterative solver performance. In unfitted methods, small cut-cell instabilities and the weak imposition of boundary conditions can significantly increase the condition number, sometimes impairing the use of iterative solvers. Fig. 26 examines the condition number behavior for the embedded square geometry shown in Fig. 23a, where the proposed approach, SBM and IBRA are used to impose Dirichlet BCs on the complete boundary. Compared to IBRA, both the formulation presented in this paper and SBM consistently produce condition numbers κ(A) that are several orders of magnitude lower. To conclude, we emphasize that our proposed method, whether combined with low order or higher order discretizations, remains stable and unaffected by small-cut cell instability, a frequent challenge in CutFEM-like approaches and other methods that require integration over trimmed elements such as IBRA. IBRA-based discretizations tend to encounter this issue, often leading to numerical difficulties that hinder the performance of iterative solvers. In our case, the resulting system remains well-conditioned, ensuring efficient and reliable solutions with both direct and iterative solvers. 6.4. Capabilities of the proposed approach In this section, we investigate two representative examples in which the proposed method is applied to domains with complex boundary geometries. Both cases involve solving the Poisson problem introduced in Eq. 1, using a Cartesian mesh composed of linear and quadratic B-Spline discretizations. To enable a controlled evaluation of the method’s accuracy, we adopt again the manufactured solution ϕ(x, y) = sin(πx) cos(πy), which is also prescribed as a Dirichlet boundary condition along both the inner and outer boundaries of the domain. In addition, a quantitative comparison with the shifted boundary method (SBM) is provided for the first example to assess the relative accuracy of the proposed approach in handling more complex immersed geometries. 27

1014

IBRA - p = 1 IBRA - p = 2 IBRA - p = 3 SBM - p = 1 SBM - p = 2 SBM - p = 3

1012

106

105

Condition Number κ(A)

1010

108

106

104

103

102

101

2 × 10−2

5 × 10−2

1 × 10−1

2 × 10−1

Proposed method - p = 1 Proposed method - p = 2 Proposed method - p = 3 Proposed method - p = 4 2 × 10−2

(a) Condition number κ(A) of IBRA and SBM system matrix

5 × 10−2

1 × 10−1

2 × 10−1

5 × 10−1

(b) Condition number κ(A) of the proposed approach system matrix

Figure 26: Comparison of the condition number behavior for IBRA (left panel) and our proposed method (right panel), considering different orders of basis functions for the embedded square case in Fig. 23a. The order of the basis functions is represented by the following symbols: triangles (p = 1), crosses (p = 2), stars (p = 3) and circles (p = 4).

The first case, shown in Fig. 27a, applies the method to the external boundary of the two-dimensional Stanford Bunny. In the figure, blue dots represent the active integration points (located within active elements), which contribute to assembling both the left- and right-hand sides of the governing equations. Red dots indicate the intersected integration points (located within intersected elements), which contribute to the stabilization term in the additional algebraic problem. The convergence results in Fig. 27b confirm that the method achieves the expected optimal convergence rate, and for low-order elements, its accuracy is comparable to that of the Shifted Boundary Method. A second, more challenging case is shown in Fig. 28a, where the domain contains multiple inner loops. The corresponding convergence study in Fig. 28b again confirms the robustness of the method and its ability to maintain optimal convergence properties. While all examples in the present work are two-dimensional, the proposed formulation is dimension-independent and can be directly extended to three-dimensional problems without modification to its mathematical structure. In 3D, the additional challenges are mainly related to the naturally higher computational effort required for element classification and gradient reconstruction, which stem from the increased problem size rather than from any limitation of the method itself. 6.5. Explicit transient heat diffusion This section evaluates the performance of the proposed method for the explicit simulation of transient heat diffusion, governed by the following partial differential equation (PDE) ∂ϕ = ∇ · (k∇ϕ) + Q(x, t) (42) ∂t where ρ is the density of the material, c p is the specific heat capacity, k is the thermal conductivity, and Q(x, t) represents a heat source that varies with position x and time t. For simplicity, we set ρ = c p = k = 1 in all subsequent derivations. To numerically solve this equation, we employ an explicit time integration scheme, specifically the explicit Euler method [65], combined with a finite element spatial discretization with linear triangular elements. Since explicit methods require the inversion of the consistent mass matrix M, a common simplification is to use a lumped mass matrix approximation, replacing M −1 with a diagonal matrix to improve computational efficiency. However, as with any explicit scheme, the stability of the method is conditionally restricted by the Courant–Friedrichs–Lewy ρc p

Proposed method (p=1) - Quadratic MLS gradient approximation

Active Elements Intersected Elements

1.4

SBM (p=1) Proposed method (p=2) - Quadratic MLS gradient approximation

10−2

SBM (p=2)

2 1

10−3

kφ − φh kL2 (Ω)

1.2

1.0

3 1 −4

0.8

10−5

0.6

10−6

0.6

0.8

1.0 x

1.2

10−2

1.4

2 × 10−2

3 × 10−2 4 × 10−2

6 × 10−2

(b) Convergence study for the L2 -norm error and comparison with SBM

(a) Computational set-up, active and intersected elements

Figure 27: Capabilities of the proposed approach. a) the two-dimensional Stanford Bunny is used as the external embedded boundary, and the blue and red dots are the active and intersected integration points respectively. b) the convergence study for linear and quadratic B-Spline discretizations, demonstrating optimal convergence and a comparison with the SBM, demonstrating similar accuracy for low-order elements study for linear quadrilateral elements

Proposed method (p=1) - Quadratic MLS gradient approximation

Active Integration Points

10−2

2.0

kφ − φh kL2 (Ω)

1.5

1.0

0.5

2 1

10−3

0.0

0.5

1.0 x

1.5

10−2

2.0

2 × 10−2

3 × 10−2

4 × 10−2

6 × 10−2

(b) Convergence study for the L2 -norm error

(a) Computational set-up and active integration points

Figure 28: Capabilities of the proposed approach. a) a complex example composed of an external and three internal embedded boundaries treated with the proposed method b) the convergence study for linear quadrilateral elements, demonstrating optimal convergence study for linear quadrilateral elements

(CFL) condition, which imposes an upper bound on the time step ∆t to prevent numerical instability (any explicit time scheme is inherently conditionally stable). With this example, we aim to demonstrate that for explicit dynamic simulations with smooth and non-smooth solution fields, good accuracy can be achieved without the need of iterating at each time step. This is because, in 29

4.0

Solution with iterations Solution without iterations Exact Solution

Number of iterations

Number of iterations for convergence N

0.75

3.5

Solution Value

0.50

3.0

0.25 0.00

2.5

−0.25

2.0

−0.50

1.5

−0.75 1.0 0.000

0.005

0.010

0.015

0.020 Time [s]

0.025

0.030

0.035

0.040

0.000

(a) Exact solution vs. Numerical solution with and without iterations

0.005

0.010

0.015

0.020 Time [s]

0.025

0.030

0.035

0.040

(b) Iterations vs. time for algorithm convergence

Figure 29: (a) Time-dependent numerical solution at the fixed point (x, y) = (0.5, 0.84) using a time step ∆t = 0.00005 s. The plot compares three different solutions: the solution obtained with iterations (solid red line), the solution obtained without iterations (solid green line over the red one), and the exact solution (dashed blue line).The results demonstrate excellent accuracy for this transient smooth problem, even without iterations at each time step. (b) Number of iterations required for convergence when the algorithm tolerance is set to ϵ = 0.001. The peak number of iterations occurs at the initial stage of the transient simulation and during instances where the solution’s concavity changes.

explicit simulations, time steps are typically very small, resulting in negligible gradient differences between successive steps. Consequently, at each time step, the normal gradient for the first (and only) iteration can be assumed to be the same as that from the previous time step. As an initial example, we examine a smooth manufactured solution for the transient Poisson problem (42) within the embedded circle domain defined in Fig. (9 which is discretized with a background mesh of linear triangular elements. The manufactured solution ϕ(x, y, t) is given by ϕ(x, y, t) = sin (πx) cos (πy) sin (50πt) The corresponding source term, Q(x, t) , is given by Q(x, y, t) = sin (πx) cos (πy)[50π cos (50πt) + 2π2 sin (50πt)] In Fig. 29a, we present the time-dependent numerical solution at a specific point in the domain with coordinates (x, y) = (0.5, 0.84), using a time step ∆t = 0.00005 s. The numerical solution is obtained using the forward Euler scheme with a lumped mass matrix. The results clearly demonstrate that, for this transient problem, very good accuracy is achieved without the need for iteration at each time step. Additionally, Fig. 29b illustrates the number of iterations required for convergence when the algorithm tolerance is set to ϵ = 0.001. As depicted in the plot, the maximum number of iterations required for convergence occurs primarily at the initial stage of the transient simulation and during instances where the solution’s concavity changes (around t = 0.02 s). As a second example, we consider a manufactured solution for the transient Poisson problem (42) exhibiting a shock-like pattern in time around t = 0.01 within the embedded circle domain shown in Fig. 9 which is discretized again with a background mesh of linear triangular elements. The manufactured solution ϕ(x, y, t) is given by, ϕ(x, y, t) = sin (πx) cos (πy) tanh [10000(t − 0.01)] The corresponding source term, Q(x, t) , is given by Q(x, y, t) = sin(πx) cos(πy)

10000 + 2π2 tanh(10000(t − 0.01)) cosh2 (10000(t − 0.01))

In Fig. 30a, we depict the time-dependent numerical solution at a specific location within the domain, given by the coordinates (x, y) = (0.5, 0.84), with a time step of ∆t = 0.00005 s. The numerical solution is obtained using 30

5.0

Solution with iterations Solution without iterations Exact Solution

4.5

0.50

Solution Value

Number of iterations

Number of iterations for convergence N

0.75

4.0

0.25

3.5 3.0

0.00

2.5

−0.25

2.0

−0.50

1.5 −0.75 1.0 0.0000

0.0025

0.0050

0.0075

0.0100 0.0125 Time [s]

0.0150

0.0175

0.0200

0.0000

(a) Exact solution vs. Numerical solution with and without iterations

0.0025

0.0050

0.0075

0.0100 0.0125 Time [s]

0.0150

0.0175

0.0200

(b) Iterations vs. time for algorithm convergence

Figure 30: (a) Time-dependent numerical solution at the fixed point (x, y) = (0.5, 0.84) using a time step ∆t = 0.00005 s. The plot compares three different solutions: the solution obtained with iterations (solid red line), the solution obtained without iterations (solid green line over the red one), and the exact solution (dashed blue line). The results demonstrate that even without iterations at each time step, the numerical solution maintains excellent accuracy for this transient non-smooth problem. (b) Number of iterations required for convergence when the algorithm tolerance is set to ϵ = 0.001. The peak number of iterations occurs at the initial stage of the transient simulation and near the shock position.

again the forward Euler scheme with a lumped mass matrix. The results indicate that, despite the transient and nonsmooth nature of the problem, high accuracy is maintained without requiring iterative corrections at each time step. Furthermore, Fig. 30b presents the number of iterations necessary for convergence when the algorithm’s tolerance is set to ϵ = 0.001. As illustrated in the figure, the peak number of iterations occurs predominantly at the initial stage of the transient simulation and around the shock location, approximately at t = 0.01 s. These results demonstrate that the proposed algorithm is well-suited for explicit dynamic simulations, as it enables the strong enforcement of Dirichlet boundary conditions with virtually no additional computational cost, even in the presence of non-smooth solution fields. 7. Conclusions In this work, we have presented a physics-agnostic and non-intrusive iterative approach for the strong-like imposition of Dirichlet boundary conditions in unfitted meshes. The method enforces the prescribed conditions by reformulating the problem as an L2 -norm error minimization, enhanced with a stabilization term that approximates the normal gradient within trimmed elements. This minimization defines an auxiliary algebraic problem, independent from the weak form of the governing equations. A key feature of this approach is that it preserves the original variational formulation, avoiding modifications to the physical system and maintaining its intrinsic properties. The method has been implemented using the Kratos Multiphysics API (release v10.1), and a central contribution of this work is to demonstrate that such a strategy can enable already validated, body-fitted, black-box solvers to operate in unfitted meshes. This is possible as long as four conditions are fulfilled: (i) scripting support, (ii) allowance for the imposition of Dirichlet boundary conditions at the node level, (iii) element deactivation, and (iv) access to the solution gradient in active elements. It is worth noting that the last condition can also be satisfied by externally reconstructing the gradient from nodal values and connectivity information, provided the element formulation is known, making it optional in practice. The treatment of Neumann boundary conditions in this black-box context has also been detailed earlier in the manuscript (see Rem. 6). This greatly expands the flexibility of traditional solvers, allowing them to handle complex geometries without the need for body-fitted meshes or intrusive code changes. The proposed method has been thoroughly validated through the numerical solution of the Poisson problem, using both Finite Element Method (FEM) and Isogeometric Analysis (IGA) discretizations. The results confirm optimal L2 -norm error convergence, demonstrating the effectiveness of the approach. A detailed comparison of different

gradient approximation strategies has shown that Moving Least Squares (MLS) interpolation yields the most accurate reconstruction inside intersected elements, likely due to its local minimization formulation. Additionally, the method eliminates the need for penalty parameter tuning and improves system conditioning compared to Isogeometric B-Rep Analysis (IBRA), while maintaining comparable accuracy. It is especially suitable for explicit dynamic simulations, where small time steps lead to minor changes in the gradient field between iterations. In such cases, the method enforces strong Dirichlet conditions with virtually no additional computational cost, even for non-smooth solutions. From a computational cost standpoint, although the proposed strategy is iterative and therefore entails a higher total runtime than any single-pass method, the additional cost per iteration is small compared to solving the main system. This is because the auxiliary problem for updating Dirichlet values involves only a small fraction of the total degrees of freedom, and in each iteration only the gradient-approximation term needs to be recomputed. In terms of accuracy, the proposed method has shown results comparable to established approaches such as IBRA and SBM for low-order discretizations, but for higher-order discretizations its accuracy is generally lower, primarily due to the gradient reconstruction procedure. It is therefore important to emphasize that the main strength of the proposed approach lies neither in outperforming alternative methods in accuracy nor in minimizing computational cost, but in enabling a solver originally designed for body-fitted solution strategies to be adapted, without intrusive modifications, to work in unfitted scenarios. While the method has shown promising results, further improvements in gradient approximation could enhance its accuracy, particularly in high-order discretizations. Extending the method to additional physical problems is a natural next step, including its application to embedded structures, fluid domains, and fluid–structure interaction. Integration into high-performance computing environments will also enable its use in large-scale simulations. Declarations Conflict of interest. The authors have no competing interests to declare that are relevant to the content of this article. Data availability Data will be made available on request. Code availability The complete implementation of the proposed method, as used to produce all results in this paper, is publicly available at https://github.com/KratosMultiphysics/Kratos/tree/iga/embedded extended gradient method. This repository contains the Kratos Multiphysics (v10.1) implementation of the algorithm. While the present work employs Kratos as the host solver, the formulation is solver-agnostic and can be implemented in any programming language that provides an interface to impose Dirichlet boundary conditions at the node level, deactivate elements, and calculate the solution gradient in active elements. Acknowledgments The authors gratefully acknowledge the Design for IGA-type discretization workflows (GECKO) project. The Design for IGA-type discretization workflows has received funding from the European Union’s Horizon Europe research and Innovation programme under grant agreement No. 101073106, Call: HORIZON-MSCA-2021-DN-01. References [1] T.-P. Fries, An extended finite element method with higher-order elements for curved cracks, Computational Mechanics 31 (2003) L29–L36. [2] A. Yazid, A. Nabbou, A. Hamouine, A state-of-the-art review of the X-FEM for computational fracture mechanics, Applied Mathematical Modelling 33 (12) (2009) 4269–4282. [3] S. Nagaraja, M. Elhaddad, M. Ambati, S. Kollmannsberger, L. De Lorenzis, E. Rank, Phase-field modeling of brittle fracture with multi-level hp-FEM and the finite cell method, Computational Mechanics 63 (6) (2019) 1283–1300.

[4] E. Burman, D. Elfverson, P. Hansbo, M. G. Larson, K. Larsson, Shape optimization using the cut finite element method, Computer Methods in Applied Mechanics and Engineering 328 (2018) 242–261. [5] M. S. Edke, K. H. Chang, Shape optimization for 2-D mixed-mode fracture using Extended FEM (XFEM) and Level Set Method (LSM), Structural and Multidisciplinary Optimization 42 (6) (2010) 725–738. doi:10.1007/s00158-010-0616-5. [6] M. Meßmer, R. N. Asl, S. Kollmannsberger, R. Wüchner, K.-U. Bletzinger, Shape optimization of embedded solids using implicit VertexMorphing, Computer Methods in Applied Mechanics and Engineering 426 (2024) 116999. doi:10.1016/j.cma.2024.116999. [7] R. Zorrilla, R. Rossi, R. Wüchner, E. Oñate, An embedded finite element framework for the resolution of strongly coupled fluid–structure interaction problems. Application to volumetric and membrane-like structures, Computer Methods in Applied Mechanics and Engineering 368 (2020) 113179. [8] B. Schott, C. Ager, W. A. Wall, Monolithic cut finite element–based approaches for fluid-structure interaction, International Journal for Numerical Methods in Engineering 119 (8) (2019) 757–796. [9] L. Boilevin-Kayl, M. A. Fernández, J.-F. Gerbeau, Numerical methods for immersed FSI with thin-walled structures, Computers & Fluids 179 (2019) 744–763. [10] E. Miranda Neiva, Large-scale tree-based unfitted finite elements for metal additive manufacturing, Ph.D. thesis, Universitat Politècnica de Catalunya (2020). [11] M. C. Wu, H. M. Muchowski, E. L. Johnson, M. R. Rajanna, M.-C. Hsu, Immersogeometric fluid–structure interaction modeling and simulation of transcatheter aortic valve replacement, Computer Methods in Applied Mechanics and Engineering 357 (2019) 112556. [12] F. Xu, S. Morganti, R. Zakerzadeh, D. Kamensky, F. Auricchio, A. Reali, T. J. Hughes, M. S. Sacks, M.-C. Hsu, A framework for designing patient-specific bioprosthetic heart valves using immersogeometric fluid–structure interaction analysis, International Journal for Numerical Methods in Biomedical Engineering 34 (4) (2018) e2938. [13] M.-C. Hsu, D. Kamensky, F. Xu, J. Kiendl, C. Wang, M. C. Wu, Y. Min, A. Reali, Y. Bazilevs, M. S. Sacks, T. J. Hughes, Dynamic and fluid–structure interaction simulations of bioprosthetic heart valves using parametric design with T-splines and Fung-type material models, Computational Mechanics 55 (6) (2015) 1211–1225. [14] D. Kamensky, M.-C. Hsu, D. Schillinger, J. A. Evans, A. Aggarwal, Y. Bazilevs, M. S. Sacks, T. J. Hughes, An immersogeometric variational framework for fluid–structure interaction: Application to bioprosthetic heart valves, Computational Methods in Applied Mechanics and Engineering 284 (2015) 1005–1053. [15] D. Kamensky, M.-C. Hsu, Y. Yu, J. A. Evans, M. S. Sacks, T. J. Hughes, Immersogeometric cardiovascular fluid–structure interaction analysis with divergence-conforming B-splines, Computational Methods in Applied Mechanics and Engineering 314 (2017) 408–472. [16] A. Nitti, J. Kiendl, A. Reali, M. D. de Tullio, An immersed-boundary/isogeometric method for fluid–structure interaction involving thin shells, Computational Methods in Applied Mechanics and Engineering 364 (2020) 112977. [17] R. Zorrilla, E. Soudah, An efficient procedure for the blood flow computer simulation of patient-specific aortic dissection, Computers in Biology and Medicine 179 (2024) 108832. doi:https://doi.org/10.1016/j.compbiomed.2024.108832. [18] F. de Prenter, C. V. Verhoosel, E. H. van Brummelen, M. G. Larson, S. Badia, Stability and Conditioning of Immersed Finite Element Methods: Analysis and Remedies, Archives of Computational Methods in Engineering 30 (2023) 3617–3656, open Access. doi:10.1007/ s11831-023-09851-6. [19] J. A. Nitsche, Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind, Abhandlungen aus dem Mathematischen Seminar der Universität Hamburg 36 (1971) 9–15. doi:10.1007/BF02995904. [20] M. Juntunen, R. Stenberg, Nitsche’s method for general boundary conditions, Mathematics of Computation 78 (267) (2009) 1353–1374. [21] I. Babuška, The finite element method with Lagrange multipliers, Numerische Mathematik 20 (1973) 179–192. doi:10.1137/0710027. [22] F. Brezzi, On the existence, uniqueness and approximation of saddle-point problems arising from Lagrangian multipliers, RAIRO - Analyse Numérique 8 (R2) (1974) 129–151. [23] I. Babuška, The finite element method with penalty, Numerische Mathematik 20 (1973) 179–192. doi:10.1137/0710027. [24] T. J. R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Prentice Hall, 1987. [25] K. Höllig, U. Reif, J. Wipper, Weighted Extended B-spline Approximation of Dirichlet Problems, SIAM Journal on Numerical Analysis 39 (2) (2001) 442–462. [26] K. Höllig, C. Apprich, A. Streit, Introduction to the Web-Method and its Applications, Advances in Computational Mathematics 23 (2005) 215–237. [27] R. A. K. Sanches, P. B. Bornemann, F. Cirak, Immersed B-spline (I-spline) Finite Element Method for Geometrically Complex Domains, Computer Methods in Applied Mechanics and Engineering 200 (13–16) (2011) 1432–1445. [28] J. Parvizian, A. Düster, E. Rank, Finite cell method: h- and p-extension for embedded domain problems in solid mechanics, Computational Mechanics 41 (1) (2007) 121–133. [29] A. Düster, J. Parvizian, Z. Yang, E. Rank, The finite cell method for three-dimensional problems of solid mechanics, Computer Methods in Applied Mechanics and Engineering 197 (45–48) (2008) 3768–3782. [30] E. Rank, M. Ruess, S. Kollmannsberger, D. Schillinger, A. Düster, Geometric modeling, isogeometric analysis and the finite cell method, Computer Methods in Applied Mechanics and Engineering 249–252 (2012) 104–115. [31] D. Schillinger, L. Dedè, M. A. Scott, J. A. Evans, M. J. Borden, E. Rank, T. J. Hughes, An isogeometric design-through-analysis methodology based on adaptive hierarchical refinement of NURBS, immersed boundary methods, and T-spline CAD surfaces, Computer Methods in Applied Mechanics and Engineering 249–252 (2012) 116–150. [32] D. Schillinger, M. Ruess, N. Zander, Y. Bazilevs, A. Düster, E. Rank, Small and large deformation analysis with the p- and B-spline versions of the Finite Cell Method (2012). [33] M. Ruess, D. Schillinger, T. Duenser, E. Rank, Weakly enforced essential boundary conditions for NURBS-embedded and trimmed NURBS geometries on the basis of the finite cell method, International Journal for Numerical Methods in Engineering 95 (10) (2013) 811–846. [34] L. C. Foucard, F. J. Vernerey, An X-FEM-based numerical–asymptotic expansion for simulating a Stokes flow near a sharp corner, International Journal for Numerical Methods in Engineering 102 (2) (2015) 79–98. [35] E. Burman, S. Claus, P. Hansbo, M. G. Larson, A. Massing, CutFEM: Discretizing geometry and partial differential equations, International

Journal for Numerical Methods in Engineering 104 (7) (2015) 472–501. [36] A. Massing, B. Schott, W. Wall, A stabilized Nitsche cut finite element method for the Oseen problem, Computer Methods in Applied Mechanics and Engineering 328 (2018) 262–300. doi:https://doi.org/10.1016/j.cma.2017.09.003. [37] M. Winter, B. Schott, A. Massing, W. Wall, A Nitsche cut finite element method for the Oseen problem with general Navier boundary conditions, Computer Methods in Applied Mechanics and Engineering 330 (2018) 220–252. doi:10.1016/j.cma.2017.10.023. URL http://dx.doi.org/10.1016/j.cma.2017.10.023 [38] C. S. Peskin, Numerical analysis of blood flow in the heart, Journal of computational physics 25 (3) (1977) 220–252. [39] M. Breitenberger, A. Apostolatos, P. Bucher, R. Wüchner, K.-U. Bletzinger, Analysis in computer aided design: Nonlinear isogeometric B-Rep analysis of shell structures, Computational Methods in Applied Mechanics and Engineering 284 (2015) 401–457. [40] T. Teschemacher, A. M. Bauer, T. Oberbichler, M. Breitenberger, R. Rossi, R. Wüchner, K.-U. Bletzinger, Realization of CAD-integrated shell simulation based on isogeometric B-Rep analysis, Advanced Modeling and Simulation in Engineering Sciences 5 (2018) 19. [41] T. Teschemacher, A. M. Bauer, R. Aristio, M. Meßmer, R. Wüchner, K.-U. Bletzinger, Concepts of data collection for the CAD-integrated isogeometric analysis, Engineering with Computers 38 (6) (2022) 5675–5693. [42] M. Meßmer, T. Teschemacher, L. F. Leidinger, R. Wüchner, K.-U. Bletzinger, Efficient CAD-integrated isogeometric analysis of trimmed solids, Computational Methods in Applied Mechanics and Engineering 400 (2022) 115584. [43] M. Meßmer, S. Kollmannsberger, R. Wüchner, K.-U. Bletzinger, Robust numerical integration of embedded solids described in boundary representation, Computer Methods in Applied Mechanics and Engineering 419 (2024) 116670. [44] A. Main, G. Scovazzi, The shifted boundary method for embedded domain computations. Part I: Poisson and Stokes problems, Journal of Computational Physics 372 (2018) 972–995. [45] A. Main, G. Scovazzi, The shifted boundary method for embedded domain computations. Part II: Linear advection–diffusion and incompressible Navier–Stokes equations, Journal of Computational Physics 372 (2018) 996–1026. [46] N. Antonelli, R. Aristio, A. Gorgi, R. Zorrilla, R. Rossi, G. Scovazzi, R. Wüchner, The Shifted Boundary Method in Isogeometric Analysis, Computer Methods in Applied Mechanics and Engineering 430 (2024) 117228. [47] R. Zorrilla, R. Rossi, G. Scovazzi, C. Canuto, A. Rodrı́guez-Ferran, A shifted boundary method based on extension operators, Computer Methods in Applied Mechanics and Engineering 421 (2024) 116782. [48] J. H. Collins, A. Lozinski, G. Scovazzi, A penalty-free Shifted Boundary Method of arbitrary order, Computer Methods in Applied Mechanics and Engineering 417 (2023) 116301, a Special Issue in Honor of the Lifetime Achievements of T. J. R. Hughes. [49] S. Löhnert, A stabilization technique for the regularization of nearly singular extended finite elements, Computational Mechanics 54 (2014) 523–533. doi:10.1007/s00466-014-1003-7. [50] W. Garhuom, K. Usman, A. Düster, An eigenvalue stabilization technique to increase the robustness of the finite cell method for finite strain problems, Computational Mechanics 69 (2022) 1225–1240. doi:10.1007/s00466-022-02140-7. [51] S. Eisenträger, L. Radtke, W. Garhuom, S. Löhnert, A. Düster, D. Juhre, D. Schillinger, An eigenvalue stabilization technique for immersed boundary finite element methods in explicit dynamics, Computers & Mathematics with Applications 166 (2024) 129–168. [52] E. Burman, Ghost penalty, Comptes Rendus Mathematique 348 (21–22) (2010) 1217–1220. [53] P. Dadvand, R. Rossi, E. Oñate, An Object-oriented Environment for Developing Finite Element Codes for Multi-disciplinary Applications, Archives of Computational Methods in Engineering 17 (3) (2010) 253–297. [54] P. Dadvand, R. Rossi, M. Gil, X. Martorell, J. Cotela, E. Juanpere, S. Idelsohn, E. Oñate, Migration of a generic multi-physics framework to HPC environments, Computers & Fluids 80 (2013) 301–309. [55] R. Sartorti, A. Düster, Data transfer within a finite cell remeshing approach applied to large deformation problems, Computational Mechanics (2024). [56] M. J. D. Powell, Radial Basis Functions for Multivariable Interpolation: A Review, Vol. 8, 1992, pp. 579–606. [57] M. D. Buhmann, Radial Basis Functions, Acta Numerica 9 (2000) 1–38. [58] D. Levin, The approximation power of moving least-squares, Mathematics of Computation 67 (224) (1998) 1517–1531. doi:10.1090/ S0025-5718-98-00928-0. [59] O. Afshar, et al., Moving Least Squares Method and its Improvement: A Concise Review, Journal of Computational Methods in Engineering 40 (3) (2021) 123–140. [60] T. J. Hughes, J. A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Computational Methods in Applied Mechanics and Engineering 194 (39-41) (2005) 4135–4195. [61] M. G. Cox, The numerical evaluation of B-splines, IMA Journal of Applied Mathematics 10 (2) (1972) 134–149. [62] C. D. Boor, On calculating with B-splines, Journal of Approximation Theory 6 (1) (1972) 50–62. [63] L. Piegl, W. Tiller, The NURBS Book, Springer Science & Business Media, 2012. [64] J. H. Collins, A. Lozinski, G. Scovazzi, A penalty-free shifted boundary method of arbitrary order, Computer Methods in Applied Mechanics and Engineering 417 (2023) 116301, a Special Issue in Honor of the Lifetime Achievements of T. J. R. Hughes. [65] E. Hairer, S. P. Nørsett, G. Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems, 2nd Edition, Vol. 8 of Springer Series in Computational Mathematics, Springer, 1993.

Record · ID 919482 · SHA-256 9d7b5c0180fc4d13
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.