**1. Introduction**

Ischemic heart disease is the first cause of death globally, accounting for 27% of fatalities in 2019 [1,2], with coronary atherosclerosis being the cause of most myocardial infarctions [3]. Atherosclerotic plaques (within the coronaries) result from a complex inflammatory process starting with the accumulation and retention of low-density lipoprotein within the intima. The result is a build-up of material (cholesterol and other lipid compositions) within the wall layers, producing stenosis and blood flux reduction in the vessel. Typically, a patient presents either stable or unstable (low or high risk of rupture) plaque. The fast distinction between these two groups is crucial regarding the treatment and disposition of the patient [4]. Thus, the need for patient-specific approaches is self-evident. It is here where computation-supported decision-making processes play a crucial role. This work contributes with a new approach to modeling two-dimensional coronary sections undergoing uniform physiological internal pressure in a quasi-static regime. Holzapfel et al. [5] showed that the pertinence and accuracy of the results depend on the method used to define material properties and to acquire in vivo patient-specific geometries. Typical methods for geometry

**Citation:** Gahima, S.; Díez, P.; Stefanati, M.; Rodríguez Matas, J.F.; García-González, A. An Unfitted Method with Elastic Bed Boundary Conditions for the Analysis of Heterogeneous Arterial Sections. *Mathematics* **2023**, *11*, 1748. https:// doi.org/10.3390/math11071748

Academic Editor: Fernando Simoes

Received: 19 February 2023 Revised: 22 March 2023 Accepted: 30 March 2023 Published: 6 April 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

acquisition include Magnetic Resonance Imaging [6,7], Computer Tomography [8,9], Optical Coherence Tomography [10,11], and Intravascular Ultrasound [12,13], among others. Hyperelastic [14,15], piece-wise homogeneous [16,17], incompressible [18,19], and plane strain [20,21] hypotheses characterize the two-dimensional models. In this work, we have used for the mechanical properties of all plaque components (i.e., normal vessel wall, loose matrix, calcification, and lipid core) a linear approximation of the stress–strain curve up to around 10% deformation reported by [22,23]. In particular, the rational behind using linear properties was to test the proposed approach, since it allows reducing the calculation times. However, recent studies use linear mechanical properties of the arterial tissue to perform clinical predictions based on geometrical and biomechanical markers obtained from finite element simulation [24]. In addition, the proposed methodology can directly accommodate a nonlinear hyper-elastic behavior of the tissues composing the coronary plaque.

One of the main problems when solving a finite element problem is to properly constrain the structure to prevent rigid body motion. The simplest approach to fixing the singularity and suppressing rigid-body modes consists in blocking three degrees of freedom (an isostatic condition with vertical displacements in two arbitrary nodes and horizontal displacement in another, for example). Thus, it is about setting essential (Dirichlet type) boundary conditions (BC) in arbitrarily selected nodes. In addition, blocking axial displacement and allowing free radial expansion [25,26] is another possibility. In addition, fixing two adjacent points [27,28], or the entire external boundary [29], or creating a soft and compressible layer surrounding the section with a zero-displacement BC constraining the latter [30,31] are examples of BCs used to suppress rigid body motions. All these methods struggle, to different degrees, to consider what surrounds the coronary section. Some of them depend on arbitrary choices (e.g., choosing the nodes where to impose BC). Few works have attempted to account for the influence of the surrounding tissue on the artery. In [32], an elastic bed boundary condition was applied along the coronary artery to simulate the effect of the myocardial tissue. In [33], the artery is surrounded by the cardiac wall for half of its circumference to simulate the coronary artery embedded into the myocardium, with the cardiac wall modeled using finite elements. We propose to assume the section to be surrounded by a material along its external boundary. This embedding matrix produces a linear elastic reaction, and it is described with an elastic bed coefficient *α* to be assessed depending on the stiffness of the surrounding medium.

