This paper pushes further the intrinsic capabilities of the global-local approach  introduced initially in [7].
We develop a distributed computing approach using MPI (Message Passing Interface) both for the global and local problems.
Regarding local problems, a specific scheduling strategy is introduced.
Then, to measure correctly the convergence of the iterative process, we introduce a reference solution that
revisits the product of classical and enriched functions.
As a consequence, we are able to propose a purely matrix-based implementation of the global-local problem.
The distributed approach is then compared to other parallel solvers
either direct or iterative with domain decomposition.
The comparison addresses the scalability as well as the elapsed time.
Numerical examples deal with linear elastic problems: a polynomial exact solution problem, a complex micro-structure, and, finally, a pull-out test (with different crack extent).

§ INTRODUCTION

Computational methods are addressing more and more complex  problems.
This paper is concerned with elasticity problems having roughly two-scales.
A global scale at the structural level and a local scale representing a micro-structure, a crack or a reinforcement for example.
Note that the local scale  appears over the whole domain
in the case of a micro-structure and only in narrow zones for cracks and reinforcements.
A direct finite element approach for these problems is to create a unique mesh and solve.
This approach does not take advantage of the two-scale nature.
This leads to difficulties: the matrix can have a huge size and can be badly conditioned.
Iterative domain decomposition resolution is one way to face huge matrix size (see [30]) but does not solve the conditioning issue unless a specific preconditioner is used.
A direct solver can cope with poor conditioning, but matrix size is still an issue.

In order to take into account the two-scale nature of a problem,  the first existing  global to local method (see for example [23]) use the raw  global scale field as a boundary condition for local problems (refined encapsulated meshes), which provide the desired accuracy to the simulation.
But since the raw global field does not take into account the local behavior, it does not provide a good boundary condition to the local problems.
The method, first introduced in [7], adds the idea of local to global interaction via enrichment functions.
More precisely,  local problems, which can overlap, yield solutions used to build enrichment functions.
These functions are  used through the partition of unity method to compute at global scale an improved solution.
The overlapping aspect of the local problems and their intermixing use in the  enriched global problem permit to partially address the boundary condition error of the raw global problem field.
In [6] this problem of boundary condition has been addressed by increasing the size of the local problems (the boundary condition at the local scale is pushed far enough) or by increasing the polynomial interpolation order  of the  global scale problem.
But it is in [24] that a more practical and robust solution to this boundary problem has been proposed by introducing iterations between scales.
From the raw results, boundary conditions are imposed on the local problems, whose solutions are used as enrichment functions for a new enriched global problem, whose solution is again used as boundary conditions for local problems and finally their solutions reused for a last new enriched global problem.

Later in [26] the idea arises that in evolutionary phenomena such as crack propagation, this scale loop can be intertwined with a propagation loop.
But it was only in [14] that the looping strategy was further investigated for a larger number of scale iterations as a remedy for the boundary problem, in a context of a crack propagation step with a single local problem.
This work shows that in a few iterations, the method provides a significant error reduction at both scales.
But, as in many other publication, the way this error is obtained is problematic because it is a relative error based on a known solution for this problem (analytical or numerical fine reference).
Another approach to deal with this boundary problem has been studied in more detail in [20] by examining how boundary conditions are imposed at the local scale: Dirichlet, Neumann and Cauchy boundary conditions have been compared, the last one providing a better convergence rate in the treated case (fracture mechanics problems). But in [27] it turned out that the Dirichlet boundary condition was the best in the given context.

The question of how the local problem solution is used to construct an enrichment function is another topic widely studied in the literature.
In [4] the SGFEM technique proposes to remove  from the fine scale field (obtained here with more conventional analytic functions) its piecewise linear interpolant on the coarse[In this article, the terms  "fine" or "micro" are sometimes used instead  of the original "local"   designation in the method.
And  "coarse" or "macro" are sometimes used instead of  "global".]-scale field.
This operation ensures a  matrix conditioning number for the enriched system equivalent to that of the non enriched system.
In [13, 12] the SGFEM technique, mentioned to potentially boost the method, was adapted to fracture mechanics with an  analytical enrichment  function improving both matrix conditioning and result quality.

The method has also  been successfully evaluated  in a nonlinear context [17, 29], in other mechanical domains [28, 11] and even integrated in a commercial software [22].
But since its inception until today, it is its intrinsic parallel capacities that have been put forward by its creators.
In [18], since  the local problems (one per patch) are independent of each other, they can be solved (with a direct solver) in parallel (in multithreaded  paradigm, each thread computes a set of patches).
In [27], the authors successfully reused this parallel implementation and addressed, among many aspects, the issue of memory consumption related to local problems.
Especially when the scale loop is used to cope with the boundary condition issue, storing the factorization of local problems reduces the cost of iterations, at the price of a significant increase in the memory footprint.
In [21] the local problems themselves have been treated with a parallel solver  keeping the idea of computing each of them in parallel.
This dual level of parallelism  was also studied at the same time in [31] (in a message passing context).

In [18, 27, 21] the local problems called "sub-local" domains were constructed from one or a few non-overlapping master local domains,  in order to apply the same discretization for all these "sub-local" domains.
This artifact greatly simplify the integration process of the enriched macro-scale problem since all its enrichment function share the same discretization  in this case, thus  the same integration points.

In this paper, the proposed method retains the idea of a loop between scales to deal with boundary condition problems, but adds the calculation of  a new error criterion  to terminate this loop at a given accuracy.
This criterion is based on the notion of reference solution which introduces a new way to formulate the enrichment function and to create the global problem.
Otherwise, to avoid any limitation imposed by computer memory, the proposed method  uses the parallel message passing interface (MPI) paradigm  to access distributed memory resources.
All the meshes are distributed over the processors, which induces an original scheduling algorithm, described in  section  <ref>, to optimize the load balancing of the processors at the fine scale.
The distributed nature of the mesh at the global scale  also leads to a parallel resolution at this scale.
Compared to the "sub-local" domain approach, the  proposed method  uses the more systematic approach of a local problem per patch while keeping the notion of a single discretization at the fine scale in order to simplify the parallel implementation, the definition of the enrichment functions and the construction of the problem at the global scale.
And finally, the chosen boundary conditions imposed on these local problems  are  of Dirichlet type.

The paper is organized as follows.
The next section recalls the classical  global-local formulation.
Then in section <ref>, we introduce a reference solution and reformulate the two-scale approach using a pure matrix format.
Section <ref>  studies the issue of scale difference between the global and local level and its impact on performance.
Section <ref> deals with the parallel global-local approach in which both the global
and local problems are distributed.
Section <ref> deals with the three numerical experiments and the paper ends with conclusion and possible future works.

§ INGREDIENTS

§.§ Continuous mechanical problem

In this work, we will only consider mechanical elasticity but other fields of application are possible.
In the  domain of interest $\Omega$,  a  region of $\mathbb{R}^3$, the strong form of the equilibrium  equation of the mechanical problem is defined as follows:

\begin{equation}
\nabla . \tens{\sigma} +\vm{f}= 0
\label{equilibrium}
\end{equation}

with  $\vm{f}$ being a prescribed volume loading and the stress tensor is  $\tens{\sigma}=\tens{C}:\tens{\epsilon}$, $\tens{C}$ being the Hook's tensor and $\tens{\epsilon}$ being the strain tensor.
The Neumann and Dirichlet conditions applied on the  $\partial\Omega$ boundary  of the $\Omega$ domain are:

\begin{equation}
\tens{\sigma}.\vm{n} = \vm{t}~ \text{on}~\partial\Omega^N, \vm{u}=\overline{\vm{u}} ~\text{on} ~\partial\Omega^D, \partial\Omega^N\cap \partial\Omega^D=\oslash
\label{BC}
\end{equation}

where $\vm{t}$ is the prescribed tensile load on the   $\partial\Omega^N$ part of $\partial\Omega$, $\vm{n}$ the outgoing normal vector of $\Omega$, $\overline{\vm{u}}$ the prescribed displacements on the  $\partial\Omega^D$ part of $\partial\Omega$ and $\vm{u}$ the displacement field.
Let $\mathcal{M}^{\mathsf{C}}$  be the continuous space of the problem defined on $\Omega$ and compatible with the Dirichlet boundary conditions ($\mathcal{M}_0^{\mathsf{C}}$ being $\mathcal{M}^{\mathsf{C}}$ with null Dirichlet boundary conditions).
The solution $\vm{u}^C \in \mathcal{M}^{\mathsf{C}}$ of the continuous problem defined by (<ref>) and (<ref>) in their weak form, is given by solving:

\begin{equation}
\forall \vm{v}^* \in \mathcal{M}_0^{\mathsf{C}},~~ A\left( \vm{u}^C,\vm{v}^* \right)_{\Omega} = B\left(\vm{v}^*\right)_{\Omega},~~\vm{u}^C=\overline{\vm{u}} ~\text{on} ~\partial\Omega^D
\label{equilib}
\end{equation}

where  the bi-linear and linear  forms are:

\begin{equation}
\begin{array}{l}
A\left( \vm{u},\vm{v} \right)_{\Omega} =\stretchint{5ex}_{\Omega}
\tens{\epsilon} \left( \vm{u} \right):\tens{C}:\tens{\epsilon}\left( \vm{v}\right)\mathrm{d}V \\
B\left( \vm{v} \right)_{\Omega} =\stretchint{5ex}_{\partial\Omega^N}
\vm{t}\cdot \vm{v} \mathrm{d}S-\stretchint{5ex}_{\Omega}
\vm{f}\cdot \vm{v} \mathrm{d}V
\end{array}
\text{with}~ \vm{u}\in \mathcal{M}^{\mathsf{C}} ~\text{and} ~\vm{v}\in \mathcal{M}_0^{\mathsf{C}}
\label{bilinform}
\end{equation}

and the strain is given by the  kinematic equation
$\tens{\epsilon}\left( \vm{u} \right)=\frac{1}{2}\left(\nabla\vm{u} +\left(\nabla \vm{u}\right)^t\right)$.

§.§ Global problem

To solve (<ref>),  the approach  discretizes the continuous problem into two steps.
It starts by discretizing $\Omega$ with a mesh composed of  macro (or global or coarse) elements capable of representing the general behavior of the problem.
On this global mesh,  where $I^g$ denotes its set of  nodes (the superscript g stands for "global scale"), we consider the standard first-order finite element approximation functions $N^i:\vm{x}\in \mathbb{R}^3\rightarrow N^i(\vm{x}) \in \left[0,1\right]\subset\mathbb{R}$, $i\in I^g$, associated to a coarse scale node $\vm{x}^i$, with a support  given by the union of all  finite elements sharing the node $\vm{x}^i$ ( this group of macro-elements is called "patch" hereafter).
These $N^i$ are used to linearly interpolate the degrees of freedom (dof), also called values in this work, which represent the classical values of the discretized displacement field at the  mesh nodes.

Some nodes of this global mesh are then enriched when their support covers certain local behaviors that must be finely described.
Let $I_{e}^g\subseteq I^g$ be the set of enriched global nodes ($card\left( I_{e}^g\right)\leqslant card\left( I^g\right)$ and the subscript e stands for "enriched").
For each node $p\in I_{e}^g$  a enrichment function $\vm{F}^p(\vm{x})$ (described in the next section) enriches the classical basis  ($N^i$) of the discretization space.
The  kinematic equation used to describe the discrete global field $\vm{U}$  at a point $\vm{x}$ is then:

\begin{equation}
\vm{U}(\vm{x})=\displaystyle\sum_{i\in I^g}  N^i(\vm{x})~ \vm{U}^{i} + \displaystyle\sum_{p\in I_{e}^g}   N^p(\vm{x})~\vm{E}^{p}*\vm{F}^p(\vm{x})
\label{kinematic}
\end{equation}

* $\vm{U}^{i}$ is the vector of classical values  for the node $\vm{x}^i$

* $\vm{E}^{i}$ is the vector of enriched values  for the node $\vm{x}^i$

* $*$ operator is the component  by component  multiplication of two vectors

Thus (<ref>) becomes to find $\vm{U}^G\in \mathcal{M}^{\mathsf{G}}\subset \mathcal{M}^{\mathsf{C}}$ such that:

\begin{equation}
\forall \vm{v}^* \in \mathcal{M}_0^{\mathsf{G}},~~ A\left( \vm{U}^G,\vm{v}^* \right)_{\Omega} = B\left(\vm{v}^*\right)_{\Omega},~~\vm{U}^G=\overline{\vm{u}} ~\text{on} ~\partial\Omega^D
\label{equilibg}
\end{equation}

with $
\mathcal{M}^{\mathsf{G}}=\left\{\vm{U}(\vm{x}):\vm{U}(\vm{x})=\displaystyle\sum_{i\in I^g}  N^i(\vm{x})~ \vm{U}^{i} + \displaystyle\sum_{p\in I_{e}^g}   N^p(\vm{x})~\vm{E}^{p}*\vm{F}^p(\vm{x}), \vm{U}(\vm{x})=\overline{\vm{u}}(\vm{x}) ~\text{for}~ \vm{x}\in \partial\Omega^D  \right\}$.
The integration of (<ref>) using appropriate Gauss  quadrature points (see [18] for example) leads to the resolution of a linear system computed with a direct solver.
This system is constructed  by eliminating  the rows of the elementary matrices related to the Dirichlet values and  moving the coupling term (columns related to the Dirichlet values) to the right side.

§.§ Local problems and enrichment functions

For a node $p\in I_{e}^g$, the enrichment function $\vm{F}^p(\vm{x})$ is constructed numerically using the solution of a micro (or local or fine) problem.
This local problem  follows the same equations as (<ref>) and (<ref>), and is discretized by a micro-mesh corresponding to the refinement (in a nested manner) of the macro elements of the  $p$ patch.
This is the second discretization step of the method.
Considering $\omega_p$ the region corresponding to the patch $p$, the solution of the local problem $p$ is obtained by finding $\vm{u}^{Q^p}\in \mathcal{M}^{{\mathsf{Q}}^p}\subset \mathcal{M}^{\mathsf{C}}$ such that:

\begin{equation}
\forall \vm{v}^* \in \mathcal{M}_0^{{\mathsf{Q}}^p},~~ A\left( \vm{u}^{Q^p},\vm{v}^* \right)_{\omega_p} = B\left(\vm{v}^*\right)_{\omega_p},
%~~\vm{u}^{Q^p}=\overline{\vm{u}} ~\text{on} ~\partial\omega_p\cap \partial\Omega^D,
~~	\vm{u}^{Q^p}=\vm{U}^G~ \text{on}~\partial\omega_p%\setminus \partial\Omega^D
\label{equilibl}
\end{equation}

* $\vm{U}^G$ is the solution of (<ref>)

* $\mathcal{M}^{{\mathsf{Q}}^p}=\left\{\vm{u}(\vm{x}):\vm{u}(\vm{x}) =\displaystyle\sum_{m\in I_p^l} n^m(\vm{x})~\vm{u}^{m^p},
%~ \vm{u}(\vm{x})=\overline{\vm{u}}(\vm{x}) ~\text{when}~ \vm{x}\in \partial\omega_p\cap \partial\Omega^D ,
~  \vm{u}(\vm{x})=\vm{U}^G(\vm{x})~ \text{for}~\vm{x}\in \partial\omega_p%\setminus \partial\Omega^D
\right\}$

* $I_p^l$ is the set of nodes associated with the fine mesh discretization of $\omega_p$ (superscript l stands for "local scale")

* $n^m:\vm{x}\in \mathbb{R}^3\rightarrow n^m(\vm{x}) \in \left[0,1\right]\subset\mathbb{R}$ is the first order finite element shape function associated to the fine scale node $\vm{x}^m$

* $\vm{u}^{m^p}$ is the vector of classical dofs  for the node $\vm{x}^m$ corresponding to the fine scale problem $p$

* $A\left( \vm{u},\vm{v} \right)_{\omega_p} =\stretchint{5ex}_{\omega_p}
\tens{\epsilon} \left( \vm{u} \right):\tens{C}:\tens{\epsilon}\left( \vm{v}\right)\mathrm{d}V$ with $\vm{u}\in \mathcal{M}^{{\mathsf{Q}}^p}$ and $\vm{v}\in \mathcal{M}_0^{{\mathsf{Q}}^p}$

* $B\left( \vm{v} \right)_{\omega_p} =\stretchint{5ex}_{\partial\omega_p\cap\partial\Omega^N}
\vm{t}\cdot \vm{v} \mathrm{d}S-\stretchint{5ex}_{\omega_p}
\vm{f}\cdot \vm{v} \mathrm{d}V$ with $\vm{v}\in \mathcal{M}_0^{{\mathsf{Q}}^p}$

Regarding the boundary condition of the fine-scale problems, different approaches have been tested (in particular in [20]), all of them using the global-scale solution for  Dirichlet, Neumann or  a mix of both.
The idea is to impose the global behavior of the problem on the boundary of the fine-scale problems so that the local behaviors that must be finely described can be revealed by the solutions of the fine-scale problems.
In this work, the global-scale solution is only imposed by the Dirichlet boundary condition as presented in (<ref>).
Note that the boundary conditions at the global scale are inherited at the fine scale. In particular the coarse-scale Dirichlet boundary condition  are imposed via $\vm{U}^G$ (since $\vm{U}^G$ is the solution of (<ref>), $\vm{U}^G=\overline{\vm{u}} ~\text{on} ~\partial\Omega^D$ and  $\vm{u}^{Q^p}=\vm{U}^G=\overline{\vm{u}} ~\text{on} ~\partial\omega_p\cap \partial\Omega^D$).

