A ROM-based BDDC solver for unfitted p-FEM level-set-based two-dimensional lattice structures

The starting point for the geometric definition of an individual cell is the untrimmed parametric domain, $\hat{\Pi}$, which is defined as the unit square:

$$\hat{\Pi}=[0,1]\times[0,1].$$

(1)

This domain represents the set where the curvilinear coordinates $(\xi_{1},\xi_{2})$ take values. These coordinates are used to explicitly define the map $\mathcal{F}:\hat{\Pi}\to\mathbb{R}^{2}$ as:

$$\boldsymbol{\mathcal{F}}\left(\xi_{1},\xi_{2}\right)=\mathcal{F}_{i}\left(\xi_{1},\xi_{2}\right)\boldsymbol{e}_{i}\quad\text{with}\quad i=1,2,$$

(2)

where $\boldsymbol{e}_{i}$ represents the vectors of the standard Euclidean basis. The four edges of the parametric domain are denoted as: the left vertical edge $\partial\hat{\Pi}_{1b}=\{(0,\xi)\}$, the right vertical edge $\partial\hat{\Pi}_{1t}=\{(1,\xi)\}$, the bottom horizontal edge $\partial\hat{\Pi}_{2b}=\{(\xi,0)\}$, and the top horizontal edge $\partial\hat{\Pi}_{2t}=\{(\xi,1)\}$, where $\xi\in[0,1]$. The untrimmed physical domain in $\mathbb{R}^{2}$ is the image of the unit square $\hat{\Pi}$ under this map:

$${\Pi}=\boldsymbol{\mathcal{F}}(\hat{\Pi}).$$

(3)

Figure 1 shows an example of such construction. At this stage, no particular assumptions are made concerning the function $\boldsymbol{\mathcal{F}}$. However, a discussion regarding its required regularity is provided in the remainder of this section.

The mapping introduced above is combined with an immersed boundary approach to define the active portion of the cell. Specifically, the trimmed parametric domain is defined as:

$$\hat{\Omega}=\left\{\left(\xi_{1},\xi_{2}\right)\in\hat{\Pi}:\phi_{0}\left(\xi_{1},\xi_{2}\right)<\mu\left(\xi_{1},\xi_{2}\right)\right\}\;,$$

(4)

where $\phi_{0}$ and $\mu$ are the level-set and threshold functions, respectively. Or equivalently, by introducing the function $\phi\left(\xi_{1},\xi_{2}\right)=\phi_{0}\left(\xi_{1},\xi_{2}\right)-\mu\left(\xi_{1},\xi_{2}\right)$ as

$$\hat{\Omega}=\left\{\left(\xi_{1},\xi_{2}\right)\in\hat{\Pi}:\phi\left(\xi_{1},\xi_{2}\right)<0\right\}\;.$$

(5)

While $\phi_{0}$ can be any arbitrary implicit function, for illustrative purposes, in this work it is selected from the class of triply periodic minimal surfaces (TPMS) adapted for 2D applications, which does not constitute a limitation for the proposed method. Common examples are defined in Table 1 and depicted in Figure 2.

Typically, the threshold is kept uniform, therefore determining which isoline of $\phi_{0}$ is adopted as domain boundary. However, to introduce more flexibility in the design of the cells, this work allows $\mu$ to vary linearly as follows:

$$\mu\left(\xi_{1},\xi_{2}\right)=\sum_{i=1}^{4}{l}_{i}\left(\xi_{1},\xi_{2}\right)\mu_{i}\;,$$

(6)

where ${l}_{i}$ are the classical first-degree Lagrangian shape functions, and the threshold parameters $\mu_{i}$ correspond to the nodal values at the bottom-left, bottom-right, top-left, and top-right corners of $\hat{\Pi}$. As for its untrimmed counterpart, the trimmed physical domain is obtained through the map:

$$\Omega=\boldsymbol{\mathcal{F}}(\hat{\Omega})\;.$$

(7)