Additionally, to improve the computational efficiency of the realistic calculation and analysis process, the numerical methodology proposed here implements the aforementioned elastic bed coefficient in an Immersed Boundary (IB) [34] framework with a generic description of the domain based on level sets. IB (combined with elastic bed BC) bases its simulations on a unique (background) mesh supporting the solutions corresponding to different configurations (different patients). It allows comparing solutions and opens the door to reduced-order modeling leading to fast simulations for different patient-specific geometries from initial medical images (with the same degrees of freedom), avoiding individual meshing and preprocessing steps. This methodology is motivated by its potential applicability with voxelized data [35], such as medical images. Via segmentation [36–39], it is possible to identify the different components' contours, and the IB performs the stress analysis on a refined voxelized background mesh to increase accuracy. In general, an IB approach allows seamless integration of structural analysis in a medical image processing pipeline.

The remainder of the paper is structured as follows: Section 2 describes the problem statement (Section 2.1), presents the level set approach (Sections 2.2 and 2.3), and includes the description of the IB framework (Section 2.4), emphasizing the details of the mathematical formulation required in the weak form of the problem (Section 2.5). Section 3 shows the results of the proposed methodology, finishing with a discussion and the main conclusions in Section 4.

### **2. Materials and Methods**

With a patient-specific application in mind, the discretization method to be used has to handle different configurations corresponding to various patients effortlessly. Moreover, input (the diversity of arterial cross-section geometries) and output (the solution in terms of deformation, strains, and stresses) data are to be expressed in homogeneous formats to ease the analysis and the possible application of reduced-order models. Here, level sets defined on a background mesh (discretizing a background domain ΩB) describe the diversity of the geometric configurations, that is, all the possible instances of the actual computational domain, Ω. It comes naturally to solve the problem with an unfitted approach. Specifically, it uses the background mesh not only to describe the geometry (actually the same background mesh for all the possible geometries) but also to compute the solution, following an IB methodology. Thus, one may prescind the conformal meshes adapted to the geometry that change from case to case. Note that to solve the problem with conformal finite elements, the mesh must be such that it tallies with Ω, matching the boundary *∂*Ω. Such an approach requires ad-hoc meshing algorithms, especially for convoluted geometries, and complicates comparing different configurations and their solutions.

### *2.1. Problem Statement*

Let the section occupy a region Ω ⊂ R<sup>2</sup> with boundary *∂*Ω. The intrinsic heterogeneity of arterial cross-sections is described by dividing Ω into different subdomains Ω1, Ω2,..., corresponding to homogeneous regions having different material properties i.e., normal vessel wall, loose matrix, calcification, and lipid core (see Figure 1). Without body forces, the equilibrium is governed by

$$\nabla \cdot \sigma(\mathbf{u}) = \mathbf{0} \quad \text{in } \Omega,\tag{1}$$

with boundary conditions

$$
\sigma(\mathbf{u}) \cdot \mathbf{\hat{n}} = \mathbf{t} \tag{2}
$$

$$
\sigma(\mathbf{u}) \cdot \mathbf{\hat{n}} = \kappa \mathbf{u} \quad (\kappa < 0) \tag{3}
$$

where *σ* is the Cauchy stress tensor and **u** is the displacement field; **t** is the surface traction, *α* is the elastic bed coefficient, and **n**ˆ is the outward unit normal to the boundary. Equation (3) represents the Robin boundary condition, physically corresponding to an elastic bed condition, simulating the surrounding tissue of the artery. The Neumann and elastic bed boundaries cover the whole boundary, i.e., *∂*Ω = Γ*N* ∪ Γ*R*.

**Figure 1.** Schematic description of Problem (1) in the Euclidean space. In particular, Ω = 24 *k*=1 Ω*k* with *∂*Ω = Γ*N* ∪ Γ*R* where Γ*N* and Γ*R* are depicted in red and cyan, respectively.

The weak form of Problem 1 (physically corresponding to the principle of virtual work) reads: find **u** ∈ H<sup>1</sup>(Ω) 2 such that