The integration of (<ref>) for all $p$, using the standard Gauss  quadrature points, lead to the construction of $card\left( I_{e}^g\right)$ linear systems (reduced at assembly by elimination of Dirichlet values).
Their resolution can be done in different ways, but,  as already indicated in the introduction, in [18] they are solved in parallel using threads.
The distribution of patches to threads uses dynamic scheduling: each thread takes a local problem from a list and processes it (with a direct solver) until the list is empty.
Kim et al. proved that local problems should be sorted in this list by decreasing cost (calculated  a priori with some metric) in order to maintain a good load balance.
They also showed that the most expensive local problem  limits the number of processes that can be used for a given problem.
Beyond this limit, parallel  scaling efficiency drops.
Section <ref> will explain how we go from a multi-thread patch resolution to a full MPI implementation.

When $\vm{u}^{Q^p}$ are available, the enrichment functions $\vm{F}^p$ can be constructed from them.
A first crud approach is to consider that:

\begin{equation}
\forall p\in I_{e}^g,~\forall \vm{x}\in \omega_p ~\vm{F}^p(\vm{x})=\vm{u}^{Q^p}(\vm{x})~ \text{and}~ \forall \vm{x}\notin \omega_p~\vm{F}^p(\vm{x})=\vm{0}
\label{crudenrich}
\end{equation}

This solution can works if the terms of $N^p(\vm{x})~\vm{F}^p(\vm{x})$  are not to close to $N^p(\vm{x})$.
Otherwise,  the created system is badly conditioned because the discretization basis is almost redundant.
But in the literature,  many authors use the SGFEM technique [4]  which removes  from $\vm{u}^{Q^p}$ its piecewise linear interpolant over the coarse-scale field:

\begin{equation}
\forall p\in I_{e}^g,~ \forall \vm{x}\in \omega_p ~\vm{F}^p(\vm{x})=\vm{u}^{Q^p}(\vm{x})-\displaystyle\sum_{j\in I_{p}^{g}}  N^j(\vm{x})\vm{u}^{Q^p}(\vm{x}^j)~ \text{and}~ \forall \vm{x}\notin \omega_p~\vm{F}^p(\vm{x})=\vm{0}
\label{sgfem}
\end{equation}

where $I_{p}^g$ is the set of macro-scale nodes associated with  patch $p$.
Eliminating the projection first removes from the enrichment function the part that can be represented by the coarse-scale discretization and thus improves the matrix conditioning of the coarse-scale problem.
Second, it forces the enrichment function to be zero at coarse-scale  nodes location, which helps solve the blending [Blending appears in a coarse level element when  all  its classical dofs are not enriched.] element problem and also greatly simplifies[With a non-zero enrichment function, imposing a Dirichlet condition on an enriched node adds the complexity of setting up an equation linking the classical and enriched dofs.] the application of the global Dirichlet boundary conditions, if any, on those enriched nodes.
In this work, a rather different approach is proposed, as presented in the section <ref>.

§.§ The scale loop

The main idea of the solver is to  iteratively solve the two types of problems:
the global problem (<ref>) and
the set of micro  problems (<ref>) over patches solved independently.
The global and fine problems are linked together.
The fine problems over the patches deliver the enrichment functions which are used in the global problem.
The global problem delivers the boundary conditions for the patches.

The initial step of the loop is a key element of the method.
In general, a global-scale solution is used as the starting point of the loop (it provides a boundary condition for fine scale problems).
This solution can be the result of a global-scale problem without enrichment or the  result of a previous  computation.
In particular,  when evolutionary phenomena are involved, the solver can effectively take advantage of using the last solution found (from the previous evolutionary step)  as the starting solution of the  loop used for the current evolution step.
Some test cases in the section <ref> illustrate this point.

In any case, all the fine-scale problems are mainly handled during this initialization phase.
The integration of (<ref>)  yields a set of linear system which  have  constant matrices  and  partially constant right-hand side vectors during iterations.
Only the Dirichlet coupling term changes right-hand side vector during the loop as $\vm{U}^G$ varies.
Thus, the  linear systems matrices and constant right-hand side vectors are  created only once during the initialization.
And by choosing a direct solver to solve systems, we can also compute the matrix factorization only once.
If preserved (which has an impact on memory, as shown in [27]), these factorizations can be reused for backward/forward resolution of fine-scale problems in the remaining steps of the loop.

§ THE PROPOSED METHOD: MATRIX-BASED RESOLUTION WITH A REFERENCE TARGET

§.§ Overview

Based on the ingredients of the section <ref>, the figure <ref> illustrates in 2D, but the same principle applies in 3D, the  general concept of the proposed method with a fictitious problem on a $\Omega$ region.

[][Global-scale problem]

[][The 25 enriched nodes with 4 of their supports colored]

[][: support union  ]

[][Fine scale discretisation]

[][Reference field discretization]

[][Fine scale problems (2 out of 25)]

(0.57702894,0.21777677)(0,0)[lt]0t]lDirichlet ($\vm{u}_d^p$)
(0.27095745,0.145)(0,0)[lt]0t]lVisual effect to show the location of the local
(0.27095745,0.115)(0,0)[lt]0t]l behavior to be  captured
(0.03399547,0.30154357)(0,0)[lt]0t]lBoundary conditions:
(0.27095745,0.0574839)(0,0)[lt]0t]lLinear relation
(0.27479608,0.00599376)(0,0)[lt]0t]lHanging nodes

(0.08816039,0.26898148)(0,0)[lt]0t]lGlobal scale enriched node ($\in I_{e}^g$) with a partially refined patch
(0.08816039,0.33843793)(0,0)[lt]0t]lGlobal scale enriched node ($\in I_{e}^g$) with a fully refined patch
(0.17760806,0.15915917)(0,0)[lt]0t]lGlobal-scale mesh
(0.17760806,0.04575297)(0,0)[lt]0t]lReference mesh
(0.72510056,0.15262199)(0,0)[lt]0t]l: refined
(0.72510056,0.02882823)(0,0)[lt]0t]l: not refined

method presented in a fictitious 2D problem (for sake of visibility fine-scale mesh is not refined as it should)

Figure <ref> shows the discretization of the global scale problem  with its boundary conditions.
The magenta area is a simple visual effect to show the position of the local behavior to be captured.
This area is just a visual artifact where the mesh needs to be refined  to appropriately simulate a crack, a specific material inclusion, a damaged material, a rivet, ...
All macro nodes for which a macro element of their support has been refined  are enriched (figure <ref> where red and green dots form the $I_e^g$ set introduced in the section <ref>).
Their enrichment functions are obtained from the solutions of the fine-scale problems  associated with their patches (figure <ref> where global and local boundary condition  are imposed, $\vm{u}_d^p$ being  introduced in the section <ref> and defined precisely  in <ref>).
The union of all the elements of the enriched patches (yellow + orange in the figure <ref>) is called the , noted SP thereafter.
Conversely, the set of macro-element that are not part of the SP (shaded in the figure <ref>)  will be noted as NSP in what follows.
One of the goals of the is to impose a single fine mesh discretization for all fine-scale problems.
More precisely, the idea is to impose a fine discretization for all patches covering the local behavior to be captured.
Thus their union (in orange on the figure <ref>) is the starting mesh of the adaptation strategy described in  <ref> (see this appendix for details on how the elements are cut/split/refined) which guarantees that all the elements of this area are at least split once and that the local behavior is well refined (figure <ref>).
Using only  orange elements prevents the refinement of any other element, in order to respect the "one hanging node per edge" rule (avoiding the enrichment of more nodes).
The yellow part of the SP is then the set of macro elements not reﬁned and therefore used without modification  at the fine-scale level.
Thus, in the SP, at the boundary between two macro-elements, the micro-meshes are either compatible or leave hanging nodes (visible in figure <ref> and presented in section <ref>).
This unified discretization imposed in the SP, allows to process  all quantities related to fine-scale grid  (integration at Gauss point, level set, damage, mapping, ... ) only once.
In this work, the union of the SP refined mesh and the NSP macro-elements (which naturally connect  to the yellow fine scale grid elements) is called the reference mesh (figure <ref>).
It represented a new vision of what the method intends to discretize at the local and global level alternatively.
The section  <ref> will detail the problem associated with this mesh.

A rich naming convention is implemented in <ref> to clarify  the different discretizations and the status of the values involved in the different problems solved by the solver.
Thus, all subscribe letters appearing in the matrix or vector in the following algorithms or discussions come from the table <ref>.

The general scheme  of this new version is given by the algorithm <ref>.

Create $\vm{U}_g^0$ macro-scale dofs

$\left(\vm{A}_{gg}^{ini},\vm{B}_{g}^{ini},\vm{A}_{FF},\vm{B}_F,\left\|\vm{B}_r\right\|,\vm{A}_{qq}^{patches},\vm{BI}_q^{patches},\vm{D}_{qd}^{patches} \right) \gets$ INIT
Retrieve macro-scale classical dofs to initialize $\vm{U}_g$ (enriched dof set to zero)
$ \vm{u}_F \gets$ UPDATE_MICRO_DOFS$\vm{U}_{g}$

$ \vm{u}_q^{patches} \gets$MICRO-SCALE_RESOLUTION$\vm{u}_F$,  $\vm{A}_{qq}^{patches}$, $\vm{BI}_q^{patches}$, $\vm{D}_{qd}^{patches}$
$\left( \vm{A}_{gg},\vm{B}_{g}\right)\gets$UPDATE_MACRO_PRB$\vm{A}_{gg}^{ini}$, $\vm{B}_{g}^{ini}$, $\vm{u}_q^{patches}$
$\vm{U}_{g}\gets \vm{A}_{gg}^{-1}\cdot \vm{B}_{g}$  $\dagger$
$ \vm{u}_F \gets$UPDATE_MICRO_DOFS$\vm{U}_{g}$
$resi\gets$COMPUTE_RESIDUAL$\vm{u}_F$, $\vm{A}_{FF}$, $\vm{B}_F$, $\left\|\vm{B}_r\right\|$
$resi< \epsilon$

algorithm: general procedure.   Procedures INIT, UPDATE_MACRO_PRB, UPDATE_MICRO_DOFS, MICRO-SCALE_RESOLUTION and COMPUTE_RESIDUAL are  respectively depicted in algorithm <ref>, <ref> , <ref>, <ref> and <ref> of <ref>.
The $\epsilon$ value is the user-defined  target accuracy that the criterion given in the equation (<ref>) must meet to end the scaling loop. See <ref> for notation convention : subscripted letters,superscript  $patches$,$\triangleright$, ...

After an initialization step,
the alternate resolution between scales starts with the MICRO-SCALE_RESOLUTION procedure (algorithm <ref>) which computes the solutions of fine-scale problems.
These solutions are used by the  UPDATE_MACRO_PRB procedure (algorithm <ref>)  to update the linear system at the global-scale ($\vm{A}_{gg}$ and $\vm{B}_{g}$).
Then, a conventional  solver ($\dagger$ in the algorithm <ref> described in section <ref>)  solves the problem at the global-scale.
The global scale solution ($\vm{U}_{g}$) is then transferred to the fine scale by the  UPDATE_MICRO_DOFS procedure (algorithm <ref>).
Compare to the classical version, the scale loop is now controlled by the computation of a $resi$ criterion (COMPUTE_RESIDUAL procedure described by the algorithm <ref>) which stops the iterations when it is lower than a threshold ( $\epsilon$) given by the user.
This  convergence check is an important contribution  because it allows to evaluate the solver and to easily compare  different version of it.
This criterion is introduced in section <ref> and its computation is given in section <ref>.
The following subsections present in more  detail the  algorithm <ref> and the associated procedures without focusing too much on the parallelism  which will be discussed in section <ref> (only the presence of the symbol $\triangleright$, on the right of a line of the algorithm, indicates that a certain communication between the processes exists for this step.).

§.§ The reference field and its associated criterion

The reference field $\vm{u}^{R}$ is  the  solution of the  problem defined by  (<ref>) and (<ref>), and  discretized by the reference mesh (figure <ref>).
Solving this reference problem  is to find $\vm{u}^{R}\in \mathcal{M}^{\mathsf{R}}\subset \mathcal{M}^{\mathsf{C}}$ such that:

\begin{equation}
\forall \vm{v}^* \in \mathcal{M}_0^{\mathsf{R}},~~ A\left( \vm{u}^{R},\vm{v}^* \right)_{\Omega} = B\left(\vm{v}^*\right)_{\Omega},
~~\vm{u}^{R}=\overline{\vm{u}} ~\text{on} ~\partial\Omega^D,
\label{equilibR}
\end{equation}

* $\mathcal{M}^{\mathsf{R}}=\left\{\vm{u}(\vm{x}):\vm{u}(\vm{x}) =\displaystyle\sum_{k\in I^l} \mathtt{N}^k(\vm{x})~\vm{u}^{k},
~  \vm{u}(\vm{x})=\overline{\vm{u}}(\vm{x})~ \text{for}~\vm{x}\in  \partial\Omega^D
\right\}$

* $I^l$ is the set of nodes associated with the reference mesh, discretization of $\Omega$

* $\mathtt{N}^k$ is either $n^k$ or $N^k$  corresponding respectively to the shape functions defined in <ref> or in <ref> depending on whether $\vm{x}$ is in a micro or macro element.

* $\vm{u}^{k}$ is the vector of classical dofs  for the node $\vm{x}^k$

This solution $\vm{u}^R$  represents the continuous field solution $\vm{u}^C$  only polluted by the discretization error.
The standard Gauss quadrature  integration of (<ref>) gives the following linear system without Dirichlet elimination:

\begin{equation}
\vm{A}_{RR}\cdot \vm{u}_R^R=\vm{B}_R%^{\eqref{equilibR}}
~\text{with $\vm{u}_R^R$ the vector of $card(R)$ values defining $\vm{u}^R$}
\label{R_sys}
\end{equation}

As already mentioned in  the section <ref>, in the reference mesh there are hanging nodes between the refined and unrefined areas (visible in figure <ref>) that need to  be fixed to avoid  discontinuity of the $\vm{u}^{R}$ field.
A simple approach is to eliminate them from the system (<ref>) by using a linear relationship to fix their value using the following formula:

\begin{equation}
\text{For an hanging node}~\vm{x}^j,~ \vm{u}^{j} =\displaystyle\sum_{i\in I_j}N^i(\vm{x}^j)~ \vm{u}^{i}
\label{hanging_eq}
\end{equation}

* $\vm{u}^{j}$ is the vector of fixed values (L-set) for the node $\vm{x}^j$

* $\vm{u}^{i}$ is the vector of free values (f-set) for the node $\vm{x}^i$

* $I_j$ is the set of nodes corresponding to the vertices of the edge or face on which hanging node $\vm{x}^j$ is located. Thus $card\left( I_j \right) \in \left\lbrace 2,3\right\rbrace $.

* $N^i(\vm{x})$ is the standard first-order finite element approximation function associated with the  $\vm{x}^i$ node, related to the coarse element support that holds the edge or face where hanging node is present.

Eliminating the linear relations (<ref>) and the Dirichlet boundary condition of (<ref>) leads to the system:

\begin{equation}
\vm{A}_{rr}\cdot \vm{u}_r^R=\vm{B}_r
\label{r_sys}
\end{equation}

The solution given by a direct solver of this system is  $\vm{u}_r^R$ the vector of $card(r)$ dofs defining $\vm{u}^R$.
But in this work, this system is never solved directly.
Instead, the proposed solver, at each scale iteration, computes a  solution $\vm{u}^{ts}$ ($\vm{u}^{ts} \in \mathcal{M}^{\mathsf{R}}$) such that it quickly converges  to the  reference solution $\vm{u}^R$.

In this work, a $resi$ criterion is proposed to control this convergence.
The $resi$ criterion is taken as the relative residual error of  system (<ref>):

\begin{equation}
resi=\frac{\left \| \vm{A}_{rr}.\vm{u}_{r}^{ts}-\vm{B}_{r} \right \|}{\left \|\vm{B}_{r} \right \|}
\label{residual_error}
\end{equation}

where $\left \| . \right \|$ is the L2-norm and $\vm{u}_r^{ts}$ is the vector of $card(r)$ dofs   corresponding to the $\vm{u}^{ts}$ solution obtained by the solver.
For simplicity, in what follows, $\vm{u}_r$ will represent the $\vm{u}_r^{ts}$ vector of the field $\vm{u}^{ts}$ where  superscripts $ts$ is made implicit.
The question of computing (<ref>) is exposed in the following section.

To  relate this  criterion to an energy error, an energy norm is defined  as follows using (<ref>):

\begin{equation}
\left\| \vm{u} \right\|_{E_{\Omega}} =\sqrt{A\left( \vm{u},\vm{u} \right)_{\Omega}}=\sqrt{\stretchint{5ex}_{\Omega}
\tens{\epsilon} \left( \vm{u} \right):\tens{C}:\tens{\epsilon}\left( \vm{u}\right)\mathrm{d}\Omega}
\label{energynorm}
\end{equation}

The error introduced by the solver (mainly on the boundary conditions imposed on the fine scale problems) in the indirect resolution of system <ref> compare to its direct resolution can thus be expressed by $\left\| \vm{u}^{ts}-\vm{u}^R \right\|_{E_{\Omega}}$.
It is related to $resi$ by $\vm{A}_{rr}\cdot \vm{u}_{r}-\vm{B}_{r}$ vector:

\begin{equation}
\left\| \vm{u}^{ts}-\vm{u}^R \right\|_{E_{\Omega}}=\sqrt{A\left( \vm{u}^{ts}-\vm{u}^R, \vm{u}^{ts}-\vm{u}^R\right)}=\sqrt{\left( \vm{A}_{rr}\left( \vm{u}_{r}-\vm{u}_r^R\right)\right) \cdot \left( \vm{u}_{r}-\vm{u}_r^R\right)}=\sqrt{\left(\vm{A}_{rr}.\vm{u}_{r}-\vm{B}_{r}\right)\cdot \left( \vm{u}_{r}-\vm{u}_r^R\right)}
\label{lien_resi_energ}
\end{equation}

