# On-the-fly adaptive sampling strategy for data-driven computational mechanics: Applications to computational homogenisation

Felipe Rocha, Auriane Platzer, Adrien Leygue, Laurent Stainier

► To cite this version:

Felipe Rocha, Auriane Platzer, Adrien Leygue, Laurent Stainier. On-the-fly adaptive sampling strategy for data-driven computational mechanics: Applications to computational homogenisation. Mechanics of Materials, 2025, 207, 10.1016/j.mechmat.2025.105382. hal-05113418

HAL Id: hal-05113418

https://hal.science/hal-05113418v1

Submitted on 28 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1 On-the-fly Adaptive Sampling strategy for Data-driven

2 Computational Mechanics: applications to Computational
Homogenisation

4 Felipe Rocha<sup>b, a</sup>, Auriane Platzer<sup>a, c</sup>, Adrien Leygue<sup>a</sup>, Laurent Stainier<sup>a</sup>

5 Abstract

To overcome well-known drawbacks of classical phenomenological constitutive modelling of non-linear heterogeneous materials, a popular approach is the so-called computational homogenisation (CH), which relies on a combined description of constitutive behaviours and spatial morphology of smaller scales constituents. Their computational implementation usually encompasses two-level finite element numerical models (FE²), which quickly leads to prohibitive computational costs, even for problems with a modest number of degrees of freedom. A popular strategy in the literature to alleviate such computational burden is to use machine learning-based surrogate constitutive models, although fundamental drawbacks such as the absence of interpretability, limited extrapolation capability, and limited mathematical analysis, are still present. On the other hand, the so-called (model-free) data-driven computational mechanics (DDCM) paradigm (Kirchdoerfer and Ortiz, 2016) proposes the direct integration of "experimental data", completely bypassing the need for explicit constitutive laws. The main goal of this work is to show how the DDCM approach can be used in synergy with CH to bypass fully-coupled multiscale computations. A naive approach to using multiscale constitutive behaviour along with DDCM is the offline construction of a database via CH by assuming some sampling of the strain-space. This approach has little interest since the region of the strain-space covered during a simulation is problem-dependent. During the iterations of the DDCM solver, a finite set of strain-stress pairs is used as input, while another set, comprising mechanically admissible states, is dynamically generated. Such a scenario naturally unveils the most relevant goal-oriented phase-space instances to guide a material dataset enrichment iterative procedure, bridging DDCM and CH towards computationally effective two-scale simulations. On the other hand, the performance is limited if the quality of the mechanically admissible states is low, particularly when the material

<footnote><sup>a</sup>Nantes Université, Ecole Centrale de Nantes, CNRS, GeM (UMR 6183), 1 Rue de la Noë, Nantes, 44300, France</footnote>

<footnote><sup>b</sup>Univ Paris Est Creteil, Univ Gustave Eiffel, CNRS, UMR 8208, MSME, Créteil, 94010, France</footnote>

<footnote><sup>c</sup>INSA Lyon, CNRS, LaMCoS, UMR5259, Villeurbanne, 69621, France</footnote>
dataset is not insufficiently dense. To address these challenges, we propose meaningful ranking scores alongside numerical techniques to enhance the DDCM solver in sparse data scenarios. As result, we show through meaningful numerical examples that the proposed framework results in significant computational savings if compared to standard FE².

Keywords: (Model-free) Data-Driven Computational Mechanics, Adaptive Sampling, Computational Homogenisation

# 1. Introduction

To overcome some limitations associated to traditional constitutive material modelling, so-called data-driven techniques have become a popular tool in computational mechanics during recent years. Two classes of approaches can be distinguished: i) machine learning-based (ML) surrogate models and ii) model-free approaches, such as the data-driven computational mechanics (DDCM) paradigm introduced in the seminal work of Kirchdoerfer and Ortiz (2016). Another related method is based on manifold learning to approximate the local behaviour of experimental data (Ibañez et al., 2016). For more information concerning the first class of approaches, we refer the reader to the review paper of Linden et al. (2023). A comparative study of ML and DDCM approaches for hyperelastic material response has recently been proposed by Zlatić et al. (2024).

In more details, DDCM can be seen as a method to solve boundary value problems, which does not require explicit constitutive models. The sought solution is constrained to satisfy relevant conservation principles, e.g. balance of linear and angular momentum, while the interplay with experimental material data is ensured by minimising a constitutive distance (Kirchdoerfer and Ortiz, 2016). It is worth mentioning that, in our context, experimental does not only simply mean empirical testing databases, but also synthetically generated databases from virtual testing, via computational homogenisation approaches for example. Recent advances in DDCM has extended its range of application to inelastic materials (Eggersmann et al., 2019; Dandin et al., 2024), noisy data (Kirchdoerfer and Ortiz, 2017), finite strains (Nguyen and Keip, 2018; Conti et al., 2020; Platzer et al., 2021), poroelasticity (Bahmani and Sun, 2021), micromorphic continua (Ulloa et al., 2024), chemomechanics (Zhou et al., 2024), piezoelectricity (Marenić et al., 2022), fracture mechanics (Carrara et al., 2020), material identification (Leygue et al., 2018; Dalémat et al., 2019; Costecalde et al., 2023), uncertainty quantification (Prume et al., 2023), frequency domain data (Salahshoor and Ortiz, 2023), etc.

While DDCM features the outstanding characteristic of relying exclusively on material data, with no interpolation nor extrapolation biases, one can only expect accurate results if the database is rich enough to faithfully
represent the material response. Naturally, such sought characteristic raises issues: i) cost of generating such database (either experimentally or numerically); ii) lack of a priori knowledge of material response range when facing a new material/problem; iii) prohibitively high computational cost for too large databases given the combinatorial complexity of DDCM (Kanno, 2019). It is easy to see that points i) and ii) are conflicting with iii) and vice-versa, such that it is impossible to comply with all the aforementioned points with the standard DDCM procedure. Therefore, it is necessary to obtain a more efficient strategy to build tailored material database, sufficiently dense on useful regions but sparse otherwise. To reach this aim, our main contributions are two-fold:

1. Boosting strategies for DDCM solver when faced with sparse databases, which correspond to cases with little a priori knowledge on the material.

2. On-the-fly adaptive sampling strategies that optimally populate the material manifold leading to goal-oriented databases.

While the second point of contribution is our final target, the first step is essential since inaccurate solutions obtained with sparse databases can lead to non-converging solutions. In this regard, we propose applying four techniques, three of which are novel contributions of this work, to be detailed in due course. Indeed, each of these ingredients have different aims, and are not mutually exclusive, which means that they can be used isolated or combined to enhance numerical performance and robustness.

The proposed enhanced DDCM solver is developed within a particular setting, but the contributions could be applied more broadly. We restrict our analysis to: i) infinitesimal strain (possibly non-linear) elasticity in two-dimensions which allows for a reduced phase space and simpler database generation/exploration and ii) the numerical generation of material database thanks to Computational Homogenisation (CH). While infinitesimal strains in assumption i) is considered only for sake of simplicity (the extension to finite strain is reasonably straightforward), the formulation for inelastic behaviours is not trivial and therefore beyond the scope of the present work. Assumption ii) naturally yields a third derived contribution: a method to accelerate $FE^2$ computations. In this regard, the work of Xu et al. (2020) proposes the use CH computations to build an offline database. Although it is shown that a reduction of the online time is achieved, it should be noted that, given the impossibility of knowing in advance the most relevant regions of the phase-space to be sampled, such method is inherently too expensive due to the necessary offline computations. Offline sampling was also employed by Karapiperis et al. (2020) in a goal-oriented manner to provide valuable data points from a micromechanical model. Following a non-DDCM alternative, the work of Chaouch and Yvonnet (2024) proposes an alternative to the $FE^2$ method by only querying CH simulations on the
barycentric points of K-means clusters computed with the strain states, stress fields being reconstructed by linear extrapolations anchored on the available simulated points. In contrast, our work proposes to use DDCM output combined with adaptive approaches to generate tailored databases in an on-the-fly manner, thus reducing the number of CH queries in a two-scale simulation.

The manuscript is structured in three main parts, each of them containing preliminary numerical examples: i) (Model-free) Data-driven Computational Mechanics (DDCM) is applied along Computational Homogenisation (CH) to set the basis of proposed framework in Section 2; ii) DDCM boosting techniques are introduced and tested in Section 3; iii) finally, on-the-fly adaptive sampling methods are designed and validated in Section 4. A fourth part in Section 5 is devoted to investigate more deeply speed-up aspects through a last unifying numerical example.