$$\int\_{\Omega} \sigma(\mathbf{u}) : \mathfrak{e}(\mathbf{v}) \, \mathrm{d}\Omega \, - \int\_{\Gamma\_R} \mathfrak{a} \mathbf{u} \cdot \mathbf{v} \, \mathrm{d}\Gamma = \int\_{\Gamma\_N} \mathbf{t} \cdot \mathbf{v} \, \mathrm{d}\Gamma,\tag{4}$$

for all **v** ∈ H<sup>1</sup>(Ω) 2 . H<sup>1</sup>(Ω) is the Sobolev space of order 1 on Ω; refer to [40] for details. Note that the test function **v** is also seen as a virtual infinitesimal displacement (a perturbation from the equilibrium configuration of the body) consistent with the imposed boundary displacements, and *ε*(**v**) = 1/2(*∇***v** + *<sup>∇</sup>***v**). The elastic bed BC (3) is an alternative to the standard practice of suppressing rigid-body motions by prescribing displacements at some arbitrarily selected points. As shown in the following, enforcing an isostatic scheme by prescribing point displacements and elastic bed BC produce similar results. We advocate the latter because the elastic bed BC includes physical information about the surrounding medium and does not require selecting arbitrary points to prescribe displacements. This is crucial for model order reduction, where one has to perform operations with the solutions of different configurations, and hence, they need to be comparable.

#### *2.2. Level Set Description of the Domain and Subdomains*

As introduced previously, the domain Ω is divided into *n* subdomains Ω*i*, *i* = 1, ... , *n*. The *n* subdomains cover Ω, that is

$$
\Omega = \bigcup\_{k=1}^{n} \Omega\_k. \tag{5}
$$

Level set functions implicitly describe the geometry of Ω and its subdomains in a unique framework. A background domain ΩB, having a simple geometry (here rectangle or square shape), is introduced to accommodate all possible instances of Ω, resulting in Ω ⊂ ΩB; see Figure 2A.

**Figure 2.** (**A**) The background domain Ω<sup>B</sup> and (**B**) (one of its possible) mesh T*h*(Ω<sup>B</sup>) (for more details regarding accurate estimations of the displacement fields at the interfaces, see Section 2.4). Inner and cut (by *∂*Ω) elements *T* ∈ T*h* are in blue and yellow, respectively.

A standard level set to describe Ω in Ω<sup>B</sup> is a continuous function *φ* taking values in Ω<sup>B</sup> such that *φ*(**x**) > 0 for **x** ∈ Ω and negative elsewhere. Thus, *φ*(**x**) = 0 for **x** ∈ *∂*Ω. Typically, *φ* is a signed distance to *∂*Ω [41,42]. For a configuration such as the one in Figure 3, with two non-connected parts of the boundary, Γ*N* and Γ*R*, it is convenient to describe Ω using two level sets to distinguish between the two. Thus, Ω is identified with *φ*(1) and *φ*(2) such that: *φ*(1)(**x**) = 0 on Γ*N*, and *φ*(2)(**x**) = 0 on Γ*R*. Both level sets are positive in Ω; see Figure 3B,C for an illustration. Note that one may recover a standard level set for Ω by just taking *φ* = *φ*(1)*φ*(2). Then, following the ideas in [43], new level set functions are introduced to

describe the *n* subdomains. Function *φ*(3) provides the information to identify Ω1 and distinguish it from the remainder subdomains. In particular, *φ*(3)(**x**) > 0 for **x** ∈ 2*nk*=<sup>2</sup> Ω*k* and is negative elsewhere in Ω (that is in Ω1). Similarly, *φ*(4) is positive in 2*nk*=<sup>3</sup> Ω*k* and negative in the remainder, that is in 2<sup>2</sup>*k*=<sup>1</sup> Ω*k*. The last hierarchical level set needed is *φ*(*n*+<sup>1</sup>) identifying Ω*<sup>n</sup>*−<sup>1</sup> because then Ω*n* is precisely the remainder (*φ*(*n*+<sup>1</sup>) > 0). The values of *φ*(*k*) with *k* = 3, ... , *n* + 1 outside Ω are not relevant. This is consistent with the hierarchical character of this approach. A visualization of the hierarchical level sets is illustrated in the panels of Figure 3 and summarized in Table 1.