So when $\vm{A}_{rr}\cdot \vm{u}_{r}-\vm{B}_{r}$ vector tends to the zero vector, both $resi$ and $\left\| \vm{u}^{ts}-\vm{u}^R \right\|_{E_{\Omega}}$ tend to zero.

It is also interesting to  compare the solution to the continuous solution.
The associated energy error is  $\left\| \vm{u}^{ts}-\vm{u}^C \right\|_{E_{\Omega}}$  and can be expressed as:

\begin{equation}
\begin{array}{r}
\left\| \vm{u}^{ts}-\vm{u}^C \right\|_{E_{\Omega}}^2 =\stretchint{5ex}_{\Omega}
\tens{\epsilon} \left( \vm{u}^{ts} -\vm{u}^R + \vm{u}^R-\vm{u}^C \right):\tens{C}:\tens{\epsilon}\left( \vm{u}^{ts}-\vm{u}^R + \vm{u}^R-\vm{u}^C\right)\mathrm{d}\Omega =\\
\left\| \vm{u}^{ts}-\vm{u}^R \right\|_{E_{\Omega}}^2+\left\| \vm{u}^R-\vm{u}^C \right\|_{E_{\Omega}}^2+2A\left( \vm{u}^R-\vm{u}^C,\vm{u}^{ts}-\vm{u}^R \right)_{\Omega} 	\end{array}
\label{scale error}
\end{equation}

Subtracting (<ref>) to (<ref>) gives as $\mathcal{M}^{\mathsf{R}}\subset \mathcal{M}^{\mathsf{C}}$:

\begin{equation}
\forall \vm{v}^* \in \mathcal{M}_0^{\mathsf{R}}, A\left( \vm{u}^R-\vm{u}^C,\vm{v}^* \right)_{\Omega} =B\left(\vm{v}^*\right)_{\Omega}-B\left(\vm{v}^*\right)_{\Omega}=0
\label{equilibA-R}
\end{equation}

Based on (<ref>) and as $\vm{u}^{ts}-\vm{u}^R \in \mathcal{M}_0^{\mathsf{R}}$, then the extra term  of (<ref>)  vanishes:

\begin{equation}
A\left( \vm{u}^R-\vm{u}^C,\vm{u}^{ts}-\vm{u}^R \right)_{\Omega} =0
\label{Nullterm}
\end{equation}

Leading to the relation:

\begin{equation}
\left\| \vm{u}^{ts}-\vm{u}^C \right\|_{E_{\Omega}}^2 =\left\| \vm{u}^{ts}-\vm{u}^R \right\|_{E_{\Omega}}^2+\left\| \vm{u}^R-\vm{u}^C \right\|_{E_{\Omega}}^2
\label{equality}
\end{equation}

The second term on the right-hand side of (<ref>) represents the error introduced by the reference mesh discretization with respect to the continuous solution.
The first one is the error introduced by the method.
Most publications on the subject do compute $\vm{u}^R$ with another solver and verify a posteriori that $\left\| \vm{u}^{ts}-\vm{u}^C \right\|_{E_{\Omega}}$ tends to $\left\| \vm{u}^R-\vm{u}^C \right\|_{E_{\Omega}}$ (i.e. that $\left\| \vm{u}^{ts}-\vm{u}^R \right\|_{E_{\Omega}}$ tends to zero) in few iterations.
In this work $resi$ carries out this verification during the iterations.
A numerical comparison of these errors is given in section <ref>.

§.§ Specific global-scale system construction: the algebraic operator

The $\vm{u}_{r}$  (as mentioned above, the simplified notation of $\vm{u}_r^{ts}$) vector will act as a bridge between computational tasks.
It stores  classical values at the nodes of $\vm{u}^{ts}$ field which will be used to compute $resi$ and impose boundary conditions on fine-scale problems.
The $\vm{u}_{r}$ dofs of this field can be related to global-scale field dofs via the kinematic equation (<ref>).
With this equation, for any node $\vm{x}^j$ of the reference mesh, it is possible to write, without any approximation, that the solution  $\vm{u}^{ts}$ is equal to the  global enriched discrete solution at this location:

\begin{equation}
\vm{u}^{ts}(\vm{x}^j)=\vm{U}^G(\vm{x}^j)=\displaystyle\sum_{i\in I^g}  N^i(\vm{x}^j)~ \vm{U}^{i} + \displaystyle\sum_{p\in I_{e}^g}   N^p(\vm{x}^j)~\vm{E}^{p}*\vm{F}^p(\vm{x}^j)
\label{relation}
\end{equation}

This can be rewritten in algebraic form as:

\begin{equation}
\vm{u}_M=\vm{T}_{MG}\cdot \vm{U}_G
\label{u_ru_g}
\end{equation}

where $\vm{T}_{MG}$ is a linear  interpolation operator whose structure and construction are given in <ref>.
In this work, the process of eliminating linear relations (<ref>) is made transparent by choosing,  from the R-set, to  first eliminate the L-set before applying  linear interpolation.
Thus, the operator in (<ref>) interpolates from the G-set to the M-set.

Considering the matrices of (<ref>) with the  L-set eliminated, and the equation (<ref>) to apply a change of variable,  the equality between the potential energy and the work of  external forces can be written as follows:

\begin{equation}
\vm{U}_G^t\cdot \vm{T}_{MG}^t\cdot  \vm{A}_{MM} \cdot \vm{T}_{MG}\cdot \vm{U}_G=\vm{U}_G^t\cdot \vm{T}_{MG}^t\cdot \vm{B}_{M}
\label{energy_g}
\end{equation}

This gives the following system (see <ref> for the detailed construction):

\begin{equation}
\vm{A}_{gg}\cdot \vm{U}_g=\vm{B}_g
\label{g_sys}
\end{equation}

From a practical point of view, the sub-blocks of $\vm{T}_{MG}$, $\vm{A}_{MM}$ and  $\vm{B}_{M}$ are large sparse matrices, so they are never fully assembled in memory but stored in blocks using a per macro-element strategy.
This strategy, which gathers in memory all the matrices and vectors necessary for the calculation relative to a macro element, is favorable to the reduction of cache misses and to multi-threading computation.
Moreover, during the loop, only the terms  related to  $\vm{F}^p$ (related to $\vm{T}_{Fe}$ operator given in <ref> ) have to be updated at each scale iteration. Thus some sub-block of $\vm{A}_{gg}$ and $\vm{B}_g$  are constant during the iteration and are computed and assembled only once at initialization time.

The INIT procedure (algorithm <ref> detailed in <ref> ) is in charge of this initialization.
Let $\omega^{e_{macro}}$ be the region of $\Omega$ covered by a macro element $e_{macro}$.
$\omega^{e_{macro}}$ is discretized by one or more  micro-elements of the reference mesh.
With two loops on all macro elements, the INIT procedure integrates and assembles (<ref>)  over $\omega^{e_{macro}}$ with the elimination of (<ref>) but keeping Dirichlet terms.
The resulting matrices and vectors  $\vm{A}_{FF}^{e_{macro}}$ and $\vm{B}_{F}^{e_{macro}}$  are stored for all SP elements.
From this first matrix assembly of $\vm{A}_{FF}$ and $\vm{B}_{F}$ per block and the  creation  of the  sub blocks of  $\vm{T}_{MG}$, all the constant terms of $\vm{A}_{gg}$ and $\vm{B}_g$ can be computed algebraically  and assembled in two dedicated memory areas, $\vm{A}_{gg}^{ini}$ and $\vm{B}_{g}^{ini}$.

§.§ global to local transfer

In this work the idea is to obtain the $\vm{u}^{ts}$ field not by solving directly (<ref>)  but by solving (<ref>) and use (<ref>) to go from $\vm{U}^G$ to $\vm{u}_F$.
This last point is achieved by the UPDATE_MICRO_DOFS procedure (algorithm <ref> detailed in <ref>) which uses the same  per-macro-element storage strategy and stores the $\vm{u}^{ts}$ field in the $\vm{u}_F^{e_{macro}}$ vectors attached to SP macro-elements.
Note that this procedure does not loop over the NSP elements and that no information about the h-set is stored in memory because in the computational tasks related to $\vm{u}^{ts}$, this information can be ignored as we will see later.

§.§ Local problem

In this work, as in the section <ref>, the solution of the local problem  for a patch $p$ will be given by solving (<ref>) but using a different mechanism to create the linear system.
Consider that $\omega_p$ is equal to $\bigcup\limits_{e_{macro}\in J^p}\omega_{e_{macro}}$, with $J^p$ the set of macro-elements of the patch $p$.
As  the local problems are parts of the reference problem sharing with it its definition  and discretization, the integration of (<ref>) has already been partially done  when integrating (<ref>)  over $\omega^{e_{macro}}$.
Thus, creating the matrix and vector for a patch $p$ is just  assembling  the blocks $\vm{A}_{FF}^{e_{macro}}$ and $\vm{B}_F^{e_{macro}}$ with $e_{macro}\in J^p$.

As for the Dirichlet boundary conditions on $\partial\omega_p$, they are directly obtained from the $\vm{u}^{ts}$ field ($\vm{u}_F^{e_{macro}}$ dofs).
It is strictly the same as using $\vm{U}^G$ because $\vm{u}_F$ which comes from (<ref>) is computed with $\vm{T}_{MG}$ constructed from (<ref>) as the space $\mathcal{M}^{\mathsf{G}}$ constructed from same kinematic equation.
The application of these boundary conditions gives the following set of local linear systems:

\begin{equation}
\forall p\in I_e^g~\vm{A}_{qq}^p\cdot \vm{u}_q^p=\vm{B}_q^p
\label{p_sys}
\end{equation}

Again, it is the INIT procedure, with a loop over the patches, that is responsible for  creating  $\vm{A}_{qq}^p$ which is constant during the iteration.
Similarly, the constant part of the vector $\vm{B}_q^p$, called  $\vm{BI}_q^p$ and the constant coupling Dirichlet term $\vm{D}_{qd}^p$ are computed in this loop over the patches (see <ref> for more details).
Note that the use of Dirichlet boundary conditions for local problems  simplifies the software implementation and offers good performance for solving  (<ref>): smaller system size, simple constant algebraic operator ($\vm{D}_{qd}^p$) for imposing the prescribed displacements, and the matrix is not polluted by additional parameters (springs) to be set.
But the consequence of this choice on the proposed loop has not been studied and should be in a future work.

§.§ The loop

When the INIT procedure is completed, all  local structures associated with  macro-elements or patches are allocated and partially computed.
Before entering the loop,  the algorithm <ref>  initializes the field $\vm{u}^{ts}$  by calling  UPDATE_MICRO_DOFS with   the vector $\vm{U}_g$, solution of a specific macro-scale problem  (in general  the unenriched macro-scale problem but it can come from another computation).

The algorithm <ref>  then enters  the loop and  first computes the solutions $\vm{u}_q^p$ of (<ref>) by calling the  MICRO-SCALE_RESOLUTION procedure (algorithm <ref> detailed in <ref>).
This procedure creates vectors $\vm{B}_q^p$   from the constant terms and the coupling Dirichlet terms multiplied by the current subset of $\vm{u}_F$.
It then solves the systems with a direct solver ($\ddagger$ of the algorithm <ref>),  factoring only once the matrices.

Having all the $\vm{u}_q^p$, the global to local step is done and the local to global step can start by calling the  UPDATE_MACRO_PRB procedure ( algorithm <ref> detailed in <ref>).
This procedure computes the enrichment function from the solution $\vm{u}_q^p$  and using  $\vm{A}_{gg}^{ini}$ and $\vm{B}_{g}^{ini}$ finishes the algebraic computation of $\vm{A}_{gg}$ and $\vm{B}_{g}$.
The system (<ref>) is then solved ($\dagger$ of the algorithm <ref>) providing $\vm{U}_g$ solution.

A call to  UPDATE_MICRO_DOFS with this new $\vm{U}_g$ solution updates $\vm{u}^{ts}$ field so that the $resi$ criterion   (<ref>) can be computed for the current iteration.
This is done by COMPUTE_RESIDUAL procedure described in the  next section.
Depending on whether $resi$ is greater than the given  $\epsilon$, the loop must continue or not.

§.§ residual criterion computation

The numerator of the equation (<ref>) can be divided as follows:

\begin{equation}
resi=\frac{\left \| \begin{array}{c}
\vm{A}_{hr}.\vm{u}_{r}-\vm{B}_{h}\\
\vm{A}_{fr}.\vm{u}_{r}-\vm{B}_{f}
\end{array} \right \|}{\left \|\vm{B}_{r} \right \|}
\label{residual_error_split}
\end{equation}

In this work, the h-set part of the residual vector is considered as zero.
This  is justified by the fact that the discretization of the NSP elements  at both scales is the same: the h-set part of the r-set is the same as the h-set part of the g-set.
Thus, since the  $\dagger$ resolution, when performed by a direct solver, provides a "zero" residual error of the g-set system (zero at the machine accuracy), all  this residual vector is zero and in particular its h-set rows.
The computation of $resi$ is therefore simplified to:

\begin{equation}
resi=\frac{\left \|
\vm{A}_{fr}\cdot \vm{u}_{r}-\vm{B}_{f}
\right \|}{\left \|\vm{B}_{r} \right \|}
\label{residual_error_simp}
\end{equation}

The organization of the data by macro-element both simplifies and complicates the computation of  (<ref>) performed by the COMPUTE_RESIDUAL procedure  (algorithm <ref> detailed in <ref>).
The simplicity comes from the computation by $e_{macro}$  of the vector $\vm{A}_{FF}^{e_{macro}}\cdot\vm{u}_{F}^{e_{macro}}-\vm{B}_{F}^{e_{macro}}$ with all the  vectors and the matrix available as data attached  to the $e_{macro}$ element.
The complexity comes from the accumulation of the special local scalar product of these vectors to obtain term $\left \| \vm{A}_{fr}\cdot \vm{u}_{r}-\vm{B}_{f} \right \|$ .

The $\left \|\vm{B}_{r} \right \|$ term of the equation (<ref>) is computed by the procedure COMPUTE_B_NORM  described by the algorithm <ref> in <ref>.
This algorithm is very similar to the algorithm <ref> in that the  scalar product of the vector $\vm{B}_r$  by itself is transformed into  a sum of local dot product.
But as $\vm{B}_r$ remains unchanged during the  loop, this procedure is called only  once during the initialization by the  INIT procedure.
The $\left \|\vm{B}_{r} \right \|$ term is therefore supplied as an argument to the COMPUTE_RESIDUAL procedure  which computes the $resi$ numerator and just divide by this argument to obtains $resi$.

§.§ enrichment function

In this work we choose to approach blending element problem  by using a shift to cancel the enrichment function at the coarse-scale enriched node location and by forcing discretization transition to be embedded in patches, contrasting with SGFEM [4] presented in (<ref>).
This last point is somewhat similar to the idea in [9] where enrichment was added to nodes initially not enriched  in blending elements, with the use of a special ramp function to correct the original enrichment function.
Here, mixed discretization patches are added, so enrichments are added (red node in figure <ref>)  on initially unenriched  blending elements nodes (unenriched  because their support did not cover the behavior being captured).
And the treatment of hanging nodes replaces the effect of the ramp function.
In this way, the original blending element (orange element connected to red node in figure <ref>)  regain their partition of the unity property and the new blending elements (yellow element in figure <ref>) are not perturbed  thanks to the use of a null enrichment function.
Some tests not provided in this work show that the proposed enrichment give a good convergence compare to that given by the SGFEM enrichment.
In this work, the  enrichment function $\vm{F}^p(\vm{x})$ is :

\begin{equation}
\forall p\in I_{e}^g,~\forall \vm{x}\in \omega_p ~\vm{F}^p(\vm{x})=\vm{u}^{Q^p}(\vm{x})-\vm{u}^{Q^p}(\vm{x}^p)~ \text{and}~ \forall \vm{x}\notin \omega_p~\vm{F}^p(\vm{x})=\vm{0}
\label{shift_enrich}
\end{equation}

* $\vm{u}^{Q^p}(\vm{x})$ and $I_{e}^g$ are defined in section <ref>

* $\vm{x}^p$ is the location of the enriched node $p$

Another specificity of this work is the wish that $\vm{u}^{ts}$ converges to  $\vm{u}^{R}$ which can only be achieved if the enrichment function is linearly interpolated on the reference mesh following  the construction and use of $\vm{T}_{MG}$.
Thus, the  kinematic equation actually used at the global scale, corresponding to the application of the $\vm{T}_{MG}$ operator, is:

\begin{equation}
\vm{U}(\vm{x})=\displaystyle\sum_{i\in I^g}  N^i(\vm{x})~ \vm{U}^{i} + \displaystyle\sum_{p\in I_{e}^g} \vm{E}^{p} \displaystyle\sum_{m\in {\bar{I}_p^l}} n^m(\vm{x})~  N^p(\vm{x}_m)~\vm{F}^p(\vm{x}_m)
\label{kinematicnew}
\end{equation}

where $\bar{I}_p^l={I}_p^l\setminus \left\lbrace \text{hanging nodes} \right\rbrace$ since the L-set has also been eliminated in the patch problems.
The  construction of the second term of (<ref>) is first illustrated, in 1D, in the figure <ref> where  artificial elevations show the different quantities of the 8 micro-edges of a  $p$ patch (composed of 2 macro-edges not shown here).