**Remark 1 (Notations).** Let us consider the following conventions. The unit vectors {$e_i$}$_{i=1}^d$ refer to the canonical basis for $\mathbb{R}^d$, for $d = 2$ or 3. If not stated otherwise, we admit the font convention $\mathbf{a} \in \mathbb{R}^d$, $\mathbf{A} \in \mathbb{R}^{d,d}$, $\mathbf{A} \in \mathbb{R}^{d,d,d}$ to discern between vectors, second- and fourth-order tensors, respectively. Space $\mathbb{R}_{\text{sym}}^{d,d}$ denotes symmetric real-valued second-order tensors, the operator $\text{Sym}(\cdot)$ stands for the symmetric part of an object, in particular $\nabla^s(\cdot) = \text{Sym}(\nabla(\cdot))$ stands for the gradient symmetric operator. Regular typeface fonts are used to denote components (or scalars) in Cartesian coordinates, e.g., $\mathbf{a} = a_i e_i$, $\mathbf{A} = A_{ij} e_i \otimes e_j$ and so-forth (Einstein's summation convention applied). The dot ($\cdot$) designates inner products between vectors or tensors, e.g., $\mathbf{a} \cdot \mathbf{b} = a_i b_i$, $\mathbf{A} \cdot \mathbf{B} = A_{ij} B_{ij}$, etc. Standard matrix-vector and fourth-second-order tensor multiplication are assumed in the absence of any sign, i.e., $\mathbf{A}\mathbf{b} := A_{ij} b_j e_i$, $\mathbf{A}\mathbf{B} := A_{ijkl} B_{kl} e_i \otimes e_j$. As usual, $L^2(\Omega; \mathbb{K})$ stands for the square-integrable functions with domain on $\Omega \subset \mathbb{R}^d$ and image on $\mathbb{K}$, $\mathbb{K} \in \{\mathbb{R}, \mathbb{R}^d, \mathbb{R}_{\text{sym}}^{d,d}\}$ in our context. The $L^2$-inner product is denoted by $\langle \bullet, \star \rangle_{L^2(\Omega;\mathbb{K})} = \int_{\Omega}(\bullet) \cdot (\star) d\Omega$ (or simply $\langle \bullet, \star \rangle$), the induced norm by $\|(\star)\|_{L^2(\Omega;\mathbb{K})} = [\int_{\Omega}(\star) \cdot (\star) d\Omega]^{1/2}$ (or simply $\|(\star)\|_{L^2}$). Simplified versions being preferred whenever the meaning is clear by the context. Finally, it is also useful to consider the spatial average operator $\langle (\star) \rangle_{\Omega} = \frac{1}{|\Omega|} \int_{\Omega}(\star) d\Omega$.

## 2. DDCM coupled with CH

In order to motivate the two main contributions of this work, we first introduce a coupling of the data-driven method with computational homogenisation to generate databases on the fly, adapted to the macroscopic mechanical response. We take this opportunity to recall the key concepts behind these two tools, and we use an initial naive approach to bring them into dialogue, as this efficiently illustrates the idea but also demonstrates its limitations.
To this end, let us consider a problem (at macroscale) in infinitesimal strain elasticity regime, occupying a domain $\Omega \subset \mathbb{R}^d$. The partition of its boundary $\partial\Omega = \partial\Omega_N \cup \partial\Omega_D$ is such that a traction vector $\bar{t}$ and a prescribed displacement $\bar{u}$ are known, defined on the Neumann ($\partial\Omega_N$) and Dirichlet ($\partial\Omega_D$) boundaries, respectively. From the linear momentum balance principle and kinematical compatibility conditions, the Cauchy stress tensor field $\sigma$ and infinitesimal strain tensor field $\varepsilon$ should respectively satisfy

$$ - \mathrm{div} \, \sigma = \mathbf{b} \quad \text{in } \Omega, \qquad \sigma \mathbf{n} = \bar{\mathbf{t}} \quad \text{on } \partial\Omega_N \quad (1a) $$

$$ \varepsilon = \nabla^s \mathbf{u} \quad \text{in } \Omega, \qquad \mathbf{u} = \bar{\mathbf{u}} \quad \text{on } \partial\Omega_D, \quad (1b) $$

where $\mathbf{u}$ is a sufficiently regular displacement field. Note that we have not explicitly introduced a link between $\sigma$ and $\varepsilon$. Indeed, in the DDCM standard setting (Section 2.1), such a relation, in general non-linear, is replaced by a finite collection of strain-stress pairs forming a dataset which can be obtained, e.g., by CH (Section 2.2). The online coupling strategy between the two scales is presented in Section 2.3, resulting in an iterative DD-CH scheme illustrated in Section 2.4.

## 2.1. Data-driven model-free solver at the macroscale

In this section, we briefly recall the (model-free) DDCM formulation, originally proposed by Kirchdoerfer and Ortiz (2016). Our formulation follows more closely the work of Nguyen et al. (2020), which is more suited for continuum problems in a variational (weak) format. This choice allows for a more straightforward finite element implementation using the FEniCS/FEniCSx library (Alnaes et al., 2015; Baratta et al., 2023), in contrast with the original DDCM version (Kirchdoerfer and Ortiz, 2016) which privileges discretised language. Particularly, our implementation, namely ddfenics, is available in a public repository (Rocha, 2023a).

As previously stated, DDCM is a model-free approach that circumvents the constitutive model, substituting it with a discrete material dataset comprising strain-stress pairs. The formulation then distinguishes two kinds of strain-stress pairs, also referred to as states: (i) mechanically admissible states that verify (1), (ii) material states which select the material response of the body from the material data set. The problem is then formulated as the minimisation of a distance between those two states, which is subsequently detailed within the relevant mathematical framework.

Let $\mathbb{Z} = \mathbb{R}_{sym}^{d,d} \times \mathbb{R}_{sym}^{d,d}$ be the so-called phase-space composed by strain-stress pairs and let $\mathcal{Z} = L^2(\Omega; \mathbb{Z})$ its associated functional space. We adopt the following useful convention to distinguish functions from their point-wise values: $\mathbf{z} = (\varepsilon, \sigma) \in \mathcal{Z}$, while $\hat{\mathbf{z}} = (\hat{\varepsilon}, \hat{\sigma}) = \mathbf{z}(\mathbf{x}) \in \mathbb{Z}$ for some $\mathbf{x} \in \Omega$. The distance between two elements $\mathbf{z}$ and $\mathbf{y}$ in $\mathcal{Z}$ is defined by

$$ d^2(\mathbf{z}, \mathbf{y}) = \int_{\Omega} \|\mathbf{z}(\mathbf{x}) - \mathbf{y}(\mathbf{x})\|_{\mathcal{C}}^2 d\Omega, \quad (2) $$
156 where

$$ \|(\hat{\epsilon}, \hat{\sigma})\|_{\mathbb{C}}^2 = \mathbb{C}\hat{\epsilon} \cdot \hat{\epsilon} + \mathbb{C}^{-1}\hat{\sigma} \cdot \hat{\sigma}, \quad (3) $$

with $\mathbb{C}$ a tensorial algorithmic parameter that defines an energy-like norm, balancing contributions of the strain and stress parts. Finally, the DDCM problem reads as

$$ \min_{z' \in \mathcal{Z}_E} \min_{z^* \in \mathcal{Z}_D} d(z', z^*), \quad (4) $$

where $\mathcal{Z}_E = \{z \in \mathcal{Z}; (1a) \text{ and } (1b) \text{ hold}\}$ is the so-called equilibrium manifold, and $\mathcal{Z}_D = \{z \in \mathcal{Z}; z(x) \in D; \forall x \text{ control point in } \Omega\}$ is the so-called material manifold, also referred to as Data-function space (Nguyen et al., 2020), with the dataset $D = \{\hat{z}_i = (\hat{\epsilon}_i, \hat{\sigma}_i) \in \mathbb{Z}, i = 1, \dots, N_d\}$, $N_d$ being the number of material data points. By convention, thereafter we denote $z^* \in \mathcal{Z}_D$ in contrast to $z \in \mathcal{Z}_E$. In practice, an appropriate finite-dimensional subspace of $\mathcal{Z}_E$ (e.g. finite elements) should be considered such that a given function $z$ is unequivocally defined by a finite number of control points (Nguyen et al., 2020), usually associated to integration points, numbered $j \in \{1, \dots, N_g\}$ from now on. Note that $\mathcal{Z}_D$ translates the dataset $D$ into discretised functions induced by the same discretization of $\mathcal{Z}_E$. The mechanically admissible states then live in $\mathcal{Z}_E$ while the material states are taken from $\mathcal{Z}_D$.

The standard algorithm (Kirchdoerfer and Ortiz, 2016) used to deal with the double minimisation problem (4) encompasses the fixed-point iteration

$$ \mathbf{z}^{(k+1)} = P_E P_D \mathbf{z}^{(k)}, \quad (5) $$

until convergence, where

$$ \begin{gathered} P_E : \mathcal{Z}_D \to \mathcal{Z}_E \\ \mathbf{z}^* \mapsto P_E \mathbf{z}^* = \arg \min_{\mathbf{z} \in \mathcal{Z}_E} d(\mathbf{z}, \mathbf{z}^*) \quad (6a) \end{gathered} $$

is defined by the solution of FE-like linear systems followed by post-processing operations and

$$ \begin{gathered} P_D : \mathcal{Z}_E \to \mathcal{Z}_D \\ \mathbf{z} \mapsto P_D \mathbf{z} = \arg \min_{\mathbf{z}^* \in \mathcal{Z}_D} d(\mathbf{z}, \mathbf{z}^*) \quad (6b) \end{gathered} $$

is a nearest neighbour search (k-NN) in a finite dataset $D$.

Given $\mathbf{z}^* \in \mathcal{Z}_D$, it is convenient to rewrite $P_E \mathbf{z}^* \in \mathcal{Z}_E$ to a more tractable format. The projection can be seen as a constrained minimisation of the distance, subject to imposed equilibrium (1a) and kinematic compatibility (1b). By incorporating Lagrange Multipliers $\lambda$ to weakly enforce the former, and writing the functional in terms of displacements instead of strains to automatically enforce the latter, we get the unconstrained saddle-point problem
as follows

$$ (\mathbf{u}, \boldsymbol{\sigma}, \boldsymbol{\lambda}) = \arg \min_{\mathbf{u}' \in \mathcal{U}} \min_{\boldsymbol{\sigma}' \in \mathcal{S}} \max_{\boldsymbol{\lambda}' \in \mathcal{V}} \left\{ \frac{1}{2} d^2((\nabla^s \mathbf{u}, \boldsymbol{\sigma}), \mathbf{z}^*) - \langle \boldsymbol{\sigma}, \nabla^s \boldsymbol{\lambda} \rangle + \Pi_{\text{ext}}(\boldsymbol{\lambda}) \right\}, \quad (7) $$

where $S = L^2(\Omega; \mathbb{R}_{sym}^{d,d})$, $\mathcal{U} = \{\mathbf{w} \in H^1(\Omega; \mathbb{R}^d), \mathbf{w}|_{\Gamma_D} = \bar{\mathbf{u}}\}$, $\mathcal{V} = \{\mathbf{w} \in H^1(\Omega; \mathbb{R}^d); \mathbf{w}|_{\Gamma_D} = \mathbf{0}\}$, $\Pi_{\text{ext}}(\star) = \langle \bar{\mathbf{t}}, \star \rangle_{L^2(\Gamma_N; \mathbb{R}^d)} + \langle \mathbf{b}, \star \rangle$ is the external (virtual) work functional, with $\bar{\mathbf{t}}$ and $\mathbf{b}$ being admissible tractions on the Neumann boundary and body forces, respectively.

Finding the stationary conditions by taking the Gateaux derivatives on the three variables and by static condensation of stresses, we obtain the following variational equations

$$ \langle \mathbb{C} \nabla^s \mathbf{u}, \nabla^s \delta \mathbf{u} \rangle = \langle \mathbb{C} \boldsymbol{\varepsilon}^*, \nabla^s \delta \mathbf{u} \rangle \quad \forall \delta \mathbf{u} \in \mathcal{V}, \quad (8a) $$

$$ \langle \mathbb{C} \nabla^s \boldsymbol{\lambda}, \nabla^s \delta \boldsymbol{\lambda} \rangle = \Pi_{\text{ext}}(\delta \boldsymbol{\lambda}) - \langle \boldsymbol{\sigma}^*, \nabla^s \delta \boldsymbol{\lambda} \rangle \quad \forall \delta \boldsymbol{\lambda} \in \mathcal{V}, \quad (8b) $$

and the post-processing strain/stress updates

$$ \boldsymbol{\varepsilon} = \nabla^s \mathbf{u}, \quad (9a) $$

$$ \boldsymbol{\sigma} = \boldsymbol{\sigma}^* + \mathbb{C} \nabla^s \boldsymbol{\lambda}. \quad (9b) $$

Finally, the resulting mechanical state $\mathbf{z} \in \mathcal{Z}_E$ used for feeding the $P_D$ subproblem (6b) should be understood at the integration point level (once finite element spaces is taken for $\mathcal{V}$), as follows

$$ \hat{\mathbf{z}}_j = (\hat{\boldsymbol{\varepsilon}}_j := \nabla^s \mathbf{u}(\hat{\mathbf{x}}_j), \hat{\boldsymbol{\sigma}}_j := \boldsymbol{\sigma}^*(\hat{\mathbf{x}}_j) + \mathbb{C} \nabla^s \boldsymbol{\lambda}(\hat{\mathbf{x}}_j)), \quad (10) $$

for a generic integration point $\hat{\mathbf{x}}_j, j \in \{1, \dots, N_g\}$.

It is worth mentioning that updates for strain/stresses (9) are performed at the Gauss point level and that the left-hand-side of the resulting linear systems in (8) is constant for all iterations, so that the reutilization of a factorization is possible for all iterations, decreasing considerably the computational cost. As for the $P_D$ part, an efficient way to perform the k-NN search is using the KD-tree data structure. It is important to remark that $P_D$ is generally more costly than $P_E$, so the importance of using efficient implementations. In particular, ddfenics (Rocha, 2023a) is responsible for orchestrating $P_D$ and $P_E$, respectively implemented by KD-Tree using scikit-learn (Pedregosa et al., 2011) and finite elements with FEniCS/FEniCSx (Alnaes et al., 2015; Baratta et al., 2023).

All-in-all, we list in the Algorithm 1, the pseudocode for the alternate minimisation loop (5). We have used a convergence criterion based on the relative change of the global distance. Also, an initial guess $\mathbf{z}_{\text{mat}}$ in form of a $D$-admissible should be given. Usually, the null state is a suitable candidate as shown in Platzer et al. (2021). The resulting solver DDCM is from now on used as a black-box, except for a slight modification of the `solve_PE` that will be introduced in Section 3.
input : z_mech, DB, tol, maxit

output: z_mech, z_mech, disp

Standard initialisations: i, d, d0, err;

while err < tol and i < maxit do

<figure><img src="image_2.png" /></figure>

z_mech, disp ← solve_PE (z_mech);

z_mech ← solve_PD (z_mech, DB);

d ← dist (z_mech, z_mech);

err ← $\frac{d-d0}{d0}$ ;

d0 ← d;

i ← i + 1 ;

end

**Algorithm 1:** Pseudo-code for the alternate minimisation iteration loop used in Data-driven mechanics.

## 2.2. Creation of the dataset using computational homogenisation of the microscale

The DDCM solver, outlined in Algorithm 1, is inherently dependent on the fundamental input, namely the material dataset *D*, represented by the variable *DB*. One convenient method for generating a suitable dataset is to perform computations at a lower scale, where the constitutive models become simpler and thus more reliable. Furthermore, the quality of the generated material data points can be readily controlled. For instance, imposing a macroscopic strain on the boundary of the microstructure allows for the computation of homogenised stress. This section provides an overview of the fundamental principles and assumptions underlying computational homogenisation (CH) and its application in the present work. Given the extensive literature on this topic, and the fact that CH is not the primary focus of the present paper, the interested reader is directed to standard references such as Hill (1965); Mandel (1972); de Souza Neto et al. (2015); Blanco et al. (2016).

For each macroscale (integration) point $\hat{\mathbf{x}} \in \Omega$ we associate a microscale domain $\Omega_{\mu} \subset \mathbb{R}^d$. It is considered that $L_{\mu}/L \ll 1$, with $L$ and $L_{\mu}$ the characteristic lengths related to the macroscale and microscale, respectively. The characteristic size $L_{\mu}$ is chosen such that $\Omega_{\mu}$ should have meaningful numbers and type of material heterogeneities while keeping the compromise in terms of computational cost.

In terms of kinematics, the basic assumption in the CH approach is that the microscale kinematics objects correspond to the macroscale analogues in the averaged sense. In other words, given a macroscale strain $\hat{\varepsilon}$, the microscopic displacement field $\mathbf{u}_{\mu} \in H^1(\Omega_{\mu}; \mathbb{R}^d)$ is admissible if $\hat{\varepsilon} = \langle \varepsilon_{\mu}(\mathbf{u}_{\mu}) \rangle_{\Omega_{\mu}}$ holds, where $\varepsilon_{\mu}(\star) := \nabla^s(\star)$ is the symmetric gradient operator on $\Omega_{\mu}$ coordinates. It is out of the scope of this work to discuss in full generality
the definition of functional spaces enforcing these constraints, we refer the interested reader to Blanco et al. (2016); de Souza Neto et al. (2015). In this work, we assume periodic fluctuations (difference between microscale displacements and the affine displacements) on $\partial\Omega_\mu$, which is well-known to guarantee the average constraint (Blanco et al., 2016; de Souza Neto et al., 2015). As a result of the balance of virtual powers in macro- and microscale, the local equilibrium problem at the microscale level and the subsequent stress homogenisation read as follows (Blanco et al., 2016): given the macroscopic strain $\hat{\varepsilon}$, find an admissible $\mathbf{u}_\mu$ and the homogenised stress $\hat{\sigma} = \hat{\sigma}(\mathbf{u}_\mu(\hat{\varepsilon}))$ such that

$$ \left\{ \begin{array}{l} \langle \sigma_\mu(\varepsilon_\mu(\mathbf{u}_\mu)), \varepsilon_\mu(\mathbf{v}) \rangle = 0 \quad \forall \mathbf{v} \text{ admissible variations,} \\ \\ \hat{\sigma} = \langle \sigma_\mu(\varepsilon_\mu(\mathbf{u}_\mu)) \rangle_{\Omega_\mu}, \end{array} \right. \qquad (11) $$

with $\sigma_\mu = \sigma_\mu(\varepsilon_\mu)$, a constitutive law at microscale. As “admissible”, we simply mean sufficiently regular functions, satisfying the average constraint and periodicity.

In this work, CH is considered a virtual testing toolbox, or constitutive box “oracle”, which enables the generation of synthetic databases. When a macroscopic strain $\hat{\varepsilon}$ is provided, the corresponding homogenised stress $\hat{\sigma}$ is generated, forming a material data point $(\hat{\varepsilon}, \hat{\sigma}) \in D$. The aforementioned “oracle” can be queried either to generate databases in an *offline* fashion, as demonstrated by Xu et al. (2020), or in an *online* fashion, in the same spirit as introduced by Platzer (2020).

The latter approach is demonstrated to have advantages over the former in Section 2.4. The present work is concerned with the investigation of an optimal on-the-fly update of the database, which will be extensively detailed in Section 4.

### 2.3. Coupling scales: a first strategy

The work of Platzer (2020) introduced a first strategy for coupling DDCM with CH in order to build material databases on-the-fly. The proposed two-scale solver is an iterative scheme which updates the material dataset $D$ with selected microscopic evaluations in between successive runs of a macroscopic DDCM solver. The update strategy relies on a specific output of DDCM: the final local distance $d_{ij}$ between a mechanical state $\hat{z}_j = (\hat{\varepsilon}_j, \hat{\sigma}_j)$ and its material state $\hat{z}_i^* = (\hat{\varepsilon}_i^*, \hat{\sigma}_i^*) \in D$, $i \in \{1, \dots, N_d\}$ at integration point $\hat{x}_j$, $j \in \{1, \dots, N_g\}$, as defined by (3):

$$ d_{ij} = \| (\hat{\varepsilon}_j - \hat{\varepsilon}_i^*), (\hat{\sigma}_j - \hat{\sigma}_i^*) \|_{\mathbb{C}}. \qquad (12) $$

Indeed, if the k-NN distance is too large (e.g. larger than a user-defined critical value $d_c$), it is interpreted that the surrounding region of the phase space is lacking information, i.e. material data points, as illustrated in
Figure 1 (left). Such mechanical states are designated as isolated mechanical states. The set of isolated integration points is denoted by I and writes

$$I = \{j = 1, \dots, N_g \mid \min_i d_{ij} > d_c\} \subset \{1, \dots, N_g\}. \quad (13)$$

Notice that $I = I(D, z)$. Then, considering a material database $D^{(k)}$ at iteration $k$, from which results a data-driven solution $z^{(k)}$, we enrich the database for the next iteration $k + 1$ via the homogenisation of isolated mechanical states as follows:

with $I^{(k)} = I(D^{(k)}, z^{(k)})$ and $\hat{\sigma}_j^H = \hat{\sigma}_j^H(\hat{\varepsilon}_j)$ the homogenised stress resulting from the resolution of the CH problem (11). This strategy is illustrated in Figure 1.

It should be noted that when a rather sparse database is used as a starting point, the number of isolated mechanical states in set I can reach $N_g$, which in turn can result in a significant number of microscopic simulations. Consequently, we henceforth define a maximum number of points, $N_a^{\max} < N_g$ to be included in the material database update. The elements in I are then sorted in descending order of the distance, and at most the first $N_a^{\max}$ (and most critical) points are selected for addition.

<figure><img src="image_3.png" /><figcaption>Figure 1: Database update strategy from Platzer (2020). New material data points are obtained from the resolution of the CH (11), with prescribed macroscopic strains $\hat{\varepsilon}_j$, selected from the isolated mechanical states, lying too far from their associated material state $(\hat{\varepsilon}_i^*, \hat{\sigma}_i^*)$.</figcaption></figure>

The resulting iterative coupled DD-CH solver is given in the Algorithm 2, where the database update (14) is implemented in the function UpdateDB. In this context, DDCM stands for the Algorithm 1. The tolerance criterion is
based on the relative global distance, and the rest of the notations are self-explanatory. In the following sections, unless otherwise stated, the initial database (DB) is populated solely with the null state, represented by the set $D_0 = \{(0, 0)\}$.

input : DB, tolin, tolout, maxitout, maxitin
output: z_mech, z mat, disp

Standard initialisations: k, error, z mat;

while error < tolout and k < maxitout do
  z mat, z mech, disp ← DDCM (z mat, DB, tolin, maxitin);
  DB, z mat ← UpdateDB (DB, z mech, z mat); /* calls to CH
  */
  error ← distance relative (z mech, z mat);
  k ← k + 1;
end

// Optional final run of the DDCM with the last DB
z mat, z mech, disp ← DDCM (z mat, DB, tolin, maxitin);

<caption>**Algorithm 2:** Coupled DD-CH iteration loop.</caption>

## 2.4. Illustrative example

### 2.4.1. Macroscale problem

The performance of the coupled DD-CH solver is assessed with a first numerical test with the macroscopic problem described in Figure 2. The holed structure (Figure 2(a)) is clamped at the bottom surface $\Gamma_D^{\text{bot}}$ and subjected to quasi-static loading scenario through imposed displacements on the top boundary $\Gamma_D^{\text{top}}$. The prescribed displacements are given by $\mathbf{u}_D^{\text{top}}(t) = (\max(0, 5(t-3)/3), \min(5, 5t/3))$, for $t \in [0, 6s]$ in 6 increments (see Figure 2(b)). Both geometry and loading are inspired from mechanical displacement-driven experiments presented in Costecalde (2023), that were designed to obtain a rich set of mechanical states from a single test. As for the FE discretization, linear triangles have been adopted, with one Gauss point per element.

### 2.4.2. Microscale model

At microscale, we consider a soft matrix ($\Omega_{\mu}^{m}$) with hard circular inclusion ($\Omega_{\mu}^{i}$) in a square unit cell as in Fig. 1 (right), such that $\Omega_{\mu} = \Omega_{\mu}^{m} \cup \Omega_{\mu}^{i}$ with $|\Omega_{\mu}^{i}|/|\Omega_{\mu}| = 1/8$. The material model is described by the following non-linear elastic strain energy

$$ \Psi(\varepsilon) = \frac{\lambda}{2} \left( (\text{tr} \, \varepsilon)^2 + \frac{\alpha}{2} (\text{tr} \, \varepsilon)^4 \right) + \mu \left( |\varepsilon|^2 + \frac{\alpha}{2} |\varepsilon|^4 \right) \quad (15) $$
<figure><img src="image_4.png" /><figcaption>(a) Holed plate geometry (dimensions in (b) Quasi-static displacement loading history).</figcaption></figure>

<caption>Figure 2: Holed plate macroscopic problem.</caption>

where the piece-wise constant Lamé parameters are such that

$$ \lambda(\mathbf{x}) = \chi_{\Omega^m}(\mathbf{x})\lambda^{(m)} + \chi_{\Omega^i}(\mathbf{x})\gamma\lambda^{(m)} \quad (16a) $$

$$ \mu(\mathbf{x}) = \chi_{\Omega^m}(\mathbf{x})\mu^{(m)} + \chi_{\Omega^i}(\mathbf{x})\gamma\mu^{(m)} \quad (16b) $$

with $\lambda^{(m)} = 60.7047$ MPa, $\mu^{(m)} = 26.0128$ MPa, $\gamma = 50$, $\chi_{\Omega^i}$ and $\chi_{\Omega^m}$ being the characteristic function for matrix and inclusion subdomains respectively. Note that the material parameters values were chosen such that volume averaged Lamé parameters match $\bar{\lambda} = 432.099$ MPa, $\bar{\mu} = 185.185$ MPa as found in Schröder et al. (2021). Moreover, the parameter $\alpha = 10^4$ rules the nonlinear behaviour of the material, even though small strains are considered for the sake of simplicity.

### 2.4.3. Ground-truth solution generation: FE²

In order to evaluate the data-driven outcome, a comparison is made with a reference solution for the aforementioned problem, which has been obtained using the established FE² method (Feyel, 2003). It is noteworthy that in this context, the constitutive behaviour is entirely governed by the microscale simulations. Both scales are simultaneously solved using a Newton-Raphson scheme, which in particular for the macroscale, requires a homogenised tangent operator given by CH, obtained semi-analytically by solving three linear microscale problems, one for each strain direction (de Souza Neto et al., 2015). It is worth remarking that the common stiffness matrix is factorised only once by the LU decomposition and then reused by backward-substitution multiple times. Noticing that the Newton-Raphson scheme at the microscale takes in average 5 iterations to converge, the combined stress and tangent homogenisation is then only about 20% more costly than the simple stress homogenisation. We assume a convergence tolerance
of $10^{-8}$ for the relative error on the displacement euclidean norm in both scales. We employ the open-source library micmacsfenics (Rocha, 2023b; Rocha et al., 2023) to orchestrate two-scale FE computations also using FEniCS/FEniCSx.

### 2.4.4. Preliminary results

The coupled DD-CH iterative scheme (cf. Algorithm 2) was run for the 6 load steps given in Figure 2(b). At each load step, the material database is initialised from the converged material database of the previous one, whereby $D^{(0, t_n)} = D^{(k_f, t_{n-1})}$, with $D^{(0, t=0)} = D_0$, $n$ being the current time iteration index and $k_f$ the converged iteration of the Algorithm 2. The maximum number of newly added points is set to $N_a^{\max} = 30$ for the sake of simplicity. Sensibility analysis concerning the $N_a^{\max}$ choice is presented later on, in Section 5. The final resulting stress fields are compared in Figure 3. The results presented herein focus on the $\sigma_{xy}$ component, yet the conclusions drawn are consistent across all components. The data-driven scheme yields mechanical states that are closely aligned with the FE$^2$ reference results.

<figure><img src="image_5.png" /><figcaption>Figure 3: Stress fields comparison between of the first proposed coupled DD-CH scheme and ground-truth FE$^2$ for the $\sigma_{xy}$ component. Similar results (not shown) are observed for the other stress components.</figcaption></figure>

Furthermore, it is possible to observe the database that has been generated in an on-the-fly manner throughout the process. At the final stage (step 6), the database comprises 5251 material data points, which can be observed in Figure 4. These points represent the strain-stress relationship in a six-dimensional space, and for visual clarity, they have been represented in two reduced spaces, namely mean and deviatoric stress vs. strain. In this representation, instead of displaying the scattered points, their density is mapped, which provides a more intuitive understanding of the data structure.
<figure><img src="image_6.png" /><figcaption>Figure 4: 2D hexagonal binning plot of the final material database generated on-the-fly by the coupled DD-CH scheme.</figcaption></figure>

Another representation of interest is the distribution of the material data points in three-dimensional strain space ($\varepsilon_{xx}, \varepsilon_{yy}, \varepsilon_{xy}$), as illustrated in Figure 5. This representation shows the evolution of the data with loading. As anticipated, the database expands between each load step. However, it is evident that the material data points are not evenly distributed, as they would have been in an offline, regular sampling of the strain space. This is an inherent feature of the method, which adds data points selectively, favouring regions of the strain space that are covered by the macroscopic solution.

### 2.4.5. Limitations

The results presented above are highly encouraging, as the proposed approach successfully builds a tailored material database for the considered problem with a decreased number of microscopic evaluations. Nevertheless, when the proposed method was applied to another macroscopic geometry (e.g. the Cook's membrane as analysed in Section 5) in order to assess its versatility, serious robustness issues have been observed. In particular, the iterative scheme failed to converge whenever the near-empty database $D_0$ is considered as initial dataset, this without including further adjustments. Indeed, the material database update ultimately deviated from the expected outcome as the newly added points were increasingly distant from the solution. This is attributed to the inherent tendency of the DDCM alternate minimisation to converge towards local minima (Kanno, 2019). Furthermore, it can be similarly anticipated that even in a converging simulation, the current scheme is not optimised. This is because a significant number of microscale queries are still being requested at too far from the reference solution, which slows down the convergence process.
<figure><img src="image_7.png" /><figcaption>Figure 5: Strain of the final material database generated on-the-fly by the coupled DD-CH scheme.</figcaption></figure>

It can be concluded that this initial approach lacks the necessary robustness to be applied to a diverse range of problems and its convergence rates can be considerably improved. In particular, two main weaknesses have been identified:

1. In the initial iterations, when the current database is sparse and does not necessarily meaningfully cover the phase-space (in the limit of a near-empty database), the convergence of the DDCM is challenging, and the newly added points may not be optimally selected. An enhanced algorithm would be able to provide more meaningful candidates.

2. It is possible to track which material data points are utilized by the DDCM solver, revealing that a significant number of points are computed at the microscale but are not relevant for the solution at the macroscale. An improved update strategy would enable the selection of more appropriate candidates for the update.

These two identified limitations provide a foundation for the remaining contributions of this work. First, boosting strategies for the DDCM solver are proposed in Section 3. These strategies aim to enhance the solver's convergence and robustness when operating with sparse databases. Second, an on-the-fly adaptive sampling paradigm is introduced in Section 4. This
framework rules the selection of isolated mechanical states that are valuable additions to the material database.

### 3. Boosting DDCM for sparse databases

As already commented, our final goal is to reduce the number of calls to the local solver until convergence according to given criteria. In this section, we address point i) of Section 2.4.5. More specifically, we introduce three novel strategies to boost the DDCM solver, each of them targeting different weaknesses of the classical DDCM solver, with the common goal of providing more meaningful candidates to the on-the-fly adaptive sampling strategy of Section 4. For the sake of convenience they are listed below followed by a short description, which is further developed in dedicated sections:

1. Accelerated Projection-to-Equilibrium (APE): enhanced version of the classical $P_E$ operator that boosts convergence of the alternated minimisation iterations (see Section 3.1).

2. Sylvester Mean Metric (SMM): replaces the metric hyperparameter $C$ by an automatic, database-guided estimation given by the solution of a Sylvester equation (see Section 3.2).

3. Defective Restarting (DR): helps to avoid local minima by multiple slightly perturbed reinitialisations (see Section 3.3).

Moreover, we have also implemented the so-called Locally Convex Reconstruction (LCR) introduced by He and Chen (2020), which can be seen as a change on the $P_D$ operator such that it looks for the closest point belonging to a locally convex manifold (see Section 3.4). The different strategies are combined in an incremental manner and assessed through another preliminary numerical example in Section 3.5, before applying them in combination with the on-the-fly database updates, which is subject of the next section (Section 4). For a reader exclusively interested in the online sampling selection part, it should be noted that one may proceed directly to Section 4 without compromising comprehension.

#### 3.1. Accelerated Projection-to-Equilibrium

The standard alternated minimisation iteration (5) writes as $\mathbf{z}^{(k+1)} = P_E P_D \mathbf{z}^{(k)}$. It has been observed that this approach alone leads to poor convergence rates, particularly for sparse data scenarios. Below we propose a modified $P_E$ operator, denoted $\tilde{P}_E$, that encapsulates an approximation of $P_D$, denoted $\tilde{P}_D$, aiming at reaching the result of two successive iterations at the expense of only one iteration.

Let us then consider the following ansatz

$$ \tilde{\mathbf{z}}^{(k+1)} = \tilde{P}_E P_D \mathbf{z}^{(k)} = P_E \tilde{P}_D P_E P_D \mathbf{z}^{(k)} = P_E \tilde{P}_D \mathbf{z}^{(k+1)}, \quad (17) $$
where $\tilde{P}_D$ can be regarded as a continuous version of $P_D$, such that $(P_E P_D)^2 \approx P_E \tilde{P}_D P_E P_D = \tilde{P}_E P_D$. For constructing the operator $\tilde{P}_D$, we have to first look at the local level and build a tangent reconstruction around the mechanical state at a given point.

Consider a mechanical state $\hat{z} = (\hat{\varepsilon}, \hat{\sigma})$ and its material states versions $\hat{z}_\sigma^H = (\hat{\varepsilon}, \sigma^H(\hat{\varepsilon}))$ and $\hat{z}_\varepsilon^H = (\varepsilon^H(\hat{\sigma}), \hat{\sigma})$, where $\sigma^H(\cdot)$ encapsulates the homogenisation procedure for the stress defined in (11), while $\varepsilon^H(\cdot)$ can simply be seen as the inverse of the homogenised stress mapping. Note that here we use the CH operators for the sake of convenience and consistency with our work, but any other direct or inverse constitutive mapping are allowed. Let us estimate the material states versions of $\hat{z}$ using the first-order Taylor expansion around its nearest-neighbour $\hat{z}^*$ on the database with

$$ \hat{z}_{\sigma}^{H}(\hat{z}) \approx \hat{z}^{*} + (\Delta\hat{\varepsilon}, \partial_{\hat{\varepsilon}}\sigma^{H}(\hat{\varepsilon}^{*})\Delta\hat{\varepsilon}) \approx \hat{z}^{*} + (\Delta\hat{\varepsilon}, \mathbb{C}\Delta\hat{\varepsilon}), \quad (18a) $$

$$ \hat{z}_{\varepsilon}^{H}(\hat{z}) \approx \hat{z}^{*} + (\partial_{\hat{\sigma}}\varepsilon^{H}(\hat{\sigma}^{*})\Delta\hat{\sigma}, \Delta\hat{\sigma}) \approx \hat{z}^{*} + (\mathbb{C}^{-1}\Delta\hat{\sigma}, \Delta\hat{\sigma}), \quad (18b) $$

where $\Delta\hat{\varepsilon} := \hat{\varepsilon} - \hat{\varepsilon}^*$, $\Delta\hat{\sigma} := \hat{\sigma} - \hat{\sigma}^*$, and the tangent operators estimates $\partial_{\hat{\varepsilon}}\sigma^H(\hat{\varepsilon}^*) \approx \mathbb{C}$ and $\partial_{\hat{\sigma}}\varepsilon^H(\hat{\sigma}^*) \approx \mathbb{C}^{-1}$ have been assumed. The reasons for such approximations are twofold: i) the proper homogenised tangent operators computations would add unnecessary computational burden, and ii) the use of the same data-driven distance yields a straightforward, explicit update for the proposed enhanced mechanical state (see Appendix B.2 for details).

For the sake of generality, let us take a convex combination of the two estimates in (18) to define

$$ \hat{z}'(\theta) = \hat{z}^* + \mathbb{G}_{\theta}(\hat{z} - \hat{z}^*), \quad \theta \in [0, 1], \quad (19) $$

$$ \text{with } \mathbb{G}_{\theta} = \begin{bmatrix} \theta \mathbb{I} & (1-\theta)\mathbb{C}^{-1} \\ \theta\mathbb{C} & (1-\theta)\mathbb{I} \end{bmatrix}. \quad (20) $$

Note that $\hat{z}'$ is in general not contained in the database, but it is an approximation of it, as illustrated in Figure 6 (i.e. it does not lie on the underlying constitutive manifold). Finally, as there is no clear guidance for choosing $\theta$, we consider $\tilde{P}_D$ as the global version of (19) for $\theta = \frac{1}{2}$. After some tedious steps, detailed in Appendix B.2, the accelerated solution is given by

$$ \tilde{z}^{(k+1)} = \frac{1}{2}(z^{(k+1)} + z^*) + \mathbb{G}_{\frac{1}{2}}(z^{(k+1)} - P_D z^{(k)}). \quad (21) $$

**Remark 2.** In fact, $\hat{z}'$ from (19) may also serve as a more reliable candidate for the online selection phase, as it combines both material and mechanical states in such a way that it approximates the tangent space of implicitly-defined CH material manifold around the annotated point $\hat{z}^*$. The specifics of this approach, referred to as Tangent Reconstruction (TR), will be discussed in greater detail in Section 4.1. It is worth mentioning that (21) can
also be used as ansätz in finite strain problems, even though some parts of
the proof (see Section Appendix B.2) strongly rely on the small strain set-
ting. Therefore, we expect an APE method that is less performant yet still
effective.

<figure><img src="image_8.png" /><figcaption>Figure 6: Convex combination of first-order Taylor expansions (right) as compared to the “exact” CH updates (left).</figcaption></figure>

3.2. Sylvester Mean Metric

The metric $\mathbb{C}$ is regarded as a hyperparameter in most of DDCM works.
A notable alternative to this approach is to perform $\mathbb{C} := \text{sym}(\mathbf{A}_{\sigma}\mathbf{A}_{\varepsilon}^{-1})$,
where the square matrices $\mathbf{A}_{\sigma}$ and $\mathbf{A}_{\varepsilon}$ are obtained by truncating the SVD
decomposition of the dataset matrix (Eggersmann et al., 2021). Note that
SVD is not able to naturally enforce the symmetry of $\mathbb{C}$. The symmetrization
is then taken as second step, which degrades the optimal low-rank properties
SVD normally possess.

In this work, we set $\mathbb{C}$ as the symmetric matrix that minimises the mean squared sum of gaps between the data points $(\hat{\varepsilon}_i, \hat{\sigma}_i)$ and the dataset barycentre $(\bar{\varepsilon}, \bar{\sigma})$, such as

$$
\mathbb{C} = \arg \min_{S \in \mathbb{R}^{N_{\varepsilon} \times N_{\varepsilon}}} \frac{1}{N_d} \sum_{i=1}^{N_d} \| S(\hat{\varepsilon}_i - \bar{\varepsilon}) - (\hat{\sigma}_i - \bar{\sigma}) \|^2. \quad (22)
$$

We can show that the solution for this problem yields the Sylvester
equation below (see Appendix B.1 for proof)

$$
\Sigma_{\varepsilon\varepsilon} \mathbb{C} + \mathbb{C} \Sigma_{\varepsilon\varepsilon} = \Sigma_{\varepsilon\sigma} + \Sigma_{\varepsilon\sigma}^T, \quad (23)
$$

where $\Sigma_{\varepsilon\varepsilon}$ and $\Sigma_{\varepsilon\sigma}$ are $N_{\varepsilon} \times N_{\varepsilon}$ blocks extracted from the empirical co-
variance matrix of the database $D$. The symmetry of $\mathbb{C}$ holds because the
same matrix is used in both left and right multiplication of $\mathbb{C}$. Positive
definiteness is not enforced a priori, but it has been successfully checked
numerically for the examples of the present article.
### 3.3. Defective Restarting

It is well known that the alternate minimisation procedure used in DDCM has no guarantee of reaching global minima (Kanno, 2019). One could resort to mixed-integer quadratic programming approaches, although these techniques have prohibitive computational costs (Kanno, 2019). We propose a method to prevent the alternate minimisation algorithm from getting trapped in local minima, which consists in perturbing the nearest-neighbour search by considering the second-nearest point instead of the first. We recommend using this strategy in restarting cycles, i.e., that the perturbation be only performed on the material states at the end of a full standard DDCM solver run, and that it be repeated $n_{\text{restarts}}$ times. This perturbation should be performed randomly to some part of the states only, that are chosen with a probability $p_d$, typically within 5% to 10%.

It is worth remarking that our method is inspired by the Dropout regularisation used by the Deep Neural Networks community, a technique that consists in neglecting some neurons during training. The choice of neurons is done randomly according to a given probability value. Dropout promotes models that generalise better by preventing the optimisation algorithm from getting trapped in local minima, similar to the Defective Restarting in the DDCM setting.

### 3.4. Locally Convex Reconstruction

As any surrogate model, Neural-Networks tend to smooth out the noise present in the data (if not overfitted) and have the ability to interpolate and/or extrapolate. In the standard DDCM setting however, the k-NN search prevents interpolation/extrapolation. While this can be a desirable property, it turns out to be a hindrance for zones of small density or large noise in the data. The so-called locally convex embedding He and Chen (2020) proposes instead the following modification for Data-functions set

$$ \hat{\mathbf{z}}^* = \sum w_i^* \hat{\mathbf{z}}_i, \qquad (24) $$

$$ \mathbf{w}_i^*(\mathbf{z}) = \arg \min_{\mathbf{w}} \|\mathbf{z} - \sum_{i \in \mathcal{N}_k(\mathbf{z})} w_i \hat{\mathbf{z}}^i\|^2, \qquad (25) $$

$$ \text{subject to :} \quad \sum_{i \in \mathcal{N}_k(\mathbf{z})} w_i = 1, w_i \ge 0, \forall i \in \mathcal{N}_k(\mathbf{z}), \qquad (26) $$

with $\mathcal{N}_k(\mathbf{z})$ being the set of k-nearest neighbours indexes. The equality constraint is relaxed and regularised with the recommended penalty terms as in He and Chen (2020) ($\bar{\mu} = 10^{-4}$, $\bar{\xi} = 10^4$, using the reference paper notation), and the non-negativity of weights is handled by the non-negative least-squares function `scipy{nls` of the Scipy Python library. Another alternative approach for DDCM regularisation is the so-called max-entropy approach of DDCM Kirchdoerfer and Ortiz (2017).
### 3.5. Preliminary results

The performance of the different boosting strategies are assessed with a preliminary numerical example using “offline” databases, for the sake of simplicity. Even though the final goal is to allow DD-AS simulations, with on-the-fly update of the database, we first aim at demonstrating the efficiency of the proposed approaches on DDCM independently. To evaluate the contribution of each strategy, the results are compared in the following order, chosen for the sake of clarity: 1) APE vs standard DD, 2) APE-SMM vs APE, 3) APE-SMM-LCR vs APE-SMM, and finally 4) APE-SMM-LCR-DR vs APE-SMM-LCR (acronyms are provided at the beginning of Section 3). The comparison is made assessing the evolution of standard relative $L^2$ errors on the mechanical fields {$e_\epsilon, e_\sigma, e_u$} with the number of points in the database, or, in some cases, the relative error in terms of the maximum norm {$e_\epsilon^{max}, e_\sigma^{max}, e_u^{max}$} (see Appendix A for more details). It is worth noticing that, since the metric $C$ is obtained independently for each realization of the coarsened dataset by the SVD decomposition as in Eggersmann et al. (2021) or by SMM, the Data-driven distances between different runs are not comparable, hence the use of relative errors on the mechanical fields instead of the objective function as a comparison tool.

