* [13] O. Ploder, C. Auer, C. Motz, T. Paireder, O. Lang, and M. Huemer, “SICNet—Low Complexity Sample Adaptive Neural Network-Based Self-Interference Cancellation in LTE-A/5G Mobile Transceivers,” _IEEE Open Journal of the Communications Society_ , vol. 3, pp. 958–972, 2022.
  * [14] A. Balatsoukas-Stimming, “Non-Linear Digital Self-Interference Cancellation for In-Band Full-Duplex Radios Using Neural Networks,” in _2018 IEEE 19th International Workshop on Signal Processing Advances in Wireless Communications (SPAWC)_ , 2018, pp. 1–5.
  * [15] S. Baumgartner, C. Böck, M. Huemer, and A. Balatsoukas-Stimming, “Joint Self-Interference Cancellation and Data Estimation for OFDM Based Full-Duplex Communication Systems,” in _Proceedings of the 57th Asilomar Conference on Signals, Systems, and Computers_ , 2023, accepted for publication.
  * [16] M. Mohammadkarimi, M. Mehrabi, M. Ardakani, and Y. Jing, “Deep Learning-Based Sphere Decoding,” _IEEE Transactions on Wireless Communications_ , vol. 18, no. 9, pp. 4368–4378, 2019.
  * [17] K. Pratik, B. D. Rao, and M. Welling, “RE-MIMO: Recurrent and Permutation Equivariant Neural MIMO Detection,” _IEEE Transactions on Signal Processing_ , vol. 69, pp. 459–473, 2021.
  * [18] N. Samuel, T. Diskin, and A. Wiesel, “Learning to Detect,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 10, pp. 2554–2564, 2019.
  * [19] H. He, C.-K. Wen, S. Jin, and G. Y. Li, “Model-Driven Deep Learning for MIMO Detection,” _IEEE Transactions on Signal Processing_ , vol. 68, pp. 1702–1715, 2020.
  * [20] J. Liao, J. Zhao, F. Gao, and G. Y. Li, “A Model-Driven Deep Learning Method for Massive MIMO Detection,” _IEEE Communications Letters_ , vol. 24, no. 8, pp. 1724–1728, 2020.
  * [21] M. Khani, M. Alizadeh, J. Hoydis, and P. Fleming, “Adaptive Neural Signal Detection for Massive MIMO,” _IEEE Transactions on Wireless Communications_ , vol. 19, no. 8, pp. 5635–5648, 2020.
  * [22] N. Shlezinger, N. Farsad, Y. C. Eldar, and A. J. Goldsmith, “ViterbiNet: A Deep Learning Based Viterbi Algorithm for Symbol Detection,” _IEEE Transactions on Wireless Communications_ , vol. 19, no. 5, pp. 3319–3331, 2020\.
  * [23] N. Shlezinger, R. Fu, and Y. C. Eldar, “DeepSIC: Deep Soft Interference Cancellation for Multiuser MIMO Detection,” _IEEE Transactions on Wireless Communications_ , vol. 20, no. 2, pp. 1349–1362, 2021.
  * [24] T. V. Luong, N. Shlezinger, C. Xu, T. M. Hoang, Y. C. Eldar, and L. Hanzo, “Deep Learning Based Successive Interference Cancellation for the Non-Orthogonal Downlink,” _IEEE Transactions on Vehicular Technology_ , vol. 71, no. 11, pp. 1–13, 2022.
  * [25] Y. Wei, M. Zhao, M. Hong, M. Zhao, and M. Lei, “Learned Conjugate Gradient Descent Network for Massive MIMO Detection,” in _Proceedings of the 2020 IEEE International Conference on Communications (ICC 2020)_ , 2020, pp. 1–6.
  * [26] K. Hornik, “Approximation Capabilities of Multilayer Feedforward Networks,” _Neural Networks_ , vol. 4, no. 2, pp. 251–257, 1991.
  * [27] J. R. Hershey, J. Le Roux, and F. Weninger, “Deep unfolding: Model-based inspiration of novel deep architectures,” _arXiv preprint arXiv:1409.2574_ , 2014.
  * [28] F. Pancaldi, G. M. Vitetta, R. Kalbasi, N. Al-Dhahir, M. Uysal, and H. Mheidat, “Single-Carrier Frequency Domain Equalization,” _IEEE Signal Processing Magazine_ , vol. 25, no. 5, pp. 37–56, 2008.
  * [29] M. Huemer, A. Koppler, R. Weigel, and L. Reindl, “A Review of Cyclically Extended Single Carrier Transmission with Frequency Domain Equalization for Broadband Wireless Transmission,” _European Transactions on Telecommunications_ , vol. 14, no. 4, pp. 329–341, 2003.
  * [30] M. Huemer, H. Witschnig, and J. Hausner, “Unique Word Based Phase Tracking Algorithms for SC/FDE-Systems,” in _Proceedings of GLOBECOM 03 - IEEE Global Telecommunications Conference_ , 2003, pp. 70–74.
  * [31] S. Baumgartner, O. Lang, and M. Huemer, “Neural Network Based Single-Carrier Frequency Domain Equalization,” in _Computer Aided Systems Theory - EUROCAST 2022_ , ser. Lecture Notes in Computer Science (LNCS), 2023, pp. 295–302.
  * [32] W.-J. Choi, K.-W. Cheong, and J. M. Cioffi, “Iterative Soft Interference Cancellation for Multiple Antenna Systems,” in _Proceedings of the 2000 IEEE Wireless Communications and Networking Conference. Conference Record_ , 2000, pp. 304–309.
  * [33] N. Shlezinger, R. Fu, and Y. C. Eldar, “Deep Soft Interference Cancellation for MIMO Detection,” in _Proceedings of the 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2020, pp. 8881–8885.
  * [34] H. Witschnig, T. Mayer, A. Springer, L. Maurer, M. Huemer, and R. Weigel, “The Advantages of a Known Sequence Versus Cyclic Prefix in a SC/FDE System,” in _Proceedings of the 5th International Symposium on Wireless Personal Multimedia Communications_ , vol. 3, 2002, pp. 1328–1332.
  * [35] M. Huemer, C. Hofbauer, and J. B. Huber, “Unique Word Prefix in SC/FDE and OFDM: A Comparison,” in _2010 IEEE Globecom Workshops_ , 2010, pp. 1296–1301.
  * [36] S. Reinhardt, T. Buzid, and M. Huemer, “Receiver Structures for MIMO-SC/FDE Systems,” in _Proceedings of the 2006 IEEE 63rd Vehicular Technology Conference_ , 2006, pp. 1401–1405.
  * [37] M. Huemer, L. Reindl, A. Springer, and R. Weigel, “Implementation Aspects on Single Carrier Transmission with Frequency Domain Equalization,” in _Proceedings of the 4th International OFDM Workshop_ , 1999, pp. 18.1 – 18.4.
  * [38] Y. Yang, T. Ihalainen, M. Rinne, and M. Renfors, “Frequency-Domain Equalization in Single-Carrier Transmission: Filter Bank Approach,” _EURASIP Journal on Advances in Signal Processing_ , vol. 2007, pp. 1–16, 2007.
  * [39] S. Baumgartner, G. Bognár, O. Lang, and M. Huemer, “Neural Network Approaches for Data Estimation in Unique Word OFDM Systems,” _IEEE Transactions on Vehicular Technology_ , 2023, early access.
  * [40] S. M. Kay, _Fundamentals of Statistical Signal Processing: Estimation Theory_. Prentice Hall, 1993, vol. 1.
  * [41] C. Szegedy, W. Liu, Y. Jia, P. Sermanet, S. Reed, D. Anguelov, D. Erhan, V. Vanhoucke, and A. Rabinovich, “Going Deeper with Convolutions,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2015, pp. 1–9.
  * [42] H. Xu, Y. Li, B. Tan, J. Wu, and D. Hu, “A Neural Network Empowered Belief Propagation Algorithm Optimized for Short-Cycles in Tanner Graph,” _IEEE Transactions on Machine Learning in Communications and Networking_ , vol. 1, pp. 31–42, 2023.
  * [43] R. Wiesmayr, G. Marti, C. Dick, H. Song, and C. Studer, “Bit Error and Block Error Rate Training for ML-Assisted Communication,” in _Proceedings of the 2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2023, pp. 1–5.
  * [44] M. Goutay, F. Ait Aoudia, and J. Hoydis, “Deep HyperNetwork-Based MIMO Detection,” in _Proceedings of the 2020 IEEE 21st International Workshop on Signal Processing Advances in Wireless Communications (SPAWC)_ , 2020, pp. 1–5.
  * [45] Y. A. LeCun, L. Bottou, G. B. Orr, and K.-R. Müller, _Efficient BackProp_. Springer Berlin Heidelberg, 2012, pp. 9–48.
  * [46] C. M. Bishop, _Neural Networks for Pattern Recognition_. Oxford University Press, Inc., 1995.
  * [47] S. Baumgartner, G. Bognár, O. Lang, and M. Huemer, “Neural Network Based Data Estimation for Unique Word OFDM,” in _Proceedings of the 55th Asilomar Conference on Signals, Systems, and Computers_ , Pacific Grove, USA, 2021\.
  * [48] J. Fakatselis, “Criteria for 2.4 GHz PHY Comparison of Modulation,” IEEE Document, 1997, IEEE P802.11-97/157r1.
  * [49] T. Akiba, S. Sano, T. Yanase, T. Ohta, and M. Koyama, “Optuna: A Next-Generation Hyperparameter Optimization Framework,” in _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, 2019, pp. 2623–2631.
  * [50] M. Huemer and O. Lang, “On Component-Wise Conditionally Unbiased Linear Bayesian Estimation,” in _Proceedings of the 48th Asilomar Conference on Signals, Systems and Computers_ , Nov. 2014, pp. 879–885.
  * [51] M. Huemer, C. Hofbauer, and J. Huber, “The Potential of Unique Words in OFDM,” in _Proceedings of the 15th International OFDM-Workshop_ , 2010, pp. 140–144.
  * [52] M. Huemer, C. Hofbauer, and J. B. Huber, “Non-Systematic Complex Number RS Coded OFDM by Unique Word Prefix,” _IEEE Transactions on Signal Processing_ , vol. 60, no. 1, pp. 285–299, 2012.
  * [53] O. Lang, C. Böck, M. Huemer, and C. Hofbauer, “Increasing the Bandwidth Efficiency in UW-OFDM,” in _Proceedings of the 53rd Asilomar Conference on Signals, Systems, and Computers_ , 2019, pp. 1857–1861.
  * [54] G. H. Golub and C. F. Van Loan, _Matrix Computations_. Johns Hopkins University Press, 2013.
  * [55] W. Haselmayr and A. Springer, “Extrinsic LLR Computation by the SISO LMMSE Detector: Four Different Approaches,” in _Computer Aided Systems Theory – EUROCAST 2015_ , ser. Lecture Notes in Computer Science (LNCS), 2015, pp. 529–536.


# A theoretical and numerical analysis of a Dirichlet-Neumann domain
decomposition method for diffusion problems in heterogeneous media

Alex Viguerie
Division of Mathematics
Gran Sasso Science Institute
L‘Aquila, AQ 67100, Italy
<EMAIL_ADDRESS>
&Silvia Bertoluzza
CNR Imati Enrico Magenes
Pavia, PV 27100, Italy
<EMAIL_ADDRESS>
Alessandro Veneziani
Department of Mathematics
Department of Computer Science
Emory University
Atlanta, GA 30322, USA
&Ferdinando Auricchio
Dipartimento di Ingegneria Civile e Architettura
Universita degli Studi di Pavia
Pavia, PV 27100, Italy

###### Abstract

Problems with localized nonhomogeneous material properties present well-known
challenges for numerical simulations. In particular, such problems may feature
large differences in length scales, causing difficulties with meshing and
preconditioning. These difficulties are increased if the region of localized
dynamics changes in time. Overlapping domain decomposition methods, which
split the problem at the continuous level, show promise due to their ease of
implementation and computational efficiency. Accordingly, the present work
aims to further develop the mathematical theory of such methods at both the
continuous and discrete levels. For the continuous formulation of the problem,
we provide a full convergence analysis. For the discrete problem, we show how
the described method may be interpreted as a Gauss-Seidel scheme or as a
Neumann series approximation, establishing a convergence criterion in terms of
the spectral radius of the system. We then provide a spectral scaling argument
and provide numerical evidence for its justification.

##  1 Introduction

Many problems in science and engineering feature localized dynamics in which
material properties may differ across various regions of the problem domain.
Such problems are common, for example, in additive manufacturing (AM) [11, 26,
29, 30, 33, 34, 36, 46, 48, 53, 54, 56, 60], where accurate simulation of the
problem requires adequate resolution of phenomena at both the micrometer and
millimeter scale. Problems of this type are a well-known source of numerical
difficulty, as the presence of multiple length scales, irregular geometries,
and the presence of interfaces arising from abrupt changes in material
properties often complicates the meshing and simulation process. Such problems
are compounded if the regions requiring small-scale resolution change in time
[58] .