[$\displaystyle\sum_{m\in {\bar{I}_p^l}}  n_m(\vm{x})N^p(\vm{x}_m)~\vm{F}^p(\vm{x}_m)$]
Construction of an interpolated generalized approximation function  for a 1D patch $p$.
For  visualization  purposes, the quantities   here are  artificial elevations of the patch mesh.

The same type of illustration of the 2D example in figure <ref> is shown  in figure <ref>.
Again, artificial elevations  show the various quantities outside the  plane formed by a patch $p$.

[$\displaystyle\sum_{m\in {\bar{I}_p^l}} n_m(\vm{x})N^p(\vm{x}_m)~\vm{F}^p(\vm{x}_m)$]

Construction of an interpolated generalized approximation function  for a  $p$ patch (based on the example in figure  <ref> ).
For visualization purposes, the quantities transformed into scalars are   drawn here as artificial elevations outside the plane formed by the patch mesh. (*enriche_construct_a,*enriche_construct_b) and (*enriche_construct_c,*enriche_construct_d) are represented from two different viewpoints in space for better visualization.

The solutions $\vm{u}^{Q^p}(\vm{x})$  (figure <ref> and <ref>) are completely arbitrary.
The enriched function $\vm{F}^p(\vm{x})$ corresponding to the equation (<ref>) is given in  figure <ref> and <ref>.
By construction, this function is null at the location of the enriched node.
When multiplied by the standard finite element approximation function to form the generalized finite element approximation function (figure <ref> and <ref>), all the boundary values of the  patches become zero and we obtain second order curves/surfaces.
The final linear interpolation is shown in the figure <ref> and <ref>.
In the 2D example, where a mixture of refined and unrefined elements forms the patch, the interpolated generalized function is necessarily zero on the unrefined element (yellow blending  elements in the figure <ref>).
This property implies that blending elements have unenriched nodes  and enriched nodes  with zero contribution from enrichment function.
Thus, the enrichment contribution is zero for the blending element and the partition of unity is saved because only the classical terms matter.
For refined elements, $N^p(\vm{x}_m)~\vm{F}^p(\vm{x}_m)$ is not zero for any location inside the patch (if $\vm{F}^p(\vm{x}_m)\neq\vm{0}$).
For the  hanging nodes, since they are eliminated in $\bar{I}_p^l$ their values are therefore zero which is consistent with the values in the unrefined elements.

Note that if the refined elements correspond only to a terminal connection specific to the adaptation tool (element added in between  figure <ref> and figure <ref> in <ref>), then, despite  the refinement, the associated generalized finite element approximation function is zero on these elements.
To avoid this situation, which can lead to a singular problem, the refined SP elements are always divided at least once.

§ THEORETICAL SEQUENTIAL ALGEBRA PERFORMANCE

Since the method uses a resolution with sparse matrix of different patterns and densities, and since the number of iterations of the loop is not known in advance, estimating the algebraic performance of the method is difficult in a general context.
It is even more complex when using parallelism.
Thus, to obtain the  numerical complexity of the solver, we will use a simple sequential numerical example to be able to count the number of flops [flop=floating-point operation: one addition, subtraction, division or multiplication of two floating-point numbers]  in each part of the algorithm.
The  number of flops used will then be expressed as  a polynomial function of the problem dimensions keeping only the leading terms to obtain the order of magnitude ($\mathcal{O}(n)$ notation).
If we take a meshed cube with a finite hexahedral element of order one in 3D, we can, using octree refinement, simply count the coarse scale dofs, the patches and the patches' dofs.

For the linear resolution steps of the method, only direct sparse solver will be considered.
Thus, the  number of flops for the factorization, in terms of the $n$ dimension of the problem, is approximately between $n^2$  and $n^3$ (dense matrix).
And for backward and forward substitution the cost is in between $2.n$ and $2.n^2$.
Depending on the strategy of the solver (multi-frontal,left-looking,...), the reordering of the symbolic factorization (nested dissection, minimum degree, ...) and the density of the initial matrix, the number of flops varies greatly and depends on a variable number of parameters.
Thus, to estimate these quantities, a sparse ratio parameter $SR$ is introduced such that the dimension used to count the number of flops $n_f$ is:

\begin{equation}
n_f(n,SR)=\sqrt{SR. n^2}
\end{equation}

with $SR\in ]0,1]$.
Then  the number of flops for factorization ($count_{fact}$), backward/forward computation  ($count_{bf}$) and resolution ($count_{resolv}$) are, using a dense estimation:

\begin{gather}\tag{\ref{TAP:costr}}
count_{resolv}(n,SR)=count_{fact}(n,SR)+count_{bf}(n,SR) \\
\begin{align}
count_{fact}(n,SR)=(n_f(n,SR))^3 \\
\end{align}
\end{gather}

This $SR$ parameter  is completely arbitrary.
It is introduced to verify that the conclusions obtained in this section are not impacted by the lack of accurate estimation of the factorization and solving phase.

For discretization, the root of the octree  uses one element to represent the problem.
It corresponds to  zero level.
A level, noted $L$ in the following, describes the depth in the octree tree.
With $L=0$ being the starting level, to go from $L$ to $L+1$, the octree refinement adds a new node in the middle of all edges, faces and elements of level $L$.
Then, each element of  the level $L$ is replaced by  8 hexahedra encapsulated in it and connected to the old and new nodes.
The number of classical dofs for a level $L$ is easily deduced from this refinement strategy and is :

\begin{equation}nb_{dof}(L) = 3.(2^{L}+1)^{3}
\label{TAP:nbdofL}
\end{equation}

Thus, the coarse enriched problem defined at the $L_c$ level  will cost for $N_l$ iteration with a sparse ratio parameter $SR_c$ using (<ref>) and (<ref>):

\begin{equation}
\label{TAP:costc}
cost_{coarse}(L_{c},N_l,SR_c)=N_l\times count_{resolv}(2\times nb_{dof}(L_{c}),SR_c)
\end{equation}

[][ratio of floating point estimation (full rank/resolution) with $N_l=30$,$SR_c=SR_p=0.017$]

floating point estimation for different value of $L_c$,L, $SR_C$ (sparse ratio for coarse problem), $SR_F$(sparse ratio for fine problems) and $N_l$ (number of  iteration)

For a given  coarse level $L_c$, it is possible to obtain in the same way the cost of computing patches (factorization and solving) for $N_l$ iterations, a sparse ratio parameter $SR_p$ and a final refinement level $L$.
This is given by the function $cost_{patch}(L_c,L,N_l,SR_p)$ detailed in <ref>.
Using (<ref>) and (<ref>) we can state that the cost of the method is:

\begin{equation}
\label{TAP:cost_ts}
\end{equation}

Fixing  $N_l,SR_c,SR_p$  the equation (<ref>) for some values of $L_c$ and $L$ give the surface of figure <ref>.
We can see on this continuous representation of (<ref>) that the iso-$L$ curves have a minimum for a specific $L_c$ value.
It means that method has an optimal coarse level for a fixed target level that provides the best performance (i.e. the lowest number of flop consumption).
This can be understood by looking at figure <ref> for a specific target $L=8$  and fixed $N_l,SR_c,SR_p$ where the equation (<ref>) is divided into patch contribution (equation (<ref>)) and coarse-level contribution (equation (<ref>)).
In this figure, it is clear that the antagonistic  evolution of the performance of patch solving and coarse-scale problem solving  lead to having a minimal value at the intersection of the curves.
In figure <ref> the variation of the parameters $N_l,SR_c,SR_p$ does not change the fact that a minimum exists.
In this example where $L=8$, when the number of loop increases, the minimum remains the same.
When we consider the variation of the sparsity of the matrix (up to dense case), the minimum is moved around its original position depending on whether the coarse problem or the patch problems are denser or not.
With denser matrices for patches, it is naturally preferable to increase $l_c$ to have smaller patches.
And with a denser coarse matrix, it is better to decrease $l_c$ to have a smaller coarse problem to solve.
When both sparsities increase, the minimum remains the same but with a higher value.

Figure <ref> presents, as a ratio, an estimate  of the performance of a full rank solver (using the equation (<ref>)) compared to the solver.
A ratio greater than, equal to, or less than 1 indicates that the full rank solver uses more, the same, or less flop than the solver respectively.
The surface  naturally has a maximum for iso-$L$ curves with the same $L_c$ as the minimum of the figure <ref> (on iso-$L$ curves the cost of the full rank solver is constant).
In this example, the  theoretical performance of the solver compared to the full rank solver is much better for $L>2$.
For $iso-L=2$ the optimal ratio at $L_c=1$ is only of 1.76.
And for $L=1$, as in this case the patch problems size are almost equivalent to the $L_c$ level of discretization  no gain is obtained.
This illustrates the fact that the method, to be effective, must be used at a target level that allows to start from a coarse level that induces a scaling effect.

The conclusion drawn from the analysis of this octree refinement example were also observed in parallel with the different test cases studied bellow.
Parallelism, the nature of the unstructured mesh, and the cost of assembly also impact the choice of the optimal coarse mesh for a given fine-scale target discretization.
But the observation of the test case <ref> led to use the following formula to define the scale jump for most cases with regular refinement:

\begin{equation}
\label{ts_jump_level_eq}
\end{equation}

[][Weight: micro embemded element count	]

[][Global element dispatching on 4 processes partially based on weight (*ts_dist_weigth_elem)]

(0.89885423,0.88901305)(0,0)[lt]0t]lprocess 1
(0.92777561,0.46741474)(0,0)[lt]0t]lprocess 2
(0.40797012,0.15228167)(0,0)[lt]0t]lprocess 3
(-0.04530285,0.99875923)(0,0)[lt]0t]lprocess 0

[][Global element dispatching on 8 processes partially based on weight (*ts_dist_weigth_elem)]

(-0.09201796,0.74600391)(0,0)[lt]0t]lprocess 7
(0.89885423,0.88901305)(0,0)[lt]0t]lprocess 2
(0.92021202,0.51974065)(0,0)[lt]0t]lprocess 3
(0.75579117,0.14654528)(0,0)[lt]0t]lprocess 5
(-0.04530285,0.99875923)(0,0)[lt]0t]lprocess 0
(0.56634056,1.17854491)(0,0)[lt]0t]lprocess 1
(0.91647254,0.38665167)(0,0)[lt]0t]lprocess 4
(0.01965601,0.18307101)(0,0)[lt]0t]lprocess 6

(0.24956891,0.69252891)(0,0)[lt]0t]lGlobal scale node with distributed support
(0.24961644,0.59985417)(0,0)[lt]0t]lGlobal scale node with non distributed support
(0.25355642,0.30656925)(0,0)[lt]0t]lGlobal-scale mesh
(0.40199406,0.04134013)(0,0)[lt]0t]lProcess id color
(0.31214749,0.51877158)(0,0)[lt]0t]lInterprocess interface
(0.26072418,0.18327801)(0,0)[lt]0t]lFine-scale mesh
(0.31484353,0.43445197)(0,0)[lt]0t]lProcess boundary

[][Fine scale problems (2 out 25) following global element dispatching on 4 processes (*ts_dist_load_balance)]

(0.03114623,0.68638138)(0,0)[lt]0t]lGlobal comm: process 1 / Local comm: process 0
(0.07733094,0.03550769)(0,0)[lt]0t]lGlobal comm: process 3 / Local comm: process 1
(0.33973626,0.2837736)(0,0)[lt]0t]lGlobal comm: process 0 / Local comm: process 0

fictitious 2D problem of figure<ref> distributed on processes

§ PARALLEL PARADIGM

This work does not use shared-memory multithreaded parallelism and uses only the distributed-memory message passing paradigm [In this work we use the MPI-3.0 (or higher) message passing standard.].
This choice comes from the thread-safety constraint which would have generated too much implementation work for this first study (for example the in-house library used to code the  method is not thread-safe).
As a result, multi-threading is never used in this study and to be fair when comparing with other methods, their multi-threading capability was also disabled.
Thus, in this section, parallelism will always concern the message passing paradigm.

Parallelism is introduced at all scales in all frameworks, but with a quite different granularity.
At the global level, the mesh is distributed over processes (see  figures <ref> and <ref> for the example of figure <ref> distributed over 4 and 8 processes respectively), based on a macro-element partitioning,  in order to allow a good load balancing.
This partitioning (not evaluated in depth in this paper), given by ParMetis ([16]) in this work, is supposed to balance the cost of creating linear systems, assembling them and solving them at both scales.
To do this, a cost per macro-element is obtained  directly  from  the number of its encapsulated micro-elements (see figure <ref>) and the number of enriched nodes.
If the entire global-scale mesh is refined (as in the examples in the section <ref> and <ref>), all macro elements have roughly the same number of encapsulated micro elements and enriched nodes.
Thus, the cost per element can be considered as a constant.
Now, if only a specific region is involved (see the 2D example in Figure <ref> or the example in section <ref>), the number of encapsulated micro-elements and enriched nodes provides weights to the ParaMetis multi-objective partitioning optimization algorithm.
In the examples studied, this partitioning provides great scalability for the assembly task, as all processes receive roughly the same amount of micro-elements to integrate and assemble (see in particular the table <ref>  in section <ref>).
It also allows the use of a fairly well balanced parallel resolution for the linear system related to the global-scale problem ($\dagger$ in the algorithm <ref>).
The section <ref> gives more details on the choice of solver to use at this level.

At the fine scale, patches follow the distribution of the global scale mesh (see figure <ref>) and are therefore naturally distributed in most processes.
This results, as in [18], in a high level of coarse-grained parallelism because many fine-scale problems are solved independently (without communication) at the same time  by many processes.
But some patches are split between the computational units due to the global mesh distribution (see the left patch in the figure <ref>).
In this work the distributed patch problems  are solved entirely in parallel without any overlapping mechanism (phantom variables, ...).
The  assembly of the problems ($\vm{A}_{qq}^{p}$ and $\vm{BI}_{q}^{p}$ of algorithm <ref> and $\vm{B}_{q}^{p}$ of algorithm <ref>) is thus done in parallel without communication.
Then, the matrix contributions of the linear system are provided to a parallel solver ($\ddagger$ of the algorithm <ref>) and treated with a mid-grain parallelism.
But this introduces some constraints on the overall process of fine-scale parallel resolution.
The distributed patches are no longer independent.
Typically, if two distributed patches share the same process, it will be impossible to run the resolution of both in that process.
The solver is not able to treat two independent matrices at the same time.
This implies that the computation of distributed patches  must be scheduled as explained in <ref>.
These distributed patches have an ambivalent impact on performance.
On the one hand, they reduce the time to solve a rather large fine-scale problem.
And on the other hand, if their number increases too much, it expands the amount of communication and can also have an impact on sequencing.

§.§ global scale problem parallel resolution

In this work, we first chose an asynchronous parallel MPI  sparse  direct  solver (MUMPS [2, 3]) for the resolution of  the linear system associated with the global problem.
It was a safe choice (even if the matrix becomes badly conditioned) for this first implementation but by no means a requirement.
The symbolic step is performed only once at the beginning of the loop.
Factorization and solving step are repeated for each iteration with the updated $\vm{A}_{gg}$ and $\vm{B}_{g}$.
With many patches, as the fine scale resolution has good scalability (see the results presented in section <ref>), many processes can be used well beyond the scalability of the parallel direct solver which solves at the coarse-scale, a problem of rather smaller size.
To limit this effect, in all simulations tested,  a subset of the available processes (taking at most $nbp_{max}$ processes, $nbp_{max}$ arbitrarily chosen to stay within the range of good performance of the direct solver) was used to solve this coarse-scale problem (gathering $\vm{A}_{gg}$ contribution of the processes eliminated from the resolution in the $nbp_{max}$ retained processes).
This version is referred to as below.
This reduction had a positive effect, but the overall speed-up[the speed-up is the ratio between the time used by an application in single-process mode and the time used in multi-process mode.
Ideally, it should give the number of processes used. ] of the resolution was still deteriorated by the coarse-scale efficiency.
To reduce this effect, as the loop iterates over a coarse-scale problem without changing its discretization, it is possible to use a parallel preconditioned iterative solver (conjugate gradient) that starts with the solution of the previous iteration.
This choice  intrinsically reduces the coarse-scale resolution by considering that the successive solutions are not too far apart from each other, so that the number of iterations in the iterative solver should be small.
And this choice  adds scalability because the conjugate gradient scales well (to the point where the scalar product becomes a problem).
A natural choice is to use for  the preconditioner  a previous factorization performed by the direct solver.
The preconditioning step  then only costs a forward/backward resolution, performed in parallel, following the distribution of the direct solver.
If not done to many times, it is cheap compared to a full factorization.
Note that this technique does not take into account the stable nature of the classical dofs of the coarse-scale problem as  is the case in [19].
This would have added additional complexity regarding load balancing and resolution.
Nevertheless, the proposed approach offers apriori a slightly richer preconditioner  compared to the block Jacobi proposed in this quote because it incorporates all terms (including the classic x enriched coupling blocks).
The use of the iterative solver and the quality of the preconditioner  are controlled using the following arbitrary rules:

* The iterative resolution is used instead of the direct resolution as soon as the criterion "$resi$"  becomes smaller than $\epsilon\times 10000$ ($\epsilon$ being the precision to be reach by "$resi$" in algorithm <ref>  ) and as soon as 2 iterations have been performed

* When the iterative resolution takes more than 13 iterations to converge, the preconditioner is recomputed: at the next iteration, the direct solver is used again to provide a new solution and a new factorization.
The iterative solver is then reused in the next iteration.