In order to test the enhanced DDCM solver, we chose the Cook’s membrane problem due to the richness of strain/stress states it covers even for simple loading scenario. Moreover, as previously stated, the Cook’s membrane presented lack of convergence issues using the naive strategy introduced in Section 2, i.e., without enabling any of the boosting strategies proposed herein. The geometry is described in Figure 7(a), where the left side is fixed and a vertical load of magnitude of 2 MPa is applied on the right side. We consider two meshes as in Figure 7(b). The microscale and other FE details are identical to those introduced the previous holed plate example of Section 2.4.

To create the offline reference database, say $D^{ref}$, we first collect the results of a FE² simulation using the fine mesh (approximately 15K data states). If we simply feed this “static” database into the standard DD solver, in principle we are able to retrieve the same reference FE² solution at the converged solution. Such comparison is unfair regarding the dynamic DD-AS approach which does not have access to the reference solution. Therefore, it is more interesting to obtain a randomly perturbed (only on the strain part) database that approximately covers regions of strain-stress space spanned by the problem. We proceed by identifying 3000 clusters for the strains in the reference database and then we randomly sample 10 points around each cluster centroid, making up a total of 30K points in the database. More precisely, each sampled strain is drawn from a normal distribution with a standard deviation of 20% the global standard deviation in each direction and centred on the centroid strain of each cluster. The sampled strain states are then used to obtain the stress states via the usual CH
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Figure 7: Model two-scale problem: Cook's membrane.</figcaption></figure>