The boundary of the parametric domain is identified as two separate portions: the conformal portion lying on the untrimmed parametric domain boundary, and the trimmed portion that is internal to the domain:

	$\displaystyle\partial\hat{\Omega}_{c}$	$\displaystyle=\left\{(\xi_{1},\xi_{2})\in\partial\hat{\Pi}:\phi_{0}\left(\xi_{1},\xi_{2}\right)<\mu\left(\xi_{1},\xi_{2}\right)\right\}\;,$		(8a)
	$\displaystyle\partial\hat{\Omega}_{t}$	$\displaystyle=\left\{(\xi_{1},\xi_{2})\in\mathrm{int}(\hat{\Pi}):\phi_{0}\left(\xi_{1},\xi_{2}\right)=\mu\left(\xi_{1},\xi_{2}\right)\right\}\;,$		(8b)

where $\mathrm{int}(\bullet)$ denotes the internal part of the domain $\bullet$, such that $\partial\hat{\Omega}=\partial\hat{\Omega}_{c}\cup\partial\hat{\Omega}_{t}$ and $\partial\hat{\Omega}_{c}\cap\partial\hat{\Omega}_{t}=\emptyset$. In Figure 1, an example of a trimmed domain is shown in parametric and physical coordinates and its correspondent threshold parameters.

The discussion so far has focused on the geometric definition of a single cell. To describe a complete structure composed of multiple cells, we now introduce a superscript $(i_{c})$ to distinguish all quantities associated with a specific cell, where $i_{c}\in\{1,\dots,n_{c}\}$ is the cell index and $n_{c}$ is the total number of cells. Therefore, the untrimmed and trimmed physical domains of the overall structure are defined, respectively, as:

$${\Pi}=\mathrm{int}\left(\bigcup_{i_{c}=1}^{n_{c}}\bar{\Pi}^{(i_{c})}\right),\quad\mathrm{and}\quad{\Omega}=\mathrm{int}\left(\bigcup_{i_{c}=1}^{n_{c}}\bar{\Omega}^{(i_{c})}\right)\;,$$

(9)

where $\bar{\bullet}$ denotes the closure of the domain $\bullet$. While this procedure can generate complex topologies and is not restricted to tensor-product grids, the cells must satisfy certain compatibility conditions at both the untrimmed and trimmed levels.

First, for untrimmed compatibility, the physical mapping of connected edges must coincide. If the $i_{c_{1}}$-th and $i_{c_{2}}$-th cells are joined at their respective edges $\partial\hat{\Pi}_{i}^{(i_{c_{1}})}$ and $\partial\hat{\Pi}_{j}^{(i_{c_{2}})}$, then:

$$\boldsymbol{\mathcal{F}}^{(i_{c_{1}})}\left(\partial\hat{\Pi}_{i}^{(i_{c_{1}})}\right)=\boldsymbol{\mathcal{F}}^{(i_{c_{2}})}\left(\partial\hat{\Pi}_{j}^{(i_{c_{2}})}\right).$$

(10)

where $i,j\in\{1b,1t,2b,2t\}$.

Second, for trimmed compatibility, the trimming function $\phi$ must be continuous across the common interface

$$\phi^{(i_{c_{1}})}(\xi)=\phi^{(i_{c_{2}})}(\xi),$$

(11)

where $\xi\in[0,1]$ is an auxiliary curvilinear coordinate that parameterizes the interface. This condition, in turn, is satisfied if both the level-set and threshold functions are independently equal: $\phi_{0}^{(i_{c_{1}})}(\xi)=\phi_{0}^{(i_{c_{2}})}(\xi)$ and $\mu^{(i_{c_{1}})}(\xi)=\mu^{(i_{c_{2}})}(\xi)$. In the particular case of TPMS geometries, the first condition ($\phi_{0}$) is automatically satisfied by their intrinsic periodicity, provided that the edge pairs $(i,j)$ in Equation (10) are restricted to the admissible set: ${(1b,1t),(1t,1b),(2b,2t),(2t,2b)}$. The second equality ($\mu$) is satisfied by simply enforcing identical threshold parameters at the common nodes of neighboring cells.

To illustrate, in Figure 3, the two adjacent cells $\Pi^{(1)}$ and $\Pi^{(2)}$ are connected at their edges $\partial{\Pi}_{1t}^{(1)}$ and $\partial{\Pi}_{1b}^{(2)}$. Compatibility of the untrimmed interface is ensured because the maps are chosen such that the following relationship holds: $\boldsymbol{\mathcal{F}}^{(1)}(1,\xi)=\boldsymbol{\mathcal{F}}^{(2)}(0,\xi)$ with $\xi\in[0,1]$. Furthermore, compatibility of the trimmed interface is achieved by imposing $\mu^{(1)}_{2}=\mu^{(2)}_{1}$, and $\mu^{(1)}_{4}=\mu^{(2)}_{3}$.

