4. Domain decomposition methods

Domain Decomposition Methods (DDM) encompass highly efficient algorithms to obtain solutions of large-scale discrete problems on parallel super-computers. They mainly consist of partitioning the domain into various subdomains and then getting the global solution through the resolution of the subdomain problems [12, 13] often in an iterative fashion. These methods can be seen as an iterative coupling by the internal and thus unknown BCS. There is a broad literature covering these approaches, and that is why this chapter, therefore, presents a short introduction for the sake of completeness. The recommended references include Bjorstad and Widlund [14], Bramble et al. and Marini and Quarteroni [15], who considered the Dirichlet-Neumann (DN) DDM and Neumann-Neumann.

Let L be an abstract linear differential operator, such as the Laplace operator, for instance. The DN-DDM scheme implies solving a series of problems in the proper sequence that requires a coloring tool (see Figure 1). Let the Dirichlet subdomains be colored in white while the Neumann subdomains are in black. Notice that the interface between subdomains is denoted by Γ. After one provides the initial guess on the primary variables on Γ, i.e., γ<sup>k</sup> must be given, then one can solve the problem on the white subdomains (Dirichlet problems), which corresponds to step 1 in Eq. (25).

$$(1)\begin{cases} Lu\_1^{k+1} = f \text{ in } \Omega\_1\\ u\_1^{k+1} = 0 \text{ on } \partial\Omega\_1 \cap \partial\Omega\\ u\_1^{k+1} = \gamma^k \text{ on } \Gamma \end{cases} \qquad (2)\begin{cases} Lu\_2^{k+1} = f \text{ in } \Omega\_2\\ u\_2^{k+1} = 0 \text{ on } \partial\Omega\_2 \cap \partial\Omega\\ \partial\_n u\_2^{k+1} = \kappa^{k+1} \text{ on } \Gamma \end{cases} \tag{25}$$

Let the primary variable be called "displacements" and their gradient "tractions," i.e., normal derivative in the boundary. Then, the tractions on the interface Γ must be computed after first solving step 1 on the white subdomains. They are then passed through communication to solve the second step on the black subdomains, i.e., Neumann subdomains. On this latter, since the tractions are known on Γ, one can solve for unknown displacements to provide feedback

Figure 1. It depicts the DNDDM.

on the next iteration level. Both displacements and tractions are often over-relaxed to improve the convergence rate. The given relaxation parameters, referred in Eq. (26) as θ<sup>D</sup> and θ<sup>N</sup>, must lie between 0 and 1:

$$\begin{aligned} \kappa^{(k+1)} &= \left(-\theta^{N} \cdot \eth\_{n} \mu\_{1}^{k+1} + \left(1 - \theta^{N}\right) \cdot \eth\_{n} \mu\_{2}^{k}\right) \text{ on } \Gamma\\ \gamma^{(k+1)} &= \left(\theta^{D} \cdot \mu\_{2}^{k+1} + \left(1 - \theta^{D}\right) \cdot \mu\_{1}^{k}\right) \text{ on } \Gamma. \end{aligned} \tag{26}$$

referred to those references that cover the topics of computational geometry, in particular how to build these NURBS entities. Let MFEM be described for linear isotropic elasticity regarding

ð

ð

Linear Thermo-Poroelasticity and Geomechanics http://dx.doi.org/10.5772/intechopen.71873

> Ω f <sup>T</sup> � vdx

(27)

231

(28)

∂Ω<sup>N</sup> t <sup>T</sup> � vds <sup>þ</sup>

<sup>T</sup> � <sup>C</sup> � <sup>ε</sup>ð Þ <sup>u</sup> dx ; lð Þ¼ <sup>v</sup>

<sup>T</sup> � <sup>Φ</sup>ds ; ½ �¼ <sup>u</sup> <sup>u</sup>ð Þ<sup>1</sup> � <sup>u</sup>ð Þ<sup>2</sup> � �

where ϒ stands for the gluing condition among subdomain interfaces and the jump ½ � u on the