procedure. Hence, the resulting perturbed database should cover a slightly larger region of the strain-stress space, as well as being relatively denser, since the several perturbations have been considered for the strain.

As illustrated in Figure 8(a), the incorporation of the APE strategy already leads to notable improvements. We have performed tests with databases containing 25% of $D^{ref}$ points (picked randomly) until the full $D^{ref}$ to test whether the boosting potential is more pronounced in the case of sparse databases. Coarse databases are created by sampling (with no repetition) the dataset $D^{ref}$ according to a given database ratio percentage. Each coarseness level is solved for 5 different realizations for the sake of robustness. Figure 8(a) shows that gains are of the same order of magnitude in all $L^2$ metrics with respect to the reference FE$^2$ solution, and are also approximately even for all coarseness levels.

We proceed our analysis incrementally by incorporating an additional method into the same analysis. In Figure 8(b), the SVD-like evaluation of
534

535

536

537

538

539

540

C is replaced by SMM. It turns out that the evaluation of the stresses is improved at the expense of the accuracy of displacement and strain fields. Overall, we consider this to be a beneficial effect, given that stress generally exhibits higher errors. Furthermore, SMM has a clear theoretical foundation, as it is derived from the stationary point for the constitutive gap in the sense of the Gâteaux derivative by consistently taking symmetrical matrices as admissible variations.

<figure><img src="image_11.png" /><figcaption>(a) Effect of enabling APE (all others techniques deactivated).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) Effect of enabling SMM (with just APE enabled).</figcaption></figure>

<caption>Figure 8: Testing effects of APE and its combination with SMM.</caption>

541

542

543

544

545

546

547

548

In the same spirit, LCR is added to the analysis as it is expected to improve accuracy specially in highly sparse data scenarios. Coarseness levels vary from 1% to 25% of the database. Contrarily to our initial intuition, the LCR effect in roughly constant as shown in Figure 9(a), with maybe an exception for the strain errors, that are indeed improved for the coarsest databases. Results displayed in Figure 9(a) were performed with $k = 4$ (number of neighbours to perform the convex reconstruction). We have also studied the sensitivity from 3 to 6 neighbours, where results for $e_σ$, chosen
for the sake of simplicity, are displayed in Figure 9(b). Negligible effects are
observed for $k > 5$ and already little effects are seen for $k > 3$. As larger $k$
entails additional computational costs, we recommend $k = 6$ at most.

<figure><img src="image_13.png" /><figcaption>(a) Effect of enabling LCR (with APE + SMM enabled).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(b) Sensitivity to the number of neighbours in LCR.</figcaption></figure>

<caption>Figure 9: Testing effects of LCR.</caption>

Finally, we test DR in the case of the full *D*<sup>ref</sup> dataset utilization, with until 3 restarts. We can see in Figure 10(a) that accuracy in *L*<sup>2</sup>-norm is strongly improved for stress predictions, particularly after first and second restart loop. Moreover, Figure 10(b) demonstrates that the enhancing effect of DR in the maximum-norm of stress is also noticed even for the subsequent restarting loops. The effect of "dropping out" the best neighbour is clearly seen in the energy minimisation curves in Figure 10(c) (blue curve). Also note that fewer iterations are needed for the DDCM solver reach the convergence if compared with the first loop. Globally, we have a monotonically decreasing normalized DD distance (black curve), showing that DR allows for escaping local minima. To balance boosting effects while keeping relatively small computational cost we recommend using at most 2 restarting
loops.

<figure><img src="image_15.png" /><figcaption>(a) Effect of enabling DR (with APE + SMM + LCR(k=6) enabled) on L²-norm relative errors.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(b) Maximum-norm effect for DR.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(c) minimisation along DD iterations with restarting.</figcaption></figure>

Figure 10: Study effect of DR regarding the number of restarts.

Accepted Manuscript
## 4. An On-the-fly Adaptive Sampling approach for DDCM

The present section aims to set up the on-the-fly adaptive sampling framework suited for the DDCM, denoted thereafter DD-AS. Central to our exposition is the concept of *unlabelled* and *labelled* datasets and the mapping between them, defined in the next paragraph. The main goal of our method is to select instances of unlabelled data at the current DDCM iteration that are optimal in the sense of some metric. Finally, in Section 4.1, we provide a method for building the *unlabelled* dataset and in Section 4.2 we establish selection criteria for selecting unlabelled data.

More precisely, in the CH context, we can consider that the homogenised stress (see (11)) is a label for the strain, hence we define $D_L = \{(\hat{\varepsilon}, \hat{\sigma}) \in \mathbb{Z}; \hat{\sigma} = \sigma^H(\hat{\varepsilon})\}$, i.e., the graph of $\sigma^H$, as the admissible labelled states set. In other words, labelled states are those living in the constitutive manifold of Fig. 6. On the other hand, the unlabelled admissible set $D_U$ is associated with material point evaluations of the mechanical and material states (for further details see Section 4.1). For the sake of simplicity, consider $D_U = \mathbb{Z}$ without loss of generality. To distinguish notations, it is useful to set $D_L \subset D_L$ and $D_U \subset D_U$ as finite collections of points in the corresponding admissible sets, denoting labelled and unlabelled datasets, respectively. Moreover, it is useful to define the Label mapping (note that in general $\hat{\sigma} \neq \sigma^H(\hat{\varepsilon})$) to convert elements of $D_U$ into $D_L$ as follows:

$$ \text{Label} : D_U \to D_L \qquad (27) \\ (\hat{\varepsilon}, \hat{\sigma}) \mapsto (\hat{\varepsilon}, \sigma^H(\hat{\varepsilon})). $$

In practice, the Label mapping should be applied with parsimony since CH is computationally intensive. In terms of implementation of the iterative procedure, at an iteration $k$, given a pool of unlabelled data $D_U^{(k)}$, the pool-sampling problem consists in selecting $\tilde{D}_U^{(k)} \subset D_U^{(k)}$ such that they correspond to the highest values of a given criterion score $M: D_U \to \mathbb{R}$. Hence, the labelled dataset increments as follows $D_L^{(k+1)} = D_L^{(k)} \cup \text{Label}(\tilde{D}_U^{(k)})$, to be used in the next DDCM iteration. Finally, the iterative algorithm follows exactly the same steps as already described in Algorithm 2, where the pool construction and the criterion $M$ are encapsulated in the UpdateDB function, described in sections 4.1 and 4.2, respectively.

**Remark 3.** It is worth mentioning that DDCM does not easily fit into the classical framework of machine learning methods; as such, it cannot be categorized strictly as either supervised or unsupervised. Strictly speaking, with regard to the mixed use of labelled and unlabelled phase-space datasets (as defined above), the proposed method can be interpreted as a semi-supervised strategy. In this context, we should acknowledge the inspiration drawn from the so-called active learning paradigm (Goodfellow et al., 2016), particularly
in the design of some of the selection metrics of Section 4.2, which aim to rank instances from a pool of unlabelled candidates for labelling. Similarly to our case, the agent responsible for labelling (the oracle) should be called sparingly, as it involves high costs (computational, time, financial, etc.). Finally, the mapping Label (27) strictly adheres to the constitutive interpretation of labelling, regarded as a regression problem. If one wishes to adopt a classification viewpoint – specifically, in our case, using the k-NN classifier – labels can be straightforwardly assigned as integer indices corresponding to states in $D_L$, for example, $(\hat{z}_1, 1), (\hat{z}_2, 2)$, where $\hat{z}_i \in D_L, \forall i \in \{1, 2\}$.

## 4.1. Unlabelled dataset

In this section, we aim at building $D_U$ knowing mechanical and material state solutions of given DDCM iteration. A straightforward method has already been introduced in Section 2.3, purely based on mechanical states data, thereafter referred to as “naive”. A second method, the so-called Tangent Reconstruction (TR), is motivated by the convex combination of first-order Taylor expansions around material states and mechanical states (see (19)) as derived next.

First, let us make a distinction between $\hat{z}' = (\hat{\varepsilon}', \hat{\sigma}') \in D_U$, $\hat{z} = (\hat{\varepsilon}, \hat{\sigma})$ and $\hat{z}^* = (\hat{\varepsilon}^*, \hat{\sigma}^*)$, point-wise evaluations of $z \in Z_E$ and $z^* \in Z_D$, respectively. Invoking (19), we take $\hat{z}'(\theta) = \hat{z}^* + G_\theta(\hat{z} - \hat{z}^*)$ as admissible elements in $D_U$, where $\theta \in [0, 1]$ is a weighting parameter. As there is no clear guidance on the choice of the weighting parameter, $\theta$ can be considered as a random variable from a uniform probability distribution in $[0, 1]$. Hence, $D_U$ is constructed by collecting $\hat{z}'$ for each Gauss point and for a random realisation of $\theta$.

As discussed in Remark 2, the tangent reconstruction operator (19) can also be understood as a way to obtain more meaningful update candidates. Indeed, $\text{Label}(\hat{z}'(\theta)) = (\varepsilon'(\theta), \sigma^H(\varepsilon'(\theta)))$ takes into account both strain and stress via $\varepsilon'(\theta) = (1-\theta)\hat{\varepsilon}^* + \theta\hat{\varepsilon} + (1-\theta)C^{-1}(\hat{\sigma} - \hat{\sigma}^*)$, which is not the case for naive choice that simply states $\hat{\varepsilon}' = \hat{\varepsilon}'(\theta=1) = \hat{\varepsilon}$, in purely strain terms. It is worth remarking that although the Label mapping is solely function of strains, the stress part of the unlabelled state remains important since it is used to compute DDCM distances in Section 4.2. From Figure 6, also note that for intermediate values of $\theta$, $\hat{z}'(\theta)$ tends to be closer to $\text{Label}(\hat{z}'(\theta))$ than $\hat{z}$ is to the $\text{Label}(\hat{z}'(1))$, the latter being the labelled state for the naive choice. Therefore, it is reasonable to consider the TR as enhanced choice to compose the pool of unlabelled data, since the associated scores would be more representative of the gains of adding $\text{Label}(\hat{z}')$. In addition, the random nature of $\theta$ naturally adds more variability to the states being labelled, which is also beneficial to escape local minima. Finally, the numerical example of Section 4.3.1 is devoted to assess numerically such performance gains.
## 4.2. Selecting unlabelled data

In this section, we describe some alternatives for the score function $M$, namely: greedy ($G$), uncertainty ($U$), and entropy ($E$). For the sake of convenience, the strategy of Section 2.3, purely based on a local distance criteria, is rephrased here with updated notation and corresponds thereafter to the $G$-score (greedy). Naturally, for all scores, we exploit the intrinsic error estimators of DDCM, namely the global metric $d$ and the local norm $||\cdot||_C$. Other strategies, that do not easily fit into the framework described next, and also that have not led to efficient results, are presented in Appendix C for the sake of completeness and comparison.

Before defining each method, it is convenient to reset some notations. Let us define the ordered local DDCM distance as

$$d_{ij} = \text{distance between } \hat{z}_j(\in D_U) \text{ and its } i\text{-th-NN in } D_L \quad (28)$$

$\forall i \in \{1, \dots, N_d\}, \forall j \in \{1, \dots, N_g\}$. It is important to mention that for the sake of performance, only a few nearest neighbours need to be queried, which is specified at the phase of the tree-based search object instantiation using Scikit-learn (Pedregosa et al., 2011). Additionally, note that we can write the least global DDCM distance (2) in terms of the ordered distance in (28) for a given $D_L$ as

$$\sqrt{\sum_{j=1}^{N_g} w_j(d_{1j})^2}, \quad (29)$$
where $w_j$ stands for the quadrature weight times the finite element volume (or area in 2D cases) of integration point $j$.

Finally, it is useful to think of the forthcoming scores as probabilistic indexes that measure the potential that a given selected state has at decreasing the DDCM distance (global or local) throughout the next iteration. Once the scores for all integration points are computed, a number $k_{update}$ of instances of $D_U$ are selected to be appended to $D_L$.

### 4.2.1. Greedy scores ($G$ and $G_v$)

The so-called greedy methods aim for labelling points that can decrease the most the cost function in each iteration (Wu et al., 2019). Further mathematical analysis of the greedy methods from the active learning community is provided in Dasgupta (2004). Greedy techniques are also a popular technique employed for optimal sampling in the context of Model Order Reduction (Binev et al., 2011).

The most straightforward strategy corresponds to exclusively focus on decreasing local DDCM distances, equivalently to the method of Section 2.3. Hence, the $G$-score reads as

$$G_j = d_{1j}, \quad \forall j = 1, \dots, N_g. \quad (30)$$
From (29), note that Gauss points associated to larger volumes will have a larger impact on the global metric. Therefore, an enhanced score can read as

$$ (G_v)_j = \sqrt{w_j} d_{1j}, \quad \forall j = 1, \dots, N_g, \qquad (31) $$

which only differs by the factor $\sqrt{w_j}$ from $G_j$. The latter is equivalent to the former approach in the special case of constant-sized and single quadrature point elements.

### 4.2.2. Uncertainty score (U)

Regarding the selection scores as probabilistic indexes, it is also reasonable to seek unlabelled points that are less likely to be associated with only one labelled point. In other words, we can look for unlabelled points that maximises the ratio between the first and second-nearest neighbours distances, in the same spirit of the so-called uncertainty index Gweon and Yu (2021) borrowed from the active learning literature, defined as

$$ \tilde{U}_j = \frac{d_{1j}}{d_{2j}}, \quad j = 1, \dots, N_g. \qquad (32) $$

Although the aforementioned concept is important by itself, we have observed that the distance magnitude cannot be completely disregarded. To be able to account for both guiding principles, we propose the use of the following metric

$$ U_j = \tilde{U}_j^{\alpha_U} \times G_j, \quad j = 1, \dots, N_g, \qquad (33) $$

also called uncertainty for the sake of simplicity, where $\alpha_U = 0.1$ is a recommended weighting parameter, which increases $\tilde{U}$ contribution for intermediates values. Finally, it is worth noticing that the U-score and the Least Confidence approach (Li and Sethi, 2006) – also proposed in the active learning context – would promote the same samples, thus we limit our analysis to the U-score.

### 4.2.3. Entropy score (E)

With similar goal than the uncertainty score, a more general concept rooted on the information theory is the so-called Shannon entropy, thereafter called E-score. While the U-score is solely based on the two nearest neighbours, Shannon entropy accounts for all distances between unlabelled and labelled points. For the sake of simplicity and computational performance, we only consider distances up to the k-th nearest neighbour. From the definition of Shannon entropy, the E-score reads as

$$ \tilde{E}_j = - \sum_{i=1}^{k} p_{ij} \log p_{ij}, \quad j = 1, \dots, N_g, \qquad (34) $$
where $p_{ij}$ is an empirical estimation for the probability of the unlabelled point $j$ being associated with its $i$-th nearest neighbour labelled point. We propose to compute the probabilities using the so-called softmax function:

$$p_{ij} = [\mathrm{softmax}(-\bar{\mathbf{d}}_{.,j})]_i, \quad j = 1, \dots, N_g, \quad i = 1, \dots, k, \qquad (35)$$

where $(\cdot)_{.,j}$ stands for the $j$-th column vector of matrix $(\cdot)$ and

$$\bar{d}_{ij} = \frac{d_{ij} - \max_{1 \le i \le k} d_{ij}}{\max_{1 \le i \le k} d_{ij}}, \quad [\mathrm{softmax}(\mathbf{y})]_i = \frac{\exp y_i}{\sum_{l=1}^{k} \exp y_l}. \qquad (36)$$

Note that by construction the softmax function yields non-negative probabilities and is such that the $\sum_{i=1}^k p_{ij} = 1$, filling the basic requirements for probabilities. Alternatively, the work of Prune et al. (2023) also provides another viable manner to construct probability estimations (once normalised), however with additional hyperparameters. Therefore, in this work we only use (35), which is parameter-free.

Similarly to the U-score, we have observed that the G-score cannot be completely disregarded, such that we use the modified E-score as follows

$$E_j = \tilde{E}_j \times G_j, \quad j = 1, \dots, N_g. \qquad (37)$$

### 4.3. Preliminary results

Before a complete final example in Section 5, it is important to investigate the novel ingredients introduced in the present section. First, we need to validate the boosting strategies of Section 3 in a fully DD-AS strategy. For the sake convenience, only the G-score and the performance of TR strategy for obtaining the unlabelled data is assessed altogether in the dedicated example of Section 4.3.1. Secondly, we need to compare the selection scores performances between themselves as addressed in Section 4.3.2.

#### 4.3.1. Assessing the performance of boosting strategies for DD-AS

In this example, we extend the previous Section 3.5 example to the DD-AS context. For ease of interpretation, we consider using only one score, the chosen Greedy method. For convenience, we consider the coarse mesh and fixed-size database updates, by collecting 20 points per outer loop iteration for labelling. The inner tolerance of Algorithm 2 is set to $10^{-4}$, while the outer tolerance is set to a very small value, to force leaving the outer loop for $\text{maxitout} = 100$ equally for every combination of methods. Moreover, the effect of considering TR for selecting elements in the pool is also assessed.

In Figure 11, we can see the evolution of relative errors for some meaningful combination of boosting DD strategies, where the acronyms stand for the letters of the enabled methods, e.g., ATSLD stands for combining APE+TR+SMM+LCR+DR effects. It is worth noticing that in this example, we did not achieve convergence for DD-AS if neither APE nor TR was
enabled, such that the simplest strategies are *T* or *A*. We can observe that *A* has more impact than *T* alone, but their combination *AT* is always gainful. We also note that *ATS* performs equivalently to *AT* for some metrics, but sensible accuracy gains are observed in stress-related metrics. Comparing *ATS* with *ATSL* (with *k* = 6), the latter strategy is recommended regarding all metrics. Finally, allowing 2 restarts, *ATSLD* improves mostly *max*-norm metrics while having slightly poorer results for the others *L*²-norm metrics. It should be mentioned that it is not by accident that *ATSLD* shows oscillation trends since it constantly seeks to avoid local minima.

Overall, the simplest strategies tend to reach a saturation point after which very little improvement is observed, especially in the stress-related metrics. The later this point is reached, the smaller the error levels. Seemingly, the maximum-norm error for the stress is the most challenging metric. Hence, we recommend using *ATSLD* as it reaches smaller errors in this metric and displays a less apparent saturation trend. Nevertheless, the final choice should also consider computational costs and the aimed accuracy levels. For example, for a database of size 1250, *ATSLD* and *ATSL* can be considered equivalent if all metrics are taken into account.

### 4.3.2. Contrasting ranking strategies

In this example, although *ATSLD* can perform potentially better than *ATSL*, we consider the latter *ATSL* as the standard method as it adds less noisy behaviour, hence simplifying the analysis. The following scores are considered: *G*, *G*<sub>v</sub>, *U* and *E*. For the sake of clarity, only combined (geometric average) error metrics $\bar{e} := \sqrt[3]{e_u e_e e_\sigma}$ and $\bar{e}^{max} := \sqrt{e_\epsilon^{max} e_\sigma^{max}}$ are reported (see discussion in Appendix A). Two mesh refinements are considered: M2 with 3868 elements, and M3 with 15125 elements. For mesh M2, we study the influence of $k_{update}$ using 40 and 80 as values, for 200 and 100 DD-AS iterations, respectively. For M3, only $k_{update} = 40$ during 800 iterations has been simulated, for the sake of simplicity. The number of iterations has been set to keep approximately the same ratio between the database size and the number of elements ($\approx 2$) for all cases. All other settings are identical to the study of Section 4.3.1.

In Figure 12 (left-hand side), we can observe as a general trend that $G_v$ performs better than the {G, U, E}-score approaches according to the metric $\bar{e}$, with the biggest gap in the case of higher $k_{update}$ with mesh M2 in Figure 12(c). The remaining metrics perform almost equivalently. For the metric $\bar{e}^{max}$ in Figure 12 (right-hand side), the best score varies for each case, with a slight advantage for the *E*-score, which performs best in 2 out of the 3 cases. Interestingly, the *U*-score is the best in case M2 ($k_{update} = 80$), but it is the worse in the other two cases. Still for $\bar{e}^{max}$ metric, we can observe a highly oscillatory, even non-monotonic, error evolution for $G_v$, which is explained by the reluctance of updating small-sized elements, even with reasonably high *G*-score. However, overall $G_v$ recovers similar performance
<figure><img src="image_18.png" /><figcaption>(a) Strain $L^2$-norm relative error norm.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(c) Stress $L^2$-norm relative error norm.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(e) Displacement $L^2$-norm relative error norm.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(b) Strain $max$-norm relative error norm.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(d) Stress $max$-norm relative error norm.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(f) Displacement $max$-norm relative error norm.</figcaption></figure>

<caption>Figure 11: Study effect of DR regarding the number of restarts.</caption>
as the others as the iterations advance.

Regarding $k_{update}$ influence, it has been consistently observed that smaller $k_{update}$ leads to slightly smaller errors for the same database size. On the other hand, smaller updates need more database update iterations (outer iterations of Algorithm 2), which also encompass some preparation operations: metric update, tree-search structures reinitialization, and stiffness matrix assembling and factorisation. Therefore, there might be an intricate trade-off between $k_{update}$ and the number of DD-AS iterations to achieve optimal computational efficiency, which is presumably problem-dependent, and whose investigation is out the scope of this work. Nevertheless, it is worth highlighting that we have chosen to depict the database size as a privileged measure of computational complexity, rather than the number of DD-AS iterations, since it is reasonable to conjecture that most of the computational burden is due to the CH computations. Indeed, we demonstrate numerically that such assumption is true in a final example in Section 5, where speed-up comparisons are presented. Moreover, it is worth mentioning that the database update is trivially scalable since CH computations can be carried out in parallel, which is not the case of the remaining preparation operations. In this context, the number of available cores comes into play to a well-advised choice of $k_{update}$. Our current implementation is serial, therefore analysing aforementioned interplay is out of the scope of this work.

As for the $U$-score, we practically see no difference concerning $G$ for the $L^2$-norm and for the final errors in the $\max$-norm. However, we can observe a slightly better performance in the $\max$-norm for the medium-sized databases for the mesh M2. Noteworthy, $U$ curve is relatively smother than $G$ curve for $\max$-norm, which demonstrates that updates are done with a more equilibrated distribution of points along the material manifold. For the mesh M3, the $U$-score has a slightly poorer performance for the $L^2$-norm, but has an equivalent performance than $G_v$ for the $\max$-norm. Finally, it is important to remark that it achieves better accuracy for the M3 than for M2 for all scores, especially for the $\max$-norm. As a result, fixing a target accuracy level at $\max$-norm, the required database size per number of elements ratio decreases as the problem size increases.
<figure><img src="image_24.png" /><figcaption>(a) $\bar{e}$ comparison for M2 and $k_{update} = 40$.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(c) $\bar{e}$ comparison for M2 and $k_{update} = 80$.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(e) $\bar{e}$ comparison for M3 and $k_{update} = 40$.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(b) $\bar{e}^{max}$ comparison for M2 and $k_{update} = 40$.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>(d) $\bar{e}^{max}$ comparison for M2 and $k_{update} = 80$.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>(f) $\bar{e}^{max}$ comparison for M3 and $k_{update} = 40$.</figcaption></figure>

<caption>Figure 12: Study for different ranking strategies.</caption>
793
794

5. Unified numerical example: speed-up comparisons on the holed plate

795
796
797
798
799
800
801
802
803
804
805
806
807
808

So far we have deliberately chosen to validate and compare several novel ingredients between them in a step-by-step, incremental fashion. Little emphasis has been given to our first goal application, which is to show how much can accelerate FE² computations. In this section, we aim to assess more in-depth the performance gains of the best-performing algorithm combinations retained from the preliminary studies. Although all time results are presented in relative terms, which decreases the hardware-dependent influence, for the sake of completeness, we report the machine specification used in all simulations: Laptop Dell Precision 7680, 13th Gen Intel® Core™ i9-13950HX × 32, RAM 32Gb, running on Ubuntu 22.04.5 LTS with FEniCS 2019.1.0. To promote a variety of numerical examples, since Cook's membrane problem has been extensively addressed in the previous sections, we return to the holed plate example as a model problem as already described in Section 2.4.

809
810
811
812
813
814
815
816
817
818
819
820
821
822

First, it is worth to reiterate the main conclusions from the previous tests: i) disregarding any boosting strategy and using the equivalent of the G-score, the general DD-AS framework has been posed in Section 2.4 for the holed plate example, showing its viability, but also some limitations; ii) in Section 3.5, the boosting combination ASLD (= APE + SMM + LCR + DR) has outperformed other combinations in terms of accuracy for the standard DDCM scenario (static database), followed by ASL, which is computationally cheaper and presents less oscillations; iii) in Section 4.3.1, all boosting strategies have been reevaluated in the full DD-AS setting using the G-score, leading to the very same most efficient methods, now further improved with TR for the unlabelled dataset selection yielding to the ATSLD (= ASLD + TR) and ATSL (= ASL + TR) strategies; iv) in Section 4.3.2 the scores performances have been assessed (ATSL enabled), showing very subtle differences between them.