To address these challenges, many methods have been proposed . The use of
classical methods involves remeshing at each time step, which requires
difficult-to-implement and expensive refinement-and-derefinement algorithms
[47, 13, 35, 38, 4]. Expensive meshing/remeshing can be avoided by resorting
to unfitted methods, such as GFEM/XFEM [32, 6] (and, more generally, partition
of unit methods, ([43]), Immersed Finite Element [49, 41], CutFEM ([12]), or
hierarchical methods such as the finite cell method with local enrichment [31]
or $hp$-$d$ methods [21, 45, 52, 55, 61]. Alternatively, one can resort to
methods that facilitate the remeshing procedure by allowing polygonal
elements, with possibly curved edges/faces, such as the virtual element method
[20, 5, 10, 14, 1].

Combining a domain decomposition method with a fictitious domain approach,
results in the method proposed in [59], referred to herein as the Two-level DD
method. More precisely, already at the continuous level, the domain is split
as the union of two regions, corresponding to two different materials, each
one homogeneous. The problem is then tackled by a non-overlapping Dirichlet-
Neumann DD method. Assuming that one of the two regions is small, the
corresponding problem (which we refer to as the local problem) is solved by a
standard finite element method on a fitted mesh. The problem in the larger
region (which we refer to as the global problem) is solved by resorting to a
fictitious domain approach: both coefficients and right hand side are suitably
extended to the whole domain, and, in the spirit of the Fat Boundary method
[42, 9], information on the jump of the normal flux on the interior interface
is retrieved from the solution of the local problem and injected as a data in
the global problem. The resulting algorithm exhibits many of the desirable
characteristics one may expect from a domain decomposition approach, as it is
easy to implement and precondition , and allows for the employment of quasi
uniform meshes, of possibly different size, for both local and global
problems. The regularity of both meshes and unchanging topology of the global
domain, usually corresponding to large scales, make this approach particularly
well-suited for problems in which the (small) region, where local-scale
phenomena occur, evolves in time, as in [58]. In such cases, contrary to most
of the other approaches considered in the literature for the solution of
problem of this kind, the approach of [58, 59] allows to completely avoid not
only remeshing, but also recomputing the entries of the local and global
stiffness matrices, as changes in the domain configurations can be handled by
comparatively simple mesh translations for the local domain.

Both domain decomposition (DD) and fictitious domain/immersed boundary methods
have a long history and are a well-studied topic in the literature [40, 57,
39, 50, 28, 24, 23, 37]. Indeed, the aforementioned Fat-boundary and $hp$-$d$
methods can themselves be regarded as variants within the classical DD
framework [52, 7]. Domain decomposition methods have many variants, and are
generally classified according to whether the subdomains are overlapping or
not, as well as on the way in which information is transferred between the
subdomains. For non-overlapping DD, the information exchange is carried out
through suitable boundary conditions at the interface such ad Dirichlet-
Neumann or Robin-Robin, the latter being at the basis of the optimized Schwarz
method. Information between subdomains can also be exchanged by the use of
Lagrange multiplier, as in the Mortar method, or by introducing an auxiliary
variable (which might be regarded as a control variable) on the interface [50,
15, 28, 8]. Approaches of this type are natural for problems exhibiting
particular geometric characteristics, and can be used both directly as
numerical solvers or as preconditioners [22, 50, 19, 16, 18, 17].

Though domain decomposition algorithms have indeed been used for problems of
the type studied in the current work, combining them with fictitious domain
type methods to tackle problems in which the domains are separated by
differences in physical materials requires some care in the design of the
transmission conditions, leading to the formulations expressed in the current
work. Though the application of such an approach in [58, 59] shows potential
from a numerical point of view, significant theoretical questions remain.
While the Two-level DD formulation was shown in [59] to be consistent with the
original problem formulation, the convergence behavior of the Schwarz type
iterations was not proven. At the continuous level, it is known that the
differences in the material properties may affect the convergence of a DD
method, requiring the introduction of under-relaxation or similar
regularization techniques. At the discrete level, the potentially large
difference in mesh resolutions and/or discretization methods necessary across
the different domains may also cause convergence issues. In the present work,
we begin to address the aforementioned issues, starting with diffusion
problems in heterogeneous media with piece-wise constant coefficients, for
which we provide a theoretical foundation for the practical application of
Dirichlet-Neumann non-overlapping domain decomposition, coupled with a
fictitious domain method for the global problem. We seek to clearly prove and
establish appropriate convergence behavior and conditions analytically.
Numerically, we aim to investigate and better understand the effect of
different materials and mesh resolutions on such convergence behavior.

The article is outlined as follows. We first introduce the model problem and
the corresponding Two-level DD formulation shown in [59] (Sect. 2). We then
proceed to provide a full convergence analysis (Sect. 3) . Next, we shift our
attention to the discrete problem. We will show that the discrete Two-level DD
method can be equivalently interpreted as a Gauss-Seidel method or a truncated
Neumann series, implying a spectral convergence condition (Sect. 4).
Accordingly, we then provide a scaling argument for the system eigenvalues in
terms of the problem parameters, which is validated through a series of
numerical simulations on two- and three-dimensional problems (Sect. 5).
Follow-ups of the present work are drawn in Sect. 6.

Figure 1: Domain and terminology. In the colored version, $\Gamma_{N}$ is
denoted in red and $\Gamma_{D}$ in green.

##  2 Statement of the Problem

We consider a domain given by $\Omega_{+}=\Omega_{A}\cup\Omega_{B}$ where
$\Omega_{A}\cap\Omega_{B}=\emptyset$ (Fig. 1 ). Given $f\in
L^{2}(\Omega_{+}),q\in H^{-1/2}(\Gamma_{N})$ and $T_{D}\in
H^{1/2}(\Gamma_{D})$, we solve the following variational problem on
$\Omega_{+}$.

######  Problem 1.

(Original Problem) Find $T\in H^{1}(\Omega_{+})$ such that:

$\displaystyle{\int\limits_{\Omega_{+}}{\kappa\nabla T\cdot\nabla
v}d{\Omega}}=\displaystyle{\int\limits_{\Omega_{+}}{fv}d{\Omega}}+\displaystyle{\int\limits_{\Gamma_{N}}{qv}d{\Gamma}},\quad\forall
v\in H^{1}_{\Gamma_{D}}(\Omega_{+}),\quad T(\Gamma_{D})=T_{D},$ (1)

with

$\kappa=\left\\{\begin{array}[]{lc}\kappa_{A}&\mathrm{in}\
\Omega_{A}\\\\[5.0pt] \kappa_{B}&\mathrm{in}\ \Omega_{B},\end{array}\right.$

$\kappa_{A},\kappa_{B}$ positive constants, where
$\Gamma_{D}\subseteq\partial\Omega_{+}$ (respectively
$\Gamma^{N}\subset\partial\Omega_{+}$) denote the Dirichlet (resp. Neumann)
boundary, and $H^{1}_{\Gamma_{D}}(\Omega_{+})=\\{u\in H^{1}(\Omega_{+}):\
u=0\text{ on }\Gamma_{D}\\}$. Precisely, $T\in{\cal
L}+H^{1}_{\Gamma_{D}}(\Omega_{+})$, where ${\cal L}$ is an arbitrary lifting
of $T_{D}$.

In the following we let
$\gamma_{i}=\partial\Omega_{+}\setminus\partial\Omega_{B}$ denote the portion
of the boundary of $\Omega_{B}$ interior to $\Omega_{+}$, and we assume that
$\bar{\gamma}_{i}\cap\bar{\Gamma}_{N}=\emptyset$. For
$\gamma\subseteq\partial\Omega_{B}$ we let
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left\langle\cdot,\cdot\right\rangle_{\gamma}}$
denote the duality relation between
$(H^{1}_{\Gamma_{D}}(\Omega_{B})|_{\gamma})^{\prime}$ and
$(H^{1}_{\Gamma_{D}}(\Omega_{B})|_{\gamma})$. We recall that for $T\in
H^{1}(\Omega_{B})$ with $\Delta T\in L^{2}(\Omega_{B})$ we have that $\nabla
T$ has a normal trace on $\gamma_{i}$ which belongs to
$(H^{1}_{\Gamma_{D}}(\Omega_{B})|_{\gamma_{i}})^{\prime}$, and that we have
the following “integration by part” identity (which actually is, for a generic
$T$ with no extra smoothness, the definition of the trace of the normal
derivative): for all $v\in H^{1}_{\Gamma_{D}}(\Omega_{+})$

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left\langle\nabla
T\cdot\bm{n},v\right\rangle_{\gamma_{i}}}=\int_{\Omega_{B}}\nabla T\cdot\nabla
v+\int_{\Omega_{B}}\Delta T\,v.$ (2)

Following [59], we then split the Problem (1) into two coupled problem for the
unknowns $T_{+}=T$ and $T_{-}=T|_{\Omega_{-}}$: a global problem defined on
all of $\Omega_{+}$ and a local problem defined on $\Omega_{-}=\Omega_{B}$ as
follows:

######  Problem 2 (Global/local Reformulation).

Find $T_{+}\in H^{1}(\Omega_{+})$ and $T_{-}\in H^{1}(\Omega_{-})$ such that

  * $\rightarrow$

Global

$\displaystyle\displaystyle{\int\limits_{\Omega_{+}}{\kappa_{+}\nabla
T_{+}\cdot\nabla
v}d{\Omega}}-(\kappa_{+}-\kappa_{-}){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left\langle\nabla
T_{-}\cdot\bm{n},v\right\rangle_{\gamma_{i}}}=\displaystyle{\int\limits_{\Omega_{+}}{\widetilde{f}v}d{\Omega}}+\displaystyle{\int\limits_{\Gamma_{N}}{\widetilde{q}v}d{\Gamma}},\quad\forall
v\in H^{1}_{\Gamma_{D}}(\Omega_{+}),$ (3) $\displaystyle T(\Gamma_{D})=T_{D}$
(4)

with $\kappa_{+}=\kappa_{A}$ on the whole $\Omega_{+}$, and with
$\widetilde{q}$ and $\widetilde{f}$ defined by,