The second rule avoids spending too much time in the iterative solver when the  preconditioner is not of high quality.
The first rule requires  waiting for the solution to reach a minimum quality so that the variation between the iteration is small enough for the preconditioner to efficiently reduce the number of iterations of the iterative solver.
This version is hereafter referred to  as .
Also note that in this version, for the direct resolution phase,  low-rank resolution (see <ref>) can also be used instead of full rank resolution.
This is left as future work.

At the coarse scale, another type of solver ( a domain decomposition solver presented in <ref>) was also tested during this study in order to gain scalability as illustrated in section <ref>.

[][weight per patch]

[][Initial graph with arbitrary node id]

[][sequence 1]

[][sequence 2 ]

[][sequence 3]

[][sequence 4]

[][Sequence 5]

[][Sequence 6]

[][Sequence 7]

[][Sequence 8]

[][Remaining local patch]

(0.20066397,1.2462974)(0,0)[lt]0t]ledge related to process 0
(0.20066397,1.13389424)(0,0)[lt]0t]ledge related to process 1
(0.20066397,1.02149109)(0,0)[lt]0t]ledge related to process 2
(0.20066397,0.90908791)(0,0)[lt]0t]ledge related to process 3
(0.20066397,0.79668479)(0,0)[lt]0t]lfirst node selected for this sequence
(0.20066397,0.68428159)(0,0)[lt]0t]ledge folowed to find impacted node by first choice
(0.20066397,0.57187843)(0,0)[lt]0t]lblocked node by first choice
(0.20066397,0.4594752)(0,0)[lt]0t]lother node(s) selected for this sequence
(0.20066397,0.34707211)(0,0)[lt]0t]ledge folowed to find impacted node by other choice(s)
(0.20066397,0.23466895)(0,0)[lt]0t]lblocked node by other choice(s)
(0.20061921,0.1234222)(0,0)[lt]0t]lnode with distributed support
(0.20066397,0.00986271)(0,0)[lt]0t]lNode with non distributed support

Patch sequencing for figure <ref> example. At each sequence, a set $\mathcal{S}$ is chosen so that as many distributed patches as possible are solved at the same time and during the same time period.  A sequence is first built  by the lowest process which chooses  one of its distributed nodes  (red node) in the graph. Then, all the other processes that are not already blocked  by this first choice do the same (blue node) in ascending order of the process identification numbers. See <ref> for a detailed description and the table <ref> for each $\mathcal{S}$. After a sequence is constructed, its nodes and connected edges are removed from the graph to give a new sub-graph used to create the next sequence until no more distributed nodes are available.

§.§ local scale problem scheduling

The 2D example in figure <ref> will be used to illustrate the static scheduling algorithm.
As already mentioned, two distributed patches that share the same process cannot start solving their linear system at the same time in that process.
To avoid this situation, a simple idea is to force distributed patches that share the same process to be computed at different times.
To do this, a static[Avoiding dynamic scheduling is normally more efficient] scheduling is created that  requires  all distributed patches over all processes to be computed in a specific order during the loop.
The general idea is to solve as many distributed patches as possible at once, thus maximizing the use of processes and reducing the number of sequences.

To achieve this goal, the first mandatory task is to identify the dependencies between the distributed patches so that they can be sequenced.
A distributed patch depends on another if it shares a process with it.
This dependency can be extended to local patches to create an undirected graph where the vertices are all patches and the edges represent a process connection (i.e. two patches are connected by an edge if they belong even partially to the same process).
Using the mesh distribution in figure <ref> this graph concept is illustrated in figure <ref>: the vertices of this graph are patches (here, for a quick visual recognition of patches, the vertices of the graph are located at the same place as the enriched nodes of the mesh and  each of them has been assigned an arbitrary identifier), the edges of this graph represent a process connection (in this figure, the edges are colored the same way as macro-element of figure <ref> to clearly identify the process identifier given in figure <ref>).

From this undirected graph, it is possible to create a set $\mathcal{S}$ of distributed vertices such that for all vertices in $\mathcal{S}$, no edge connect them.
Thus, these vertices  can be computed at the same time because they do not share a common process (no edges connect them).
Finding this  set $\mathcal{S}$ (not unique)  corresponds in graph theory to  the "independent set problem" and will correspond to finding a sequence of the scheduling.
It is interesting to maximize the size of $\mathcal{S}$ (which corresponds to the NP-hard  "maximum independent set problem" in graph theory) because on the one hand it maximizes a priori the number of used processes  and therefore the load balancing is good. On the other hand, the distributed vertices are consumed faster and therefore the number of sequences can be minimized.
When  the set $\mathcal{S}$  is assigned to a sequence (mainly by creating a colored  MPI communicator associated with it),  its vertices (and their connecting edges) can be removed from the graph.
The maximum independent set problem is then solved on the remaining subgraph to find the next sequence.
All sequences are found this way  until all distributed patches are removed from the subgraph.
Then the remaining local patches can be processed in any order as they are independent.

[distribution of figure <ref>]
8|c|Total number of distributed patches: 11

8|c|Maximum number of distributed patches per process: 8

8|c|Minimum number of distributed patches per process: 4

8|c|Average number of distributed patches per process: 6

[origin=c]90sequence    [origin=c]90Active patches    [origin=c]90Active distributed [origin=c]90patches    [origin=c]90% Active distributed [origin=c]90 patches	    [origin=c]90Active processes
[origin=c]90% Active processes
[origin=c]90$\mathcal{S}$ 	   	 [origin=c]90Figure

1   2   1   50   4   100   18,3   <ref>

2   3   1   33.3   4   100   12,2,20   <ref>

3   2   1   50   4   100   10,21   <ref>

4   3   1   33.3   4   100   13,1,23   <ref>

5   2   2   100   4   100   17,6   <ref>

6   2   2   100   4   100   7,19   <ref>

7   3   1   33.3   4   100   8,11,22   <ref>

8   2   2   100   4   100   15,5   <ref>

average   2.4   1.4   62.5   4   100   2c

[distribution of figure  <ref>]

7|c|Total number of distributed patches: 17

7|c|Maximum number of distributed patches per process: 6

7|c|Minimum number of distributed patches per process: 4

7|c|Average number of distributed patches per process: 5.1

[origin=c]90sequence    [origin=c]90Active patches    [origin=c]90Active distributed [origin=c]90patches    [origin=c]90% Active distributed [origin=c]90patches	    [origin=c]90Active processes
[origin=c]90% Active processes

1   3   3   100   8   100   12,7,17

2   3   2   66.7   8   100   10,18,21

3   3   3   100   8   100   13,6,20

4   4   4   100   8   100   8,5,19,15

5   4   4   100   8   100   14,3,11,25

6   3   1   33.3   4   50   9,1,23

average   3.3   2.8   83.3   7.3   91.7   1c

Distributed patches sequencing with the algorithm <ref> for the   fictitious 2D problem in figure<ref> distributed over 4 and 8 processes (figure <ref>). It does not correspond  the calculation of all patches  because the average does not take into account the computation of the non-constrained local patches.

The condition of maximizing the size of $\mathcal{S}$  is not sufficient to obtain a balanced load.
The vertices of $\mathcal{S}$ must have the same "load" (in a sequence, if a vertex $v_i$ of $\mathcal{S}$ takes longer to compute than the others, most of the process will be idle while waiting for $v_i$ to finish).
Thus, with distributed patches  having different computational consumption (called weights hereafter), the choice of $\mathcal{S}$ must  be oriented by these weights in order to have approximately the same "load" for all the vertices of  $\mathcal{S}$.
In this work, these weights are simply the number of micro-elements in a patch which reflects its number of dofs and thus its computational cost (see figure <ref> computed from the weights in figure <ref>).
The NP-hard nature of the maximum independent set problem, this weighting constraint, and the a priori distributed nature of the undirected graph lead to the choice of a basic sub-optimal  heuristic to quickly find a static scheduling.

More precisely, the chosen algorithm, formalized in the algorithms <ref>, <ref> and <ref>  detailed  in <ref>, is played in the table <ref> and figure <ref> for the example of the figure <ref>.
In this 2D example, 8 sequences are required, which corresponds to the maximum number of distributed patches per process: in the process with 8 distributed patches, there cannot be less than 8 sequences.
With the example of figure <ref>, distributed on  8 processes (see figure <ref>), the table <ref> shows that in this case 6 sequences must be computed.
But in this case, the last sequence is less efficient because process 3, 4, 5 and 7 are note used.
This shows that the scalability of the distributed method will be influenced by the quality of the algorithm <ref>.
In the section <ref>, and in particular in the table <ref> of this section, we study in a real example the heuristic effect of this algorithm.

§ NUMERICAL SIMULATION

All simulation were run on the Liger cluster (see  <ref> for  details) in the exclusive node condition (no other users at the node level but  shared resources at the network and I/O level).
In many cases, the solver  is  compared  to the full rank direct solver MUMPS (denoted by "fr" hereafter), to a low-rank direct solver based on MUMPS (see <ref>, denoted by "blr" hereafter) and to a domain decomposition solver (see <ref>, denoted by "dd" hereafter).
Hereinafter, the  and versions are referred to as  "ts" and "tsi" respectively.
In terms of notation E will be the Young's modulus and $\nu$ the Poisson's ratio.

§.§ Analytic solution

This test is intended to validate the proposed work using a given cubic displacement field over a plate:

figureAnalytical test:  plate under volume and surface loading. A, B, and C are the corner points used to fix the rigid body modes.

\begin{equation}
\Vec{u}^C(x,y,z)= \frac{(\nu+1)(1-2\nu)F}{E}\left[ \begin{array}{l} (x^2(\frac{K}{2}-\frac{x}{3})+2\nu(y^2-z^2))\\
-4\nu xy\\
4\nu xz
\end{array}
\right]
\label{eqUC}
\end{equation}

where :

* $E=36.5$GPa

* $\nu=0.2$

* $F$ is a scalar corresponding to a force

* $K$ is the plate length  (figure <ref> )

The imposed loads are obtained from the equation (<ref>) (writing the equilibrium of the system) and the plate is made of a  homogeneous isotropic elastic material.
The plate is discretized with a coarse mesh (595 tetrahedrons, 188 vertices) which corresponds, hereafter, to the level 0 (L0).
Up to 4 other discretizations (L1,L2,L3,L4) derived from this level are used as the global-scale problem.
Each of them is built using the adaptation strategy described in <ref>  (considering the whole mesh as an area of interest: no hanging nodes).
Each level corresponds to the previous one with all its elements divided once (L2 is L1 divided once).
This corresponds exactly to the refinement process used for scaling.
Thus, since all nodes are enriched, the 4 meshes are identical to  any refined SP targeting the same level (L2 mesh will be the same as the refined created by refining L0 mesh 2 times).

We evaluate the performance of the proposed parallel solver  and  compare it to the solvers "fr", "blr" and "dd".
In all figures, for the curves,  "L$Y$ from L$X$" or "from L$X$" (when the fine scale-level "$Y$" is implicit) means to use a "$X$" level  for the discretization of the global level and to use  a "$Y$" level  for the discretization of the ($Y>X$).
Since the full rank direct solver provides a residual error close to the machine accuracy, a relatively fair comparison requires to adopt a value of $1.e-7$  for $\epsilon$ in the algorithm <ref>,<ref> and <ref>.
As noted in section <ref>, for comparison purposes, the multithreading capability of MUMPS and underlying BLAS is not used.
Thus, the number of cores (abscissa in many figures) directly represents  the number of MPI process dispatched on one or more nodes of the cluster.

every node/.append style=scale=0.65
(0.000,0.000) rectangle (5.625,5.206);
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (0.860,0.640)–(5.264,0.640);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (0.860,0.640)–(1.040,0.640);
[gp path] (5.264,0.640)–(5.084,0.640);
[gp node right] at (0.740,0.640) $0.1$;
[gp path] (0.860,1.334)–(0.950,1.334);
[gp path] (5.264,1.334)–(5.174,1.334);
[gp path] (0.860,1.739)–(0.950,1.739);
[gp path] (5.264,1.739)–(5.174,1.739);
[gp path] (0.860,2.027)–(0.950,2.027);
[gp path] (5.264,2.027)–(5.174,2.027);
[gp path] (0.860,2.251)–(0.950,2.251);
[gp path] (5.264,2.251)–(5.174,2.251);
[gp path] (0.860,2.433)–(0.950,2.433);
[gp path] (5.264,2.433)–(5.174,2.433);
[gp path] (0.860,2.587)–(0.950,2.587);
[gp path] (5.264,2.587)–(5.174,2.587);
[gp path] (0.860,2.721)–(0.950,2.721);
[gp path] (5.264,2.721)–(5.174,2.721);
[gp path] (0.860,2.839)–(0.950,2.839);
[gp path] (5.264,2.839)–(5.174,2.839);
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (0.860,2.944)–(5.264,2.944);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (0.860,2.944)–(1.040,2.944);
[gp path] (5.264,2.944)–(5.084,2.944);
[gp node right] at (0.740,2.944) $1$;
[gp path] (0.860,3.638)–(0.950,3.638);
[gp path] (5.264,3.638)–(5.174,3.638);
[gp path] (0.860,4.044)–(0.950,4.044);
[gp path] (5.264,4.044)–(5.174,4.044);
[gp path] (0.860,4.332)–(0.950,4.332);
[gp path] (5.264,4.332)–(5.174,4.332);
[gp path] (0.860,4.555)–(0.950,4.555);
[gp path] (5.264,4.555)–(5.174,4.555);
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (0.860,0.640)–(0.860,4.555);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (0.860,0.640)–(0.860,0.820);
[gp path] (0.860,4.555)–(0.860,4.375);
[gp node center] at (0.860,0.440) 1;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (1.489,0.640)–(1.489,4.555);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (1.489,0.640)–(1.489,0.820);
[gp path] (1.489,4.555)–(1.489,4.375);
[gp node center] at (1.489,0.440) 2;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (2.118,0.640)–(2.118,4.555);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (2.118,0.640)–(2.118,0.820);
[gp path] (2.118,4.555)–(2.118,4.375);
[gp node center] at (2.118,0.440) 4;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (2.747,0.640)–(2.747,4.555);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (2.747,0.640)–(2.747,0.820);
[gp path] (2.747,4.555)–(2.747,4.375);
[gp node center] at (2.747,0.440) 8;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (3.376,0.640)–(3.376,4.555);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (3.376,0.640)–(3.376,0.820);
[gp path] (3.376,4.555)–(3.376,4.375);
[gp node center] at (3.376,0.440) 16;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (4.005,0.640)–(4.005,4.555);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (4.005,0.640)–(4.005,0.820);
[gp path] (4.005,4.555)–(4.005,4.375);
[gp node center] at (4.005,0.440) 32;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (4.634,0.640)–(4.634,4.555);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (4.634,0.640)–(4.634,0.820);
[gp path] (4.634,4.555)–(4.634,4.375);
[gp node center] at (4.634,0.440) 64;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (5.263,0.640)–(5.263,4.555);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (5.263,0.640)–(5.263,0.820);
[gp path] (5.263,4.555)–(5.263,4.375);
[gp node center] at (5.263,0.440) 128;
[gp path] (0.860,4.555)–(0.860,0.640)–(5.264,0.640)–(5.264,4.555)–cycle;
[gp node center,rotate=-270] at (0.190,2.597) elapsed time in s;
[gp node center] at (3.062,0.140) number of cores;
[gp node right] at (1.460,4.913) fr;
rgb color=0.000,1.000,1.000
[gp path] (1.580,4.913)–(2.240,4.913);
[gp path] (0.860,3.196)–(1.489,2.707)–(2.118,2.237)–(2.747,2.234)–(3.376,2.410)
3pointgp mark 1(0.860,3.196)
3pointgp mark 1(1.489,2.707)
3pointgp mark 1(2.118,2.237)
3pointgp mark 1(2.747,2.234)
3pointgp mark 1(3.376,2.410)
3pointgp mark 1(4.005,2.538)
3pointgp mark 1(4.634,2.571)
3pointgp mark 1(5.263,2.650)
3pointgp mark 1(1.910,4.913)
color=gp lt color border
[gp node right] at (1.460,4.688) blr;
rgb color=0.000,0.000,0.545
[gp path] (1.580,4.688)–(2.240,4.688);
[gp path] (0.860,3.273)–(1.489,2.840)–(2.118,2.443)–(2.747,2.393)–(3.376,2.594)
3pointgp mark 2(0.860,3.273)
3pointgp mark 2(1.489,2.840)
3pointgp mark 2(2.118,2.443)
3pointgp mark 2(2.747,2.393)
3pointgp mark 2(3.376,2.594)
3pointgp mark 2(4.005,2.699)
3pointgp mark 2(4.634,2.776)
3pointgp mark 2(5.263,2.898)
3pointgp mark 2(1.910,4.688)
color=gp lt color border
[gp node right] at (2.960,4.913) dd;
rgb color=0.753,0.251,0.000
[gp path] (3.080,4.913)–(3.740,4.913);
[gp path] (0.860,3.196)–(1.489,2.830)–(2.118,2.373)–(2.747,2.090)–(3.376,1.886)
3pointgp mark 15(0.860,3.196)
3pointgp mark 15(1.489,2.830)
3pointgp mark 15(2.118,2.373)
3pointgp mark 15(2.747,2.090)
3pointgp mark 15(3.376,1.886)
3pointgp mark 15(4.005,1.507)
3pointgp mark 15(4.634,1.307)
3pointgp mark 15(5.263,1.233)
3pointgp mark 15(3.410,4.913)
rgb color=0.580,0.000,0.827
[gp path] (0.860,4.490)–(1.489,3.912)–(2.118,3.273)–(2.747,2.862)–(3.376,2.559)
3pointgp mark 4(0.860,4.490)
3pointgp mark 4(1.489,3.912)
3pointgp mark 4(2.118,3.273)
3pointgp mark 4(2.747,2.862)
3pointgp mark 4(3.376,2.559)
3pointgp mark 4(4.005,2.391)
3pointgp mark 4(4.634,2.495)
3pointgp mark 4(5.263,2.187)
rgb color=0.933,0.510,0.933
[gp path] (0.860,4.476)–(1.489,3.866)–(2.118,3.242)–(2.747,2.848)–(3.376,2.469)
3pointgp mark 5(0.860,4.476)
3pointgp mark 5(1.489,3.866)
3pointgp mark 5(2.118,3.242)
3pointgp mark 5(2.747,2.848)
3pointgp mark 5(3.376,2.469)
3pointgp mark 5(4.005,2.315)
3pointgp mark 5(4.634,2.423)
3pointgp mark 5(5.263,2.101)
color=gp lt color border
[gp node right] at (2.960,4.688) ideal;
rgb color=0.000,0.000,0.000
[gp path] (3.080,4.688)–(3.740,4.688);
[gp path] (0.860,4.476)–(0.904,4.427)–(0.949,4.378)–(0.993,4.329)–(1.038,4.280)
color=gp lt color border
[gp path] (0.860,4.555)–(0.860,0.640)–(5.264,0.640)–(5.264,4.555)–cycle;