Notably, a spline surface can be easily adopted as the untrimmed physical domain, as it naturally satisfies the compatibility condition in Equation (10). The only preprocessing required is a Bézier extraction for each of its elements, which in turn defines the mapping for each cell in the lattice structure.

The variational statement for two-dimensional linear elasticity reads: find the displacement field $\boldsymbol{u}\in[\mathcal{H}^{1}]^{2}_{\bar{\boldsymbol{u}}}$ such that

$$\mathfrak{L}_{int}(\boldsymbol{v},\boldsymbol{u})=\mathfrak{L}_{ext}(\boldsymbol{v})\;,\quad\quad\forall\boldsymbol{v}\in[\mathcal{H}^{1}]^{2}_{\boldsymbol{0}}\;,$$

(12)

where $[\mathcal{H}^{1}]^{2}_{\bar{\boldsymbol{u}}}$ and $[\mathcal{H}^{1}]^{2}_{\boldsymbol{0}}$ represent the subspaces of $[\mathcal{H}^{1}]^{2}$ for which Dirichlet boundary conditions are satisfied, in the first case, and where the test function $\boldsymbol{v}$ is null on the Dirichlet boundary, in the second one. It is further assumed that only pure Dirichlet or Neumann boundary conditions are considered, meaning that mixed boundary conditions, where Dirichlet is applied to one component of the displacement and Neumann to another, are excluded for simplicity. We note that supporting such conditions is a straightforward extension of the present framework.

The bilinear and the linear forms in Equation (12) are also referred to as virtual work of the internal and external forces, respectively, and are defined as:

		$\displaystyle\mathfrak{L}_{int}(\boldsymbol{v},\boldsymbol{u})=\sum_{i_{c}=1}^{n_{c}}\mathfrak{L}_{int}^{(i_{c})}(\boldsymbol{v},\boldsymbol{u}),$		(13a)
		$\displaystyle\mathfrak{L}_{ext}(\boldsymbol{v})=\sum_{i_{c}=1}^{n_{c}}\mathfrak{L}_{ext}^{(i_{c})}(\boldsymbol{v})\;,$		(13b)

where the virtual work of the internal and external forces for the $i_{c}$-th cell are defined as:

		$\displaystyle\mathfrak{L}_{int}^{(i_{c})}(\boldsymbol{v},\boldsymbol{u})=\int_{\Omega^{(i_{c})}}\boldsymbol{\sigma}(\boldsymbol{u}):\boldsymbol{\varepsilon}(\boldsymbol{v})\,\mathrm{d}\Omega\;,$		(14a)
		$\displaystyle\mathfrak{L}_{ext}^{(i_{c})}(\boldsymbol{v})=\int_{\Omega^{(i_{c})}}\boldsymbol{v}\cdot\boldsymbol{b}\,\mathrm{d}\Omega+\int_{\partial\Omega^{(i_{c})}_{N}}\boldsymbol{v}\cdot\boldsymbol{t}\,\mathrm{d}\partial\Omega\;,$		(14b)

where $\boldsymbol{b}$ and $\boldsymbol{t}$ are the distributed domain force and boundary traction, respectively. The components of the strain tensor $\boldsymbol{\varepsilon}$ are:

$$\varepsilon_{ij}(\boldsymbol{u})=\frac{1}{2}\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right)\;,$$

(15)

whereas the components of the stress tensor are obtained through the constitutive relationship:

$$\sigma_{ij}=C_{ijkl}\,\varepsilon_{kl}\;,$$

(16)

being $C_{ijkl}$ the elasticity tensor that for isotropic materials reads as:

$$C_{ijkl}=\lambda\,\delta_{ij}\delta_{kl}+\mu\,(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}),$$

(17)

where $\lambda$ and $\mu$ are the Lamé parameters of the material. The linear elasticity variational statement can be reformulated through an adequate transformation directly in the parametric domain of the cell, which will be useful later in this article for the fast assembly procedure. To do so, one shall notice that