$\widetilde{q}=\left\\{\begin{array}[]{lc}q&\mathrm{on}\
\Gamma_{N}\cap\partial\Omega_{A}\\\\[5.0pt]
\dfrac{\kappa_{A}}{\kappa_{B}}q&\mathrm{in}\
\Gamma_{N}\cap\partial\Omega_{B},\end{array}\right.\qquad\widetilde{f}=\left\\{\begin{array}[]{lc}f&\mathrm{in}\
\Omega_{A}\\\\[5.0pt] \dfrac{\kappa_{A}}{\kappa_{B}}f&\mathrm{in}\
\Omega_{B}.\end{array}\right.$ (5)

  * $\rightarrow$

Local

$\displaystyle\displaystyle{\int\limits_{\Omega_{-}}{\kappa_{-}\nabla
T_{-}\cdot\nabla
v}d{\Omega}}=\displaystyle{\int\limits_{\Omega_{-}}{fv}d{\Omega}}+\displaystyle{\int\limits_{\partial\Omega_{-}\cap\Gamma_{N}}{qv}d{\Gamma}},\qquad\forall
v\in H^{1}_{\gamma}(\Omega_{-}),$ (6) $\displaystyle
T_{-}(\partial\Omega_{-}\cap\Gamma_{D})=T_{D},\qquad
T_{-}(\gamma_{i})=T_{+}(\gamma_{i})$ (7)

where $\bm{n}$ denotes the unit normal pointing outwards from $\Omega_{-}$,
and $\gamma=\gamma_{i}\cup(\Gamma_{D}\cap\partial\Omega_{-})$ (for uniformity
of notation we set $\kappa_{-}\equiv\kappa_{B}$).

A more detailed explanation describing the derivation of Problem 2, and
establishing its consistency with Problem 1, may be found in [59] . The
following theorem establishes the consistency of the coupled formulation
(3)-(6) with Problem 1; it was proved in [59]. Such consistency was also shown
for non-constant coefficient and unsteady variants of the basic problem;
however for the purposes of this work we will restrict our attention to the
steady case with constant coefficients.

Remark.  We note that $\kappa_{+}$ and $\kappa_{-}$ are quantities defined
after discretization, and for this reason we elect to keep their definitions
distinct from $\kappa_{A}$ and $\kappa_{B}$. In particular, $\kappa_{+}$ is
defined on all of $\Omega_{+}$, while $\kappa_{A}$ is defined only in
$\Omega_{A}$. For this reason, $\kappa_{+}$ naturally takes the defintion of
$\kappa_{A}$ on $\Omega_{A}$, however, as its domain also includes
$\Omega_{B}$, it requires an extension over this region. The implications of
this extension are important from the numerical point of view and are explored
in section 5.5.

######  Theorem 2.1.

Let $\kappa_{A}=\kappa_{+}$, $\kappa_{B}=\kappa_{-}$, with $\kappa_{+}$ and
$\kappa_{-}$ positive constants, and let $\widetilde{q}$ and $\widetilde{f}$
be given by (5). If a function $T\in H^{1}(\Omega_{+})$, with $\Delta
T|_{\Omega_{-}}\in L^{2}(\Omega_{-})$ and $\Delta
T|_{\Omega_{+}\setminus\Omega_{-}}\in L^{2}(\Omega_{+}\setminus\Omega_{-})$,
is a weak solution of Problem 1, then, for $T_{+}=T$, $T_{-}=T|_{\Omega_{-}}$,
the couple $(T_{+},T_{-})$ solves Problem 2. Conversely, if $(T_{+},T_{-})\in
H^{1}(\Omega_{+})\times H^{1}(\Omega_{-})$ with $\Delta T_{+}|_{\Omega_{-}}\in
L^{2}(\Omega_{-})$ and $\Delta T_{+}|_{\Omega_{+}\setminus\Omega_{-}}\in
L^{2}(\Omega_{+}\setminus\Omega_{-})$, solves Problem 2, then $T_{-}=T_{+}$ in
$\Omega_{-}$, and $T_{+}$ solves Problem 1.

######  Remark 2.1.

Theorem 2.1 was proven in [59] under the stronger assumption that
$T|_{\Omega_{-}}\in H^{2}(\Omega_{-})$ and
$T|_{\Omega_{+}\setminus\Omega_{-}}\in H^{2}(\Omega_{+}\setminus\Omega_{-})$.
It is however not difficult not see that, under the assumption
$\bar{\gamma}_{i}\cap\Gamma_{N}=\emptyset$, thanks to the identity (2), such
strong assumption can be replaced by the weaker assumptions considered above.

Two-level DD Method. The split problem formulation (3)-(6) forms the backbone
of our algorithm for solving (1):

Step 0:

Solve

$\begin{array}[]{l}\kappa_{+}\displaystyle{\int\limits_{\Omega_{+}}{\nabla
T_{+}^{0}\cdot\nabla
v}d{\Omega}}=\displaystyle{\int\limits_{\Omega_{+}}{fv}d{\Omega}}+\displaystyle{\int\limits_{\Gamma_{N}}{\widetilde{q}v}d{\Gamma}},\qquad\forall
v\in H^{1}_{\Gamma_{D}}(\Omega_{+})\\\\[5.0pt]
T_{0}(\Gamma_{D})=T_{D}\end{array}$ (8)

Step ($k+{1/2}$):

Given $T_{+}^{k}\in H^{1}(\Omega_{+})$, compute $T_{-}^{k+1/2}$ by solving

$\begin{array}[]{l}\kappa_{-}\displaystyle{\int\limits_{\Omega_{-}}{\nabla
T_{-}^{k+1/2}\cdot\nabla
v}d{\Omega}}=\displaystyle{\int\limits_{\Omega_{-}}{fv}d{\Omega}}+\displaystyle{\int\limits_{\partial\Omega_{-}\cap\Gamma_{N}}{qv}d{\Gamma}},\qquad\forall
v\in{H^{1}_{\gamma}}(\Omega_{-})\\\\[5.0pt] \quad
T_{-}^{k+1/2}(\partial\Omega_{-}\cap\Gamma_{D})=T_{D},\qquad
T_{-}^{k+1/2}(\gamma_{i})=T_{+}^{k}(\gamma_{i})\end{array}$ (9)

where $\Gamma_{-\,D}=\Gamma_{D}\cap\partial\Omega_{-}$ and
$H^{1}_{\gamma}=\\{v\in H^{1}(\Omega_{-}):\ v=0\ \text{ on }\gamma\\}$.

Step ($k+1$):

given $T_{-}^{k+1/2}\in H^{1}(\Omega_{-})$, compute
$\widetilde{T}_{+}^{k+1}\in H^{1}(\Omega_{+})$ by solving

$\displaystyle\kappa_{+}\displaystyle{\int\limits_{\Omega_{+}}{\nabla\widetilde{T}_{+}^{k+1}\cdot\nabla
v}d{\Omega}}=\displaystyle{\int\limits_{\Omega_{+}}{fv}d{\Omega}}+(\kappa_{+}-\kappa_{-}){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left\langle\nabla
T_{-}^{k+1/2}\cdot\bm{n},v\right\rangle_{\gamma_{i}}}+\displaystyle{\int\limits_{\Gamma_{N}}{\widetilde{q}v}d{\Gamma}},\qquad\forall
v\in H^{1}_{\Gamma_{D}}(\Omega_{+}),$ (10)
$\displaystyle\widetilde{T}_{+}^{k+1}(\Gamma_{D})=T_{D}.$ (11)

Relaxation step:

Set

$T_{+}^{k+1}=\theta\widetilde{T}_{+}^{k+1}+\left(1-\theta\right)T_{+}^{k},\qquad
0<\theta\leq 1.$ (12)

Check convergence:

If convergence criteria are met, terminate, otherwise, repeat steps $k+1/2$
and $k+1$.

Notice that, in Step $k+1/2$, the Dirichlet boundary condition
$T_{-}^{k+1/2}(\gamma_{i})=T_{+}^{k}(\gamma_{i})$ can be enforced weakly (by
penalization or more sophisticated approaches). Here we stick to a traditional
essential treatment of the Dirichlet conditions (i.e. in the functional
space). Note that later, when we discuss the discrete problem, we will use a
penalization formulation to enforce the Dirichlet-type interface condition on
$\gamma_{i}$.

Remark.  Before continuing with the analysis, we would like to briefly discuss
how the above method differs from extant methods in the literature. The Two-
level DD method can be seen as a combination of two ingredients: a non-
overlapping domain decomposition method (similar to [18, 57, 52]) combined
with a fictitious domain method [3, 51]. For the domain decomposition, we
decompose the domain into the local problem 6, in which we consider a fitted
mesh for $\Omega_{B}$. For the global problem 3, we then consider the problem
on $\Omega_{B}$ in an unfitted manner, applying a fictitious domain
discretization. The transmission conditions considered here are Dirichlet-
Neumann (for the local and global problems, respectively); however, in
principle, we may use any other number of transmission conditions, including
Robin-Robin [40, 39, 57] or the optimization-based approach of [28, 27]. While
the related fictitious domain and non-overlapping domain decomposition methods
are well-studied in the literature, an approach combining these two paradigms,
as is done here, is, to the authors knowledge, novel, though it has some
strong similarities with the approach of [3], the main difference being that
in this last paper the transmission between local and global problem is
obtained via a distributed multiplier involving possibly cumbersome integrals
of the product of functions ”living” on two different meshes.

##  3 Convergence Theorem

In this section we establish the convergence of the iterative procedure
(8)-(12). This Theorem is one of the novel contributions of this paper.
Theorem:  The iterative procedure defined by (8)-(12) converges to the
solution $T$ of Problem 1, provided that $\theta$ is sufficiently small for
$\kappa_{+}/\kappa_{-}>2$.

We begin by analyzing the homogenous case in which
$f,\,\widetilde{f},\,q,\,\widetilde{q}$, and $T_{D}$ are uniformly zero. We
let $\bm{G}:H^{1}(\Omega_{-})\to H^{1}_{\Gamma_{D}}(\Omega_{+})$ and
$\bm{L}:H_{\Gamma_{D}}^{1}(\Omega_{+})\to H^{1}_{\Gamma_{-\,D}}(\Omega_{-})$
respectively denote the solution operators for the global problem (10) and for
the local problem (9) with homogeneous data. Given $T_{+}^{k}\,\in
H_{\Gamma_{D}}^{1}(\Omega_{+})$, $T_{+}^{k+1}$ is given by

$\displaystyle\begin{split}T_{+}^{k+1}&=\theta\bm{G}\circ\bm{L}\,[T_{+}^{k}]+(1-\theta)T_{+}^{k}.\end{split}$
(13)

We start by proving a stability bound.  We observe that
$v=T_{-}^{k+1/2}-T^{k}_{+}\in H^{1}_{\gamma}(\Omega_{-})$. We can then take
such a function as a test function in (9) and (for $f=0,q=0$) we get

$\displaystyle\begin{split}\kappa_{-}\int_{\Omega_{-}}\big{|}\nabla
T_{-}^{k+1/2}\big{|}^{2}-\kappa_{-}\int_{\Omega_{-}}\nabla
T_{-}^{k+1/2}\cdot\nabla T_{+}^{k}&=0.\end{split}$ (14)

We then immediately get the stability bound

$\displaystyle\begin{split}\kappa_{-}\big{|}T_{-}^{k+1/2}\big{|}_{H^{1}(\Omega_{-})}^{2}\leq\kappa_{-}\big{|}T_{-}^{k+1/2}\big{|}_{H^{1}(\Omega_{-})}\big{|}T_{+}\big{|}_{H^{1}(\Omega_{-})},\end{split}$
(15)

from which it follows that:

$\displaystyle\big{|}T_{-}^{k+1/2}\big{|}_{H^{1}(\Omega_{-})}$
$\displaystyle\leq\big{|}T_{+}^{k}\big{|}_{H^{1}(\Omega_{-})}.$ (16)

We may then follow a similar procedure for (10) in order to bound
$\widetilde{T}_{+}^{k+1}=\bm{G}\circ\bm{L}[T_{+}^{k}]$: using standard trace
bounds for harmonic functions we can write

$\displaystyle\begin{split}\kappa_{+}\big{|}\widetilde{T}_{+}^{k+1}\big{|}_{H^{1}(\Omega_{+})}^{2}&=\left(\kappa_{+}-\kappa_{-}\right){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left\langle\nabla
T_{-}^{k+1/2}\cdot\bm{n},\widetilde{T}_{+}^{k+1}\right\rangle_{\gamma_{i}}}\\\
&\leq\big{|}\kappa_{+}-\kappa_{-}\big{|}\big{|}T_{-}^{k+1/2}\big{|}_{H^{1}(\Omega_{-})}\big{|}\widetilde{T}_{+}^{k+1}\big{|}_{H^{1}(\Omega_{-})},\end{split}$
(17)

implying (recalling that $\kappa_{+}$ and $k_{-}$ are both positive):

$\displaystyle\begin{split}\big{|}\widetilde{T}_{+}^{k+1}\big{|}_{H^{1}(\Omega_{+})}&\leq\frac{\big{|}\kappa_{+}-\kappa_{-}\big{|}}{\kappa_{+}}\big{|}T_{-}^{k+1/2}\big{|}_{H^{1}(\Omega_{-})}\\\
&\leq\frac{\big{|}\kappa_{+}-\kappa_{-}\big{|}}{\kappa_{+}}\big{|}T_{+}^{k}\big{|}_{H^{1}(\Omega_{+})}.\end{split}$
(18)

Observe, at this point, that, if

$\frac{|\kappa_{+}-\kappa_{-}|}{\kappa_{+}}<1\Leftrightarrow(0<)\
\dfrac{\kappa_{-}}{2}<\kappa_{+}$ (19)

then the operator $\bm{G}\circ\bm{L}$ is a contraction, and the sequence
$T_{+}^{k}$ converges to a unique fixed point for all $\theta\in]0,1]$. In
particular, we can then take $\theta=1$ and the relaxation step (12) is not
necessary.

We then consider the case $\kappa_{+}\leq\kappa_{-}/2$. We let
$\widetilde{T}_{+}^{k+1}=\bm{G}\circ\bm{L}\,[T_{+}^{k}]$ and recall that
$f,\,\widetilde{f},\,q,\,\widetilde{q}$ and $T_{D}$ are all zero. We can write

$\displaystyle\begin{split}\kappa_{+}\int_{\Omega_{+}}\nabla\widetilde{T}_{+}^{k+1}\cdot\nabla
T_{+}^{k}&=(\kappa_{+}-\kappa_{-}){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left\langle\nabla
T_{-}^{k+1/2}\cdot\bm{n},T_{+}^{k}\right\rangle_{\gamma_{i}}}\\\
&=(\kappa_{+}-\kappa_{-}){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left\langle\nabla
T_{-}^{k+1/2}\cdot\bm{n},T_{-}^{k+1/2}\right\rangle_{\gamma_{i}}}\\\
&=(\kappa_{+}-\kappa_{-}){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left\langle\nabla
T_{-}^{k+1/2}\cdot\bm{n},T_{-}^{k+1/2}\right\rangle_{\partial\Omega_{-}}},\end{split}$
(20)

where we exploit the fact that $T^{k+1/2}_{-}=T^{k}_{+}$ on $\gamma_{i}$ and
the homogeneous boundary conditions on $\Gamma_{D},\,\Gamma_{N}$, where
$\partial\Omega_{-}=\gamma_{i}\cup(\partial\Omega_{-}\cap\partial\Omega_{+})$.
Integration by parts, $\Delta T_{-}^{k+1/2}=0$ in $\Omega_{-}$, and
$(\kappa_{+}-\kappa_{-})<0$ yield:

$\displaystyle\begin{split}\kappa_{+}\int_{\Omega_{+}}\nabla\widetilde{T}_{+}^{k+1}\cdot\nabla
T_{+}^{k}&=(\kappa_{+}-\kappa_{-})\left[\int_{\Omega_{-}}(\Delta
T_{-}^{k+1/2})T_{-}^{k+1/2}+\int_{\Omega_{-}}\big{|}\nabla
T_{-}^{k+1/2}\big{|}^{2}\right]\\\
&=(\kappa_{+}-\kappa_{-})|T_{-}^{k+1/2}|_{H^{1}(\Omega_{-})}^{2}\leq
0.\end{split}$ (21)

We then apply relaxation as in (12), giving:

$T_{+}^{k+1}=\theta\widetilde{T}_{+}^{k+1}+(1-\theta)T_{+}^{k}=T_{+}^{k}+\theta(\widetilde{T}_{+}^{k+1}-T_{+}^{k}).$
(22)

It then follows that:

$\displaystyle\begin{split}\kappa_{+}\big{|}T_{+}^{k+1}\big{|}_{H^{1}(\Omega_{+})}^{2}&=\kappa_{+}\bigg{[}\big{|}T_{+}^{k}\big{|}_{H^{1}(\Omega_{+})}^{2}+\theta^{2}\big{|}\widetilde{T}_{+}^{k+1}\big{|}_{H^{1}(\Omega_{+})}^{2}+\theta^{2}\big{|}T_{+}^{k}\big{|}_{H^{1}(\Omega_{+})}^{2}\\\
&\quad-2\theta^{2}(\nabla\widetilde{T}_{+}^{k+1},\,\nabla
T_{+}^{k})+2\theta(\nabla\widetilde{T}_{+}^{k+1},\,\nabla
T_{+}^{k})-2\theta\big{|}T_{+}^{k}\big{|}_{H^{1}(\Omega_{+})}^{2}\bigg{]}\\\
&=\kappa_{+}\left(1+\theta^{2}-2\theta\right)\big{|}T_{+}^{k}\big{|}_{H^{1}(\Omega_{+})}^{2}+\kappa_{+}\theta^{2}\big{|}\widetilde{T}_{+}^{k+1}\big{|}_{H^{1}(\Omega_{+})}^{2}\\\
&\quad+\kappa_{+}2(\theta-\theta^{2})(\nabla\widetilde{T}_{+}^{k+1},\nabla
T_{+}^{k})\\\
&\leq\kappa_{+}\left(1-\theta\right)^{2}\big{|}T_{+}^{k}\big{|}_{H^{1}(\Omega_{+})}^{2}+\theta^{2}\dfrac{|\kappa_{+}-\kappa_{-}|^{2}}{\kappa_{+}}\
\big{|}T_{+}^{k}\big{|}_{H^{1}(\Omega_{+})}^{2},\end{split}$ (23)

where the last line follows from (18) and (21). From (23), and the fact that
$|\kappa_{+}-\kappa_{-}|$ in this case is $\kappa_{-}-\kappa_{+}$:

$\displaystyle\big{|}T_{+}^{k+1}\big{|}_{H^{1}(\Omega_{+})}^{2}$
$\displaystyle\leq\left((1-\theta)^{2}+\theta^{2}\left(\dfrac{(\kappa_{-}-\kappa_{+})^{2}}{\kappa_{+}^{2}}\right)\right)\big{|}T_{+}^{k}\big{|}_{H^{1}(\Omega_{+})}^{2},$
(24)

so that

$\displaystyle\big{|}T_{+}^{k+1}\big{|}_{H^{1}(\Omega_{+})}^{2}$
$\displaystyle\leq\left(\left(\left(\dfrac{\kappa_{-}-\kappa_{+}}{\kappa_{+}}\right)^{2}+1\right)\theta^{2}-2\theta+1\right)\big{|}T_{+}^{k}\big{|}_{H^{1}(\Omega_{+})}^{2}$
(25)

The parabola
$\left(\left(\dfrac{\kappa_{-}-\kappa_{+}}{\kappa_{+}}\right)^{2}+1\right)\theta^{2}-2\theta+1$
has value 1 and slope negative for $\theta=0$ (the slope being -2), minimum in
$\theta_{opt}=\left(\left(\dfrac{\kappa_{-}-\kappa_{+}}{\kappa_{+}}\right)^{2}+1\right)^{-1}$
with value $1-\theta_{opt}<1$. So, there exists an interval $(0,\bar{\theta}$)
where it takes values $<1$, proving that the map is a contraction also in this
case.

We now consider the non-homogeneous case, with $f,\,\widetilde{f}\in
L^{2}(\Omega_{+})$ and $q,\,\widetilde{q}\in H^{-1/2}(\Gamma_{N})$ and
$T_{D}\in H^{1/2}(\Gamma_{D})$. Let $T_{+}^{*},\,\widehat{T}_{+}^{*}\in
H^{1}(\Omega_{+})$ and $T_{-}^{*}\in H^{1}(\Omega_{-})$ be auxiliary functions
defined as the solutions to the following problems, the first one being a
standard elliptic boundary value problem with given data, the second and third
being also standard elliptic boundary value problems with data depending on
the solution of the first and second problem, respectively:

$\displaystyle\begin{split}\int_{\Omega_{+}}\kappa_{+}\nabla
T_{+}^{*}\cdot\nabla
v&=\int_{\Omega_{+}}\widetilde{f}v+\int_{\Gamma_{N}}\widetilde{q}v,\quad\forall
v\in H_{\Gamma_{D}}^{1}(\Omega_{+}),\qquad
T_{+}^{*}(\Gamma_{D})=T_{D};\end{split}$ (26)

$\displaystyle\begin{split}&\int_{\Omega_{-}}\kappa_{-}\nabla
T_{-}^{*}\cdot\nabla
v=\int_{\Omega_{-}}fv+\int_{\Gamma_{N}\cap\partial\Omega_{-}}qv,\qquad\forall
v\in H_{\gamma}^{1}(\Omega_{-})\\\
&T_{-}^{*}(\Gamma_{D}\cap\partial\Omega_{-})=T_{D},\qquad
T_{-}^{*}(\gamma_{i})=T_{+}^{*})(\gamma_{i});\end{split}$ (27)

$\displaystyle\begin{split}\int_{\Omega_{+}}\kappa_{+}\nabla\widehat{T}_{+}^{*}\cdot\nabla
v&=(\kappa_{+}-\kappa_{-}){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left\langle\nabla
T_{-}^{*}\cdot\bm{n},v\right\rangle_{\gamma_{i}}},\qquad\forall v\in
H_{\Gamma_{D}}^{1}(\Omega_{+}),\qquad\widehat{T}_{+}^{*}(\Gamma_{D})=0.\end{split}$
(28)

To study the convergence of the sequence $T_{+}^{k}$ of the iteration
(8)–(12), we investigate the existence of a limit to the sequence

$w_{+}^{k}=T_{+}^{k}-T_{+}^{*}.$

To this aim we show that the operator mapping $w_{+}^{k}$ to $w_{+}^{k+1}$ is
a contraction. We start by rewriting (9) and (10) in terms of $w_{+}^{k}$ and
of the new local unknown

$w_{-}^{k+1/2}=T_{-}^{k+1/2}-T_{-}^{*}.$

This gives us the following equation for $w_{-}^{k+1/2}$:

$\begin{array}[]{l}\kappa_{-}\displaystyle{\int\limits_{\Omega_{-}}{\nabla(w_{-}^{k+1/2}+T_{-}^{*})\cdot\nabla
v}d{\Omega}}=\displaystyle{\int\limits_{\Omega_{-}}{fv}d{\Omega}}+\displaystyle{\int\limits_{\partial\Omega_{-}\cap\Gamma_{N}}{qv}d{\Gamma}},\qquad\forall
v\in H^{1}_{\Gamma_{D}\cup\gamma_{i}}(\Omega_{-})\\\\[5.0pt]
\quad(w_{-}^{k+1/2}+T_{-}^{*})(\partial\Omega_{-}\cap\Gamma_{D})=T_{D},\quad(w_{-}^{k+1/2}+T_{-}^{*})(\gamma_{i})=(w_{+}^{k}+T_{+}^{*})(\gamma_{i}).\end{array}$
(29)

Substituting (LABEL:TMinusStar) into (29), it follows promptly that:

$\displaystyle w_{-}^{k+1/2}$ $\displaystyle=\bm{L}w_{+}^{k},$ (30)

where $\bm{L}$ is the homogeneous local-problem solution operator studied
above.

Following a similar procedure for (10), we see that, setting
$\widetilde{w}_{+}^{k+1}=\widetilde{T}_{+}^{k+1}-T_{+}^{*}$, we can write, for
all $v$ in $H^{1}_{\Gamma_{D}}(\Omega_{+})$,

$\displaystyle\begin{split}\int_{\Omega_{+}}\kappa_{+}\nabla\left(\widetilde{w}_{+}^{k+1}+T_{+}^{*}\right)\cdot\nabla
v&=\int_{\Omega_{+}}\widetilde{f}v+\int_{\Gamma_{N}}\widetilde{q}v+(\kappa_{+}-\kappa_{-}){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left\langle\nabla
T_{-}^{k+1/2}\cdot\bm{n},v\right\rangle_{\gamma_{i}}},\qquad\forall v\in
H^{1}_{\Gamma_{D}}(\Omega_{+})\\\
&\widetilde{w}_{+}^{k+1}(\Gamma_{D})+T_{+}^{*}(\Gamma_{D})=T_{D}\end{split}$
(31)

whence, using (26) and (28)

$\displaystyle\begin{split}\int_{\Omega_{+}}\kappa_{+}\nabla\widetilde{w}_{+}^{k+1}\cdot\nabla
v&=(\kappa_{+}-\kappa_{-}){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left\langle\nabla\left(w_{-}^{k+1/2}+T_{-}^{*}\right)\cdot\bm{n},v\right\rangle_{\gamma_{i}}}\\\
&=(\kappa_{+}-\kappa_{-}){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left\langle\nabla\left(w_{-}^{k+1/2}\right)\cdot\bm{n},v\right\rangle_{\gamma_{i}}}+\int_{\Omega_{+}}\kappa_{+}\nabla\widehat{T}_{+}^{*}\cdot\nabla
v.\end{split}$ (32)

Then, $\widetilde{w}_{+}^{k+1}-\widehat{T}_{+}^{*}=\bm{G}w_{-}^{k+1/2}$, and,
therefore,

$\widetilde{w}_{+}^{k+1}=\bm{G}[w_{-}^{k+1/2}]+\widehat{T}_{+}^{*}=\bm{G}\circ\bm{L}[w_{+}^{k}]+\widehat{T}_{+}^{*},$
(33)

Applying relaxation, one obtains:

$\displaystyle\begin{split}w_{+}^{k+1}&=\theta\widetilde{w}_{+}^{k+1}+(1-\theta)w_{+}^{k}\\\
&=\theta(\bm{G}\circ\bm{L}[w_{+}^{k}]+\widehat{T}_{+}^{*})+(1-\theta)w_{+}^{k}\\\
&=\left(\theta\bm{G}\circ\bm{L}[w_{+}^{k}]+(1-\theta)w_{+}^{k}\right)+\theta\widehat{T}_{+}^{*}.\end{split}$
(34)

From (34), the map $w_{+}^{k+1}={\cal M}[w_{+}^{k+1}]$ is contractive.

It follows that $T_{+}^{k+1}=T_{+}^{*}+w_{+}^{k+1}$, obtained directly by
solving the nonhomogeneous problem, converges to the solution of the coupled
problem (3)-(6).

######  Remark 3.1.

Observe that, in order for the iterative procedure to converge, we did not
require that $q$ and $\widetilde{q}$ (resp. $f$ and $\widetilde{f}$) satisfy
(5). Of course, such relations are needed if we want the solution of the
coupled problem to coincide with the solution of the original problem, as
stated in Theorem 2.1.  However, it could be interesting to exploit the
freedom in the choice of the extension of $f$ in the design of the method. In
particular, in the spirit of [2], one could look for an extension
$\widetilde{f}$ such that the jump along $\gamma_{i}$ of $\nabla
T^{+}\cdot\bm{n}$ vanishes, thus allowing for optimal convergence rates.

##  4 Discrete Problem

Having formally established the convergence of the continuous problem in the
preceding section, in the present section we consider the discrete version of
Problem (1) and of its corresponding Two-level DD formulation (8)-(10). We
first introduce a discrete monolithic problem. Let $\mathcal{T}_{+}$ and
$\mathcal{T}_{-}$ denote appropriate discretizations of $\Omega_{+}$ and
$\Omega_{-}$, with $\mathcal{E}_{i}$ denoting the discretization of
$\gamma_{i}$ induced by $\mathcal{T}_{-}$, and let $X_{\mathcal{T}_{+}}\subset
H^{1}(\Omega_{+})$, $X_{\mathcal{T}_{-}}\subset H^{1}(\Omega_{-})$ denote the
corresponding finite dimensional approximation spaces. Using a penalization
method with parameter $\alpha>0$ for the local problem to impose the
Dirichlet-type coupling condition on $\gamma_{i}$, we consider the following
discrete equations: find $(T_{+}^{h},T_{-}^{h})$ in $X_{\mathcal{T}_{+}}\times
X_{\mathcal{T}_{-}}$ such that for all $v_{h}$ in $X_{\mathcal{T}_{+}}\cap
H^{1}_{\Gamma_{D}}(\Omega_{+})$ and $w_{h}$ in $X_{\mathcal{T}_{-}}\cap
H^{1}_{\Gamma_{-\,D}}(\Omega_{-})$:

$\displaystyle\begin{split}\int_{\Omega_{+}}\kappa_{+}\nabla
T_{+}^{h}\cdot\nabla
v_{h}-{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum_{e\in\mathcal{E}_{i}}\int_{e}\left(\kappa_{+}-\kappa_{-}\right)\left(\nabla
T_{-}^{h}\cdot\bm{n}\right)v_{h}}&=\int_{\Omega_{+}}\widetilde{f}\,v_{h}+\int_{\Gamma_{N}\cap\partial\Omega_{-}}\widetilde{q}\,v_{h}\\\
\int_{\Omega_{-}}\kappa_{-}\nabla T_{-}^{h}\cdot\nabla
w_{h}+\alpha\int_{\gamma_{i}}T_{-}^{h}w_{h}-\alpha\int_{\gamma_{i}}T_{+}^{h}w_{h}&=\int_{\Omega_{-}}f\,w_{h}+\int_{\Gamma_{N}\cap\partial\Omega_{-}}q\,w_{h}\\\
T_{+}^{h}&=T_{0}\text{ on }\Gamma_{D}.\end{split}$ (35)