823
824
825
826
827
828
829
830
831
832
833

For the remainder of this section, we retain the *G*-score and *E*-score with ATSL, to be compared with the corresponding FE² simulations. Comparisons are regarding the ratio of CH calls and the ratio of overall runtime (speed-up), FE² over DD-AS in both cases. As stopping criteria, we choose a fixed value 0.1% for the tolerance *tolout* for the Algorithm 2, which allows to focus the comparison only in terms of computational cost. For simplicity, we have fixed $k_{update} \in \{20, 30, 50\}$ for meshes M1, M2 and M3 (increases with the mesh refinement), respectively, to keep similar number of DD-AS updates across the different meshes, although it has been shown in the previous examples that the influence of $k_{update}$ is quite marginal, the database size being the most relevant indicator for both performance and accuracy.

834
835

Regarding the ratio of CH calls, we can see in Figure 13(a) (solid lines) that the gains of using DD-AS with the *G*-score approach over FE² increases
<figure><img src="image_30.png" /><figcaption>(a) Ratio of microscale solver calls.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>(b) Runtime ratio.</figcaption></figure>

<caption>Figure 13: Assessing speed-up gains for per load step.</caption>

with the mesh refinement. In contrast, it has been observed that the number of Newton-Raphson calls needed to converge $FE^2$ is practically independent of the mesh size, which means that the number of CH calls scales practically linearly with problem size for the ground-truth case. Also, DD-AS performs even better for the last load steps, given the reuse of previous states on the dataset. Note that the loading campaign for this problem is such that the imposed displacement changes direction after the third step while keeping one the components magnitude (see Figure 2(b)), which explains the aforementioned behaviour. As for the runtime speed-up, Figure 13(b) (solid lines) follows the very same trends as for the CH calls ratio, however, with decreased runtime ratios. This proves both DD-AS and $FE^2$ are dominated by CH calls costs, with the remaining overhead computational burden being a little higher in DD-AS than in $FE^2$.

As for the *E*-score results, we can see in Figure 13(a) (dashed lines) that the performances are in general superior to the *G*-score for finer meshes,
with the best ratios of more than 100× for the former against about 70× for the latter. However, when compared in terms of runtime in Figure 13(b) (dashed lines), speed-ups are only slightly better compared to the *G*-score. We attribute this to increased number of internal DDCM loops, which encompass more nearest-neighbour searches, given that outer DD-AS loops are reduced as the CH call ratio increases.

<figure><img src="image_32.png" /><figcaption>Figure 14: Strain components scatter plots of “useful” states (related to $z^* \in Z_D$) and “unused” states (related to $D_L$), in blue and grey, respectively.</figcaption></figure>

One reasonable way of assessing the performance of a given score is to check the “unused” CH calls, i.e., material states that have been appended to the dataset, possibly used during iterations, but that, at the end of the simulation, are not the nearest neighbour of any of the mechanical states in the converged solution. In Figure 14, we can see that by comparing the scatter plots views of the sampled strain cloud, the *G*-score and *E*-score has similar of “unused” points. Note that the results are plotted for the last load step, hence, regions exhibiting low strain magnitudes display nearly identical densities since they have been useful for the previous steps.

Concerning the spatial variation of mechanical points that have led to
material updates, Figure 15 shows in red the elements that have led to at least one CH call. Clearly, the *E*-score needed fewer points than *G*-score, although it has covered essentially the same regions. Deepening the analysis, Figure 16 provides the number of CH calls per element. In both cases, regions bordering the holes have clearly requested more CH calls, which is natural given the stress concentrations. In general, the *G*-score queried more points: 63× against 44× (E-entropy) for the most evaluated element. It should be noticed that although these extreme points showed slightly higher evaluations than the FE² case, which is 38× for each Gauss points in this example, the vast majority of points needed less than 5 CH calls, if any. All in all, the results of Figures 15 and 16 reinforces the differences already observed on the CH calls ratios of Figures 13(a).

Finally, as expected from the previous case, this example also demonstrates the tight superiority of entropy over the greedy method, which is also translated to runtime speed-up.

<figure><img src="image_33.png" /><figcaption>Figure 15: DD-AS update flag per element: i) (red) at least one DD-AS update; ii) (grey) none DD-AS update.</figcaption></figure>
<figure><img src="image_34.png" /><figcaption>Figure 16: DD-AS update counter per element.</figcaption></figure>

## 6. Concluding Remarks

In this work, we have proposed an on-the-fly adaptive sampling approach for the (constitutive model-free) data-driven computational mechanics, named DD-AS method. As an application of such framework, we have addressed multiscale solid mechanics, classically solved by the FE². As a result, we have demonstrated through several numerical examples that the DD-AS outperforms FE². Particularly for the holed plate case with finer mesh, the gains range from 20× to 100× considering the CH calls ratio, and by 15× to 60× regarding the runtime, depending on the load step and specific ranking score. It is worth highlighting that higher speed-ups are expected as the problem size increases.

More specifically, to obtain the aforementioned results, several techniques have been proposed to sharpen the numerical performances. These strategies can be split into two groups: i) DDCM boosting and ii) adaptive selection via ranking scores. For the first group, we have introduced three novel techniques, namely: Accelerated Projection-to-Equilibrium (APE), Sylvester Mean Metric (SMM) and Defective Restarting (DR). In addition, the so-called Locally Convex Reconstruction method He and Chen (2020) has been employed in synergy with the other proposed techniques. In particular, as shown through extensive numerical tests, combining all methods leads to better performances than isolated techniques for the DDCM, and more importantly, they boost the algorithms used in the selection step. For the second group of techniques, we have introduced one novel technique to obtain more efficient unlabelled data, namely Tangential Reconstruction (TR), and also four novel DDCM-distance-based scores: Greedy (G), volumetric Greedy ($G_v$), Uncertainty (U), Entropy (E). Although formulated in very distinct manners, reported performances are not drastically different at first sight. With a more refined analysis, our recommendation favours
the E-score with little margin. In practice, depending on the specific example and the error metric regarded, the four scores can be ranked differently, but we expect little difference margins. Notwithstanding, the choice of any of the four proposed scores drastically outperforms other possible classical choices such as the random selection and K-means clustering (see Appendix Appendix C).

It is worth mentioning that although special attention has been taken to provide as accurate as possible runtime speed-up values (same finite element solvers, same machine, etc), results reported are still implementation dependent. In this sense, the ratio of CH calls value is a more viable performance metric, even though the user should be conscious of the overhead computational costs of DDCM, which are minor but non-negligible. In this context, parallel computing can be used to further exploit the beneficial embarrassingly parallel properties of the CH problem. Moreover, we have not explored further selection scores, as its adaptation from the active learning literature, although interesting, it is not always evident.

Finally, the proposed method can be applied out of the box to more challenging contexts as finite strains and three-dimensional problems, which are subject of future investigation. To decrease the problem-dependence and score-dependence effects on the speed-up levels, the automatic choice of ranking scores and its combinations could also be further explored to obtain a fine-tuning for the DD-AS performance. Extension to inelastic material is not trivial, although of utmost practical importance. The difficulty arises from the history dependence, requiring either to store microscopic states along with material data (Chaouch and Yvonnet, 2024), or to apply the method to full histories (Gorgogianni et al., 2023).

## 7. Acknowledgments

F. Rocha, A. Leygue and L. Stainier gratefully acknowledge financial support from ANR through project D3MecA (ANR-19-CE46-0012). This work was partially funded by the "Support for research for newly appointed Associate Professors" (Faculté des Sciences et Technologie, Université Paris-Est Créteil, France).

## Appendix A. Error definitions

The relative errors used to compare DD-AS results and FE² results throughout this work are exhaustively defined as follows.

* L²-based relative errors:

$$e_X = \frac{\|X - X^{FE^2}\|_{L^2}}{\|X^{FE^2}\|_{L^2}} \quad \text{for } X \in \{\varepsilon, \sigma, \mathbf{u}\}. \quad (A.1)$$

Such errors are coined to give a global error measure.
• Maximum-norm relative errors:

$$e_X^{\max} = \frac{\|X - X^{\text{FE}2}\|_{\text{max}}}{\|X^{\text{FE}2}\|_{\text{max}}} \quad \text{for } X \in \{\varepsilon, \sigma\}, \qquad (A.2)$$

where $(\|\star\|_{\text{max}} = \max_j \|(\hat{\star})_i\|_2)$, with $j$ running over all integration points. Such errors are coined to give a local error measure. It is worth noticing that we do not define max-norm for displacements, since this error, even evaluated in a pointwise manner, would provide global information.

• Geometric average of errors:

$$\bar{e} = \sqrt[3]{e_u e_\sigma e_\varepsilon}, \qquad (A.3a)$$

$$\bar{e}^{\max} = \sqrt[2]{e_{\varepsilon}^{\max} e_{\sigma}^{\max}}, \qquad (A.3b)$$

When comparing individual errors we may see opposite trends (e.g. stress errors decreases while strain errors increases). To combine errors, we have assumed a geometric average, since errors for each variable can have different orders of magnitude. Therefore, averaged errors yield to a unified, easier-to-compare indexes.

## Appendix B. Miscellaneous proofs for strategies of Section 3

### Appendix B.1. Sylvester Mean Metric proof

It is straightforward to see that the original problem (22) proposed in Section 3.2 can be equivalently rewritten as

$$\min_{S \in \mathbb{R}_{sym}^{N_{\varepsilon} \times N_{\varepsilon}}} \left\{ \mathcal{J}(S) := \frac{1}{2} \|\tilde{\mathbf{D}}_{\varepsilon}S - \tilde{\mathbf{D}}_{\sigma}\|_{F}^{2} \right\}, \qquad (B.1)$$