gp plot 10.860cm0.640cm5.264cm4.555cm

every node/.append style=scale=0.65
(0.000,0.000) rectangle (5.625,4.987);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (0.540,0.640)–(0.630,0.640);
[gp path] (5.264,0.640)–(5.174,0.640);
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (0.540,0.743)–(5.264,0.743);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (0.540,0.743)–(0.720,0.743);
[gp path] (5.264,0.743)–(5.084,0.743);
[gp node right] at (0.420,0.743) $1$;
[gp path] (0.540,1.419)–(0.630,1.419);
[gp path] (5.264,1.419)–(5.174,1.419);
[gp path] (0.540,1.815)–(0.630,1.815);
[gp path] (5.264,1.815)–(5.174,1.815);
[gp path] (0.540,2.096)–(0.630,2.096);
[gp path] (5.264,2.096)–(5.174,2.096);
[gp path] (0.540,2.314)–(0.630,2.314);
[gp path] (5.264,2.314)–(5.174,2.314);
[gp path] (0.540,2.492)–(0.630,2.492);
[gp path] (5.264,2.492)–(5.174,2.492);
[gp path] (0.540,2.642)–(0.630,2.642);
[gp path] (5.264,2.642)–(5.174,2.642);
[gp path] (0.540,2.772)–(0.630,2.772);
[gp path] (5.264,2.772)–(5.174,2.772);
[gp path] (0.540,2.887)–(0.630,2.887);
[gp path] (5.264,2.887)–(5.174,2.887);
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (0.540,2.990)–(5.264,2.990);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (0.540,2.990)–(0.720,2.990);
[gp path] (5.264,2.990)–(5.084,2.990);
[gp node right] at (0.420,2.990) $10$;
[gp path] (0.540,3.667)–(0.630,3.667);
[gp path] (5.264,3.667)–(5.174,3.667);
[gp path] (0.540,4.062)–(0.630,4.062);
[gp path] (5.264,4.062)–(5.174,4.062);
[gp path] (0.540,4.343)–(0.630,4.343);
[gp path] (5.264,4.343)–(5.174,4.343);
[gp path] (0.540,4.561)–(0.630,4.561);
[gp path] (5.264,4.561)–(5.174,4.561);
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (0.540,0.640)–(0.540,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (0.540,0.640)–(0.540,0.820);
[gp path] (0.540,4.561)–(0.540,4.381);
[gp node center] at (0.540,0.440) 1;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (1.130,0.640)–(1.130,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (1.130,0.640)–(1.130,0.820);
[gp path] (1.130,4.561)–(1.130,4.381);
[gp node center] at (1.130,0.440) 2;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (1.721,0.640)–(1.721,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (1.721,0.640)–(1.721,0.820);
[gp path] (1.721,4.561)–(1.721,4.381);
[gp node center] at (1.721,0.440) 4;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (2.311,0.640)–(2.311,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (2.311,0.640)–(2.311,0.820);
[gp path] (2.311,4.561)–(2.311,4.381);
[gp node center] at (2.311,0.440) 8;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (2.902,0.640)–(2.902,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (2.902,0.640)–(2.902,0.820);
[gp path] (2.902,4.561)–(2.902,4.381);
[gp node center] at (2.902,0.440) 16;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (3.492,0.640)–(3.492,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (3.492,0.640)–(3.492,0.820);
[gp path] (3.492,4.561)–(3.492,4.381);
[gp node center] at (3.492,0.440) 32;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (4.083,0.640)–(4.083,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (4.083,0.640)–(4.083,0.820);
[gp path] (4.083,4.561)–(4.083,4.381);
[gp node center] at (4.083,0.440) 64;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (4.673,0.640)–(4.673,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (4.673,0.640)–(4.673,0.820);
[gp path] (4.673,4.561)–(4.673,4.381);
[gp node center] at (4.673,0.440) 128;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (5.264,0.640)–(5.264,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (5.264,0.640)–(5.264,0.820);
[gp path] (5.264,4.561)–(5.264,4.381);
[gp node center] at (5.264,0.440) 256;
[gp path] (0.540,4.561)–(0.540,0.640)–(5.264,0.640)–(5.264,4.561)–cycle;
[gp node center] at (2.902,0.140) number of cores;
rgb color=0.000,1.000,1.000
[gp path] (0.540,3.614)–(1.130,3.073)–(1.721,2.506)–(2.311,2.112)–(2.902,1.660)
3pointgp mark 1(0.540,3.614)
3pointgp mark 1(1.130,3.073)
3pointgp mark 1(1.721,2.506)
3pointgp mark 1(2.311,2.112)
3pointgp mark 1(2.902,1.660)
3pointgp mark 1(3.492,1.625)
3pointgp mark 1(4.083,2.003)
3pointgp mark 1(4.673,2.390)
rgb color=0.000,0.000,0.545
[gp path] (0.540,3.734)–(1.130,3.200)–(1.721,2.579)–(2.311,2.148)–(2.902,1.757)
3pointgp mark 2(0.540,3.734)
3pointgp mark 2(1.130,3.200)
3pointgp mark 2(1.721,2.579)
3pointgp mark 2(2.311,2.148)
3pointgp mark 2(2.902,1.757)
3pointgp mark 2(3.492,1.795)
3pointgp mark 2(4.083,2.131)
3pointgp mark 2(4.673,2.455)
rgb color=0.753,0.251,0.000
[gp path] (0.540,3.608)–(1.130,3.336)–(1.721,2.928)–(2.311,2.957)–(2.902,2.567)
3pointgp mark 15(0.540,3.608)
3pointgp mark 15(1.130,3.336)
3pointgp mark 15(1.721,2.928)
3pointgp mark 15(2.311,2.957)
3pointgp mark 15(2.902,2.567)
3pointgp mark 15(3.492,1.937)
3pointgp mark 15(4.083,1.703)
3pointgp mark 15(4.673,1.253)
3pointgp mark 15(5.264,0.760)
color=gp lt color border
[gp node right] at (1.860,4.694) ts from L0;
rgb color=0.580,0.000,0.827
[gp path] (1.980,4.694)–(2.640,4.694);
[gp path] (0.540,4.046)–(1.130,3.452)–(1.721,2.841)–(2.311,2.379)–(2.902,1.992)
3pointgp mark 4(0.540,4.046)
3pointgp mark 4(1.130,3.452)
3pointgp mark 4(1.721,2.841)
3pointgp mark 4(2.311,2.379)
3pointgp mark 4(2.902,1.992)
3pointgp mark 4(3.492,1.684)
3pointgp mark 4(4.083,1.626)
3pointgp mark 4(4.673,1.249)
3pointgp mark 4(2.310,4.694)
color=gp lt color border
[gp node right] at (4.080,4.694) tsi from L0;
rgb color=0.933,0.510,0.933
[gp path] (4.200,4.694)–(4.860,4.694);
[gp path] (0.540,4.049)–(1.130,3.451)–(1.721,2.838)–(2.311,2.382)–(2.902,1.983)
3pointgp mark 5(0.540,4.049)
3pointgp mark 5(1.130,3.451)
3pointgp mark 5(1.721,2.838)
3pointgp mark 5(2.311,2.382)
3pointgp mark 5(2.902,1.983)
3pointgp mark 5(3.492,1.677)
3pointgp mark 5(4.083,1.634)
3pointgp mark 5(4.673,1.213)
3pointgp mark 5(4.530,4.694)
rgb color=0.784,0.784,0.000
[gp path] (0.540,4.487)–(1.130,3.870)–(1.721,3.202)–(2.311,2.694)–(2.902,2.237)
3pointgp mark 12(0.540,4.487)
3pointgp mark 12(1.130,3.870)
3pointgp mark 12(1.721,3.202)
3pointgp mark 12(2.311,2.694)
3pointgp mark 12(2.902,2.237)
3pointgp mark 12(3.492,1.825)
3pointgp mark 12(4.083,1.499)
3pointgp mark 12(4.673,1.349)
3pointgp mark 12(5.264,1.241)
rgb color=1.000,1.000,0.000
[gp path] (0.540,4.387)–(1.130,3.801)–(1.721,3.179)–(2.311,2.629)–(2.902,2.150)
3pointgp mark 13(0.540,4.387)
3pointgp mark 13(1.130,3.801)
3pointgp mark 13(1.721,3.179)
3pointgp mark 13(2.311,2.629)
3pointgp mark 13(2.902,2.150)
3pointgp mark 13(3.492,1.682)
3pointgp mark 13(4.083,1.260)
3pointgp mark 13(4.673,1.026)
3pointgp mark 13(5.264,0.880)
rgb color=0.000,0.000,0.000
[gp path] (0.540,4.049)–(0.582,4.001)–(0.623,3.953)–(0.665,3.905)–(0.707,3.858)
color=gp lt color border
[gp path] (0.540,4.561)–(0.540,0.640)–(5.264,0.640)–(5.264,4.561)–cycle;

gp plot 10.540cm0.640cm5.264cm4.561cm

every node/.append style=scale=0.65
(0.000,0.000) rectangle (5.625,4.987);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (0.660,0.640)–(0.750,0.640);
[gp path] (5.264,0.640)–(5.174,0.640);
[gp path] (0.660,0.860)–(0.750,0.860);
[gp path] (5.264,0.860)–(5.174,0.860);
[gp path] (0.660,1.032)–(0.750,1.032);
[gp path] (5.264,1.032)–(5.174,1.032);
[gp path] (0.660,1.171)–(0.750,1.171);
[gp path] (5.264,1.171)–(5.174,1.171);
[gp path] (0.660,1.289)–(0.750,1.289);
[gp path] (5.264,1.289)–(5.174,1.289);
[gp path] (0.660,1.392)–(0.750,1.392);
[gp path] (5.264,1.392)–(5.174,1.392);
[gp path] (0.660,1.482)–(0.750,1.482);
[gp path] (5.264,1.482)–(5.174,1.482);
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (0.660,1.563)–(5.264,1.563);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (0.660,1.563)–(0.840,1.563);
[gp path] (5.264,1.563)–(5.084,1.563);
[gp node right] at (0.540,1.563) $10$;
[gp path] (0.660,2.094)–(0.750,2.094);
[gp path] (5.264,2.094)–(5.174,2.094);
[gp path] (0.660,2.405)–(0.750,2.405);
[gp path] (5.264,2.405)–(5.174,2.405);
[gp path] (0.660,2.625)–(0.750,2.625);
[gp path] (5.264,2.625)–(5.174,2.625);
[gp path] (0.660,2.796)–(0.750,2.796);
[gp path] (5.264,2.796)–(5.174,2.796);
[gp path] (0.660,2.936)–(0.750,2.936);
[gp path] (5.264,2.936)–(5.174,2.936);
[gp path] (0.660,3.054)–(0.750,3.054);
[gp path] (5.264,3.054)–(5.174,3.054);
[gp path] (0.660,3.156)–(0.750,3.156);
[gp path] (5.264,3.156)–(5.174,3.156);
[gp path] (0.660,3.247)–(0.750,3.247);
[gp path] (5.264,3.247)–(5.174,3.247);
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (0.660,3.327)–(5.264,3.327);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (0.660,3.327)–(0.840,3.327);
[gp path] (5.264,3.327)–(5.084,3.327);
[gp node right] at (0.540,3.327) $100$;
[gp path] (0.660,3.859)–(0.750,3.859);
[gp path] (5.264,3.859)–(5.174,3.859);
[gp path] (0.660,4.169)–(0.750,4.169);
[gp path] (5.264,4.169)–(5.174,4.169);
[gp path] (0.660,4.390)–(0.750,4.390);
[gp path] (5.264,4.390)–(5.174,4.390);
[gp path] (0.660,4.561)–(0.750,4.561);
[gp path] (5.264,4.561)–(5.174,4.561);
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (0.660,0.640)–(0.660,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (0.660,0.640)–(0.660,0.820);
[gp path] (0.660,4.561)–(0.660,4.381);
[gp node center] at (0.660,0.440) 1;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (1.235,0.640)–(1.235,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (1.235,0.640)–(1.235,0.820);
[gp path] (1.235,4.561)–(1.235,4.381);
[gp node center] at (1.235,0.440) 2;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (1.811,0.640)–(1.811,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (1.811,0.640)–(1.811,0.820);
[gp path] (1.811,4.561)–(1.811,4.381);
[gp node center] at (1.811,0.440) 4;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (2.386,0.640)–(2.386,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (2.386,0.640)–(2.386,0.820);
[gp path] (2.386,4.561)–(2.386,4.381);
[gp node center] at (2.386,0.440) 8;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (2.962,0.640)–(2.962,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (2.962,0.640)–(2.962,0.820);
[gp path] (2.962,4.561)–(2.962,4.381);
[gp node center] at (2.962,0.440) 16;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (3.537,0.640)–(3.537,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (3.537,0.640)–(3.537,0.820);
[gp path] (3.537,4.561)–(3.537,4.381);
[gp node center] at (3.537,0.440) 32;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (4.113,0.640)–(4.113,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (4.113,0.640)–(4.113,0.820);
[gp path] (4.113,4.561)–(4.113,4.381);
[gp node center] at (4.113,0.440) 64;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (4.688,0.640)–(4.688,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (4.688,0.640)–(4.688,0.820);
[gp path] (4.688,4.561)–(4.688,4.381);
[gp node center] at (4.688,0.440) 128;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (5.264,0.640)–(5.264,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (5.264,0.640)–(5.264,0.820);
[gp path] (5.264,4.561)–(5.264,4.381);
[gp node center] at (5.264,0.440) 256;
[gp path] (0.660,4.561)–(0.660,0.640)–(5.264,0.640)–(5.264,4.561)–cycle;
[gp node center] at (2.962,0.140) number of cores;
rgb color=0.000,1.000,1.000
[gp path] (0.660,4.546)–(1.235,4.047)–(1.811,3.742)–(2.386,3.196)–(2.962,2.851)
3pointgp mark 1(0.660,4.546)
3pointgp mark 1(1.235,4.047)
3pointgp mark 1(1.811,3.742)
3pointgp mark 1(2.386,3.196)
3pointgp mark 1(2.962,2.851)
3pointgp mark 1(3.537,2.459)
3pointgp mark 1(4.113,2.295)
3pointgp mark 1(4.688,2.209)
3pointgp mark 1(5.264,2.335)
rgb color=0.000,0.000,0.545
[gp path] (0.660,4.267)–(1.235,3.799)–(1.811,3.332)–(2.386,2.937)–(2.962,2.590)
3pointgp mark 2(0.660,4.267)
3pointgp mark 2(1.235,3.799)
3pointgp mark 2(1.811,3.332)
3pointgp mark 2(2.386,2.937)
3pointgp mark 2(2.962,2.590)
3pointgp mark 2(3.537,2.302)
3pointgp mark 2(4.113,2.172)
3pointgp mark 2(4.688,2.158)
3pointgp mark 2(5.264,2.313)
rgb color=0.753,0.251,0.000
[gp path] (0.660,4.546)–(1.235,4.194)–(1.811,3.842)–(2.386,3.712)–(2.962,3.348)
3pointgp mark 15(0.660,4.546)
3pointgp mark 15(1.235,4.194)
3pointgp mark 15(1.811,3.842)
3pointgp mark 15(2.386,3.712)
3pointgp mark 15(2.962,3.348)
3pointgp mark 15(3.537,2.816)
3pointgp mark 15(4.113,2.689)
3pointgp mark 15(4.688,2.164)
3pointgp mark 15(5.264,1.778)
rgb color=0.580,0.000,0.827
[gp path] (0.660,4.348)–(1.235,3.868)–(1.811,3.368)–(2.386,3.066)–(2.962,2.752)
3pointgp mark 4(0.660,4.348)
3pointgp mark 4(1.235,3.868)
3pointgp mark 4(1.811,3.368)
3pointgp mark 4(2.386,3.066)
3pointgp mark 4(2.962,2.752)
3pointgp mark 4(3.537,2.482)
3pointgp mark 4(4.113,2.366)
3pointgp mark 4(4.688,1.902)
rgb color=0.933,0.510,0.933
[gp path] (0.660,4.331)–(1.235,3.844)–(1.811,3.362)–(2.386,3.067)–(2.962,2.752)
3pointgp mark 5(0.660,4.331)
3pointgp mark 5(1.235,3.844)
3pointgp mark 5(1.811,3.362)
3pointgp mark 5(2.386,3.067)
3pointgp mark 5(2.962,2.752)
3pointgp mark 5(3.537,2.483)
3pointgp mark 5(4.113,2.356)
3pointgp mark 5(4.688,1.915)
color=gp lt color border
[gp node right] at (1.980,4.694) ts from L1;
rgb color=0.784,0.784,0.000
[gp path] (2.100,4.694)–(2.760,4.694);
[gp path] (0.660,4.135)–(1.235,3.628)–(1.811,3.090)–(2.386,2.724)–(2.962,2.274)
3pointgp mark 12(0.660,4.135)
3pointgp mark 12(1.235,3.628)
3pointgp mark 12(1.811,3.090)
3pointgp mark 12(2.386,2.724)
3pointgp mark 12(2.962,2.274)
3pointgp mark 12(3.537,1.904)
3pointgp mark 12(4.113,1.533)
3pointgp mark 12(4.688,1.207)
3pointgp mark 12(5.264,0.941)
3pointgp mark 12(2.430,4.694)
color=gp lt color border
[gp node right] at (4.200,4.694) tsi from L1;
rgb color=1.000,1.000,0.000
[gp path] (4.320,4.694)–(4.980,4.694);
[gp path] (0.660,4.056)–(1.235,3.541)–(1.811,3.056)–(2.386,2.624)–(2.962,2.227)
3pointgp mark 13(0.660,4.056)
3pointgp mark 13(1.235,3.541)
3pointgp mark 13(1.811,3.056)
3pointgp mark 13(2.386,2.624)
3pointgp mark 13(2.962,2.227)
3pointgp mark 13(3.537,1.866)
3pointgp mark 13(4.113,1.470)
3pointgp mark 13(4.688,1.150)
3pointgp mark 13(5.264,0.904)
3pointgp mark 13(4.650,4.694)
rgb color=0.000,0.000,0.000
[gp path] (0.660,4.056)–(0.701,4.019)–(0.741,3.981)–(0.782,3.944)–(0.823,3.906)
color=gp lt color border
[gp path] (0.660,4.561)–(0.660,0.640)–(5.264,0.640)–(5.264,4.561)–cycle;

gp plot 10.660cm0.640cm5.264cm4.561cm

every node/.append style=scale=0.65
(0.000,0.000) rectangle (5.625,4.987);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (0.980,0.640)–(1.070,0.640);
[gp path] (5.264,0.640)–(5.174,0.640);
[gp path] (0.980,0.726)–(1.070,0.726);
[gp path] (5.264,0.726)–(5.174,0.726);
[gp path] (0.980,0.803)–(1.070,0.803);
[gp path] (5.264,0.803)–(5.174,0.803);
[gp path] (0.980,0.871)–(1.160,0.871);
[gp path] (5.264,0.871)–(5.084,0.871);
[gp node right] at (0.860,0.871) $10$;
[gp path] (0.980,1.319)–(1.070,1.319);
[gp path] (5.264,1.319)–(5.174,1.319);
[gp path] (0.980,1.582)–(1.070,1.582);
[gp path] (5.264,1.582)–(5.174,1.582);
[gp path] (0.980,1.768)–(1.070,1.768);
[gp path] (5.264,1.768)–(5.174,1.768);
[gp path] (0.980,1.912)–(1.070,1.912);
[gp path] (5.264,1.912)–(5.174,1.912);
[gp path] (0.980,2.030)–(1.070,2.030);
[gp path] (5.264,2.030)–(5.174,2.030);
[gp path] (0.980,2.130)–(1.070,2.130);
[gp path] (5.264,2.130)–(5.174,2.130);
[gp path] (0.980,2.216)–(1.070,2.216);
[gp path] (5.264,2.216)–(5.174,2.216);
[gp path] (0.980,2.292)–(1.070,2.292);
[gp path] (5.264,2.292)–(5.174,2.292);
[gp path] (0.980,2.360)–(1.160,2.360);
[gp path] (5.264,2.360)–(5.084,2.360);
[gp node right] at (0.860,2.360) $100$;
[gp path] (0.980,2.809)–(1.070,2.809);
[gp path] (5.264,2.809)–(5.174,2.809);
[gp path] (0.980,3.071)–(1.070,3.071);
[gp path] (5.264,3.071)–(5.174,3.071);
[gp path] (0.980,3.257)–(1.070,3.257);
[gp path] (5.264,3.257)–(5.174,3.257);
[gp path] (0.980,3.402)–(1.070,3.402);
[gp path] (5.264,3.402)–(5.174,3.402);
[gp path] (0.980,3.520)–(1.070,3.520);
[gp path] (5.264,3.520)–(5.174,3.520);
[gp path] (0.980,3.619)–(1.070,3.619);
[gp path] (5.264,3.619)–(5.174,3.619);
[gp path] (0.980,3.706)–(1.070,3.706);
[gp path] (5.264,3.706)–(5.174,3.706);
[gp path] (0.980,3.782)–(1.070,3.782);
[gp path] (5.264,3.782)–(5.174,3.782);
[gp path] (0.980,3.850)–(1.160,3.850);
[gp path] (5.264,3.850)–(5.084,3.850);
[gp node right] at (0.860,3.850) $1000$;
[gp path] (0.980,4.299)–(1.070,4.299);
[gp path] (5.264,4.299)–(5.174,4.299);
[gp path] (0.980,4.561)–(1.070,4.561);
[gp path] (5.264,4.561)–(5.174,4.561);
[gp path] (0.980,0.640)–(0.980,0.820);
[gp path] (0.980,4.561)–(0.980,4.381);
[gp node center,font=8.0pt9.6pt] at (0.980,0.440) 2;
[gp path] (1.408,0.640)–(1.408,0.820);
[gp path] (1.408,4.561)–(1.408,4.381);
[gp node center,font=8.0pt9.6pt] at (1.408,0.440) 4;
[gp path] (1.837,0.640)–(1.837,0.820);
[gp path] (1.837,4.561)–(1.837,4.381);
[gp node center,font=8.0pt9.6pt] at (1.837,0.440) 8;
[gp path] (2.265,0.640)–(2.265,0.820);
[gp path] (2.265,4.561)–(2.265,4.381);
[gp node center,font=8.0pt9.6pt] at (2.265,0.440) 16;
[gp path] (2.694,0.640)–(2.694,0.820);
[gp path] (2.694,4.561)–(2.694,4.381);
[gp node center,font=8.0pt9.6pt] at (2.694,0.440) 32;
[gp path] (3.122,0.640)–(3.122,0.820);
[gp path] (3.122,4.561)–(3.122,4.381);
[gp node center,font=8.0pt9.6pt] at (3.122,0.440) 64;
[gp path] (3.550,0.640)–(3.550,0.820);
[gp path] (3.550,4.561)–(3.550,4.381);
[gp node center,font=8.0pt9.6pt] at (3.550,0.440) 128;
[gp path] (3.979,0.640)–(3.979,0.820);
[gp path] (3.979,4.561)–(3.979,4.381);
[gp node center,font=8.0pt9.6pt] at (3.979,0.440) 256;
[gp path] (4.407,0.640)–(4.407,0.820);
[gp path] (4.407,4.561)–(4.407,4.381);
[gp node center,font=8.0pt9.6pt] at (4.407,0.440) 512;
[gp path] (4.836,0.640)–(4.836,0.820);
[gp path] (4.836,4.561)–(4.836,4.381);
[gp node center,font=8.0pt9.6pt] at (4.836,0.440) 1024;
[gp path] (5.264,0.640)–(5.264,0.820);
[gp path] (5.264,4.561)–(5.264,4.381);
[gp node center,font=8.0pt9.6pt] at (5.264,0.440) 2048;
[gp path] (0.980,4.561)–(0.980,0.640)–(5.264,0.640)–(5.264,4.561)–cycle;
[gp node center,rotate=-270] at (0.206,2.600) elapsed time in s;
[gp node center] at (3.122,0.086) number of cores;
rgb color=0.000,1.000,1.000
[gp path] (2.265,4.433)–(2.694,4.096)–(3.122,3.822)–(3.550,3.518)–(3.979,3.303)
3pointgp mark 1(2.265,4.433)
3pointgp mark 1(2.694,4.096)
3pointgp mark 1(3.122,3.822)
3pointgp mark 1(3.550,3.518)
3pointgp mark 1(3.979,3.303)
3pointgp mark 1(4.407,3.157)
rgb color=0.000,0.000,0.545
[gp path] (1.408,4.202)–(1.837,3.849)–(2.265,3.502)–(2.694,3.219)–(3.122,3.081)
3pointgp mark 2(1.408,4.202)
3pointgp mark 2(1.837,3.849)
3pointgp mark 2(2.265,3.502)
3pointgp mark 2(2.694,3.219)
3pointgp mark 2(3.122,3.081)
3pointgp mark 2(3.550,2.961)
3pointgp mark 2(3.979,2.844)
3pointgp mark 2(4.407,2.806)
rgb color=0.753,0.251,0.000
[gp path] (2.694,3.885)–(3.122,3.480)–(3.550,2.927)–(3.979,2.615)–(4.407,2.275)
3pointgp mark 15(2.694,3.885)
3pointgp mark 15(3.122,3.480)
3pointgp mark 15(3.550,2.927)
3pointgp mark 15(3.979,2.615)
3pointgp mark 15(4.407,2.275)
3pointgp mark 15(4.836,1.948)
3pointgp mark 15(5.264,1.772)
rgb color=0.784,0.784,0.000
[gp path] (1.408,3.868)–(1.837,3.472)–(2.265,3.108)–(2.694,2.782)–(3.122,2.437)
3pointgp mark 12(1.408,3.868)
3pointgp mark 12(1.837,3.472)
3pointgp mark 12(2.265,3.108)
3pointgp mark 12(2.694,2.782)
3pointgp mark 12(3.122,2.437)
3pointgp mark 12(3.550,2.137)
3pointgp mark 12(3.979,1.891)
3pointgp mark 12(4.407,1.802)
rgb color=1.000,1.000,0.000
[gp path] (1.408,3.868)–(1.837,3.466)–(2.265,3.103)–(2.694,2.780)–(3.122,2.416)
3pointgp mark 13(1.408,3.868)
3pointgp mark 13(1.837,3.466)
3pointgp mark 13(2.265,3.103)
3pointgp mark 13(2.694,2.780)
3pointgp mark 13(3.122,2.416)
3pointgp mark 13(3.550,2.128)
3pointgp mark 13(3.979,1.875)
3pointgp mark 13(4.407,1.813)
color=gp lt color border
[gp node right] at (2.300,4.694) ts from L2;
rgb color=0.000,0.392,0.000
[gp path] (2.420,4.694)–(3.080,4.694);
[gp path] (0.980,4.096)–(1.408,3.645)–(1.837,3.181)–(2.265,2.759)–(2.694,2.450)
3pointgp mark 8(0.980,4.096)
3pointgp mark 8(1.408,3.645)
3pointgp mark 8(1.837,3.181)
3pointgp mark 8(2.265,2.759)
3pointgp mark 8(2.694,2.450)
3pointgp mark 8(3.122,2.076)
3pointgp mark 8(3.550,1.740)
3pointgp mark 8(3.979,1.426)
3pointgp mark 8(4.407,1.201)
3pointgp mark 8(4.836,1.046)
3pointgp mark 8(5.264,1.015)
3pointgp mark 8(2.750,4.694)
color=gp lt color border
[gp node right] at (4.520,4.694) tsi from L2;
rgb color=0.000,1.000,0.000
[gp path] (4.640,4.694)–(5.300,4.694);
[gp path] (0.980,4.046)–(1.408,3.558)–(1.837,3.108)–(2.265,2.707)–(2.694,2.372)
3pointgp mark 9(0.980,4.046)
3pointgp mark 9(1.408,3.558)
3pointgp mark 9(1.837,3.108)
3pointgp mark 9(2.265,2.707)
3pointgp mark 9(2.694,2.372)
3pointgp mark 9(3.122,1.969)
3pointgp mark 9(3.550,1.608)
3pointgp mark 9(3.979,1.271)
3pointgp mark 9(4.407,0.980)
3pointgp mark 9(4.836,0.738)
3pointgp mark 9(5.264,0.764)
3pointgp mark 9(4.970,4.694)
rgb color=0.545,0.000,0.000
[gp path] (1.408,4.070)–(1.837,3.657)–(2.265,3.250)–(2.694,2.851)–(3.122,2.569)
3pointgp mark 10(1.408,4.070)
3pointgp mark 10(1.837,3.657)
3pointgp mark 10(2.265,3.250)
3pointgp mark 10(2.694,2.851)
3pointgp mark 10(3.122,2.569)
3pointgp mark 10(3.550,2.315)
3pointgp mark 10(3.979,2.164)
3pointgp mark 10(4.407,2.109)
3pointgp mark 10(4.836,2.057)
3pointgp mark 10(5.264,2.169)
rgb color=1.000,0.000,0.000
[gp path] (1.408,3.784)–(1.837,3.357)–(2.265,2.955)–(2.694,2.535)–(3.122,2.178)
3pointgp mark 11(1.408,3.784)
3pointgp mark 11(1.837,3.357)
3pointgp mark 11(2.265,2.955)
3pointgp mark 11(2.694,2.535)
3pointgp mark 11(3.122,2.178)
3pointgp mark 11(3.550,1.899)
3pointgp mark 11(3.979,1.550)
3pointgp mark 11(4.407,1.327)
3pointgp mark 11(4.836,1.180)
3pointgp mark 11(5.264,1.161)
rgb color=0.000,0.000,0.000
[gp path] (0.980,4.046)–(1.015,4.010)–(1.049,3.974)–(1.084,3.937)–(1.118,3.901)
color=gp lt color border
[gp path] (0.980,4.561)–(0.980,0.640)–(5.264,0.640)–(5.264,4.561)–cycle;

gp plot 10.980cm0.640cm5.264cm4.561cm

every node/.append style=scale=0.65
(0.000,0.000) rectangle (5.625,4.987);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (0.780,0.640)–(0.870,0.640);
[gp path] (5.264,0.640)–(5.174,0.640);
[gp path] (0.780,0.985)–(0.870,0.985);
[gp path] (5.264,0.985)–(5.174,0.985);
[gp path] (0.780,1.230)–(0.870,1.230);
[gp path] (5.264,1.230)–(5.174,1.230);
[gp path] (0.780,1.420)–(0.870,1.420);
[gp path] (5.264,1.420)–(5.174,1.420);
[gp path] (0.780,1.575)–(0.870,1.575);
[gp path] (5.264,1.575)–(5.174,1.575);
[gp path] (0.780,1.707)–(0.870,1.707);
[gp path] (5.264,1.707)–(5.174,1.707);
[gp path] (0.780,1.820)–(0.870,1.820);
[gp path] (5.264,1.820)–(5.174,1.820);
[gp path] (0.780,1.921)–(0.870,1.921);
[gp path] (5.264,1.921)–(5.174,1.921);
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (0.780,2.010)–(5.264,2.010);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (0.780,2.010)–(0.960,2.010);
[gp path] (5.264,2.010)–(5.084,2.010);
[gp node right] at (0.660,2.010) $100$;
[gp path] (0.780,2.601)–(0.870,2.601);
[gp path] (5.264,2.601)–(5.174,2.601);
[gp path] (0.780,2.946)–(0.870,2.946);
[gp path] (5.264,2.946)–(5.174,2.946);
[gp path] (0.780,3.191)–(0.870,3.191);
[gp path] (5.264,3.191)–(5.174,3.191);
[gp path] (0.780,3.381)–(0.870,3.381);
[gp path] (5.264,3.381)–(5.174,3.381);
[gp path] (0.780,3.536)–(0.870,3.536);
[gp path] (5.264,3.536)–(5.174,3.536);
[gp path] (0.780,3.667)–(0.870,3.667);
[gp path] (5.264,3.667)–(5.174,3.667);
[gp path] (0.780,3.781)–(0.870,3.781);
[gp path] (5.264,3.781)–(5.174,3.781);
[gp path] (0.780,3.881)–(0.870,3.881);
[gp path] (5.264,3.881)–(5.174,3.881);
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (0.780,3.971)–(5.264,3.971);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (0.780,3.971)–(0.960,3.971);
[gp path] (5.264,3.971)–(5.084,3.971);
[gp node right] at (0.660,3.971) $1000$;
[gp path] (0.780,4.561)–(0.870,4.561);
[gp path] (5.264,4.561)–(5.174,4.561);
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (0.780,0.640)–(0.780,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (0.780,0.640)–(0.780,0.820);
[gp path] (0.780,4.561)–(0.780,4.381);
[gp node center,font=8.0pt9.6pt] at (0.780,0.440) 16;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (1.356,0.640)–(1.356,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (1.356,0.640)–(1.356,0.820);
[gp path] (1.356,4.561)–(1.356,4.381);
[gp node center,font=8.0pt9.6pt] at (1.356,0.440) 32;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (1.932,0.640)–(1.932,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (1.932,0.640)–(1.932,0.820);
[gp path] (1.932,4.561)–(1.932,4.381);
[gp node center,font=8.0pt9.6pt] at (1.932,0.440) 64;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (2.508,0.640)–(2.508,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (2.508,0.640)–(2.508,0.820);
[gp path] (2.508,4.561)–(2.508,4.381);
[gp node center,font=8.0pt9.6pt] at (2.508,0.440) 128;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (3.084,0.640)–(3.084,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (3.084,0.640)–(3.084,0.820);
[gp path] (3.084,4.561)–(3.084,4.381);
[gp node center,font=8.0pt9.6pt] at (3.084,0.440) 256;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (3.660,0.640)–(3.660,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (3.660,0.640)–(3.660,0.820);
[gp path] (3.660,4.561)–(3.660,4.381);
[gp node center,font=8.0pt9.6pt] at (3.660,0.440) 512;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (4.236,0.640)–(4.236,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (4.236,0.640)–(4.236,0.820);
[gp path] (4.236,4.561)–(4.236,4.381);
[gp node center,font=8.0pt9.6pt] at (4.236,0.440) 1024;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (4.812,0.640)–(4.812,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (4.812,0.640)–(4.812,0.820);
[gp path] (4.812,4.561)–(4.812,4.381);
[gp node center,font=8.0pt9.6pt] at (4.812,0.440) 2048;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (5.264,0.640)–(5.264,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (5.264,0.640)–(5.264,0.820);
[gp path] (5.264,4.561)–(5.264,4.381);
[gp node center,font=8.0pt9.6pt] at (5.264,0.440) 3528;
[gp path] (0.780,4.561)–(0.780,0.640)–(5.264,0.640)–(5.264,4.561)–cycle;
[gp node center] at (3.022,0.086) number of cores;
rgb color=0.753,0.251,0.000
[gp path] (3.660,4.182)–(4.236,3.692)–(4.812,3.292)–(5.264,2.893);
3pointgp mark 15(3.660,4.182)
3pointgp mark 15(4.236,3.692)
3pointgp mark 15(4.812,3.292)
3pointgp mark 15(5.264,2.893)
color=gp lt color border
[gp node right] at (2.100,4.694) ts from L3;
rgb color=0.545,0.000,0.000
[gp path] (2.220,4.694)–(2.880,4.694);
[gp path] (0.780,4.288)–(1.356,3.709)–(1.932,3.232)–(2.508,2.766)–(3.084,2.387)
3pointgp mark 10(0.780,4.288)
3pointgp mark 10(1.356,3.709)
3pointgp mark 10(1.932,3.232)
3pointgp mark 10(2.508,2.766)
3pointgp mark 10(3.084,2.387)
3pointgp mark 10(3.660,2.084)
3pointgp mark 10(4.236,1.903)
3pointgp mark 10(4.812,1.936)
3pointgp mark 10(5.264,2.062)
3pointgp mark 10(2.550,4.694)
color=gp lt color border
[gp node right] at (4.320,4.694) tsi from L3;
rgb color=1.000,0.000,0.000
[gp path] (4.440,4.694)–(5.100,4.694);
[gp path] (0.780,4.204)–(1.356,3.614)–(1.932,3.091)–(2.508,2.591)–(3.084,2.092)
3pointgp mark 11(0.780,4.204)
3pointgp mark 11(1.356,3.614)
3pointgp mark 11(1.932,3.091)
3pointgp mark 11(2.508,2.591)
3pointgp mark 11(3.084,2.092)
3pointgp mark 11(3.660,1.673)
3pointgp mark 11(4.236,1.293)
3pointgp mark 11(4.812,0.973)
3pointgp mark 11(5.264,0.901)
3pointgp mark 11(4.770,4.694)
rgb color=0.000,1.000,0.000
[gp path] (3.660,2.343)–(4.236,1.999)–(4.812,1.651);
3pointgp mark 9(3.660,2.343)
3pointgp mark 9(4.236,1.999)
3pointgp mark 9(4.812,1.651)
rgb color=0.000,0.000,0.000
[gp path] (0.780,4.204)–(0.815,4.168)–(0.850,4.132)–(0.885,4.096)–(0.920,4.061)
color=gp lt color border
[gp path] (0.780,4.561)–(0.780,0.640)–(5.264,0.640)–(5.264,4.561)–cycle;

gp plot 10.780cm0.640cm5.264cm4.561cm

every node/.append style=scale=0.65
(0.000,0.000) rectangle (5.625,4.987);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (0.780,0.928)–(0.870,0.928);
[gp path] (5.264,0.928)–(5.174,0.928);
[gp path] (0.780,1.382)–(0.870,1.382);
[gp path] (5.264,1.382)–(5.174,1.382);
[gp path] (0.780,1.734)–(0.870,1.734);
[gp path] (5.264,1.734)–(5.174,1.734);
[gp path] (0.780,2.021)–(0.870,2.021);
[gp path] (5.264,2.021)–(5.174,2.021);
[gp path] (0.780,2.265)–(0.870,2.265);
[gp path] (5.264,2.265)–(5.174,2.265);
[gp path] (0.780,2.475)–(0.870,2.475);
[gp path] (5.264,2.475)–(5.174,2.475);
[gp path] (0.780,2.661)–(0.870,2.661);
[gp path] (5.264,2.661)–(5.174,2.661);
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (0.780,2.827)–(5.264,2.827);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (0.780,2.827)–(0.960,2.827);
[gp path] (5.264,2.827)–(5.084,2.827);
[gp node right] at (0.660,2.827) $1000$;
[gp path] (0.780,3.921)–(0.870,3.921);
[gp path] (5.264,3.921)–(5.174,3.921);
[gp path] (0.780,4.561)–(0.870,4.561);
[gp path] (5.264,4.561)–(5.174,4.561);
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (0.780,0.640)–(0.780,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (0.780,0.640)–(0.780,0.820);
[gp path] (0.780,4.561)–(0.780,4.381);
[gp node center] at (0.780,0.440) 128;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (1.719,0.640)–(1.719,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (1.719,0.640)–(1.719,0.820);
[gp path] (1.719,4.561)–(1.719,4.381);
[gp node center] at (1.719,0.440) 256;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (2.658,0.640)–(2.658,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (2.658,0.640)–(2.658,0.820);
[gp path] (2.658,4.561)–(2.658,4.381);
[gp node center] at (2.658,0.440) 512;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (3.597,0.640)–(3.597,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (3.597,0.640)–(3.597,0.820);
[gp path] (3.597,4.561)–(3.597,4.381);
[gp node center] at (3.597,0.440) 1024;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (4.536,0.640)–(4.536,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (4.536,0.640)–(4.536,0.820);
[gp path] (4.536,4.561)–(4.536,4.381);
[gp node center] at (4.536,0.440) 2048;
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (5.264,0.640)–(5.264,4.561);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (5.264,0.640)–(5.264,0.820);
[gp path] (5.264,4.561)–(5.264,4.381);
[gp node center] at (5.264,0.440) 3504;
[gp path] (0.780,4.561)–(0.780,0.640)–(5.264,0.640)–(5.264,4.561)–cycle;
[gp node center] at (3.022,0.140) number of cores;
[gp node right] at (2.100,4.694) ts from L4;
rgb color=0.000,0.000,0.545
[gp path] (2.220,4.694)–(2.880,4.694);
[gp path] (3.597,4.303)–(4.536,3.879);
3pointgp mark 6(3.597,4.303)
3pointgp mark 6(4.536,3.879)
3pointgp mark 6(2.550,4.694)
color=gp lt color border
[gp node right] at (4.320,4.694) tsi from L4;
rgb color=0.000,0.000,1.000
[gp path] (4.440,4.694)–(5.100,4.694);
[gp path] (0.780,3.891)–(1.719,2.956)–(2.658,2.271)–(3.597,1.755)–(4.536,1.353)
3pointgp mark 7(0.780,3.891)
3pointgp mark 7(1.719,2.956)
3pointgp mark 7(2.658,2.271)
3pointgp mark 7(3.597,1.755)
3pointgp mark 7(4.536,1.353)
3pointgp mark 7(5.264,1.013)
3pointgp mark 7(4.770,4.694)
rgb color=0.000,0.000,0.000
[gp path] (0.780,3.891)–(0.818,3.847)–(0.856,3.803)–(0.894,3.758)–(0.932,3.714)
color=gp lt color border
[gp path] (0.780,4.561)–(0.780,0.640)–(5.264,0.640)–(5.264,4.561)–cycle;

gp plot 10.780cm0.640cm5.264cm4.561cm

Cubic field: elapsed time in second versus the number of cores. Log/log scale. The results for a fine-scale refinement at  L2 , L3, L4, L5, L6  and L7 level are presented  in
*UOL_curve_times_L6 and *UOL_curve_times_L7 respectively.
The ideal curve shows the  slope when the elapsed time of a single process is perfectly divided by the number of processes used for the calculation (this corresponds to the ideal speed-up).
This curve is always shifted so that it starts at the same point as the best solution for the smallest number of cores.

Figure <ref> presents, for different SP refinement, an analysis of the elapsed time used by all solvers.
The points not present on these curves are mainly due to the fact that, with this cluster, it was impossible to compute them: not enough memory, too much computation time, too many processes for the solver or integer overflow.
These limitations only apply to the actual cluster installation and are in no way an absolute judgment on the software used, especially in a non-multi-threaded context.

As discussed  in section <ref>, the solver has an ideal scale ratio.
Thus, for the refinement of levels L3, L4, L5 and L6 different coarse-scale levels were tested.
The observed  optimal scale  ratio  is 3 (best elapsed times curves are "tsi from L1", "tsi from L2" and  "tsi from L3" respectively for the refinement of levels L4 (<ref>), L5 (<ref>) and L6 (<ref>): 1+3=4, 2+3=5 and 3+3=6), which is in good agreement with the analysis in section <ref>  and gives  equation  (<ref>).
For L7, only this ratio was used to save computing resources.
Regarding the number of dofs, the figure <ref> shows that this ratio corresponds to a fixed number of dofs per patch appearing visually as a horizontal line (circle between 3435 and 4928).

every node/.append style=scale=0.60
(0.000,0.000) rectangle (8.125,5.250);
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (1.010,0.592)–(7.794,0.592);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (1.010,0.592)–(1.190,0.592);
[gp path] (7.794,0.592)–(7.614,0.592);
[gp node right] at (0.900,0.592) 1e+02;
[gp path] (1.010,0.746)–(1.100,0.746);
[gp path] (7.794,0.746)–(7.704,0.746);
[gp path] (1.010,0.835)–(1.100,0.835);
[gp path] (7.794,0.835)–(7.704,0.835);
[gp path] (1.010,0.899)–(1.100,0.899);
[gp path] (7.794,0.899)–(7.704,0.899);
[gp path] (1.010,0.949)–(1.100,0.949);
[gp path] (7.794,0.949)–(7.704,0.949);
[gp path] (1.010,0.989)–(1.100,0.989);
[gp path] (7.794,0.989)–(7.704,0.989);
[gp path] (1.010,1.023)–(1.100,1.023);
[gp path] (7.794,1.023)–(7.704,1.023);
[gp path] (1.010,1.053)–(1.100,1.053);
[gp path] (7.794,1.053)–(7.704,1.053);
[gp path] (1.010,1.079)–(1.100,1.079);
[gp path] (7.794,1.079)–(7.704,1.079);
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (1.010,1.102)–(7.794,1.102);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (1.010,1.102)–(1.190,1.102);
[gp path] (7.794,1.102)–(7.614,1.102);
[gp node right] at (0.900,1.102) 1e+03;
[gp path] (1.010,1.256)–(1.100,1.256);
[gp path] (7.794,1.256)–(7.704,1.256);
[gp path] (1.010,1.346)–(1.100,1.346);
[gp path] (7.794,1.346)–(7.704,1.346);
[gp path] (1.010,1.410)–(1.100,1.410);
[gp path] (7.794,1.410)–(7.704,1.410);
[gp path] (1.010,1.459)–(1.100,1.459);
[gp path] (7.794,1.459)–(7.704,1.459);
[gp path] (1.010,1.499)–(1.100,1.499);
[gp path] (7.794,1.499)–(7.704,1.499);
[gp path] (1.010,1.534)–(1.100,1.534);
[gp path] (7.794,1.534)–(7.704,1.534);
[gp path] (1.010,1.563)–(1.100,1.563);
[gp path] (7.794,1.563)–(7.704,1.563);
[gp path] (1.010,1.589)–(1.100,1.589);
[gp path] (7.794,1.589)–(7.704,1.589);
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (1.010,1.613)–(7.794,1.613);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (1.010,1.613)–(1.190,1.613);
[gp path] (7.794,1.613)–(7.614,1.613);
[gp node right] at (0.900,1.613) 1e+04;
[gp path] (1.010,1.766)–(1.100,1.766);
[gp path] (7.794,1.766)–(7.704,1.766);
[gp path] (1.010,1.856)–(1.100,1.856);
[gp path] (7.794,1.856)–(7.704,1.856);
[gp path] (1.010,1.920)–(1.100,1.920);
[gp path] (7.794,1.920)–(7.704,1.920);
[gp path] (1.010,1.969)–(1.100,1.969);
[gp path] (7.794,1.969)–(7.704,1.969);
[gp path] (1.010,2.010)–(1.100,2.010);
[gp path] (7.794,2.010)–(7.704,2.010);
[gp path] (1.010,2.044)–(1.100,2.044);
[gp path] (7.794,2.044)–(7.704,2.044);
[gp path] (1.010,2.073)–(1.100,2.073);
[gp path] (7.794,2.073)–(7.704,2.073);
[gp path] (1.010,2.100)–(1.100,2.100);
[gp path] (7.794,2.100)–(7.704,2.100);
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (1.010,2.123)–(7.794,2.123);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (1.010,2.123)–(1.190,2.123);
[gp path] (7.794,2.123)–(7.614,2.123);
[gp node right] at (0.900,2.123) 1e+05;
[gp path] (1.010,2.276)–(1.100,2.276);
[gp path] (7.794,2.276)–(7.704,2.276);
[gp path] (1.010,2.366)–(1.100,2.366);
[gp path] (7.794,2.366)–(7.704,2.366);
[gp path] (1.010,2.430)–(1.100,2.430);
[gp path] (7.794,2.430)–(7.704,2.430);
[gp path] (1.010,2.480)–(1.100,2.480);
[gp path] (7.794,2.480)–(7.704,2.480);
[gp path] (1.010,2.520)–(1.100,2.520);
[gp path] (7.794,2.520)–(7.704,2.520);
[gp path] (1.010,2.554)–(1.100,2.554);
[gp path] (7.794,2.554)–(7.704,2.554);
[gp path] (1.010,2.584)–(1.100,2.584);
[gp path] (7.794,2.584)–(7.704,2.584);
[gp path] (1.010,2.610)–(1.100,2.610);
[gp path] (7.794,2.610)–(7.704,2.610);
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (1.010,2.633)–(7.794,2.633);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (1.010,2.633)–(1.190,2.633);
[gp path] (7.794,2.633)–(7.614,2.633);
[gp node right] at (0.900,2.633) 1e+06;
[gp path] (1.010,2.787)–(1.100,2.787);
[gp path] (7.794,2.787)–(7.704,2.787);
[gp path] (1.010,2.877)–(1.100,2.877);
[gp path] (7.794,2.877)–(7.704,2.877);
[gp path] (1.010,2.940)–(1.100,2.940);
[gp path] (7.794,2.940)–(7.704,2.940);
[gp path] (1.010,2.990)–(1.100,2.990);
[gp path] (7.794,2.990)–(7.704,2.990);
[gp path] (1.010,3.030)–(1.100,3.030);
[gp path] (7.794,3.030)–(7.704,3.030);
[gp path] (1.010,3.064)–(1.100,3.064);
[gp path] (7.794,3.064)–(7.704,3.064);
[gp path] (1.010,3.094)–(1.100,3.094);
[gp path] (7.794,3.094)–(7.704,3.094);
[gp path] (1.010,3.120)–(1.100,3.120);
[gp path] (7.794,3.120)–(7.704,3.120);
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (1.010,3.143)–(7.794,3.143);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (1.010,3.143)–(1.190,3.143);
[gp path] (7.794,3.143)–(7.614,3.143);
[gp node right] at (0.900,3.143) 1e+07;
[gp path] (1.010,3.297)–(1.100,3.297);
[gp path] (7.794,3.297)–(7.704,3.297);
[gp path] (1.010,3.387)–(1.100,3.387);
[gp path] (7.794,3.387)–(7.704,3.387);
[gp path] (1.010,3.451)–(1.100,3.451);
[gp path] (7.794,3.451)–(7.704,3.451);
[gp path] (1.010,3.500)–(1.100,3.500);
[gp path] (7.794,3.500)–(7.704,3.500);
[gp path] (1.010,3.541)–(1.100,3.541);
[gp path] (7.794,3.541)–(7.704,3.541);
[gp path] (1.010,3.575)–(1.100,3.575);
[gp path] (7.794,3.575)–(7.704,3.575);
[gp path] (1.010,3.604)–(1.100,3.604);
[gp path] (7.794,3.604)–(7.704,3.604);
[gp path] (1.010,3.630)–(1.100,3.630);
[gp path] (7.794,3.630)–(7.704,3.630);
color=gp lt color axes
gp lt axes
gp dt axes
[gp path] (1.010,3.654)–(7.794,3.654);
color=gp lt color border
gp lt border
gp dt solid
[gp path] (1.010,3.654)–(1.190,3.654);
[gp path] (7.794,3.654)–(7.614,3.654);
[gp node right] at (0.900,3.654) 1e+08;
[gp path] (1.010,3.807)–(1.100,3.807);
[gp path] (7.794,3.807)–(7.704,3.807);
[gp path] (1.010,3.897)–(1.100,3.897);
[gp path] (7.794,3.897)–(7.704,3.897);
[gp path] (1.010,3.961)–(1.100,3.961);
[gp path] (7.794,3.961)–(7.704,3.961);
[gp path] (1.010,4.010)–(1.100,4.010);
[gp path] (7.794,4.010)–(7.704,4.010);
[gp path] (1.010,4.051)–(1.100,4.051);
[gp path] (7.794,4.051)–(7.704,4.051);
[gp path] (1.010,4.085)–(1.100,4.085);
[gp path] (7.794,4.085)–(7.704,4.085);
[gp path] (1.010,4.115)–(1.100,4.115);