respectively, with corresponding discrete function spaces
$X_{\mathcal{T}_{+}}$ and $X_{\mathcal{T}_{-}}$. We denote discrete functions
with the subscript $h$.

We define the following matrices resulting from the discretization of the
bilinear forms in (35):

$\displaystyle K_{+}$ $\displaystyle\ \text{for}\
\int_{\Omega_{+}}\kappa_{+}\nabla T_{+}^{h}\cdot\nabla v_{h}$ (36)
$\displaystyle S$ $\displaystyle\ \text{for}\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum_{e\in\mathcal{E}_{i}}\int_{e}\left(\kappa_{+}-\kappa_{-}\right)\left(\nabla
T_{-}^{h}\cdot\bm{n}\right)v_{h}}$ (37) $\displaystyle D$ $\displaystyle\
\text{for}\ -\alpha\int_{\gamma_{i}}T_{+}^{h}w_{h}$ (38) $\displaystyle K_{-}$
$\displaystyle\ \text{for}\ \int_{\Omega_{-}}\kappa_{-}\nabla
T_{-}^{h}\cdot\nabla w_{h}+\alpha\int_{\gamma_{i}}T_{-}^{h}w_{h},$ (39)

using which the algebraic form of (35) reads

$\displaystyle\begin{bmatrix}K_{+}&S\\\
D&K_{-}\end{bmatrix}\begin{bmatrix}\bm{T}_{+}^{h}\\\
\bm{T}_{-}^{h}\end{bmatrix}$ $\displaystyle=\begin{bmatrix}\bm{f}_{+}\\\
\bm{f}_{-}\end{bmatrix}.$ (40)

We now offer two useful and equivalent algebraic interpretations of Algorithm
in (8)-(12).

###  4.1 Gauss-Seidel-type interpretation

We note that (35) can be written equivalently as:

$\displaystyle\begin{bmatrix}K_{+}&S\\\
0&K_{-}\end{bmatrix}\begin{bmatrix}\bm{T}_{+}^{h}\\\
\bm{T}_{-}^{h}\end{bmatrix}$ $\displaystyle=\begin{bmatrix}\bm{f}_{+}\\\
\bm{f}_{-}\end{bmatrix}+\begin{bmatrix}0&0\\\
-D&0\end{bmatrix}\begin{bmatrix}\bm{T}_{+}^{h}\\\
\bm{T}_{-}^{h}\end{bmatrix}.$ (41)

The above splitting (41) can then be employed to solve (40) iteratively as an
inverted block Gauss-Seidel-type method111The reason why we call it “inverted”
is because it works with an upper triangular matrix as opposed to a lower
triangular one, as in the standard GS.\- given
$[\bm{T}_{+}^{h,0},\,\bm{T}_{-}^{h,0}]^{T}$, one solves until convergence:

$\displaystyle\begin{bmatrix}K_{+}&S\\\
0&K_{-}\end{bmatrix}\begin{bmatrix}\bm{T}_{+}^{h,k+1}\\\
\bm{T}_{-}^{h,k+1}\end{bmatrix}$ $\displaystyle=\begin{bmatrix}\bm{f}_{+}\\\
\bm{f}_{-}\end{bmatrix}+\begin{bmatrix}0&0\\\
-D&0\end{bmatrix}\begin{bmatrix}\bm{T}_{+}^{h,k}\\\
\bm{T}_{-}^{h,k}\end{bmatrix}.$ (42)

One may modify (42) to incorporate an under-relaxation parameter $\theta$,

$\displaystyle\begin{bmatrix}\frac{1}{\theta}K_{+}&S\\\
0&K_{-}\end{bmatrix}\begin{bmatrix}\bm{T}_{+}^{h,k+1}\\\
\bm{T}_{-}^{h,k+1}\end{bmatrix}$ $\displaystyle=\begin{bmatrix}\bm{f}_{+}\\\
\bm{f}_{-}\end{bmatrix}+\begin{bmatrix}\frac{1-\theta}{\theta}K_{+}&0\\\
-D&0\end{bmatrix}\begin{bmatrix}\bm{T}_{+}^{h,k}\\\
\bm{T}_{-}^{h,k}\end{bmatrix}.$ (43)

The scheme (42) leads to solving the sequence of problems:

$\displaystyle\begin{split}K_{-}\bm{T}_{-}^{h,k+1}&=\bm{f}_{-}-D\bm{T}_{+}^{h,k}\,,\end{split}$
(44)
$\displaystyle\begin{split}K_{+}\bm{T}_{+}^{h,k+1}&=\bm{f}_{+}-S\bm{T}_{-}^{h,k+1}\,,\end{split}$
(45)

where (44) and (45) are the algebraic versions of step $k+1/2$ and $k+1$,
respectively, in (8)-(10). Incorporating the under-relaxation step (12) is
equivalent to solving the modified system (43):

$\displaystyle\begin{split}K_{-}\bm{T}_{-}^{h,k+1}&=\bm{f}_{-}-D\bm{T}_{+}^{h,k}\,,\end{split}$
(46)
$\displaystyle\begin{split}\frac{1}{\theta}K_{+}\bm{T}_{+}^{h,k+1}&=\bm{f}_{+}+\frac{1-\theta}{\theta}K_{+}\bm{T}_{+}^{h,k}-S\bm{T}_{-}^{h,k+1}\,.\end{split}$
(47)

It is well-known that the block Gauss-Seidel iteration (43) converges to the
solution of (40) provided that the spectral radius of the system:

$\displaystyle\begin{bmatrix}\frac{1}{\theta}K_{+}&S\\\
0&K_{-}\end{bmatrix}^{-1}\begin{bmatrix}\frac{1-\theta}{\theta}K_{+}&0\\\
-D&0\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\left(1-\theta\right)I_{+}+\theta
K_{+}^{-1}SK_{-}^{-1}D&0\\\ -K_{-}^{-1}D&0\end{bmatrix}$ (48)

is less than one (see e.g. [25]), where $I_{+}$ denotes the identity matrix of
appropriate size for the function space $X_{\mathcal{T}_{+}}$ ($I_{-}$ will be
used similarly). As (48) features a zero block on the diagonal, the
convergence criterion reduces to the condition:

$\displaystyle\rho\left(\left(1-\theta\right)I_{+}+\theta
K_{+}^{-1}SK_{-}^{-1}D\right)<1.$ (49)

This is the algebraic counterpart of the condition (19); both rely on the
relative values of $\kappa_{-}$ and $\kappa_{+}$ (as we will confirm with the
numerical experiments).

###  4.2 Neumann series interpretation ($\theta=1$)

Following classical arguments related to the Schur complement [25], we may
factor (40) in the following block-UL system:

$\displaystyle\begin{bmatrix}K_{+}-SK_{-}^{-1}D&S\\\
0&K_{-}\end{bmatrix}\begin{bmatrix}I_{+}&0\\\
K_{-}^{-1}D&I_{-}\end{bmatrix}\begin{bmatrix}\bm{T}_{+}^{h}\\\
\bm{T}_{-}^{h}\end{bmatrix}$ $\displaystyle=\begin{bmatrix}\bm{f}_{+}\\\
\bm{f}_{-}\end{bmatrix}.$ (50)

Solving the system yields:

$\displaystyle\begin{split}\bm{T}_{+}^{h}&=\left[K_{+}-SK_{-}^{-1}D\right]^{-1}\left(\bm{f}_{+}-SK_{-}^{-1}\bm{f}_{-}\right)\\\
\bm{T}_{-}^{h}&=K_{-}^{-1}\bm{f}_{-}-K_{-}^{-1}D\bm{T}_{+}^{h}\\\
&=K_{-}^{-1}\bm{f}_{-}-K_{-}^{-1}D\left[K_{+}-SK_{-}^{-1}D\right]^{-1}\left(\bm{f}_{+}-SK_{-}^{-1}\bm{f}_{-}\right).\end{split}$
(51)

Similarly, for an iteration $k+1$ of (44) -(45) (assuming $\theta=1$):

$\displaystyle\begin{split}\bm{T}_{+}^{h,k+1}&=K_{+}^{-1}\bm{f}_{+}-K_{+}^{-1}S\bm{T}_{-}^{h,k+1}\\\
&=K_{+}^{-1}\bm{f}_{+}-K_{+}^{-1}S\left[K_{-}^{-1}\bm{f}_{-}-K_{-}^{-1}D\bm{T}_{+}^{h,k}\right]\\\
&=K_{+}^{-1}\left(\bm{f}_{+}-SK_{-}^{-1}\bm{f}_{-}\right)+K_{+}^{-1}SK_{-}^{-1}D\bm{T}_{+}^{h,k}\\\
\bm{T}_{-}^{h,k+1}&=K_{-}^{-1}\bm{f}_{-}-K_{-}^{-1}D\bm{T}_{+}^{h,k}.\end{split}$
(52)

######  Proposition 1.

Given an initial guess $\bm{T}_{+}^{h,0}$, at an iteration $k$ of (44)-(45) we
may write $\bm{T}_{+}^{h,k}$ as:

$\bm{T}_{+}^{h,k}=\left[\sum_{j=0}^{k-1}\left(K_{+}^{-1}SK_{-}^{-1}D\right)^{j}\right]K_{+}^{-1}\left(\bm{f}_{+}-SK_{-}^{-1}\bm{f}_{-}\right)+\left(K_{+}^{-1}SK_{-}^{-1}D\right)^{k}\bm{T}_{+}^{h,0}.$

Further, if (49) holds, then the iteration given by (44)-(45) converges to the
solution of (40).

Proof.  Proceed by induction. The case $k=1$ is trivially true. Assume that
the hypothesis holds for $k$. At $k+1$:

$\displaystyle\begin{split}\bm{T}_{+}^{h,k+1}&=K_{+}^{-1}\left(\bm{f}_{+}-SK_{-}^{-1}\bm{f}_{-}\right)+K_{+}^{-1}SK_{-}^{-1}D\bm{T}_{+}^{h,k}\\\
&=K_{+}^{-1}\left(\bm{f}_{+}-SK_{-}^{-1}\bm{f}_{-}\right)+\\\ &\quad
K_{+}^{-1}SK_{-}^{-1}D\left[\left[\sum_{j=0}^{k-1}\left(K_{+}^{-1}SK_{-}^{-1}D\right)^{j}\right]K_{+}^{-1}\left(\bm{f}_{+}-SK_{-}^{-1}\bm{f}_{-}\right)+\left(K_{+}^{-1}SK_{-}^{-1}D\right)^{k}\bm{T}_{+}^{h,0}\right]\\\
&=K_{+}^{-1}\left(\bm{f}_{+}-SK_{-}^{-1}\bm{f}_{-}\right)+\\\
&\quad\left[\sum_{j=1}^{k}\left(K_{+}^{-1}SK_{-}^{-1}D\right)^{j}\right]K_{+}^{-1}\left(\bm{f}_{+}-SK_{-}^{-1}\bm{f}_{-}\right)+\left(K_{+}^{-1}SK_{-}^{-1}D\right)^{k+1}\bm{T}_{+}^{h,0}\\\
&=\left[\sum_{j=0}^{k}\left(K_{+}^{-1}SK_{-}^{-1}D\right)^{j}\right]K_{+}^{-1}\left(\bm{f}_{+}-SK_{-}^{-1}\bm{f}_{-}\right)+\left(K_{+}^{-1}SK_{-}^{-1}D\right)^{k+1}\bm{T}_{+}^{h,0},\end{split}$
(53)

completing the first part of the proof.

It remains to show that if (49) for $\theta=1$ holds, the solutions (51) and
(52) are equivalent as $k\to\infty$. Note that for $\theta=1$ (49) implies the
convergence of the Neumann series

$\displaystyle\begin{split}\left(K_{+}-SK_{-}^{-1}D\right)^{-1}&=\left[K_{+}\left(I_{+}-K_{+}^{-1}SK_{-}^{-1}D\right)\right]^{-1}\\\
&=\left(I_{+}-K_{+}^{-1}SK_{-}^{-1}D\right)^{-1}K_{+}^{-1}\\\
&=\left[\sum_{j=0}^{\infty}\left(K_{+}^{-1}SK_{-}^{-1}D\right)^{j}\right]K_{+}^{-1},\end{split}$
(54)

hence from (51):

$\displaystyle\bm{T}_{+}^{h}$
$\displaystyle=\left[\sum_{j=0}^{\infty}\left(K_{+}^{-1}SK_{-}^{-1}D\right)^{j}\right]K_{+}^{-1}\left(\bm{f}_{+}-SK_{-}^{-1}\bm{f}_{-}\right).$
(55)

Condition (49) for $\theta=1$ implies additionally that:

$\displaystyle\lim_{k\to\infty}\left(K_{+}^{-1}SK_{-}^{-1}D\right)^{k}$
$\displaystyle=0,$ (56)

with $0$ denoting the zero matrix of appropriate dimension. Taken together,
(53), (55), and (56) imply the result.