$$\frac{\partial u_{i}}{\partial x_{j}}=\frac{\partial u_{i}}{\partial\xi_{k}}\frac{\partial\xi_{k}}{\partial x_{j}}\;,\quad\mathrm{d}\Omega=j_{\Omega}\mathrm{d}\hat{\Omega}\;,\quad\mathrm{d}\partial\Omega=j_{\partial\Omega}\mathrm{d}\partial\hat{\Omega}\;,\quad\mathrm{and}\quad\varepsilon_{ij}C_{ijkl}\varepsilon_{kl}=\frac{\partial u_{i}}{\partial x_{j}}C_{ijkl}\frac{\partial u_{k}}{\partial x_{l}}\,$$

(18)

where the last equivalence comes from the symmetries of the elasticity tensor, $j_{\Omega}$ is the determinant of the Jacobian of the map in Equation (2), and $j_{\partial\Omega}$ is the curve Jacobian associated with the boundary. Substituting into Equation (14), the virtual works of the internal and external forces become:

		$\displaystyle\mathfrak{L}_{int}^{(i_{c})}(\boldsymbol{v},\boldsymbol{u})=\int_{\hat{\Omega}^{(i_{c})}}\frac{\partial u_{i}}{\partial\xi_{j}}\hat{C}_{ijkl}\frac{\partial u_{k}}{\partial\xi_{l}}\,\mathrm{d}\hat{\Omega}\;,$		(19a)
		$\displaystyle\mathfrak{L}_{ext}^{(i_{c})}(\boldsymbol{v})=\int_{\hat{\Omega}^{(i_{c})}}\boldsymbol{v}\cdot\hat{\boldsymbol{b}}\,\mathrm{d}\hat{\Omega}+\int_{\partial\hat{\Omega}^{(i_{c})}_{N}}\boldsymbol{v}\cdot\hat{\boldsymbol{t}}\,\mathrm{d}\partial\hat{\Omega}\;,$		(19b)

where the elasticity tensor and external forces have been modified as:

		$\displaystyle\hat{C}_{ijkl}=j_{\Omega}\;\frac{\partial\xi_{m}}{\partial x_{j}}C_{imkn}\frac{\partial\xi_{n}}{\partial x_{l}}\;,$		(20a)
		$\displaystyle\hat{\boldsymbol{b}}=j_{\Omega}\;\boldsymbol{b}\;,$		(20b)
		$\displaystyle\hat{\boldsymbol{t}}=j_{\partial\Omega}\;\boldsymbol{t}\;.$		(20c)

It shall be mentioned that in the previous equations, the term relative to the differential geometry of the map should be enriched with a superscript $(i_{c})$ since this is specific to the cell. However, to enhance readability and simplify the notation the superscript has not been used here.

The approximation space employed in this work is based on Lagrange polynomials defined on Gauss-Lobatto-Legendre (GLL) nodes, constructed over the untrimmed parametric domain $\hat{\Pi}^{(i_{c})}$. Accordingly, the polynomials are expressed as functions of the curvilinear coordinates $\xi_{1}$ and $\xi_{2}$. Figure 4 illustrates the univariate polynomials of degree 6 in both directions, along with the associated tensor-product GLL nodes. The corresponding bivariate polynomials are defined as:

$${B}_{ij}(\xi_{1},\xi_{2})={b}_{i}(\xi_{1})\cdot{b}_{j}(\xi_{2}),\quad\mathrm{with}\quad i,j\in\{0,\dots,p\},$$

(21)

where ${b}_{i}(\xi_{1})$ and ${b}_{j}(\xi_{2})$ denote the $i$-th and $j$-th univariate GLL polynomials in the $\xi_{1}$ and $\xi_{2}$ directions, respectively, and $p$ represents the degree of the approximation space. Formally, the approximation space for the $i_{c}$-th cell is given by

$$\mathcal{S}_{h}=\mathrm{span}\{{B}_{ij}\circ\boldsymbol{\mathcal{F}}^{-1}:i\in\{0,\dots,p\},\;j\in\{0,\dots,p\}\},$$

(22)