**Table 1.** Level set-based criteria to classify a point **x** in Ω and its subdomains.


**Figure 3.** (**A**) Ω embedded in ΩB. Level set (**B**) *φ*(1) = 0 describes Γ*N*, and (**C**) *φ*(2) = 0 describes Γ*R*. (**D**) *φ*(3) = 0 describes the interfaces between the subset 2<sup>4</sup>*k*=<sup>2</sup> Ω*k* and the background domain. (**E**) *φ*(4) = 0 describes the interfaces between Ω3 ∪ Ω4 and Ω<sup>B</sup> \ Ω3 ∪ Ω4 and finally, (**F**) *φ*(5) = 0 describes the interface between Ω4 and the rest.

The approach described above is similar to the front-tracking method used to simulate multiphase flow with a fixed grid for the flow [44], with the difference that the front does not change with time, and therefore, the level set is computed only once at the beginning of the analysis.

#### *2.3. Discretization of the Level Set Functions in a Background Mesh*

With a finite element (FE) discretization of Ω<sup>B</sup> (see Figure 2), the level set approach is implemented. A tessellation T*h* of Ω<sup>B</sup> consisting of *n*e elements *Te*, *e* = 1, 2, ... , *n*e (*h* stands for the characteristic size of the elements) is introduced, such that Ω<sup>B</sup> = 2*<sup>n</sup>*e *<sup>e</sup>*=1 *Te*. The number of nodes in the mesh is denoted by *<sup>n</sup>*p and the corresponding shape functions are denoted by *Ni*(**x**), for *i* = 1, 2, ... , *<sup>n</sup>*p c[cyan]. Thus, each level set *φ*(*k*), *k* = 1, ... , *n* + 1, is represented in the background mesh as

$$\phi^{(k)}(\mathbf{x}) \approx \sum\_{i=1}^{n\_{\text{P}}} [\Phi^k]\_i N\_i(\mathbf{x}) \,\tag{6}$$

with **Φ***<sup>k</sup>* ∈ R*<sup>n</sup>*p being the vector of nodal values describing *φ*(*k*). In this framework, *n* + 1 vectors in R*<sup>n</sup>*p describe any geometrical configuration. This standardized representation allows for dimensionality reduction given the variance in the population of input samples (each corresponding to a different patient). To ease the task of the machine learning algorithms to be used in dimensionality reduction, standard geometric normalizations are performed previously to store the information in the discrete level set format. For instance, all samples are centered (their barycenter is translated to the origin of coordinates) and rotated such that the principal axes of inertia are parallel to the coordinate axes.

#### *2.4. Unfitted Approach: Solving the Problem in the Background Mesh*

The framework for approximating the level set over T*h*(Ω<sup>B</sup>) just described is used to solve the original problem (4) using an unfitted approach based on the ideas of the Immersed Boundary Method (IBM). Thus, the displacement field **u**(**x**) is approximated in the background mesh using a standard FE approximation, namely

$$\mathbf{u}(\mathbf{x}) \approx \sum\_{i=1}^{n\_{\mathrm{P}}} \mathbf{U}\_{i} N\_{i}(\mathbf{x}),\tag{7}$$

with **U***i* ∈ R<sup>2</sup> being the displacement vector in node *i*. All vectors **U***i*, *i* = 1, 2, ... , *<sup>n</sup>*p, are collected in the standard vector of nodal displacements **U** ∈ R2*n*p . Using the Galerkin strategy to solve Equation (4) results in a linear system of equations for **U**:

$$[\mathbf{K} + \mathbf{M}]\mathbf{U} = \mathbf{F},\tag{8}$$

where matrices **K** and **M** in R2*n*p×2*n*p are the discrete counterparts of the two bilinear forms in the left-hand side of Equation (4) and **F** ∈ R2*n*p is the discretization of the linear form in the right-hand side.

Note that a node *i* in the mesh is represented by the degrees of freedom = 2(*i* − 1) + 1 and +1 in **U**, and some other node *j* is represented by ˜ = 2(*j* −<sup>1</sup>) +1 and ˜ +1. Assuming these relations, some illustrative examples of the expressions for the corresponding entries in the matrices and the right-hand-side vector are given below

$$[\mathbf{K}]\_{\ell\ell} = \int\_{\Omega} \sigma \left( \begin{bmatrix} \mathbf{N}\_{\mathrm{i}}(\mathbf{x}) \\ 0 \end{bmatrix} \right) : \mathfrak{e} \left( \begin{bmatrix} \mathbf{N}\_{\mathrm{j}}(\mathbf{x}) \\ 0 \end{bmatrix} \right) \mathrm{d}\Omega \text{ } ; \ [\mathbf{K}]\_{\ell\bar{\ell}+1} = \int\_{\Omega} \mathfrak{e} \left( \begin{bmatrix} \mathbf{N}\_{\mathrm{i}}(\mathbf{x}) \\ 0 \end{bmatrix} \right) : \mathfrak{e} \left( \begin{bmatrix} 0 \\ \mathbf{N}\_{\mathrm{j}}(\mathbf{x}) \end{bmatrix} \right) \mathrm{d}\Omega$$

$$[\mathbf{M}]\_{\ell\bar{\ell}} = -\int\_{\Gamma\_{\mathrm{R}}} \mathfrak{e} \begin{bmatrix} \mathbf{N}\_{\mathrm{i}}(\mathbf{x}) \\ 0 \end{bmatrix} \cdot \begin{bmatrix} \mathbf{N}\_{\mathrm{j}}(\mathbf{x}) \\ 0 \end{bmatrix} \mathrm{d}\Gamma \text{ } ; \ [\mathbf{M}]\_{\ell\bar{\ell}+1} = -\int\_{\Gamma\_{\mathrm{R}}} \mathfrak{e} \begin{bmatrix} \mathbf{N}\_{\mathrm{i}}(\mathbf{x}) \\ 0 \end{bmatrix} \cdot \begin{bmatrix} 0 \\ \mathbf{N}\_{\mathrm{j}}(\mathbf{x}) \end{bmatrix} \mathrm{d}\Gamma$$

$$\left[\mathbf{F}\right]\_{\ell} = \int\_{\Gamma\_N} \mathbf{t} \cdot \begin{bmatrix} N\_i(\mathbf{x}) \\ 0 \end{bmatrix} \mathrm{d}\Gamma \; ; \; \left[\mathbf{F}\right]\_{\ell+1} = \int\_{\Gamma\_N} \mathbf{t} \cdot \begin{bmatrix} 0 \\ N\_i(\mathbf{x}) \end{bmatrix} \mathrm{d}\Gamma.$$

Note that all the integrals in the expressions above are defined in Ω, Γ*R*, and Γ*N*, and not in the background domain Ω<sup>B</sup> where the FE functions *Ni*(**x**) are supported. In particular, evaluating the local contributions (the integrals are restricted to some element *Te*) requires identifying whether an element intersects Γ*R* or Γ*N*. Thus, the main implementation challenge of the unfitted approach is classifying the elements T*h* of Ω<sup>B</sup> inside Ω, those outside, and those crossed by the interfaces. For a given configuration, the geometrical information is encoded in the level sets, as described in Section 2.2. This allows elaborating a list of the elements in T*h* that are completely inside Ω, namely IΩ such that if *e* ∈ IΩ, then *Te* ⊂ Ω. Similarly, lists <sup>I</sup>Γ*R* and <sup>I</sup>Γ*N* are such that if *e* ∈ <sup>I</sup>Γ*R* then *Te* 4 Γ*R* = ∅, and if *e* ∈ <sup>I</sup>Γ*N* then *Te* 4 Γ*N* = ∅. Figure 4 shows an example of such classification. The elements indexed in these three lists are *active*, meaning that they play a role in the solution for the configuration described by the level sets. Thus, *Te* is said to be active if *e* ∈ IΩ 2 <sup>I</sup>Γ*R* 2 <sup>I</sup>Γ*N* . Accordingly, all the nodes belonging to active elements are active nodes since the corresponding degrees of freedom are the unknowns of (8) (the non-active nodes are to be eliminated from the system).