Proposition 1 provides powerful intuition regarding the performance of the
discrete Algorithm (8)-(12), as it can be seen as a finite approximation of
the Neumann series expansion. We therefore expect that the convergence rate
will depend on $\rho\left(K_{+}^{-1}SK_{-}^{-1}D\right)$.

###  4.3 Eigenvalue scaling

Both the Gauss-Seidel (for $\theta=1$) and Neumann series interpretations of
Algorithm in (8)-(12) require the same condition for convergence:
$\rho\left(K_{+}^{-1}SK_{-}^{-1}D\right)<1$. Although a fully rigorous
spectral analysis of this system is a worthy subject of future work, we will
provide a more heuristic approach here which we will validate with numerical
tests.

Looking at the operators (36)-(39), one may argue that $\rho$ scales in the
following way:

$\displaystyle\begin{split}\rho\left(K_{+}^{-1}SK_{-}^{-1}D\right)&\sim\frac{1}{\kappa_{+}}\left|-\kappa_{+}+\kappa_{-}\right|\frac{1}{\kappa_{-}+\alpha}\alpha\\\
&=\frac{|\kappa_{-}-\kappa_{+}|}{\kappa_{+}}\frac{\alpha}{\kappa_{-}+\alpha}\\\
&=\frac{\left|\kappa_{-}-\kappa_{+}\right|\alpha}{\kappa_{+}\left(\frac{\kappa_{-}}{\alpha}+1\right)\alpha}\\\
&=\frac{|\kappa_{-}-\kappa_{+}|}{\kappa_{+}\left(\frac{\kappa_{-}}{\alpha}+1\right)}.\end{split}$
(57)

We recall that $\alpha$ is a penalization parameter and is large in general,
implying that $\kappa_{-}/\alpha\approx 0$, and therefore that:

$\displaystyle\begin{split}\rho\left(K_{+}^{-1}SK_{-}^{-1}D\right)&\sim\frac{|\kappa_{-}-\kappa_{+}|}{\kappa_{+}}=\left|\frac{\kappa_{-}}{\kappa_{+}}-1\right|.\end{split}$
(58)

This analysis is consistent with the results of the unrelaxed scheme in Sect.
3, and specifically with (19), for $\bm{G}\circ\bm{L}$ to be a contraction.

We therefore postulate the following scaling behavior:

$\displaystyle\rho\left(\frac{\kappa_{-}}{\kappa_{+}}\right)$
$\displaystyle\sim
C_{h_{+},h_{-},\gamma_{i},\Omega_{+},\Omega_{-}}\left|\frac{\kappa_{-}}{\kappa_{+}}-1\right|,$
(59)

where $C_{h_{+},h_{-},\gamma,\Omega_{+},\Omega_{-}}$ is a constant depending
on the geometric and mesh parameters and will be denoted by $\widetilde{C}$
hereafter.

##  5 Numerical experiments

In this section, we will perform a series of numerical experiments to confirm
the theory discussed in previous sections. We wish to confirm/examine the
following:

  1. 1.

That the convergence of the Two-level DD method is dependent on
$\rho(K_{+}^{-1}SK_{-}^{-1}D)$ as postulated;

  2. 2.

That the scaling behavior (59) holds;

  3. 3.

How $\rho(K_{+}^{-1}SK_{-}^{-1}D)$ is affected by changes in global mesh level
$h_{+}$, local mesh level $h_{-}$, spatial dimension, and polynomial degree.

To answer these questions, we will run series of 2D and 3D tests, organizing
the analysis of the results analysis into two distinct portions:

  1. 1.

Convergence analysis. This will focus directly on the relationship between
$\rho(K_{+}^{-1}SK_{-}^{-1}D)$ and convergence of the Two-level DD method. For
a given case, we will use the collected data examine the impact of
$\rho(K_{+}^{-1}SK_{-}^{-1}D)$ on convergence. We then verify the point at
which the Two-level DD method no longer converges, where
$\rho(K_{+}^{-1}SK_{-}^{-1}D)>1$.

  2. 2.

Spectral growth analysis. After confirming the importance of
$\rho(K_{+}^{-1}SK_{-}^{-1}D)$ with the convergence analysis, we will examine
the full range of results in 2D and 3D in order to understand how differences
in global and local mesh resolution, degree of polynomial approximation, and
spatial dimension affect the value of $\rho(K_{+}^{-1}SK_{-}^{-1}D)$. We hope
to observe that $\rho(K_{+}^{-1}SK_{-}^{-1}D)\sim\mathcal{O}\left(\log
h_{+}/h_{-}\right)$, for all polynomial degrees in both 2D and 3D, which will
indicate that refinement of the local mesh does not lead to a deterioration of
convergence behavior.

The underlying data and problem setup is identical for both analyses. For this
reason, we outline the setup and experimental procedure below for each case
before proceeding to the discussion of results.

###  5.1 Problem setup

Figure 2: Basic test setup for the problem solved throughout this work. This
is the two-dimensional version; the three-dimensional problem features the
same general setup with a cubic domain rather than a square domain.

In order to control the experimental setting as much as possible, we define a
standardized problem setup. Referring to notation shown in Fig 2, we solve the
2D problem in a square with $L=1/40$ and $H=1/40$. We define $H_{-}=1/160$ and
set $T_{0}=293.15$. The heat flux profile is assigned along the top of the
domain and defined by:

$\displaystyle q_{2D}(x)$
$\displaystyle=.4\text{e}5\exp\left(\frac{-\left(L/2-x\right)^{4}}{1\text{e-}12}\right).$
(60)

We fix the global mesh level $h_{+}=1/160$ and $\kappa_{+}=1.0$.

The 3D is analogous, considering instead a cube of dimension $1/40\times
1/40\times 1/40$, and modifying the heat flux profile accordingly to:

$\displaystyle q_{3D}(x,y)$
$\displaystyle=.4\text{e}5\exp\left(\frac{-\left(\left(L/2-x\right)^{4}+\left(L/2-y\right)^{4}\right)}{1\text{e-}12}\right).$
(61)

###  5.2 Experimental procedure

Using the standardized problem setup defined previously, each experimental
case is identified by the following characteristics:

  * •

Spatial dimension: $d$ = 2, 3;

  * •

Degree of polynomial approximation $\mathcal{P}^{m}$: $m$= 1, 2;

  * •

Local mesh resolution $h_{-}$: $1/320$, $1/640$, $1/1280$, $1/2560$. As global
mesh resolution $h_{+}$ is fixed, these are identified by the ratio
$h_{+}/h_{-}=$ $2$, $4$, $8$, $16$ respectively.

For brevity and clarity, we will refer to each case with this terminology:
Case $d$=2, $m$=1, $h_{+}/h_{-}$=4 is understood as the two-dimensional case
with linear polynomial approximation and $h_{-}$=$1/640$.

Each case then consists of 8 simulations, with $\kappa_{-}=.5^{l}$ for
$l=1,\,2,\,...,\,8$.222Due to memory constraints, the case $d$=3, $m$=2,
$h_{+}/h_{-}$=16 was not considered. We then compute linear and quadratic
least-squares polynomial fittings between $\kappa_{-}/\kappa_{+}$ and
$\rho(K_{+}^{-1}SK_{-}^{-1}D)$:

$\displaystyle\begin{split}\hat{\rho}_{L}\left(\frac{\kappa_{-}}{\kappa_{+}}\right)&\approx\left|a_{1}\left(\frac{\kappa_{-}}{\kappa_{+}}\right)+a_{0},\right|\\\
\hat{\rho}_{Q}\left(\frac{\kappa_{-}}{\kappa_{+}}\right)&\approx\left|b_{2}\left(\frac{\kappa_{-}}{\kappa_{+}}\right)^{2}+b_{1}\left(\frac{\kappa_{-}}{\kappa_{+}}\right)+b_{0}\right|.\end{split}$
(62)

If (59) holds, we expect that $|b_{2}|\approx 0$ and $a_{1}\approx-a_{0}$,
yielding:

$\displaystyle\hat{\rho}\left(\frac{\kappa_{-}}{\kappa_{+}}\right)$
$\displaystyle\approx\left|-a_{1}\left(\frac{\kappa_{-}}{\kappa_{+}}-1\right)\right|,$
(63)

implying that:

$\displaystyle\widetilde{C}$ $\displaystyle\approx-a_{1},$ (64)

giving us an estimate of $\widetilde{C}$ for each case.

To summarize clearly, we will

  * •

For each spatial dimension $d$, polynomial degree $m$, and local mesh
resolution $h_{+}/h_{-}$, simulate the model problem with
$\kappa_{-}=\frac{1}{2}^{k}$ for $k$=$1,\,2,\,...,\,8$. ;

  * •

The results of these simulations will be used to compute (62) and ideally
estimate $\widetilde{C}$.

We start considering the unrelaxed case ($\theta=1$), for which we know from
(19) and (58) that for $\kappa_{-}/\kappa_{+}<2$ convergence of the continuous
scheme (8)-(12) is guaranteed. Successively, we consider the impact of the
relaxation parameter $\theta$, where we find that the optimal relaxation
parameter in the case of $\kappa_{-}/\kappa_{+}\geq 2$ is given by
$\kappa_{+}/\kappa_{-}$.

###  5.3 Results

####  5.3.1 Convergence analysis

We first seek to verify that $\rho(K_{+}^{-1}SK_{-}^{-1}D)$ does indeed
determine convergence, and that the scaling law (59) is valid. We examine the
case $d$=2, $m$=2, $h_{+}/h_{-}$=8, in detail, preferring a two-dimensional
example in order to both run simulations and compute
$\rho(K_{+}^{-1}SK_{-}^{-1}D)$ quickly.

In Fig. 3, we plot the computed $\rho(K_{+}^{-1}SK_{-}^{-1}D)$ compared to
$\kappa_{-}/\kappa_{+}$ from the simulation pool. The results of (62) give
$a_{1}$=.4637, $a_{0}$=-.4637, and $b_{2}$=$8.84\text{e}-7$, confirming (59)
and implying that $\widetilde{C}\approx.4637$.

As pointed out, we should observe convergence for all
$\kappa_{-}/2<\kappa_{+}$. The critical value for the divergence should be for
$\widetilde{C}|\kappa_{-}/\kappa_{+}-1|\approx 1$, that means in this case
$\kappa_{-}/\kappa_{+}=1/\widetilde{C}+1\approx 3.16$. In Table 1, we report
the number of iterations necessary for convergence, computed
$\rho(K_{+}^{-1}SK_{-}^{-1}D)$, and predicted $\rho(K_{+}^{-1}SK_{-}^{-1}D)$
for a range of $\kappa_{-}/\kappa_{+}$. We define convergence as when relative
difference in $L^{2}$ norm between consecutive iterations of
$\bm{T}_{+}^{h,k}$ drops below 1e-8.

Referring to the data, we indeed observe the predicted behavior. The predicted
and computed $\rho(K_{+}^{-1}SK_{-}^{-1}D)$ are in perfect agreement,
providing strong evidence for (59).

This evidence strongly validates the theory detailed in previous sections.
Note as well that, in general, $\kappa_{-}<\kappa_{+}$ for our preferred
application of additive manufacturing, implying that for cases of this type
convergence should not be a problem [59].

##### Effects of relaxation

We note that the negative slope of the principal eigenvalue shown in Fig. 3
(note absolute value is shown, the signed value of the largest eigenvalue
become negative) implies that the relaxation scheme (49) will be effective for
this problem. In Fig 4, we show the impact of $\theta$ on (49). We therefore
expect that, even for larger gaps for which Fig. 3 indicates divergence will
occur, in Fig. 4, appropriate $\theta$ can resolve this problem and ensuring
convergence. We verify this numerically and, indeed, with $\theta$=.67 for
$\kappa_{-}/\kappa_{+}$=3.16, the simulation converged in 12 iterations. We
note also that Fig. 4 suggests this parameter can be optimized depending on
the problem, something that may be explored in future work.

Figure 3: Observed and predicted values of $\rho(K_{+}^{-1}SK_{-}^{-1}D)$ for $d$=2, $m$=2, $h_{+}/h_{-}$=8. Divergence predicted for $\kappa_{-}/\kappa_{+}\approx 3.16$. Figure 4: Predicted values of $\rho\left(\left(1-\theta\right)I_{+}+\theta K_{+}^{-1}SK_{-}^{-1}D\right)$ for $d$=2, $m$=2, $h_{+}/h_{-}$=8, various different values of $\theta$. $\kappa_{-}/\kappa_{+}$ | Predicted $\rho$ | Actual $\rho$ | Num. iterations
---|---|---|---
1.5 | .2319 | .2319 | 8
2.0 | .4637 | .4637 | 14
2.5 | .6956 | .6956 | 30
3.0 | .9274 | .9274 | 143
3.1 | .9738 | .9738 | 408
3.15 | .997 | .997 | 3988
3.16 | 1.001 | 1.001 | No convergence
Table 1: Testing the relationship between $\kappa_{-}/\kappa_{+}$ and
convergence. The predicted and computed $\rho(K_{+}^{-1}SK_{-}^{-1}D)$ are in
perfect agreement, and divergence occurs at the point predicted by
$\rho(K_{+}^{-1}SK_{-}^{-1}D)$.

####  5.3.2 Spectral growth analysis