It is worth mentioning that in the p-FEM context, the mesh is not refined beyond the level of the individual cell, thereby making the approximation element and the geometric cell equivalent. Despite this seemingly coarse discretization of the individual cell, the high polynomial order ensures substantial approximation accuracy as it will be shown in Section 5.

Moreover, unlike other high-order bases, the GLL basis features nodes located at the domain boundaries, as illustrated in Figure 4 for $p=6$. Consequently, degrees of freedom at cell interfaces could, in principle, be strongly enforced. However, as will be detailed in Section 3, the domain decomposition strategy adopted in this work does not impose a strong coupling of all interface degrees of freedom, but still exploits the coincidence of their position.

The polynomial basis and the location of the nodes are initially defined on the untrimmed parametric domain $\hat{\Pi}^{(i_{c})}$. After the trimming operation, however, the support of each basis function is restricted to its active portion corresponding to the trimmed parametric domain of the cell, $\hat{\Omega}^{(i_{c})}$. Importantly, regardless of the shape of the trimming, each basis function remains active, even if the associated node lies outside the active domain, leading to fully dense stiffness matrices.

One might argue that, for a comparable approximation accuracy, using lower-degree polynomials with mesh refinement could provide the advantage of increased sparsity. Nevertheless, the property that the basis functions remain active irrespective of the trimming geometry is particularly advantageous for building a surrogate model of the stiffness, as discussed in Section 4, and is a primary motivation for the adoption of p-FEM elements.

Unfitted discretizations naturally introduce three main challenges: quadrature, application of Dirichlet boundary conditions, and conditioning.

Coarse DoFs are selected skeleton DoFs, or combinations thereof, used to capture the global behavior of the structure. In this work, three types of coarse DoFs are considered: cross-point values, edge averages, and edge moments (see Figure 7). This choice is inspired by [68, 69], although not identical. Cross points are nodes located at the intersection of mesh lines. The coarse DoFs associated with the $i_{x}$-th cross point are defined as

	$\displaystyle c_{i_{x}}^{u_{1}}$	$\displaystyle=u_{1}|_{\boldsymbol{x}_{i_{x}}}\,,$		(33a)
	$\displaystyle c_{i_{x}}^{u_{2}}$	$\displaystyle=u_{2}|_{\boldsymbol{x}_{i_{x}}}\,,$		(33b)

where $u_{1}$ and $u_{2}$ denote the components of the displacement vector, and $\boldsymbol{x}_{i_{x}}$ are the coordinates of the $i_{x}$-th cross point. The remaining coarse DoFs are associated to the $i_{e}$-th edge. Two average coarse DoFs are introduced, one for each displacement component:

	$\displaystyle c_{i_{e}}^{a_{1}}$	$\displaystyle=\int_{\partial\Omega_{i_{e}}}u_{1}\,\text{d}\partial\Omega\,,$		(34a)
	$\displaystyle c_{i_{e}}^{a_{2}}$	$\displaystyle=\int_{\partial\Omega_{i_{e}}}u_{2}\,\text{d}\partial\Omega\,,$		(34b)

where $\partial\Omega_{i_{e}}$ denotes the domain of the $i_{e}$-th edge. In addition, one moment-based coarse DoF is defined as

$$c_{i_{e}}^{m}=\int_{\partial\Omega_{i_{e}}}\boldsymbol{f}\cdot\boldsymbol{u}\,\text{d}\partial\Omega\,,$$

(35)

where $\boldsymbol{f}=\begin{pmatrix}x_{2}-\bar{x}_{2}\\ \bar{x}_{1}-x_{1}\end{pmatrix}$ and $(\bar{x}_{1},\bar{x}_{2})$ are the coordinates of the center of the edge. Collecting all coarse DoFs into the vector $c$ of size $n_{C}$, the following relationship holds:

$$c=Qu\;,$$

(36)

where $Q\in\mathbb{R}^{n_{C}\times n_{U}}$ is the matrix that enforces the above definitions through adequate integration and evaluation of basis functions. Figure 7 shows the different types of coarse DoFs for a single cell and for a two-by-two cells structure. As for the skeleton DoFs, the coarse DoFs vector on the $i$-th cell, denoted by $c^{(i)}$ and of size $n_{C}^{(i)}$, is obtained through an extension/restriction operator $\mathcal{R}_{C}^{(i)}\in\mathbb{R}^{n_{C}^{(i)}\times n_{C}}$,