**Figure 4.** (**A**) Elements *Te* for *e* ∈ IΩ 2 <sup>I</sup>Γ*R* 2 <sup>I</sup>Γ*N* , are colored in violet (*e* ∈ IΩ), magenta (*e* ∈ <sup>I</sup>Γ*N* ), and cyan (*e* ∈ <sup>I</sup>Γ*R* ), being Γ*N* and Γ*R* the black lines. The square background domain Ω<sup>B</sup> (2.5 × 2.5 mm2) is meshed with *<sup>n</sup>*p = 100<sup>2</sup> nodes and *n*e = 2 × 99<sup>2</sup> elements. Close-ups for better illustration in panels (**C**), and (**D**). Panel (**B**) illustrates that in the elements crossed by the boundary, the quadrature is enriched to avoid having no integration points in the part of the element outside Ω. This suggests using in these elements closed quadratures (as the third-degree closed Newton–Cotes quadrature [45]).

The computation of the elementary contributions to the stiffness matrix **K** is standard for the elements completely inside *Te* for *e* ∈ IΩ (violet elements in Figure 4A). The only particular feature to be accounted for is that the material properties of each Gauss point

in the numerical quadrature belong to a subdomain Ω*k*. With the values of the level sets interpolated at the Gauss point, following the classification described in Table 1, the material properties are quickly recovered. Note that for the example shown in Figure 4A, only two level sets, *φ*(1) and *φ*(2), are required. In the elements crossed by the interfaces Γ*R* and Γ*N*, the integration has to exclude the part of the domain outside Ω. There, a more refined quadrature is used, and null material properties are assigned to the integration points outside Ω. A closed quadrature is preferred to avoid accounting for integration points inside elements with a small portion inside Ω. These geometric checks are performed by setting a tolerance and considering that the distance to the interface is zero when it is below this value. Computing elementary contributions to matrix **M** and vector **F** requires integrating within the portion of Γ*R* or Γ*N* in the element *Te*. Thus, for *e* ∈ <sup>I</sup>Γ*N* , *Te* intersects Γ*N* and contributes to **F**. Analogously, for *e* ∈ <sup>I</sup>Γ*R* , *Te* intersects Γ*R* and contributes to **M**.

For *e* ∈ <sup>I</sup>Γ*N*, the elementary contribution from element *Te* to **F** requires computing

$$\int\_{\Gamma\_N \cap \prod^T \mathbf{t}} [\mathbf{t}]\_1 N\_{\mathbf{i}}(\mathbf{x}) d\Gamma \text{ and } \int\_{\Gamma\_N \cap \prod^T \mathbf{t}} [\mathbf{t}]\_2 N\_{\mathbf{i}}(\mathbf{x}) d\Gamma\_{\mathbf{i}} \tag{9}$$