where $\|\cdot\|_F$ denotes the Frobenius norm of matrices and $\tilde{\mathbf{D}}_\varepsilon, \tilde{\mathbf{D}}_\sigma \in \mathbb{R}^{N_d \times N_\varepsilon}$ are matrices built stacking row-wise vectors $(\hat{\varepsilon}_i - \bar{\varepsilon})$ and $(\hat{\sigma}_i - \bar{\sigma})$, respectively, for $i = 1, \dots, N_d$, given a dataset $D$ with $N_d$ instances. Evaluating the stationary points of the functional $\mathcal{J}$ via the Gâteaux derivative, we get

$$\delta \mathcal{J}(S; \delta S) = \frac{d \mathcal{J}(S + \tau \delta S)}{d \tau} \bigg|_{\tau=0} = -(\tilde{\mathbf{D}}_{\varepsilon} \delta S) \cdot \tilde{\mathbf{D}}_{\sigma} + (\tilde{\mathbf{D}}_{\varepsilon} \delta S) \cdot (\tilde{\mathbf{D}}_{\varepsilon} S) \\ = \delta S \cdot (\tilde{\mathbf{D}}_{\varepsilon}^{T} \tilde{\mathbf{D}}_{\varepsilon} S - \tilde{\mathbf{D}}_{\varepsilon}^{T} \tilde{\mathbf{D}}_{\sigma}) = 0, \quad \forall \delta S \in \mathbb{R}_{sym}^{N_{\varepsilon} \times N_{\varepsilon}}. \qquad (B.2)$$

Hence, $\text{Sym}(\tilde{\mathbf{D}}_{\varepsilon}^{T} \tilde{\mathbf{D}}_{\varepsilon} S - \tilde{\mathbf{D}}_{\varepsilon}^{T} \tilde{\mathbf{D}}_{\sigma}) = 0$ should hold, yielding to

$$S \underbrace{\left(\frac{1}{N_d} \tilde{\mathbf{D}}_{\varepsilon}^{T} \tilde{\mathbf{D}}_{\varepsilon}\right)}_{:= \Sigma_{\varepsilon\varepsilon}} + \left(\frac{1}{N_d} \tilde{\mathbf{D}}_{\varepsilon}^{T} \tilde{\mathbf{D}}_{\varepsilon}\right) S = \frac{1}{N_d} \tilde{\mathbf{D}}_{\sigma}^{T} \tilde{\mathbf{D}}_{\varepsilon} + \underbrace{\frac{1}{N_d} \tilde{\mathbf{D}}_{\varepsilon}^{T} \tilde{\mathbf{D}}_{\sigma}}_{:= \Sigma_{\varepsilon\sigma}}, \qquad (B.3)$$
where the $N_d$ factor has been incorporated for the sake of convenience such that the above defined matrices are the blocks of the empirical covariance matrix $\Sigma = [\Sigma_{\epsilon\epsilon} \Sigma_{\epsilon\sigma} \Sigma_{\sigma\sigma}]$ for the dataset D. Finally, we get $\mathbb{C}$ as the solution of the Sylvester equation in (B.3).

Appendix B.2. Accelerated Projection-to-Equilibrium proof
Let $\tilde{P}_D$ the global version of (19) for $\theta = 1/2$ be defined as below

$$ \tilde{P}_D : \mathcal{Z}^E \to \mathcal{Z} \\ \mathbf{z} \mapsto P_E^{-1}\mathbf{z} + \mathbb{G}_{1/2}(\mathbf{z} - P_E^{-1}\mathbf{z}). \qquad (B.4) $$

Note that $P_E^{-1}$ should be understood in the sense $P_E^{-1}\mathbf{z}^{(k+1)} = P_D\mathbf{z}^{(k)} := \mathbf{z}^*$. Explicitly expressing $\tilde{\mathbf{z}}^* := \tilde{P}_D\mathbf{z}^{(k+1)}$, we get

$$ \tilde{\mathbf{z}}^* = P_E^{-1}\mathbf{z}^{(k+1)} + \mathbb{G}_{1/2}(\mathbf{z}^{(k+1)} - P_E^{-1}\mathbf{z}^{(k+1)}) \\ = \mathbf{z}^* + \mathbb{G}_{1/2}(\mathbf{z}^{(k+1)} - \mathbf{z}^*). \qquad (B.5) $$

Using (9) and decoupling the strain/stress components of $\tilde{\mathbf{z}}^*$, we have

$$ \tilde{\varepsilon}^* = \varepsilon^* + \frac{1}{2}(\varepsilon^{(k+1)} - \varepsilon^*) + \frac{1}{2}\mathbb{C}^{-1}(\boldsymbol{\sigma}^{(k+1)} - \boldsymbol{\sigma}^*) \\ = \frac{1}{2}(\nabla^s \mathbf{u}^{(k+1)} + \varepsilon^* + \nabla^s \boldsymbol{\lambda}^{(k+1)}), \qquad (B.6) $$

$$ \tilde{\boldsymbol{\sigma}}^* = \boldsymbol{\sigma}^* + \frac{1}{2}(\boldsymbol{\sigma}^{(k+1)} - \boldsymbol{\sigma}^*) + \frac{1}{2}\mathbb{C}(\varepsilon^{(k+1)} - \varepsilon^*) \\ = \frac{1}{2}(\boldsymbol{\sigma}^{(k+1)} + \boldsymbol{\sigma}^* + \mathbb{C}(\varepsilon^{(k+1)} - \varepsilon^*)) \\ = \boldsymbol{\sigma}^* + \frac{1}{2}\mathbb{C}(\varepsilon^{(k+1)} - \varepsilon^* + \nabla^s \boldsymbol{\lambda}^{(k+1)}) \qquad (B.7) $$

Recalling (17), $\tilde{\mathbf{z}}^{(k+1)} = P_E\tilde{\mathbf{z}}^*$, we need to solve the system (8) and perform the update (9) replacing $\mathbf{z}^*$ by $\tilde{\mathbf{z}}^*$. From (8a), we get

$$ \begin{aligned} 2\langle\mathbb{C}\nabla^s \tilde{\mathbf{u}}^{(k+1)}, \nabla^s \delta\mathbf{u}\rangle &= 2\langle\tilde{\mathbb{C}}\varepsilon^*, \nabla^s \delta\mathbf{u}\rangle \\ &= \langle\mathbb{C}(\nabla^s \mathbf{u}^{(k+1)} + \varepsilon^* + \boldsymbol{\lambda}^{(k+1)}), \nabla^s \delta\mathbf{u}\rangle \\ &= \langle\mathbb{C}(2\mathbf{u}^{(k+1)} + \boldsymbol{\lambda}^{(k+1)}), \nabla^s \delta\mathbf{u}\rangle \quad \forall \delta\mathbf{u} \in \mathcal{V}. \end{aligned} \qquad (B.8) $$

As result, $\tilde{\mathbf{u}}^{(k+1)} = \mathbf{u}^{(k+1)} + \frac{1}{2}\boldsymbol{\lambda}^{(k+1)} + \mathbf{c}$, where $\mathbf{c}$ is a constant field. As $\tilde{\mathbf{u}}^{(k+1)}, \mathbf{u}^{(k+1)} \in \mathcal{U}$ and $\boldsymbol{\lambda}^{(k+1)} \in \mathcal{V}$, we conclude that $\mathbf{c} \in \mathcal{V}$. Taking into account the homogeneous Dirichlet boundary conditions, $\mathbf{c}$ vanishes, yielding to $\tilde{\mathbf{u}}^{(k+1)} = \mathbf{u}^{(k+1)} + \frac{1}{2}\boldsymbol{\lambda}^{(k+1)}$, uniquely determined.
On the other hand, equating the external virtual powers of (8b) for both $\mathbf{z}^*$ and $\tilde{\mathbf{z}}^*$, we get

$$ 
\begin{aligned}
\langle \mathbb{C} \nabla^s \tilde{\lambda}^{(k+1)}, \nabla^s \delta \lambda \rangle + \langle \tilde{\sigma}^*, \nabla^s \delta \lambda \rangle &= \langle \mathbb{C} \nabla^s \lambda^{(k+1)}, \nabla^s \delta \lambda \rangle + \langle \sigma^*, \nabla^s \delta \lambda \rangle \\
\Rightarrow \langle \mathbb{C} \nabla^s \tilde{\lambda}^{(k+1)}, \nabla^s \delta \lambda \rangle &= \langle \sigma^* + \mathbb{C} \nabla^s \lambda^{(k+1)} - \tilde{\sigma}^*, \nabla^s \delta \lambda \rangle \\
&= \frac{1}{2} \langle \mathbb{C}(\nabla^s \lambda^{(k+1)} + \varepsilon^* - \varepsilon^{(k+1)}), \nabla^s \delta \lambda \rangle \\
&= \frac{1}{2} \langle \mathbb{C} \nabla^s \lambda^{(k+1)}, \nabla^s \delta \lambda \rangle, \quad \forall \delta \lambda \in \mathcal{V},
\end{aligned}
\quad (B.9) $$

where (8a) has been used to see that the term $\varepsilon^* - \varepsilon^{(k+1)}$ vanishes in the sense of the variational equation. As result, $\tilde{\lambda}^{(k+1)} = \frac{1}{2}\lambda^{(k+1)}$ (the same reasoning done in (B.8) is used to see that the constant field vanishes). Finally, using (9), the improved updates read as

$$ 
\begin{aligned}
\tilde{\varepsilon}^{(k+1)} &= \nabla^s \tilde{\mathbf{u}}^{(k+1)} = \nabla^s \left( \mathbf{u}^{(k+1)} + \frac{1}{2} \lambda^{(k+1)} \right) \\
&= \varepsilon^{(k+1)} + \frac{1}{2} \mathbb{C}^{-1} (\boldsymbol{\sigma}^{(k+1)} - \boldsymbol{\sigma}^*),
\end{aligned}
\quad (B.10) $$

$$ 
\begin{aligned}
\tilde{\boldsymbol{\sigma}}^{(k+1)} &= \tilde{\boldsymbol{\sigma}}^* + \mathbb{C} \nabla^s \tilde{\boldsymbol{\lambda}}^{(k+1)} \\
&= \frac{1}{2}(\boldsymbol{\sigma}^{(k+1)} + \boldsymbol{\sigma}^*) + \frac{1}{2}\mathbb{C}(\varepsilon^{(k+1)} - \varepsilon^*) + \frac{1}{2}\mathbb{C} \nabla^s \boldsymbol{\lambda}^{(k+1)} \\
&= \boldsymbol{\sigma}^{(k+1)} + \frac{1}{2}\mathbb{C}(\varepsilon^{(k+1)} - \varepsilon^*).
\end{aligned}
\quad (B.11) $$

Noticing that $\mathbb{G}_{1/2} = \frac{1}{2} \begin{pmatrix} \mathbb{O} & \mathbb{C}^{-1} \\ \mathbb{C} & \mathbb{O} \end{pmatrix} + \begin{pmatrix} \mathbb{I} & \mathbb{O} \\ \mathbb{O} & \mathbb{I} \end{pmatrix}$, the matrix forms of (B.10) and (B.11) read

$$ 
\begin{aligned}
\tilde{\mathbf{z}}^{(k+1)} &= \begin{bmatrix} \varepsilon^{(k+1)} \\ \boldsymbol{\sigma}^{(k+1)} \end{bmatrix} + \frac{1}{2} \begin{bmatrix} \mathbb{O} & \mathbb{C}^{-1} \\ \mathbb{C} & \mathbb{O} \end{bmatrix} \begin{bmatrix} \varepsilon^{(k+1)} - \varepsilon^* \\ \boldsymbol{\sigma}^{(k+1)} - \boldsymbol{\sigma}^* \end{bmatrix} \\
&= \mathbf{z}^{(k+1)} + \frac{1}{2}(2\mathbb{G}_{1/2} - \mathbb{I})(\mathbf{z}^{(k+1)} - \mathbf{z}^*) \\
&= \frac{1}{2}(\mathbf{z}^{(k+1)} + \mathbf{z}^*) + \mathbb{G}_{1/2}(\mathbf{z}^{(k+1)} - \mathbf{z}^*),
\end{aligned}
\quad (B.12) $$

which corresponds to the update formula presented in (21).

## Appendix C. Sanity-check: alternative non distance-based scores

In this last Appendix, we provide alternatives to the ranking scores proposed in the paper. Indeed, it is perfectly possible to propose DD-AS approaches that disregard the DDCM distance for selecting the updates. Two of them were also studied:
1. **Random (R):** The simplest possible strategy is to select $\tilde{D}_U \subset D_U$ by randomly picking a specified number of states of the pool.

2. **K-means (K):** Another natural choice is to create $\tilde{D}_U$ by the centroids of a K-means clustering of $D_U$. Note that, for the K-means case, $\tilde{D}_U \not\subset D_U$, but only constructed from $D_U$. We test two variants:

(a) $K_\epsilon$: Since only strains are needed for labelling, the clustering can only rely on strains.

(b) $K_z$: a clustering using the full strain-stress states is also possible by considering the Euclidean mapping $(\hat{\epsilon}, \hat{\sigma}) \mapsto (\mathbb{C}^{1/2}\hat{\epsilon}, \mathbb{C}^{-1/2}\hat{\sigma})$, which scales magnitudes of strain and stress (Platzer, 2020). The obtained centroids should be transformed back to the original strain-stress space, where only the strain part is considered for labelling.

In this sanity-check numerical example, we aim at contrasting the methods defined above with some already studied strategies. For the sake of simplicity, we only depict the $G$-score that follows approximately the same trend as the {$G_v, U, E$}-score (see Figure 12). It is worth mentioning that the random strategy is clearly non-optimal, but it is normally used as a point of comparison to assess the performance of other methods in the active learning community. On the other hand, the K-means approach is indeed a reasonable strategy for optimally exploring of the strain(-stress) space, as it has been already noticed in the literature Chaouch and Yvonnet (2024).

In Figure C.17, we can see the evolution for the errors $\bar{e}$ and $\bar{e}^{max}$. As expected, approaches $R$, $K_\epsilon$ and $K_z$ are by far less performing than $G$, $G_v$ and $U$, excepting the very first updates regarding error $\bar{e}$. This result reiterates our initial hypotheses, that using DDCM-distance-based selection scores would naturally yield to more efficient methods by taking advantage of the underlying DDCM formulation.

As for the K-means method, we observed that differences between $K_\epsilon$ and $K_z$ are not substantial, such that one can use them interchangeably. The initial enhanced observed trend can be explained by the beneficial effect of exploring more regions of the strain-stress space rather than concentrating in some zones at early stages. Naturally, combination of clustering methods and DDCM-distance-based method can also be employed, but their application is out of the scope of this contribution.
<figure><img src="image_35.png" /></figure>

<figure><img src="image_36.png" /><figcaption>Figure C.17: Sanity-check study comparing Random and K-means strategies.</figcaption></figure>
1011 References

1012 Alnaes, M.S., Blechta, J., Hake, J., Johansson, A., Kehlet, B., Logg, A., Richardson, C.N., Ring, J., Rognes, M.E., Wells, G.N., 2015. The FEniCS project version 1.5. Archive of Numerical Software 3. doi:10.11588/ans.2015.100.20553.