$$c^{(i)}=\mathcal{R}_{C}^{(i)}c.$$

(37)

Similarly to $\mathcal{R}_{U}^{(i)}$, the matrix $\mathcal{R}_{C}^{(i)}$ consists of properly placed zeros and ones according to the global and local orderings of the coarse DoFs. By stacking the local coarse DoFs, the vector of redundant coarse DoFs $\tilde{c}$ of size $\tilde{n}_{C}=\sum_{i=1}^{n_{c}}n_{C}^{(i)}$ is defined as

$$\tilde{c}=\mathrm{vstack}\left(c^{(1)},\dots,c^{(n_{c})}\right)\;,$$

(38)

and the associated matrix $\tilde{\mathcal{R}}_{C}\in\mathbb{R}^{\tilde{n}_{C}\times n_{C}}$

$$\tilde{\mathcal{R}}_{C}=\mathrm{vstack}\left(\mathcal{R}_{C}^{(1)},\dots,\mathcal{R}_{C}^{(n_{c})}\right)$$

performs the operation

$$\tilde{c}=\tilde{\mathcal{R}}_{C}\,c\;.$$

Additionally, the constraint matrices $C^{(i)}\in\mathbb{R}^{n_{C}^{(i)}\times n_{U}^{(i)}}$ are introduced to extract the coarse DoFs of cell $i$ from its skeleton DoFs:

$$C^{(i)}=\mathcal{R}_{C}^{(i)}Q\mathcal{R}_{U}^{(i)\top}\;,$$

associated to the relation

$$c^{(i)}=C^{(i)}u^{(i)}\;.$$

The global version of the precedent relation is expressed as:

$$\tilde{c}=\tilde{C}\tilde{u}\;,$$

where matrix $\tilde{C}\in\mathbb{R}^{\tilde{n}_{C}\times\tilde{n}_{U}}$ is defined as

$$\tilde{C}=\mathrm{blockdiag}\left(C^{(1)},\dots,C^{(n_{c})}\right)\;.$$

Associated with the coarse DoFs $c$, the matrix $\tilde{\Psi}$ of size $\tilde{n}_{U}\times n_{C}$ is introduced, whose columns represent a coarse basis for the global correction. The matrix $\tilde{\Psi}$ is defined through

$$\begin{bmatrix}\tilde{S}&\tilde{C}^{\top}\\ \tilde{C}&0\end{bmatrix}\begin{bmatrix}\tilde{\Psi}\\ \tilde{\Lambda}\end{bmatrix}=\begin{bmatrix}0\\ \mathcal{\tilde{R}}_{C}\end{bmatrix}\;,$$

which implies that the $i$-th column of the matrix $\tilde{\Psi}$ is obtained as the solution associated with the $i$-th column of the matrix $\tilde{\mathcal{R}}_{C}$. The matrix $\tilde{\Lambda}$ contains the Lagrange multipliers associated with the coarse DoF constraints. Additionally, due to the block structures of $\tilde{S}$, $\tilde{C}$, and $\tilde{\mathcal{R}}_{C}$, the previous problem can be decomposed, and therefore parallelized, into local problems as

$$\begin{bmatrix}S^{(i)}&C^{(i)\top}\\ C^{(i)}&0\end{bmatrix}\begin{bmatrix}\Psi^{(i)}\\ \Lambda^{(i)}\end{bmatrix}=\begin{bmatrix}0\\ \mathcal{R}_{C}^{(i)}\end{bmatrix}\;,$$

(39)

where $\Lambda^{(i)}$ contains the Lagrange multipliers associated to coarse DoFs at the cell level. Furthermore, it should be noted that while $\tilde{\mathcal{R}}_{C}$ has at least two nonzero elements in every column, this is not necessarily the case for the local matrix $\mathcal{R}_{C}^{(i)}$ in cell $i$. Therefore, only a subset of the problems (39) needs to be solved for each cell. The coarse basis functions $\tilde{\Psi}$ are then assembled as

$$\tilde{\Psi}=\mathrm{vstack}\left(\Psi^{(1)},\dots,\Psi^{(n_{c})}\right)\;.$$