for all the nodes *i* in element *Te*. If the load corresponds to a pressure *p* applied in the internal wall, then **t** = −*p***n**ˆ, recalling that **n**ˆ = [*<sup>n</sup>*1 *<sup>n</sup>*2] is the outward unit normal. Thus, [**t**]1 = −*pn*<sup>1</sup> and [**t**]2 = −*pn*2. This integral, as it is standard in the FE practice, is computed in a reference element (for linear triangles, it is handy using the triangle with vertices (0, <sup>0</sup>), (1, 0) and (0, <sup>1</sup>), see Figure 5), where the shape functions are defined (and available in their analytical expressions) in terms of the reference coordinates (*ξ*, *η*), namely *<sup>N</sup>*<sup>ˆ</sup>(*i*)(*ξ*, *η*), for (*i*) = 1, 2, 3. Mesh connectivity provides the link between the local numbering of the node inside the element, (*i*) (from 1 to 3 in the case of linear triangles), and the global numbering *i* (from 1 to *n*p). Since *Te* is crossed by Γ*N*, it is important to identify the entry and exit points in the element, that is the points {*PI*, *PI I*} = Γ*N* 4 *∂Te*; see Figure 5. This task is performed while identifying the elements in <sup>I</sup>Γ*N* , and it is straightforward after the nodal values of *φ*(1). Recall that *φ*(1)(**x**) = 0 for **x** ∈ Γ*N*. Same rationale works for Γ*R*, using *φ*(2). A quadrature is required to integrate along the segmen<sup>t</sup> *PIPI I* (a portion of Γ*N*). Here, a Simpson quadrature is adopted and involves computing the values of the function to be integrated on the endpoints of the interval and in the midpoint, *P*m; see Figure 5. The general expression for Simpson quadrature to approximate the integral of some function *ψ* reads:

$$\int\_{P\_I}^{P\_{II}} \psi \,\mathrm{d}\Gamma \approx \frac{|P\_I P\_{II}|}{6} (\psi(P\_I) + 4\psi(P\_{\mathrm{m}}) + \psi(P\_{II})),\tag{10}$$

where |*PIPI I*| is the length of the interval *PIPI I*. Thus, computing the terms in (9) requires obtaining the values of *Ni* in the three points *PI*, *PI I* and *P*m. These values are easily obtained after their coordinates in the reference element, (*ξ<sup>I</sup>*, *ηI*), (*ξI I*, *ηI I*) and (*ξ*m, *η*m). Then, it suffices using the quadrature given in (10) for *ψ*(**x**) = −*pn*1*Ni*(**x**) and *ψ*(**x**) = −*pn*2*Ni*(**x**) to obtain the horizontal and vertical components of the nodal forces (on node *i* from element *e*).

**Figure 5.** Element *Te* (right) described in the Cartesian coordinate system (*<sup>x</sup>*1, *<sup>x</sup>*2) is mapped into reference element (left) described by (*ξ*, *η*) coordinates. Outward normal **n**ˆ to the portion of Γ*N* (respectively, Γ*R*) in *Te* is to be determined in the Cartesian framework. The two points where the interface meets the boundary of *Te* (entry and exit points, here denoted by *I* and *I I*) are necessary to numerically integrate the coefficients of **F** (respectively **M**) along Γ*N* (respectively, Γ*R*).

Similarly, for *e* ∈ <sup>I</sup>Γ*R* , the elementary contribution from *Te* to **M** requires computing terms of the form

$$\int\_{\Gamma\_R \cap \prod T\_\varepsilon} a N\_i(\mathbf{x}) N\_j(\mathbf{x}) d\Gamma. \tag{11}$$

This is achieved by taking *ψ*(**x**) = *αNi*(**x**)*Nj*(**x**) and using the quadrature (10) accordingly. In the unfitted solution, the degrees of freedom corresponding to nodes that do not belong to any of the active elements (those with index *e* in 2 <sup>I</sup>Γ*R* 2 <sup>I</sup>Γ*N* ) must be removed from system (8).

The number of active nodes (that is the number of nodes in the active elements) is denoted by *n*act and indicates the measure of the size of the system to be solved. Note that in conformal FE, *n*act is the number of nodes in the mesh. On the other hand, in an unfitted approach, *n*act < *<sup>n</sup>*p.