1016 Bahmani, B., Sun, W., 2021. A kd-tree-accelerated hybrid data-driven/model-based approach for poroelasticity problems with multi-fidelity multi-physics data. Computer Methods in Applied Mechanics and Engineering 382, 113868. doi:10.1016/j.cma.2021.113868.

1020 Baratta, I.A., Dean, J.P., Dokken, J.S., Habera, M., Hale, J.S., Richardson, C.N., Rognes, M.E., Scroggs, M.W., Sime, N., Wells, G.N., 2023. Dolfinx: The next generation fenics problem solving environment. URL: https://doi.org/10.5281/zenodo.10447666, doi:10.5281/zenodo.10447666.

1024 Binev, P., Cohen, A., Dahmen, W., DeVore, R., Petrova, G., Wojtaszczyk, P., 2011. Convergence rates for greedy algorithms in reduced basis methods. SIAM Journal on Mathematical Analysis 43, 1457-1472. doi:10.1137/100795772.

1028 Blanco, P., Sánchez, P., de Souza Neto, E., Feijóo, R., 2016. Variational foundations and generalized unified theory of RVE-based multiscale models. Archives of Computational Methods in Engineering 23, 191-253. URL: http://dx.doi.org/10.1007/s11831-014-9137-5, doi:10.1007/s11831-014-9137-5.

1033 Carrara, P., De Lorenzis, L., Stainier, L., Ortiz, M., 2020. Data-driven fracture mechanics. Computer Methods in Applied Mechanics and Engineering 372, 113390. URL: https://www.sciencedirect.com/science/article/pii/S0045782520305752, doi:https://doi.org/10.1016/j.cma.2020.113390.

1038 Chaouch, S., Yvonnet, J., 2024. Unsupervised machine learning classification for accelerating fe2 multiscale fracture simulations. Computer Methods in Applied Mechanics and Engineering 432, 117278. URL: https://www.sciencedirect.com/science/article/pii/S0045782524005346, doi:https://doi.org/10.1016/j.cma.2024.117278.

1044 Conti, S., Müller, S., Ortiz, M., 2020. Data-driven finite elasticity. Archive for Rational Mechanics and Analysis 237, 1-33. doi:10.1007/s00205-020-01490-x.

1047 Costecalde, L., 2023. Data-Driven Model Identification for Hyperelasticity: Mapping the Strain Energy throughout Multiaxial Experiments. Ph.D. thesis. École Centrale de Nantes.
1050

1051

1052

Costecalde, L., Leygue, A., Coret, M., Verron, E., 2023. Data-Driven Identification of hyperelastic models by measuring the strain energy density field. Rubber Chemistry and Technology doi:10.5254/rct-23.386903.

1053

1054

1055

Dalémat, M., Coret, M., Leygue, A., Verron, E., 2019. Measuring stress field without constitutive equation. Mechanics of Materials 136, 103087. doi:10.1016/j.mechmat.2019.103087.

1056

1057

1058

1059

Dandin, H., Leygue, A., Stainier, L., 2024. Graph-based representation of history-dependent material response in the Data-Driven Computational Mechanics framework. Computer Methods in Applied Mechanics and Engineering 419, 116694. doi:10.1016/j.cma.2023.116694.

1060

1061

1062

1063

1064

Dasgupta, S., 2004. Analysis of a greedy active learning strategy, in: Saul, L., Weiss, Y., Bottou, L. (Eds.), Advances in Neural Information Processing Systems, MIT Press. URL: https://proceedings.neurips.cc/paper_files/paper/2004/file/c61fbef63df5ff317aecdc3670094472-Paper.pdf.

1065

1066

1067

1068

1069

1070

Eggersmann, R., Kirchdoerfer, T., Reese, S., Stainier, L., Ortiz, M., 2019. Model-free data-driven inelasticity. Computer Methods in Applied Mechanics and Engineering 350, 81-99. URL: https://www.sciencedirect.com/science/article/pii/S0045782519300878, doi:https://doi.org/10.1016/j.cma.2019.02.016.

1071

1072

1073

1074

1075

1076

Eggersmann, R., Stainier, L., Ortiz, M., Reese, S., 2021. Model-free data-driven computational mechanics enhanced by tensor voting. Computer Methods in Applied Mechanics and Engineering 373, 113499. URL: https://www.sciencedirect.com/science/article/pii/S0045782520306848, doi:https://doi.org/10.1016/j.cma.2020.113499.

1077

1078

1079

1080

Feyel, F., 2003. A multilevel finite element method (fe2) to describe the response of highly non-linear structures using generalized continua. Computer Methods in Applied Mechanics and Engineering 192, 3233-3244. doi:10.1016/S0045-7825(03)00348-7.

1081

1082

Goodfellow, I., Bengio, Y., Courville, A., 2016. Deep Learning. MIT Press. http://www.deeplearningbook.org.

1083

1084

1085

1086

Gorgogianni, A., Karapiperis, K., Stainier, L., Ortiz, M., Andrade, J., 2023. Adaptive goal-oriented data sampling in data-driven computational mechanics. Computer Methods in Applied Mechanics and Engineering 409. doi:10.1016/j.cma.2023.115949.
1087

1088

1089

1090

1091

Gweon, H., Yu, H., 2021. A nearest neighbor-based active learning method and its application to time series classification. Pattern Recognition Letters 146, 230-236. URL: https://www.sciencedirect.com/science/article/pii/S0167865521001033, doi:https://doi.org/10.1016/j.patrec.2021.03.016.

1092

1093

1094

1095

1096

1097

He, Q., Chen, J.S., 2020. A physics-constrained data-driven approach based on locally convex reconstruction for noisy database. Computer Methods in Applied Mechanics and Engineering 363, 112791. URL: https://www.sciencedirect.com/science/article/pii/S0045782519306838, doi:https://doi.org/10.1016/j.cma.2019.112791.

1098

1099

1100

1101

1102

Hill, R., 1965. A self-consistent mechanics of composite materials. Journal of the Mechanics and Physics of Solids 13, 213-222. URL: http://www.sciencedirect.com/science/article/pii/0022509665900104, doi:http://dx.doi.org/10.1016/0022-5096(65)90010-4.

1103

1104

1105

1106

1107

Ibañez, R., Abisset-Chavanne, E., Aguado, J.V., Gonzalez, D., Cueto, E., Chinesta, F., 2016. A Manifold Learning Approach to Data-Driven Computational Elasticity and Inelasticity. Archives of Computational Methods in Engineering 25, 47-57. URL: https://hal.science/hal-03653378, doi:10.1007/s11831-016-9197-9.

1108

1109

1110

1111

Kanno, Y., 2019. Mixed-integer programming formulation of a data-driven solver in computational elasticity. Optimization Letters 13, 1505-1514. URL: http://dx.doi.org/10.1007/s11590-019-01409-w, doi:10.1007/s11590-019-01409-w.

1112

1113

1114

Karapiperis, K., Stainier, L., Ortiz, M., Andrade, J.E., 2020. Data-driven multiscale modeling in mechanics. Journal of the Mechanics and Physics of Solids, 104239doi:10.1016/j.jmps.2020.104239.

1115

1116

1117

1118

1119

Kirchdoerfer, T., Ortiz, M., 2016. Data-driven computational mechanics. Computer Methods in Applied Mechanics and Engineering 304, 81-101. URL: https://www.sciencedirect.com/science/article/pii/S0045782516300238, doi:https://doi.org/10.1016/j.cma.2016.02.001.

1120

1121

1122

1123

1124

Kirchdoerfer, T., Ortiz, M., 2017. Data driven computing with noisy material data sets. Computer Methods in Applied Mechanics and Engineering 326, 622-641. URL: https://www.sciencedirect.com/science/article/pii/S0045782517304012, doi:https://doi.org/10.1016/j.cma.2017.07.039.
1125

1126

1127

1128

1129

Leygue, A., Coret, M., Réthoré, J., Stainier, L., Verron, E., 2018. Data-based derivation of material response. Computer Methods in Applied Mechanics and Engineering 331, 184–196. URL: https://www.sciencedirect.com/science/article/pii/S0045782517307156, doi:https://doi.org/10.1016/j.cma.2017.11.013.

1130

1131

1132

Li, M., Sethi, I., 2006. Confidence-based active learning. IEEE Transactions on Pattern Analysis and Machine Intelligence 28, 1251–1261. doi:10.1109/TPAMI.2006.156.

1133

1134

1135

1136

Linden, L., Klein, D.K., Kalina, K.A., Brummund, J., Weeger, O., Kästner, M., 2023. Neural networks meet hyperelasticity: A guide to enforcing physics. Journal of the Mechanics and Physics of Solids 179, 105363. doi:10.1016/j.jmps.2023.105363.

1137

Mandel, J., 1972. Plasticité classique et viscoplasticité. Springer-Verlag.

1138

1139

1140

1141

1142

1143

Marenić, E., Seychal, G., Passieux, J.C., 2022. Data driven approach in multiphysics framework: Application to coupled electro-mechanical problems. Computer Methods in Applied Mechanics and Engineering 395, 114959. URL: https://www.sciencedirect.com/science/article/pii/S0045782522002195, doi:https://doi.org/10.1016/j.cma.2022.114959.

1144

1145

1146

Nguyen, L.T.K., Keip, M.A., 2018. A data-driven approach to nonlinear elasticity. Computers & Structures 194, 97–115. doi:10.1016/j.compstruc.2017.07.031.

1147

1148

1149

1150

1151

1152

Nguyen, L.T.K., Rambausek, M., Keip, M.A., 2020. Variational framework for distance-minimizing method in data-driven computational mechanics. Computer Methods in Applied Mechanics and Engineering 365, 112898. URL: https://www.sciencedirect.com/science/article/pii/S0045782520300815, doi:https://doi.org/10.1016/j.cma.2020.112898.

1153

1154

1155

1156

1157

Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., Duchesnay, E., 2011. Scikit-learn: Machine learning in Python. Journal of Machine Learning Research 12, 2825–2830.

1158

1159

1160

1161

Platzer, A., 2020. Finite strain data-driven computational mechanics: From tailored data to adaptive solvers for multiscale simulations. Theses. École centrale de Nantes. URL: https://tel.archives-ouvertes.fr/tel-03164128.
1162

1163

1164

Platzer, A., Leygue, A., Stainier, L., Ortiz, M., 2021. Finite element solver for data-driven finite strain elasticity. Computer Methods in Applied Mechanics and Engineering.

1165

1166

1167

1168

1169

Prume, E., Reese, S., Ortiz, M., 2023. Model-free data-driven inference in computational mechanics. Computer Methods in Applied Mechanics and Engineering 403, 115704. URL: https://www.sciencedirect.com/science/article/pii/S0045782522006594, doi:https://doi.org/10.1016/j.cma.2022.115704.

1170

1171

1172

Rocha, F., 2023a. ddfenics: a fenics-based (model-free) data-driven computational mechanics implementation. URL: https://doi.org/10.5281/zenodo.7646226, doi:10.5281/zenodo.7646226.

1173

1174

1175

1176

Rocha, F., 2023b. micmacsfenics: a fenics-based implementation of two-level finite element simulations (fe2) using computational homogenization. URL: https://doi.org/10.5281/zenodo.7643777, doi:10.5281/zenodo.7643777.

1177

1178

1179

1180

1181

Rocha, F., Deparis, S., Antolin, P., Buffa, A., 2023. Deepbnd: A machine learning approach to enhance multiscale solid mechanics. Journal of Computational Physics 479, 111996. URL: https://www.sciencedirect.com/science/article/pii/S0021999123000918, doi:https://doi.org/10.1016/j.jcp.2023.111996.

1182

1183

1184

1185

1186

Salahshoor, H., Ortiz, M., 2023. Model-free data-driven viscoelasticity in the frequency domain. Computer Methods in Applied Mechanics and Engineering 403, 115657. URL: https://www.sciencedirect.com/science/article/pii/S0045782522006120, doi:https://doi.org/10.1016/j.cma.2022.115657.

1187

1188

1189

1190

1191

1192

1193

Schröder, J., Wick, T., Reese, S., Wriggers, P., Müller, R., Kollmannsberger, S., Kästner, M., Schwarz, A., Igelbüscher, M., Viebahn, N., Bayat, H.R., Wulfinghoff, S., Mang, K., Rank, E., Bog, T., D'Angella, D., Elhaddad, M., Hennig, P., Düster, A., Garhuom, W., Hubrich, S., Walloth, M., Wollner, W., Kuhn, C., Heister, T., 2021. A selection of benchmark problems in solid mechanics and applied mathematics. Archives of Computational Methods in Engineering 28, 713-751.

1194

1195

1196

1197

1198

de Souza Neto, E., Blanco, P., Sánchez, P., Feijóo, R., 2015. An RVE-based multiscale theory of solids with micro-scale inertia and body force effects. Mechanics of Materials 80, 136-144. URL: http://www.sciencedirect.com/science/article/pii/S0167663614001872, doi:http://dx.doi.org/10.1016/j.mechmat.2014.10.007.
1199

1200

1201

1202

1203

1204

Ulloa, J., Stainier, L., Ortiz, M., Andrade, J.E., 2024. Data-driven micromorphic mechanics for materials with strain localization. Computer Methods in Applied Mechanics and Engineering 429, 117180. URL: https://www.sciencedirect.com/science/article/pii/S0045782524004365, doi:https://doi.org/10.1016/j.cma.2024.117180.

1205

1206

1207

1208

1209

Wu, D., Lin, C.T., Huang, J., 2019. Active learning for regression using greedy sampling. Information Sciences 474, 90-105. URL: https://www.sciencedirect.com/science/article/pii/S0020025518307680, doi:https://doi.org/10.1016/j.ins.2018.09.060.

1210

1211

1212

1213

1214

1215

1216

Xu, R., Yang, J., Yan, W., Huang, Q., Giunta, G., Belouet-tar, S., Zahrouni, H., Zineb, T.B., Hu, H., 2020. Data-driven multiscale finite element method: From concurrence to separation. Computer Methods in Applied Mechanics and Engineering 363, 112893. URL: https://www.sciencedirect.com/science/article/pii/S004578252030075X, doi:https://doi.org/10.1016/j.cma.2020.112893.

1217

1218

1219

1220

1221

1222

Zhou, L., Liu, W., Mao, Y., Hou, S., 2024. Model-free chemomechanical interfaces: History-dependent damage under transient mass diffusion. Computer Methods in Applied Mechanics and Engineering 431, 117286. URL: https://www.sciencedirect.com/science/article/pii/S0045782524005425, doi:https://doi.org/10.1016/j.cma.2024.117286.

1223

1224

1225

1226

1227

1228

1229

Zlatić, M., Rocha, F., Stainier, L., Canadija, M., 2024. Data-driven methods for computational mechanics: A fair comparison between neural networks based and model-free approaches. Computer Methods in Applied Mechanics and Engineering 431, 117289. URL: https://www.sciencedirect.com/science/article/pii/S0045782524005450, doi:https://doi.org/10.1016/j.cma.2024.117289.