Having confirmed with the previous that $\rho(K_{+}^{-1}SK_{-}^{-1}D)$
determines the convergence and scales according to (59), we now wish to
analyze the influence of $d$, $m$, and $h_{+}/h_{-}$ on
$\rho(K_{+}^{-1}SK_{-}^{-1}D)$. As we have strong evidence that (59) holds, we
expect the influence of the polynomial degree, spatial dimension and mesh
resolution on $\rho(K_{+}^{-1}SK_{-}^{-1}D)$ to be reflected in the value of
$\widetilde{C}$, predicted by (64). In particular, we are interested mostly in
how this value grows as we refine $h_{-}$. Ideally we would like to observe
that:

$\displaystyle\widetilde{C}$
$\displaystyle\sim\mathcal{O}\left(\log\frac{h_{+}}{h_{-}}\right).$ (65)

In view of (59), this indicates that as the local mesh is refined, the
convergence dose not deteriorate, something of great practical importance for
the application of the Two-level DD method.

In Fig. 5 we display the relationship between $\rho(K_{+}^{-1}SK_{-}^{-1}D)$
over the different $d$, $m$, and $h_{+}/h_{-}$. The plots appear indicate
that, for a given $d$ and $m$, $\rho(K_{+}^{-1}SK_{-}^{-1}D)$ displays a
logarithmic growth trend, as hoped.

Figure 5: Relationship between $\kappa_{-}/\kappa_{+}$ for different $d$, $m$
and $h_{+}/h_{-}$. The top row corresponds to 2D cases, the bottom to 3D; the
left row corresponds to linear polynomial approximations, the right to
quadratic. For each $d$ and $m$, the increase in
$\rho(K_{+}^{-1}SK_{-}^{-1}D)$ with $h_{+}/h_{-}$ appears logarithmic.  Figure
6: Growth of $\widetilde{C}$ as $h_{-}$ is refined across different $d$ and
$m$. For all cases, we observe a logarithmic growth trend, suggesting that the
performance of the algorithm is robust to refinement of the local mesh.

We confirm this logarithmic scaling behavior further in Fig. 6, which plots
the computed value of $\widetilde{C}$ against different $h_{+}/h_{-}$ levels
for each $d$ and $m$. The evidence here suggests that the degree of polynomial
approximation influences the spectral behavior far more than spatial
dimension, at least for this problem. In particular, we note that when $m$=2,
the value of $\widetilde{C}$ is larger than when $m$=1, yet it grows very
little.

The numerical experiments performed suggest several conclusions. First and
foremost, we observe strong validation of our theoretical prediction; namely,
that the convergence of the Two-level DD method is determined by
$\rho\left(K_{+}^{-1}SK_{-}^{-1}D\right)$. Further, this quantity depends
primarily on $\kappa_{+}$ and $\kappa_{-}$. Trivially, for
$\kappa_{+}\approx\kappa_{-}$ the convergence is fast, as this cases
approaches the no-jump case in the coefficients, perfectly approximated by a
monolithic approach. Looking over the entire dataset, we observed the
dependence of $\rho(K_{+}^{-1}SK_{-}^{-1}D)$ on the local mesh resolution to
be logarithmic, a critical observation, as we do not expect rapid performance
deterioration with mesh refinement.

###  5.4 Comparison with monolithic solver on a fitted mesh

In this section, we now compare the performance of the Two-level DD method
with the use of a monolithic solver in the entire domain. One of the clear
advantages of the Two-level DD method is in its ability to employ two uniform
meshes, each with homogeneous values of $\kappa$, potentially leading to
better conditioning properties. Such improvements in conditioning, when
sufficiently large, may render the use of a DD method competitive, even
despite the need to iterate.

We will proceed by considering the same basic two-dimensional setup shown in
the previous sections, now using Jacobi-preconditioned GMRES for all linear
solves. For the purposes of comparison, we consider the solution time and
necessary GMRES iterations for the Two-level solves and for the monolithic
solution over a range of $\kappa_{-}/\kappa_{+}$ and $h_{+}/h_{-}$, fixing
$\kappa_{+}$ as 1.0,$m=2$. The monolithic solution is performed on a fitted
mesh, such that the mesh features a size of $h_{-}$ in $\Omega_{-}$ and
$h_{+}$ in $\Omega_{+}$, where all nodes along the boundary between
$\Omega_{-}$ and $\Omega_{+}$ are conformal to the separation in the physical
domain

$\kappa_{-}/\kappa_{+}$ | $h_{+}/h_{-}$ | Local GMRES Its. | Global GMRES Its. | Monolithic GMRES Its. | Time (Two-lvl) | Time (monolithic)
---|---|---|---|---|---|---
1 | 8 | 73 | 290 | 3608 | .69s | 7.9s
.1 | 8 | 77 | 310 | 1557 | 1.11s | 3.56s
.01 | 8 | 99 | 382 | 1854 | 1.34s | 4.16s
1 | 16 | 176 | 290 | 6302 | 2.88s | 26.72s
.1 | 16 | 137 | 309 | 1916 | 4.61s | 8.59s
.01 | 16 | 222 | 381 | No convergence | 5.74s | NA
Table 2: Comparison with monolithic solver: observing the relationship between
$\kappa_{-}/\kappa_{+}$, $h_{+}/h_{-}$, and numerical solution behavior in
terms of linear solve difficulty and solution time.

We show the results in Table 2. We generally observe faster solution times and
better behavior with respect to changes in $\kappa_{-}/\kappa_{+}$ for the
Two-level DD solution. Notably, for $h_{+}/h_{-}$=16 and
$\kappa_{-}/\kappa_{+}$=.01, the Jacobi-preconditioned GMRES fails to converge
for the monolithic solution, while converging in a reasonable amount of
iterations for the Two-level DD case. The behavior with respect to mesh size
is a general increase in necessary iterations for each method, as expected.
However, for $\kappa_{-}/\kappa_{+}$ the effect is less clear; it is clear,
however, that the Two-level DD method appears only mildly sensitive in this
respect, with the number of necessary GMRES iterations remaining around the
same order of magnitude in all cases. In contrast, the monolithic method shows
highly variable and unpredictable behavior, likely a result of conditioning
issues brought on by the large heterogeneities present in the problem.

###  5.5 Extension to the nonlinear case

We will now consider the same general problem setup, but will introduce
nonlinear problems in which the thermal conductivity depends on the
temperature, with the nature of this dependence different in the local and
global domains. In [59], it was shown that, through some modification of the
transmission terms, the Two-level method is also consistent with the original
problem formulation in such a case. We refer the reader to this work for the
additional details. We note in this instance the necessity to iterate no
longer represents a serious issue, as the nonlinear nature of the problem
requires that some sort of iterative method be employed, regardless of whether
one solves in a monolithic fashion or using a Two-level DD method. While this
test is, strictly speaking, outside the analysis shown within the present
work, we nonetheless feel its inclusion is important, given the eventual
desired application of the proposed method to more complex problems.

In terms of setup, the problem is identical to those shown in the previous
sections, with the important distinction now that the thermal conductivity in
the local domain is considered as the powder form of stainless steel 316L (see
e.g. [58, 44]). For the global problem, we define $\kappa_{+}$ as the solid
form of stainless steel 316L in $\Omega_{A}$. This now introduces an
additional challenge, as we must extend the definition of $\kappa_{+}$ in
$\Omega_{B}$ in some way. We denote this extension as $\kappa_{+,B}$. For ease
of computation, it would be convenient to consider $\kappa_{+,B}$ as constant,
and we will examine the viability of this choice over a range of different
constants.

In each case, we solve a nonlinear problem using a Picard-type linearization
and each linear system using GMRES with a Jacobi-style preconditioner. We
compare the Two-level DD solution to a monolithic solution on a conformal
nonuniform mesh in terms of GMRES iterations, nonlinear iterations, and
overall solution time.

$\kappa_{+,B}$ | Avg. Global GMRES Iterations | Num. nonlinear iterations | Time (s)
---|---|---|---
.1 | 643 | 13 | 26.8
.2 | 621 | 11 | 26.1
.3 | 707 | 10 | 25.6
.4 | 728 | 12 | 28.15
.5 | 725 | 11 | 28.01
.6 | 721 | 12 | 28.87
.7 | 779 | 14 | 30.42
.8 | 875 | 15 | 32.21
Table 3: Nonlinear problem, observing the relationship between the constant
extension $\kappa_{+,B}$ of $\kappa_{B}$ and the difficulty of the linear and
nonlinear solves. We see that increasing $\kappa_{+,B}$ generally leads to
more difficult linear solves, while a $\kappa_{+,B}$ closer to .3 minimizes
the number of necessary nonlinear iterations.

In Tab. 3, we display the results of these tests. Solving the monolithic
reference problem, we note that the mean value of $\kappa_{+}$ in near
$\Omega_{B}$ is approximately .306, and hence we expect values for
$\kappa_{+,B}$ close to .306 to provide superior numerical performance. In
accordance with our expectation, we indeed find that extending $\kappa_{+,B}$
as a constant such that the jump between $\kappa_{+,A\setminus B}$ and
$\kappa_{+,B}$ is small reduces the necessary number of nonlinear iterations
for convergence. As for the difficulty of linear solves, the behavior is less
clear; lower values of $\kappa_{+,B}$ seem to result in easier linear solves,
though not necessarily fewer nonlinear iterations. For the purposes of
comparison, the monolithic reference solution required 8 nonlinear iterations,
each such iteration requiring an average of 6354 GMRES iterations, resulting
in an overall solve time of 126.3 seconds. In this test, we observe both
superior conditioning and solver performance when compared to the monolithic
approach, as well as an important numerical confirmation that extending
$\kappa_{B}$ as a constant inside $\Omega_{B}$ may be a reasonable choice.
Based on the results of these tests, if one does not have a good estimate for
what constant value to use, in general lower values seem to provide better
linear and nonlinear conditioning behavior.

##  6 Conclusions and Future Work

With this work, we present a theoretical analysis of the Two Level DD method
for the solution of heterogeneous material problems, obtained by combining a
non-overlapping Dirichlet-Neumann DD technique with a fictitious domain
approach. Under the simplifying assumption that the coefficients are piece-
wise constants, we have established that, at the continuous level, the
considered formulation converges to the desired solution. At the discrete
level, we showed that one may interpret an algorithm of this type as a block
Gauss-Seidel type iteration or as a truncated Neumann series, from which we
were able to postulate a convergence criterion and its scaling behavior. We
then performed several two- and three-dimensional simulations which validated
our predictions. Importantly, our numerical tests show that the convergence of
the method is expected to be robust with respect to local refinement, which is
important information of great practical interest for the application of such
algorithms. To demonstrate the potential of the approach, we also performed
some test in the nonlinear case, and made two important conclusions: that the
superior conditioning afforded by the DD method provides a numerical advantage
over applying a fitted non uniform mesh with a monolithic approach, and that
one may extend $\kappa_{+}$ as a constant inside $\Omega_{B}$ while
maintaining good numerical performance, provided such a value is chosen
carefully.

While this work is a first step in understanding the theoretical and numerical
features of the approach considered, many important questions remain. While we
numerically examined the behavior of the method on nonlinear problem, in the
theoretical analysis we only consider the steady problem and we rely on the
simplifying assumption that the coefficients are piece-wise constants. The
techniques shown here should extend promptly to unsteady problems with
constant coefficients, however, for non-constant coefficients and for
nonlinear problems, additional care must be taken in both the design of the
method and convergence analysis. This extends to the discrete problem as well,
as the Gauss-Seidel/Neumann series interpretations shown here will require
significant adaptation. Finally, the approach here considered can be extended
in different directions. To achieve higher efficiency, other domain
decomposition strategies may be used in combination with fictitious domains,
instead of the simple Dirichlet-Neumann iterations here considered. To achieve
optimal error, one might consider the possibility of giving up the requirement
that $T^{+}=T$ on $\Omega_{-}$ and look for an extension $f$ for which $T^{+}$
is smoother than $T$ (which has a jump in the normal derivative along
$\gamma_{i}$, thus limiting the $H^{1}(\Omega_{+})$ error to $h^{1/2}$).

##  7 Acknowledgments

This work was partially supported by the Italian Minister of University and
Research through the project ”A BRIDGE TO THE FUTURE: Computational methods,
innovative applications, experimental validations of new materials and
technologies” (No. 2017L7X3CS) within the PRIN 2017 program, as well as
Regione Lombardia through the project ”MADE4LO - Metal ADditivE for LOmbardy”
(No. 240963) within the POR FESR 2014-2020 program.