� � <sup>þ</sup> <sup>ϒ</sup> vh; <sup>Λ</sup><sup>h</sup> ð Þ¼ <sup>l</sup> vh ð Þ

the parameters in Eq. (28) are as follows: Φ<sup>h</sup> represents the mortar space while vh corresponds to the weighting space and Λ<sup>h</sup> is the Lagrange multiplier space, i.e., the linear combination of mortar functions, often polynomial functions, and Lagrange multiplier degrees of freedom. Let

<sup>M</sup> be a conforming partition of the so-called parametric space, <sup>Ω</sup>, whose image serves as the mortar's geometrical entity, i.e., curve or surface, composed of line-segments (n ¼ 2Þ or quadrilaterals (n ¼ 3). One takes the mortar space as a finite-dimensional subspace of the continu-

bilinear forms, a and ϒ defined in Eq. (27) below [2, 3],

ð

Figure 2. Ω<sup>1</sup> is in the top, Ω<sup>2</sup> is in the bottom, and the interface Γ is the bold curve.

Ω εð Þv

ð

Γ ½ � u

displacements is required to vanish in an integral or "weak" sense, thus:

(

a uh; vh

ϒ uh; Φ<sup>h</sup> � � <sup>¼</sup> <sup>0</sup>

a u; v � � <sup>¼</sup>

ϒ u; Φ � � <sup>¼</sup>

T <sup>h</sup>

ous Sobolev spaces, that is:

It happens that this approach requires at least a two-entry coloring tool or even more, i.e., there may be subdomains with mixed interfaces, colored as gray [12]. There is a lack of parallelism in the sense that black subdomains must wait for the white ones to communicate their tractions. An initial guess for tractions should be prescribed to mitigate this issue, but this latter is not feasible in most cases. A straightforward way to obtain an initial estimate for the multiplier γ<sup>k</sup> is by computing the so-called coarse-run that implies solving the same problem in a coarser mesh and interpolating over Γ by using the smaller's problem FE space. The reader may refer to the literature [16, 17] for further reading and proof of convergence and also revise [2] for a more detailed description that includes implementation details, which this chapter omits for the sake of brevity.

### 5. The mortar FE method (MFEM)

The primary goal here is to extend MFEM to glue curved interfaces such as the one shown in Figure 2 where MFEM treats non-matching interfaces. The section first introduces a brief description of non-uniform rational B-Splines curves and surfaces (NURBS) in [2, 3, 18]. The reader is

Figure 2. Ω<sup>1</sup> is in the top, Ω<sup>2</sup> is in the bottom, and the interface Γ is the bold curve.

on the next iteration level. Both displacements and tractions are often over-relaxed to improve the convergence rate. The given relaxation parameters, referred in Eq. (26) as θ<sup>D</sup> and θ<sup>N</sup>, must

> <sup>k</sup>þ<sup>1</sup> <sup>þ</sup> <sup>1</sup> � <sup>θ</sup><sup>D</sup> � <sup>u</sup><sup>1</sup> <sup>k</sup> on Γ:

It happens that this approach requires at least a two-entry coloring tool or even more, i.e., there may be subdomains with mixed interfaces, colored as gray [12]. There is a lack of parallelism in the sense that black subdomains must wait for the white ones to communicate their tractions. An initial guess for tractions should be prescribed to mitigate this issue, but this latter is not feasible in most cases. A straightforward way to obtain an initial estimate for the multiplier γ<sup>k</sup> is by computing the so-called coarse-run that implies solving the same problem in a coarser mesh and interpolating over Γ by using the smaller's problem FE space. The reader may refer to the literature [16, 17] for further reading and proof of convergence and also revise [2] for a more detailed description that includes implementation details, which this chapter omits for

The primary goal here is to extend MFEM to glue curved interfaces such as the one shown in Figure 2 where MFEM treats non-matching interfaces. The section first introduces a brief description of non-uniform rational B-Splines curves and surfaces (NURBS) in [2, 3, 18]. The reader is

<sup>k</sup>þ<sup>1</sup> <sup>þ</sup> <sup>1</sup> � <sup>θ</sup><sup>N</sup> � <sup>∂</sup>nu<sup>2</sup> <sup>k</sup> on Γ

(26)

<sup>κ</sup>ð Þ <sup>k</sup>þ<sup>1</sup> ¼ �θ<sup>N</sup> � <sup>∂</sup>nu<sup>1</sup>

<sup>γ</sup>ð Þ <sup>k</sup>þ<sup>1</sup> <sup>¼</sup> <sup>θ</sup><sup>D</sup> � <sup>u</sup><sup>2</sup>

230 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

lie between 0 and 1:

Figure 1. It depicts the DNDDM.

the sake of brevity.

5. The mortar FE method (MFEM)

referred to those references that cover the topics of computational geometry, in particular how to build these NURBS entities. Let MFEM be described for linear isotropic elasticity regarding bilinear forms, a and ϒ defined in Eq. (27) below [2, 3],

$$\begin{aligned} a\left(\underline{\underline{u},v}\right) &= \int\_{\Omega} \underline{\underline{c}(\underline{\underline{u}})}^T \cdot \underline{\underline{C}} \cdot \underline{\underline{c}(\underline{\underline{u}})} d\mathfrak{x} \; ; \; l(\underline{\underline{v}}) = \int\_{\partial\Omega^{\mathbb{N}}} \underline{\underline{t}}^T \cdot \underline{\underline{v}} d\mathbf{s} + \int\_{\Omega} \underline{\underline{t}}^T \cdot \underline{\underline{v}} d\mathbf{x} \\\ \Upsilon\left(\underline{\underline{u}}, \underline{\underline{\underline{\underline{D}}}}\right) &= \int\_{\Gamma} \underline{[\underline{u}]}^T \cdot \underline{\underline{\underline{D}}} d\mathbf{s} \; ; \; \; \left[\underline{\underline{u}}\right] = \left(\underline{\underline{u}}^{(1)} - \underline{\underline{u}}^{(2)}\right) \end{aligned} \tag{27}$$

where ϒ stands for the gluing condition among subdomain interfaces and the jump ½ � u on the displacements is required to vanish in an integral or "weak" sense, thus:

$$\begin{cases} a(\underline{u}\_h, \underline{v}\_h) + \Upsilon(\underline{v}\_h, \underline{\Delta}\_h) = l(\underline{v}\_h) \\ \Upsilon(\underline{u}\_h, \underline{\Delta}\_h) = 0 \end{cases} \tag{28}$$

the parameters in Eq. (28) are as follows: Φ<sup>h</sup> represents the mortar space while vh corresponds to the weighting space and Λ<sup>h</sup> is the Lagrange multiplier space, i.e., the linear combination of mortar functions, often polynomial functions, and Lagrange multiplier degrees of freedom. Let T <sup>h</sup> <sup>M</sup> be a conforming partition of the so-called parametric space, <sup>Ω</sup>, whose image serves as the mortar's geometrical entity, i.e., curve or surface, composed of line-segments (n ¼ 2Þ or quadrilaterals (n ¼ 3). One takes the mortar space as a finite-dimensional subspace of the continuous Sobolev spaces, that is:

$$\mathcal{L}\_k\left(\overline{T}\_h^{\ M}\right) = \left\{\Phi \in L^2\left(\overline{\Omega}\right) \,:\, \forall \mathbf{e}^{\mathcal{M}} \in \overline{T}\_h^{\ M}, \Phi|\_{\boldsymbol{\varepsilon}^{\mathcal{M}}} \in \mathbb{P}\_k\left(\boldsymbol{e}^{\mathcal{M}}\right)\right\} \tag{29}$$

discretization in each subdomain, <sup>P</sup>kð Þ<sup>e</sup> , as well as mortar space <sup>P</sup><sup>k</sup> <sup>e</sup><sup>M</sup> . It also utilizes piecewise linear polynomials for the space discretizations in all examples herein that were run on a MacBook Pro laptop equipped with an Intel(R) Quad-Core(TM) i7-4870HQ CPU @ 2.5 GHz and 16 GB of RAM. The author chose this laptop for the sake of convenience, in particular, the availability of debugging tools free of charge, such as the Microsoft Visual Studio Community.

The example is a manufactured problem where the solution is a priori chosen. Then, one substitutes the given pressure field in the governing partial differential equation to obtain the source term, i.e., loading, that reproduces the input field. The problem in strong form looks like:

where the domain of interest corresponds to the unitary square and Dirichlet BCS are enforced.

Figure 3 shows the pressure field that corresponds to the problem 6.1 whose discretization encompasses three subdomains: two of them (the top and bottom ones) consist of triangular meshes while the one in the middle was discretized by a regular Cartesian quadrilateral mesh.

The pressure field is on the right-top corner, and its horizontal derivative is in the bottom-left corner, while the discrepancy between the numerical and exact solutions, i.e., the absolute

<sup>¼</sup> <sup>f</sup> in <sup>Ω</sup> ; p <sup>¼</sup> <sup>p</sup><sup>0</sup> on <sup>Γ</sup><sup>D</sup> <sup>¼</sup> <sup>Γ</sup>, (33)

Linear Thermo-Poroelasticity and Geomechanics http://dx.doi.org/10.5772/intechopen.71873 233

p xð Þ¼ ; <sup>y</sup> xy � ð Þ� <sup>x</sup> � <sup>1</sup> ð Þ� <sup>y</sup> � <sup>1</sup> exp � <sup>x</sup><sup>2</sup> <sup>þ</sup> <sup>y</sup><sup>2</sup> ; <sup>K</sup> <sup>¼</sup> <sup>I</sup>: (34)

Aside, one can achieve some level of parallelism due to the multi-core technology.

6.1. Example 1: Two-dimensional steady state single-phase flow

�∇ � K∇p

The top-left corner of the figure shows the mesh that is employed.

The input pressure field is given by:

Figure 3. The MFEM solution to problem 6.1.

herein P<sup>k</sup> e<sup>M</sup> � � stands for the space of polynomials of total degree less than or equal to k while C<sup>k</sup> T <sup>h</sup> <sup>M</sup> � � represents test functions that are continuous along the edges of <sup>e</sup><sup>M</sup>.

One can write in a matrix or algebraic form, Eq. (28) as:

$$
\begin{bmatrix}
\begin{bmatrix} k^{(1)} \\ & \end{bmatrix} & \begin{bmatrix} \mathbf{Y}^{(1)} \end{bmatrix}^{T} \\
\begin{bmatrix} \mathbf{0} \end{bmatrix} & \begin{bmatrix} k^{(2)} \\ \end{bmatrix} & -\begin{bmatrix} \mathbf{Y}^{(2)} \end{bmatrix}^{T} \\
\begin{bmatrix} \mathbf{Y}^{(1)} \end{bmatrix} & -\begin{bmatrix} \mathbf{Y}^{(2)} \end{bmatrix} & \begin{bmatrix} \underline{\mathbf{u}}^{(1)} \\ \underline{\mathbf{u}}^{(2)} \end{bmatrix} \\
\end{bmatrix} \cdot \begin{bmatrix} \underline{\mathbf{u}}^{(1)} \\ \end{bmatrix} = \begin{bmatrix} \underline{\mathbf{l}}^{(1)} \\ \underline{\mathbf{l}}^{(2)} \\ \end{bmatrix} \cdot \tag{30}
$$

The equation above corresponds to the so-called "saddle-point problem (SPP)." Notice that subdomains are only connected using the Lagrange multiplier Λ if they happen to be known (it is well-known that for elasticity, the multipliers are the unknown tractions on the interface), then one can decouple the system in Eq. (30) and then one just needs to perform subdomain solves. For the SPP (30), one may match displacements or tractions in the interface. The Dirichlet-Neumann scheme that the section presents is only a particular case of the most general Robin-Robin domain decomposition scheme [2, 3]. The rectangular matrices <sup>ϒ</sup>ð Þ<sup>i</sup> � �, i <sup>¼</sup> <sup>1</sup>…2, are denoted as projectors since they permit to map to and from the given mortar space [2, 3].

The following line integral defines the projector, for 2-D problems, as:

$$\Upsilon\_{ij}^{(k)} = \underset{\overline{\Omega}}{\int\_{\overline{\Omega}} \varphi\_j^{(k)}(\xi) \Phi\_i(\xi) \cdot \|d\_\xi \underline{\mathcal{L}}(\xi)\| \, d\xi}{} \tag{31}$$

where wð Þ<sup>k</sup> <sup>j</sup> represents the global non-mortar side interpolation functions and Φ<sup>i</sup> are the mortar space basis functions, while k k dξC is the length of the tangent vector associated to the B-Spline or NURBS curve. Similarly, 3-D problems imply:

$$\Upsilon\_{ij}^{(k)} = \iint \varphi\_j^{(k)} \left( \underline{\xi} \right) \Phi\_i(\underline{\xi}) \cdot \left| \left| \partial\_{\xi} \underline{\xi}(\underline{\xi}) \times \partial\_{\eta} \underline{\xi}(\underline{\xi}) \right| \right| d\xi \, d\eta \tag{32}$$

where <sup>∂</sup>ξ<sup>S</sup> � <sup>∂</sup>η<sup>S</sup> � � � � is the norm of the surface's normal vector. Particular quadrature rules to compute these integrals must be developed. The reader should refer to [3] for a detailed explanation including the proper algorithm in pseudo-code.

### 6. Numerical examples

The author implemented these FE models in the Integrated Parallel Finite Element Analysis program (IPFA) that is a C++ application whose main characteristics are described in [2, 12]. IPFA employs standard continuous Lagrange polynomials as shape functions for the space discretization in each subdomain, <sup>P</sup>kð Þ<sup>e</sup> , as well as mortar space <sup>P</sup><sup>k</sup> <sup>e</sup><sup>M</sup> . It also utilizes piecewise linear polynomials for the space discretizations in all examples herein that were run on a MacBook Pro laptop equipped with an Intel(R) Quad-Core(TM) i7-4870HQ CPU @ 2.5 GHz and 16 GB of RAM. The author chose this laptop for the sake of convenience, in particular, the availability of debugging tools free of charge, such as the Microsoft Visual Studio Community. Aside, one can achieve some level of parallelism due to the multi-core technology.

### 6.1. Example 1: Two-dimensional steady state single-phase flow

C<sup>k</sup> T <sup>h</sup> <sup>M</sup> � �

One can write in a matrix or algebraic form, Eq. (28) as:

232 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

½ � 0 k

ð Þ<sup>2</sup> h i

<sup>ϒ</sup>ð Þ<sup>1</sup> � � � <sup>ϒ</sup>ð Þ<sup>2</sup> � � ½ � <sup>0</sup>

<sup>k</sup>ð Þ<sup>1</sup> h i

C<sup>k</sup> T <sup>h</sup> <sup>M</sup> � �

where wð Þ<sup>k</sup>

where <sup>∂</sup>ξ<sup>S</sup> � <sup>∂</sup>η<sup>S</sup> � � �

6. Numerical examples

<sup>¼</sup> <sup>Φ</sup> <sup>∈</sup>L<sup>2</sup> <sup>Ω</sup>

� � : ∀e<sup>M</sup> ∈T <sup>h</sup>

herein P<sup>k</sup> e<sup>M</sup> � � stands for the space of polynomials of total degree less than or equal to k while

� <sup>ϒ</sup>ð Þ<sup>2</sup> � �<sup>T</sup>

The equation above corresponds to the so-called "saddle-point problem (SPP)." Notice that subdomains are only connected using the Lagrange multiplier Λ if they happen to be known (it is well-known that for elasticity, the multipliers are the unknown tractions on the interface), then one can decouple the system in Eq. (30) and then one just needs to perform subdomain solves. For the SPP (30), one may match displacements or tractions in the interface. The Dirichlet-Neumann scheme that the section presents is only a particular case of the most general Robin-Robin domain decomposition scheme [2, 3]. The rectangular matrices <sup>ϒ</sup>ð Þ<sup>i</sup> � �, i <sup>¼</sup> <sup>1</sup>…2, are denoted

<sup>j</sup> represents the global non-mortar side interpolation functions and Φ<sup>i</sup> are the mortar

� � � <sup>∂</sup>η<sup>S</sup> <sup>ξ</sup> � � � � �

� is the norm of the surface's normal vector. Particular quadrature rules to

space basis functions, while k k dξC is the length of the tangent vector associated to the B-Spline

� � � <sup>∂</sup>ξ<sup>S</sup> <sup>ξ</sup>

compute these integrals must be developed. The reader should refer to [3] for a detailed

The author implemented these FE models in the Integrated Parallel Finite Element Analysis program (IPFA) that is a C++ application whose main characteristics are described in [2, 12]. IPFA employs standard continuous Lagrange polynomials as shape functions for the space

represents test functions that are continuous along the edges of e<sup>M</sup>.

½ � <sup>0</sup> <sup>ϒ</sup>ð Þ<sup>1</sup> � �<sup>T</sup>

as projectors since they permit to map to and from the given mortar space [2, 3].

Ω

wð Þ<sup>k</sup>

The following line integral defines the projector, for 2-D problems, as:

ϒð Þ<sup>k</sup> ij ¼ ð

Ω

explanation including the proper algorithm in pseudo-code.

wð Þ<sup>k</sup> <sup>j</sup> ξ � �Φ<sup>i</sup> ξ

or NURBS curve. Similarly, 3-D problems imply:

ϒð Þ<sup>k</sup> ij ¼ ðð <sup>M</sup>; <sup>Φ</sup><sup>j</sup>

uð Þ<sup>1</sup> uð Þ<sup>2</sup> Λ

l ð Þ1 l ð Þ2 0

<sup>j</sup> ð Þ ξ Φið Þ� ξ k k dξCð Þ ξ dξ (31)

3 7

<sup>5</sup>� (30)

�dξdη (32)

2 6 4

2 6 4

eM ∈ P<sup>k</sup> e <sup>M</sup> n o � � (29)

> The example is a manufactured problem where the solution is a priori chosen. Then, one substitutes the given pressure field in the governing partial differential equation to obtain the source term, i.e., loading, that reproduces the input field. The problem in strong form looks like:

$$-\nabla \cdot \left(\underline{\underline{K}} \nabla p\right) = f \quad \text{in } \Omega \; ; \; p = p\_0 \text{ on } \Gamma\_\mathcal{D} = \Gamma,\tag{33}$$

where the domain of interest corresponds to the unitary square and Dirichlet BCS are enforced. The input pressure field is given by:

$$p(\mathbf{x}, y) = \mathbf{x}y \cdot (\mathbf{x} - \mathbf{1}) \cdot (y - \mathbf{1}) \cdot \exp\left[-\left(\mathbf{x}^2 + y^2\right)\right]; \underline{\mathbf{K}} = \underline{\mathbf{I}}.\tag{34}$$

Figure 3 shows the pressure field that corresponds to the problem 6.1 whose discretization encompasses three subdomains: two of them (the top and bottom ones) consist of triangular meshes while the one in the middle was discretized by a regular Cartesian quadrilateral mesh. The top-left corner of the figure shows the mesh that is employed.

The pressure field is on the right-top corner, and its horizontal derivative is in the bottom-left corner, while the discrepancy between the numerical and exact solutions, i.e., the absolute

Figure 3. The MFEM solution to problem 6.1.


This slope agrees with the theory that predicts a rate of <sup>O</sup> <sup>h</sup><sup>3</sup>=<sup>2</sup> [2, 3]. However, the resulting slope is slightly lower because of numerical errors, such as quadrature and linear system

Linear Thermo-Poroelasticity and Geomechanics http://dx.doi.org/10.5772/intechopen.71873 235

Finally, Figure 5 shows pressure snapshots that represent four different Dirichlet-Neumann iteration levels evolving from left-to-right and top-to-bottom. The fact that no initial guess for pressure was provided explains the mismatch in the first snapshot. That is why one needs to eliminate discrepancies by running the process to match up those subdomains, i.e., the traction residual in the interface must vanish, which for this case occurs in just a handful of iterations. The stop criterion precisely involves the residual in the tractions in the interface that is required to fall below the given tolerance. For this particular problem, the iterative process spent six

This example analyzes a coupled flow and mechanics simulation in a reconstructed reservoir (RS) model with different meshes for the flow and mechanics physics [18]. The author proposed such a reconstruction workflow in [18] which permits this latter feature by computing a projection operator to mapping pressures from the original flow mesh into the so-called reference mechanics mesh. Toward that end, the example employs the slightly compressible flow formulation loosely combined with the mechanics model as shown in Eq. (15). The objective is to show a realistic field level RS compaction and subsidence coupled computation. The goal is thus working three different cases for the mechanics part in which one only changes the resolution of the reconstructed mechanics mesh in the pay-zone while preserving the mechanical properties constant as well as the geometry, BCS, and the depletion scenario. The exercise admits the actual static properties as being in the pay-zone such as porosity f and permeability for the isotropic case Kx ¼ Kz ¼ Ky as shown in Figure 6, whose depiction is three

solving errors.

iterations to achieve a residual lower than 10�<sup>6</sup> [3].

6.2. Example 2: Coupled flow and mechanics

Figure 5. The numerical L<sup>2</sup> convergence rate for problem 6.1.

Table 1. Meshes for example 6.1.

error, was rendered in the right-bottom corner. Table 1 represents the number of elements and points of each mesh from top to bottom. The mortars as geometrical entities correspond to two B-Splines interpolants (NURBS with all weights equal 1) that were constructed by interpolating a sinusoidal wave as the figure shows (see [3] for details). Thirty-two quadratic mortar elements per curve were utilized to glue these three subdomains. A direct frontal solver was used to solve the global SPP in Eq. (30) [3]. The results that are summarized on Figure 3 are in good agreement with the analytical solution. The absolute error against the correct answer is also displayed. The discrepancy is of the order of 10�<sup>4</sup> . Notice that besides the example only matched the displacements on the interface, a good accordance is also obtained for the horizontal derivative.

Whether or not one utilizes the SPP approach, the local problems are completely disconnected. This fact can be exploited to reduce the computational time significantly. Indeed, these sub problems can be handled in separate threads using a shared memory approach, i.e., multithreading assembling. A convergence analysis was also performed, by successively running refined meshes [3] and by keeping a refinement ratio of 2:1 between subdomains. The exercise used a piecewise quadratic mortar space where the number of mortar elements equals the number of coarse edges in the non-mortar sides. It tackled meshes of size 8, 16, 32, 64, 128 and 256 respectively. Figure 4 displays the resulting convergence rate in a log � log plot. The slope of the least-squares straight line is 1.44143, where the coefficient of determination is <sup>R</sup><sup>2</sup> <sup>¼</sup> <sup>84</sup>%.

Figure 4. Snapshots showing the evolution of the DN-DDM applied to problem 6.1.

This slope agrees with the theory that predicts a rate of <sup>O</sup> <sup>h</sup><sup>3</sup>=<sup>2</sup> [2, 3]. However, the resulting slope is slightly lower because of numerical errors, such as quadrature and linear system solving errors.

Finally, Figure 5 shows pressure snapshots that represent four different Dirichlet-Neumann iteration levels evolving from left-to-right and top-to-bottom. The fact that no initial guess for pressure was provided explains the mismatch in the first snapshot. That is why one needs to eliminate discrepancies by running the process to match up those subdomains, i.e., the traction residual in the interface must vanish, which for this case occurs in just a handful of iterations. The stop criterion precisely involves the residual in the tractions in the interface that is required to fall below the given tolerance. For this particular problem, the iterative process spent six iterations to achieve a residual lower than 10�<sup>6</sup> [3].