From a physical point of view, the $i$-th basis function (i.e., the $i$-th column of $\tilde{\Psi}$) represents the displacement field obtained by solving in each cell a local Dirichlet problem where, in addition to homogeneous Dirichlet conditions on the external boundary DoFs $u_{D}^{(i)}$, one coarse DoF is set equal to one while all the others are set to zero.

In agreement with the PCG scheme, the input to the BDDC preconditioner is the residual force vector on the skeleton. As shown in the next subsection, the first step of the preconditioner is to distribute this vector to each cell by means of the extension/restriction operator $\tilde{\mathcal{R}}_{U}$. In doing so, it is important that the sum of the extended contributions equals one and that the forces are distributed proportionally to the local stiffness associated with their DoFs. For this purpose, the diagonal matrix $\tilde{D}$ is introduced, constructed as

$$\tilde{D}=\mathrm{blockdiag}\left(D^{(1)},\dots,D^{(n_{c})}\right),$$

(40)

where, for the $i$-th cell, the $j\,j$-th entry of the diagonal matrix $D^{(i)}$, corresponding to the node $X_{l}$, is defined as

$$D^{(i)}_{[j\,j]}=\frac{K_{[j\,j]}^{(i)}}{K_{l}^{tot}},$$

(41)

where $K_{l}^{tot}$ denotes the sum of all stiffness contributions associated with node $X_{l}$ from the cells sharing that node.

The BDDC preconditioner, denoted by $S_{BDDC}^{-1}$, provides an approximation of the inverse of the global Schur complement. Its action decomposes the correction of the residual into two contributions acting on distinct $\tilde{S}$-orthogonal subspaces. The first contribution is a fine correction: independent local Neumann-type problems are solved in each subdomain. Because such local Neumann problems are in general not invertible, the correction is computed in the subspace of functions that vanish at the coarse DoFs (i.e., $\ker(C^{(i)})$), which is enforced through local saddle-point systems (see Algorithm 2). The second contribution is a coarse correction: the residual is projected onto the subspace spanned by the coarse basis functions $\tilde{\Psi}$, and a global coarse system involving $S_{C}=\tilde{\Psi}^{\top}\tilde{S}\tilde{\Psi}^{\top}$ is solved. This coarse solve balances the local fine corrections, hence the name of the method. The constraints thus serve a dual purpose: they ensure the invertibility of local problems and define a coarse basis that allows the solver to scale effectively as the problem size grows. While the exact solution is globally continuous, the approximate correction obtained from the direct sum of these two subspaces is only guaranteed to be continuous at the coarse DoFs. To restore global continuity across the entire interface skeleton, the partially discontinuous correction is projected back onto the continuous space using the restriction operator $\tilde{\mathcal{R}}_{U}$ and the diagonal scaling matrix $\tilde{D}$. This weighting is essential for ensuring the robustness and convergence of the method.

The formal expression of the preconditioner reads

$$S_{BDDC}^{-1}=\tilde{\mathcal{R}}_{U}^{\top}\tilde{D}^{\top}\tilde{\Psi}S_{C}^{-1}\tilde{\Psi}^{\top}\tilde{D}\tilde{\mathcal{R}}_{U}+\sum_{i=1}^{n_{c}}\begin{bmatrix}{D^{(i)}}{\mathcal{R}}_{U}^{(i)}\\ 0\end{bmatrix}^{\top}\begin{bmatrix}S^{(i)}&C^{(i)^{\top}}\\ C^{(i)}&0\end{bmatrix}^{-1}\begin{bmatrix}{D^{(i)}}{\mathcal{R}}_{U}^{(i)}\\ 0\end{bmatrix}\;.$$

(42)

In this expression, the first term corresponds to the coarse correction and the second term to the fine correction. The application of $S_{BDDC}^{-1}$ is carried out through a sequence of matrix-vector multiplications and local solves, as detailed in Algorithm 1. While the local problems are fully decoupled and can be solved in parallel, the coarse problem requires the solution of a global system. Since $S_{C}$ is explicitly assembled, the corresponding linear systems are solved using a direct method based on a parallel Cholesky factorization (specifically, the MUMPS solver [70], accessed through PETSc).