## References

  * [1] F. Aldakheel, B. Hudobivnik, E. Artioli, L. B. da Veiga, and P. Wriggers.  Curvilinear virtual elements for contact mechanics.  COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING, 372, DEC 1 2020.
  * [2] C. Atamian, Q. Dinh, R. Glowinski, J. He, and J. Periaux.  CONTROL APPROACH TO FICTITIOUS-DOMAIN METHODS APPLICATION TO FLUID-DYNAMICS AND ELECTRO-MAGNETICS.  In Glowinski, R and Kuznetsov, YA and Meurant, G and Periaux, J and Widlund, OB, editor, Fourth International Symposium on Domain Decomposition Methods for Partial Differential Equations, pages 275–309. ACAD SCI USSR; SOC IND & APPL MATH; SOC MATH APPL & IND, 1991.
  * [3] F. Auricchio, D. Boffi, L. Gastaldi, A. Lefieux, and A. Reali.  On a fictitious domain method with distributed Lagrange multiplier for interface problems.  APPLIED NUMERICAL MATHEMATICS, 95(SI):36–50, 2015.
  * [4] J. Baiges, M. Chiumenti, C. A. Moreira, M. Cervera, and R. Codina.  An adaptive finite element strategy for the numerical simulation of additive manufacturing processes.  Additive Manufacturing, page 101650, 2020.
  * [5] L. Beirão da Veiga, F. Brezzi, L. D. Marini, and A. Russo.  The hitchhiker’s guide to the virtual element method.  Mathematical models and methods in applied sciences, 24(08):1541–1573, 2014.
  * [6] T. Belytschko, R. Gracie, and G. Ventura.  A review of extended/generalized finite element methods for material modeling.  17(4):043001, apr 2009.
  * [7] C. Bertoglio, P. Moireau, and J. F. Gerbeau.  Sequential parameter estimation for fluid-structure problems. Application to hemodynamics.  International Journal for Numerical Methods in Biomedical Engineering, 28(4):434–455, 2012.
  * [8] S. Bertoluzza, F. Brezzi, and G. Sangalli.  The method of mothers for non-overlapping non-matching DDM.  NUMERISCHE MATHEMATIK, 107(3):397–431, SEP 2007.
  * [9] S. Bertoluzza, M. Ismail, and B. Maury.  Analysis of the fully discrete fat boundary method.  Numerische Mathematik, 118:49–77, 2011.
  * [10] F. Brezzi, K. Lipnikov, and V. Simoncini.  A family of mimetic finite difference methods on polygonal and polyhedral meshes.  Mathematical Models and Methods in Applied Sciences, 15(10):1533–1551, 2005.
  * [11] C. Bruna-Russo, A. Gökhan Demir, and B. Previtali.  Selective laser melting finite element modeling: validation with high–speed imaging and lack of fusion defects prediction.  Materials and Design, 156:143–153, 2018.
  * [12] E. Burman, D. Elfverson, P. Hansbo, M. G. Larson, and K. Larsson.  HYBRIDIZED CutFEM FOR ELLIPTIC INTERFACE PROBLEMS.  SIAM JOURNAL ON SCIENTIFIC COMPUTING, 41(5):A3354–A3380, 2019.
  * [13] M. Carraturo, C. Giannelli, A. Reali, and R. Vázquez.  Suitably graded thb-spline refinement and coarsening: Towards an adaptive isogeometric analysis of additive manufacturing processes.  Computer Methods in Applied Mechanics and Engineering, 348:660 – 679, 2019.
  * [14] L. Chen, H. Wei, and M. Wen.  An interface-fitted mesh generator and virtual element methods for elliptic interface problems.  JOURNAL OF COMPUTATIONAL PHYSICS, 334:327–348, APR 1 2017\.
  * [15] J. Côté, M. Gander, L. Laayouni, and S. Loisel.  Comparison of the dirichlet-neumann and optimal schwarz method on the sphere.  In Domain decomposition methods in science and engineering, pages 235–242. Springer, 2005.
  * [16] L. B. Da Veiga, C. Chinosi, C. Lovadina, and L. F. Pavarino.  Robust bddc preconditioners for reissner–mindlin plate bending problems and mitc elements.  SIAM journal on numerical analysis, 47(6):4214–4238, 2010.
  * [17] L. B. Da Veiga, D. Cho, L. Pavarino, and S. Scacchi.  Isogeometric schwarz preconditioners for linear elasticity systems.  Computer Methods in Applied Mechanics and Engineering, 253:439–454, 2013.
  * [18] L. B. Da Veiga, D. Cho, L. F. Pavarino, and S. Scacchi.  Overlapping schwarz methods for isogeometric analysis.  SIAM Journal on Numerical Analysis, 50(3):1394–1416, 2012.
  * [19] L. B. Da Veiga, C. Lovadina, and L. F. Pavarino.  Positive definite balancing neumann–neumann preconditioners for nearly incompressible elasticity.  Numerische Mathematik, 104(3):271–296, 2006.
  * [20] B. A. de Dios, K. Lipnikov, and G. Manzini.  The nonconforming virtual element method.  ESAIM: Mathematical Modelling and Numerical Analysis, 50(3):879–904, 2016.
  * [21] A. Duster, A. Niggl, and E. Rank.  Applying the hp-d version of the FEM to locally enhance dimensionally reduced models.  Computer Methods in Applied Mechanics and Engineering, 196:3524–3533, 2007.
  * [22] L. Giraud and R. Tuminaro.  Algebraic domain decomposition preconditioners.  Mesh partitioning techniques and domain decomposition methods, pages 187–216, 2006.
  * [23] R. Glowinski, T. Pan, T. Hesla, and D. Joseph.  A distributed Lagrange multiplier fictitious domain method for particulate flows.  INTERNATIONAL JOURNAL OF MULTIPHASE FLOW, 25(5):755–794, 1999.
  * [24] R. Glowinski, T. Pan, and J. Periaux.  A Fictitious Domain Method for Dirichlet Problem and Applications.  COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING, 111(3-4):283–303, 1994.
  * [25] G. H. Golub and C. F. Van Loan.  Matrix Computations.  The Johns Hopkins University Press, fourth edition, 2013.
  * [26] D. Gu, W. Meiners, K. Wissenbach, and R. Poprawe.  Laser additive manufacturing of metallic components: materials, processes, and mechanisms.  International Materials Reviews, 57(3):133–164, 2012.
  * [27] M. Gunzburger, M. Heinkenschloss, and H. Lee.  Solution of elliptic partial differential equations by an optimization-based domain decomposition method.  APPLIED MATHEMATICS AND COMPUTATION, 113(2-3):111–139, 2000\.
  * [28] M. Gunzburger, J. Peterson, and H. Kwon.  An optimization based domain decomposition method for partial differential equations.  Computers & Mathematics with Applications, 37(10):77–93, 1999\.
  * [29] Y. Huang, L. Yang, X. Du, and Y. Yang.  Finite element analysis of thermal behavior of metal powder during selective laser melting.  International Journal of Thermal Sciences, 104:146–157, 2016.
  * [30] J. Irwin and P. Michaleris.  A line heat input model for additive manufacturing.  Journal of Manufacturing Science and Engineering, 138(11), 2016\.
  * [31] M. Joulaian and A. Duester.  Local enrichment of the finite cell method for problems with material interfaces.  COMPUTATIONAL MECHANICS, 52(4):741–762, OCT 2013.
  * [32] K. Kergrene, I. Babuska, and U. Banerjee.  Stable Generalized Finite Element Method and associated iterative schemes; application to interface problems.  COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING, 305:1–36, JUN 15 2016.
  * [33] W. E. King, A. T. Anderson, R. M. Ferencz, N. E. Hodge, C. Kamath, S. A. Khairallah, and A. M. Rubenchik.  Laser powder bed fusion additive manufacturing of metals; physics, computational, and materials challenges.  Applied Physics Reviews, 2(4), 2015.
  * [34] S. Kollmannsberger, A. Ozcan, M. Carraturo, N. Zander, and E. Rank.  A hierarchical computational model for moving thermal loads and phase changes with applications to selective laser melting.  Computers and Mathematics with Applications, 75(5):1483–1497, 2018\.
  * [35] S. Kollmannsberger, A. Özcan, M. Carraturo, N. Zander, and E. Rank.  A hierarchical computational model for moving thermal loads and phase changes with applications to selective laser melting.  Comput. Math. Appl., 75(5):1483–1497, Mar. 2018.
  * [36] N. Labonnote, A. Ronnquist, B. Manum, and P. Ruther.  Additive construction: State-of-the-art, challenges and opportunities.  Automation in Construction, 72(3):347–366, 2016.
  * [37] R. Leveque and Z. Li.  The immersed interface method for elliptic-equations with discontinuous coefficients and singular sources.  SIAM JOURNAL ON NUMERICAL ANALYSIS, 31(4):1019–1044, 1994\.
  * [38] C. Li, E. R. Denlinger, M. F. Gouge, J. E. Irwin, and P. Michaleris.  Numerical verification of an octree mesh coarsening strategy for simulating additive manufacturing processes.  Additive Manufacturing, 30:100903, 2019.
  * [39] P.-L. Lions.  On the schwarz alternating method. iii: a variant for nonoverlapping subdomains.  In Third international symposium on domain decomposition methods for partial differential equations, volume 6, pages 202–223. SIAM Philadelphia, PA, 1990.
  * [40] P.-L. Lions et al.  On the schwarz alternating method. i.  In First international symposium on domain decomposition methods for partial differential equations, volume 1, page 42. Paris, France, 1988.
  * [41] W. K. Liu, D. W. Kim, and S. Tang.  Mathematical foundations of the immersed finite element method.  Computational Mechanics, 39(3):211–222, 2007.
  * [42] B. Maury.  A Fat boundary method for the Poisson problem in a domain with holes.  Journal of Scientific Computing, 16(3):319–339, 2001.
  * [43] J. Melenk and I. Babuska.  The partition of unity finite element method: Basic theory and applications.  COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING, 139(1-4):289–314, DEC 15 1996.
  * [44] K. Mills.  Recommended values of thermophysical properties for selected commercial alloys.  Woodhead Publishing Limited, 2002.
  * [45] L. Parussini and V. Pediroda.  Fictitious domain approach with $hp$-finite element approximation for incompressible fluid flow.  Journal of Computational Physics, 228:3891–3910, 2009.
  * [46] N. Patil, D. Pal, H. Khalid Rafi, K. Zeng, A. Moreland, A. Hicks, D. Beeler, and B. Stucker.  A generalized feed forward dynamic adaptive mesh refinement and derefinement finite element framework for metal laser sintering–Part I: Formulation and algorithm development.  Journal of Manufacturing Science and Engineering, 137(4), 2015.
  * [47] N. Patil, D. Pal, H. Khalid Rafi, K. Zeng, A. Moreland, A. Hicks, D. Beeler, and B. Stucker.  A Generalized Feed Forward Dynamic Adaptive Mesh Refinement and Derefinement Finite Element Framework for Metal Laser Sintering—Part I: Formulation and Algorithm Development.  Journal of Manufacturing Science and Engineering, 137(4):041001, Aug. 2015.
  * [48] N. Patil, D. Pal, K. Kutty, K. Zeng, A. Moreland, A. Hicks, D. Beeler, and B. Stucker.  A generalized feed forward dynamic adaptive mesh refinement and derefinement finite element framework for metal laser sintering–Part II: Nonlinear thermal simulations and validations.  Journal of Manufacturing Science and Engineering, 138(6), 2016.
  * [49] C. S. Peskin.  Numerical analysis of blood flow in the heart.  Journal of computational physics, 25(3):220–252, 1977.
  * [50] A. Quarteroni and A. Valli.  Domain Decomposition Methods for Partial Differential Equations.  Oxford University Press, Oxford, UK, 1999.
  * [51] I. Ramiere, P. Angot, and M. Belliard.  A fictitious domain approach with spread interface for elliptic problems with general boundary conditions.  Computer Methods in Applied Mechanics and Engineering, 196:766–781, 2007.
  * [52] E. Rank.  Adaptive remeshing and $hp$ domain decomposition.  Computer Methods in Applied Mechanics and Engineering, 101(1–3):299–313, 1992.
  * [53] D. Riedlbauer, P. Steinmann, and J. Mergheim.  Thermomechanical finite element simulations of selective electron beam melting processes: Performance considerations.  Computational Mechanics, 54(1):109–122, 2014.
  * [54] M. Russell, A. Suoto-Iglesias, and T. Zohdi.  Numerical simualtion of Laser Fusion Additive Manufacturing processes using the SPH method.  Computer Methods in Applied Mechanics and Engineering, 341:163–187, 2018.
  * [55] D. Schillinger, A. Duster, and E. Rank.  The $hp$-$d$ adaptive finite cell method for geometrically nonlinear problems of solid mechanics.  International Journal for Numerical Methods in Engineering, 89:1171–1202, 2012.
  * [56] S. Tammas-Williams and I. Todd.  Design for additive manufacturing with site-specific properties in metals and alloys.  Scripta Materiala, 135:105–110, 2017.
  * [57] A. Toselli and O. Widlund.  Domain decomposition methods-algorithms and theory, volume 34.  Springer Science & Business Media, 2004.
  * [58] A. Viguerie and F. Auricchio.  Numerical solution of additive manufacturing problems using a two-level method.  International Journal for Numerical Methods in Engineering, 2021\.
  * [59] A. Viguerie, S. Bertoluzza, and F. Auricchio.  A fat boundary-type method for localized nonhomogeneous material problems.  Computer Methods in Applied Mechanics and Engineering, 364:112983, 2020.
  * [60] Q. Wang, J. Li, M. Gouge, A. Nassar, P. Michaleris, and E. Reutzel.  Physics-based multivariable modeling and feedback linearization control of melt-pool geometry and temperature in directed energy deposition.  Journal of Manufacturing Science and Engineering, 139(2), 2017.
  * [61] N. Zander, T. Bog, M. Elhaddad, F. Frischmann, S. Kollmannsberger, and E. Rank.  The multi-level $hp$ method for three-dimensional problems: Dynamically changing high-order mesh refinement with arbitrary hanging nodes.  Computer Methods in Applied Mechanics and Engineering, 310:252–277, 2016.