The proposed methodology has been entirely implemented in Matlab R2022b, The MathWorks Inc, and executed in a 3.2 GHz Apple M1 with 8 GB RAM.

### *2.5. Validating the Methodolgy*

The solution of an infinite linear elastic solid with a cylindrical cavity of radius *rint* subjected to internal pressure *p* (Figure 6A) is considered.

**Figure 6.** Verification problem. (**A**) Infinite solid with a cylindrical cavity of radius *rint* subjected to internal pressure *p*; (**B**) Infinite cylinder of inner radius *rint* and external radius *rext* subjected to an internal pressure *p* and elastic bed boundary conditions at *rext*.

Assuming cylindrical coordinates, the analytical solution of the problem is [46]

$$
\mu\_r = p \frac{r\_{int}^2}{2\mu r}, \quad \mu\_\theta = \mu\_z = 0,\tag{12}
$$

where *μ* = *μ*(*<sup>E</sup>*, *v*) is the shear modulus of the material, *E* and *ν* are its Young's modulus and Poisson ratio, respectively, and *r* ≥ *rint* is the radial coordinate. The strain and stress fields are obtained with the constitutive equations for a linear elastic solid as

$$
\varepsilon\_{rr} = \frac{\partial u\_r}{\partial r} = -\frac{p}{2\mu} \left(\frac{r\_{int}}{r}\right)^2,\tag{13}
$$

$$
\varepsilon\_{\theta\theta} = \frac{u\_r}{r} = \frac{p}{2\mu} \left(\frac{r\_{int}}{r}\right)^2,\tag{14}
$$

$$
\varepsilon\_{zz} = 0,\tag{15}
$$

$$
\sigma\_{rr} = -p \left(\frac{r\_{int}}{r}\right)^2,\tag{16}
$$

$$
\sigma\_{\theta\theta} = p \left( \frac{r\_{int}}{r} \right)^2,\tag{17}
$$

$$
\sigma\_{\varpi} = 0.\tag{18}
$$

This problem is equivalent to that of an infinite cylinder of internal radius *rint* and external radius, *rext*, subjected to internal pressure, *p* (applied at *r* = *rint*), and elastic bed BC on *r* = *rext* (see Figure 6B), with the *ballast* coefficient *α* given as

$$
\alpha = -\frac{2\mu}{r\_{ext}},
\tag{19}
$$

obtained from (3) together with (12), and (16)–(18).

With this problem at hand, the accuracy of the methodology is quantified in terms of local and global quantities. The displacement field is a local quantity of accuracy, and the total deformation energy (TDE) is used as a global metric to assess the convergence for the numerical solution. The TDE for the infinite cylinder in Figure 6B is given by

$$\begin{split} \text{TDE} &= \frac{1}{2} \int\_{0}^{2\pi} \int\_{r\_{int}}^{r\_{ext}} (\sigma\_{rr} \varepsilon\_{rr} + \sigma\_{\theta\theta} \varepsilon\_{\theta\theta}) r \, \text{d}r \, \text{d}\theta = \\ &= \pi \int\_{r\_{int}}^{r\_{ext}} \left[ \left( \frac{p^2 r\_{int}^4}{2\mu r^4} \right) + \left( \frac{p^2 r\_{int}^4}{2\mu r^4} \right) \right] r \, \text{d}r = \frac{\pi p^2 r\_{int}^4}{2\mu} \left( \frac{1}{r\_{int}^2} - \frac{1}{r\_{ext}^2} \right) > 0. \end{split} \tag{20}$$

Note that (20) only accounts for the elastic energy stored in the cylinder and not in the elastic bed. Numerically, the TDE is calculated as

$$\text{TDE}\_{\text{num}} = \frac{1}{2} \mathbf{U}^{\top} \mathbf{K} \,\mathbf{U},\tag{21}$$

where **U** corresponds to the displacement vector of the active nodes in the background mesh, and **K** is the stiffness matrix associated with the active elements in the background mesh.
