# A k-medoids-based partitioned method for computational homogenization of heterogeneous materials

Modesar Shakoor

► To cite this version:

Modesar Shakoor. A k-medoids-based partitioned method for computational homogenization of heterogeneous materials. Computers & Structures, 2025, 316, pp.107875. 10.1016/j.compstruc.2025.107875. hal-04861315v2

HAL Id: hal-04861315

https://hal.science/hal-04861315v2

Submitted on 24 Jun 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
<figure><img src="image_2.png" /></figure>

# A k-medoids-based partitioned method for computational homogenization of heterogeneous materials☆

Modesar Shakoor

IMT Nord Europe, Institut Mines-Télécom, Univ. Lille, Centre for Materials and Processes, F-59000 Lille, France

## HIGHLIGHTS

* Proposition of a k-medoids-based partitioned finite element squared approach.

* Linearization with cross contributions within clusters and analysis of their effects.

* Solution of fine scale problems after convergence to update internal variables.

* Simplification of the nonlinear solution algorithm.

## ARTICLE INFO

**Keywords:**

Computational homogenization
FE²
Clustering
K-medoids
Reduced-order modeling

## ABSTRACT

Although they are very interesting for structural simulations involving heterogeneous materials, Finite Element (FE) squared approaches, often coined FE², are well-known to require great computational resources. The main challenge is that, for each integration point of the *coarse* structure, a so-called *fine* scale problem should be solved. In this work, a k-medoids-based partitioned FE² approach is proposed to directly tackle this challenge by effectively reducing the number of fine scale solves. At each coarse scale nonlinear iteration, coarse scale integration points are partitioned *a priori* based on the current coarse scale displacement gradient and fine scale internal variables using the k-medoids-based clustering algorithm. Stresses and tangent moduli are computed only for cluster medoids, and are then extended to the remaining non-medoid coarse scale integration points. Results with nonlinear material behavior such as hyperelasticity and elasto-plasticity show that the proposed method is a promising candidate for reducing the computational cost of FE² simulations.

## 1. Introduction

Numerical modeling of heterogeneous materials is well-known to be challenging because the simultaneous discretization of the heterogeneity at the scale of the simulation domain is often computationally intractable. For instance, fiber-reinforced composite structures typically contain millions of fibers. The idea of separating the fine scale of the heterogeneity from the coarse scale of the domain or structure through homogenization theory has been investigated in various studies [11,17]. In this approach, as shown in Fig. 1, the heterogeneity is modeled through infinitely small unit cells or Representative Volume Elements (RVEs), which alleviates the computational burden substantially. Boundary conditions for fine scale problems are provided by point-wise strains from the coarse scale, while stresses and tangent

moduli for the coarse scale problem are provided by homogenized stresses and tangent moduli from the fine scale [11,17].

Even with this scale separation, however, the computational cost remains an issue because a fine scale problem must be solved at each point of the coarse scale domain. For nonlinear problems, in addition, all fine scale computations must happen on-the-fly during the coarse scale simulation with two-way coupling between all solvers (coarse and fine). This is the reason why various approaches can be found in the literature to solve the fine scale problems, even though the coarse scale problem is often discretized and solved using the Finite Element (FE) method [11,17]. In the FE² method [8], fine scale problems are also solved, in their full complexity, with the FE method. Alternative full-field methods such as the fast Fourier transform based numerical method have also been proposed in the literature to solve fine scale problems [12]. Simplifying

<footnote>☆ This work is licensed under a Creative Commons Attribution 4.0 International License.
Email address: modesar.shakoor@imt-nord-europe.fr.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Example of FE² result showing the equivalent stress field: (a) at the coarse scale, (b) at the fine scale for a coarse scale integration point at the top left corner, (c) at the fine scale for a coarse scale integration point between holes B and C.</figcaption></figure>

assumptions have been employed in many works to reduce the computational cost of fine scale problems, for instance using mean-field models such as the Hashin–Shtrikman model [14] or the Mori–Tanaka model [18].

In the recent literature, many researchers have attempted to reduce the computational cost without such assumptions, mainly using reduced-order models such as proper orthogonal decomposition [26] or k-means clustering [10]. Recent advances in the field of machine learning and deep learning have also been exploited [15]. Generally, these techniques consist in generating a database of results computed using a conventional numerical method such as the FE method. A reduced-order model is then trained on this database, so that it can be used to accelerate multiscale simulations. The reduced-order model can be a material genome database [25], a radial numerically explicit potential [9], a thermodynamics-based artificial neural network [16] or a recurrent neural network [24].

An interesting and promising alternative path has recently been proposed in Refs. [2,4]. Instead of attempting to simplify or reduce fine scale problems, this alternative approach called Partitioned FE² (PFE²) in the following consisted in applying model-order reduction directly at the coarse scale to effectively reduce the number of fine scale problems to solve. At each nonlinear iteration of the coarse scale solver, coarse scale

integration points were partitioned using k-means clustering with some *a priori* criterion, and only one fine scale problem was solved per cluster. Stresses were then computed for all coarse scale integration points by using an approximation which avoided solving all fine scale problems. This approach has been shown to provide accurate results with reduced computation times as compared to the FE² method, first for hyperelastic, viscoelastic and elasto-plastic material behaviors [4], and then for problems involving localization and damage [5].

Although a number of issues with the initial version have been addressed in Refs. [4,5], a few challenges remain with the PFE² method and are the object of the present work:

* K-means clustering has been used in Refs. [2,4] and the fine scale problem solved for each cluster was a virtual fine scale problem whose loading conditions corresponded to the average of loading conditions for all fine scale problems within the cluster. This approach required to create and maintain a virtual model with its virtual internal variables for each cluster.

* At each clustering change, virtual internal variables were reassigned between old and new clusters through a mapping algorithm. For each new cluster, this algorithm consisted of looking for the closest old cluster and copying its internal variables.
* Challenging convergence issues were revealed in Ref. [4] and an elaborate nonlinear solution algorithm relying on clustering freezing and iteration restarts was proposed to address them.

In the present work, the following original contributions are proposed to the PFE² method:

* K-means clustering is replaced by k-medoids clustering. This partitioning technique defines as cluster centroid one of the fine scale problems within the cluster (e.g., the cluster's medoid). Consequently, stresses and tangent moduli are computed from actual fine scale problems, which effectively avoids the introduction of a virtual fine scale problem.

* The approximation proposed in Ref. [4] is still employed, but it is shown that a second derivative tensor appears as well as cross contributions between all integration points within a cluster. The influence of these cross contributions on convergence is investigated.

* Internal variables are systematically updated for all fine scale problems, in addition to the ones associated with cluster medoids. At any step, due to clustering changes, any integration point can become a cluster medoid and be solved.

* A simplified nonlinear solution algorithm based on clustering freezing, time-step reduction and fall back to the FE² scheme is proposed to tackle convergence issues. Except for the requirement to freeze clustering, this algorithm is very similar to the one used for the FE² scheme for stiff or strongly non linear problems.

Computational homogenization is introduced in Section 2. The reference FE² method is briefly summarized in Section 3. The new k-medoids-based PFE² method is then detailed in Section 4. Numerical results demonstrating the capabilities of this new method are presented in Section 5. Section 6 finally concludes the paper.

## 2. Governing equations

Under the assumption of scale separation, homogenization consists of introducing at each point X of the coarse scale domain $\Omega_0^M \subset \mathbb{R}^3$ a fine scale domain $\Omega_0^m$. The latter is either a repetitive unit cell or an RVE which represents the heterogeneity only for a small window. For fiber-reinforced composites for instance, this domain might contain a single fiber or up to hundred fibers with some particular distribution carefully calibrated with experimental measurements (e.g., microscopy or tomography).

In the following, the coarse scale problem is first stated, and then the fine scale problem. The symbol X is used for points of the coarse scale, while the symbol Y is used for points of the fine scale. This description is given under the assumption of finite strains, and an example of material model with no internal variable for the fine scale problem is also given to complete the problem statement. Finally, simplified equations are formulated under the assumption of small strains with an elastoplastic material model which features an internal variable at the fine scale.

### 2.1. Coarse scale problem

The coarse scale domain is assumed to be homogeneous, and the unknown displacement $u^M$ is defined as the difference between the deformed coordinate x at some instant $t \in [0, T]$, $T > 0$ and the initial coordinate X. Neglecting body forces, this coarse scale displacement field can be obtained by solving the balance equation:

$$ \nabla_X \cdot \mathbf{P}^M (\mathbf{G}(\mathbf{X})) = \mathbf{0}, \forall \mathbf{X} \in \Omega_0^M, \quad (1) $$

with $G = \nabla_X u^M$, and completing this equation with appropriate boundary conditions.

### 2.2. Fine scale problem

As shown explicitly in Eq. (1), at each point $X \in \Omega_0^M$, the first Piola-Kirchhoff stress tensor $\mathbf{P}^M$ depends on the displacement gradient $G(X)$. In computational homogenization, it is actually computed by solving the following fine scale problem:

$$ \left\{ \begin{array}{l} \nabla_Y \cdot \mathbf{P}^m (\nabla_Y u^m(Y)) = \alpha, \forall Y \in \Omega_0^m, \\ \mathbf{P}^m (\nabla_Y u^m(Y)) \cdot \mathbf{N}^m(Y) = \beta \cdot \mathbf{N}^m(Y), \forall Y \in \partial \Omega_0^m, \\ \frac{1}{|\Omega_0^m|} \int_{\Omega_0^m} u^m(Y) dY = u^M(X), \\ \frac{1}{|\Omega_0^m|} \int_{\Omega_0^m} \nabla_Y u^m(Y) dY = G(X), \end{array} \right. \quad (2) $$

where $\mathbf{P}^m$ is the first Piola-Kirchhoff stress tensor field at the fine scale and $\mathbf{N}^m$ is the outgoing normal vector at the fine scale domain boundary. Details on the derivation of this fine scale problem can be found in the literature (e.g., [8,11,17]).

This fine scale problem is derived from the principle of multiscale virtual power [3], which can rely on periodic boundary conditions as most works on computational homogenization, or on volume constraints. As previously developed for unsteady flows in porous media [22], volume constraints are chosen and are imposed at the fine scale through Lagrange multipliers over the whole fine scale domain. In Eq. (2), indeed, averages of the fine scale displacement field and its gradient are tied to be equal, respectively, to the coarse scale displacement and its gradient through Lagrange multipliers $\alpha \in \mathbb{R}^3$ and $\beta \in \mathbb{R}^{3 \times 3}$. The fine scale problem in Eq. (2) should hence be simultaneously solved to find both Lagrange multipliers and the fine scale displacement field $u^m$.

On the one hand, as acceleration and body forces are neglected, $\alpha$ has solely the role of suppressing rigid body translations at the fine scale. On the other hand, the coarse scale stress tensor can be computed through $\mathbf{P}^M(X) = \beta$. In order to solve the coarse scale problem in Eq. (1), linearization is typically used, which requires to compute $\frac{\partial \mathbf{P}^M}{\partial G}(X) = \frac{\partial \beta}{\partial G}$.

### 2.3. Example of finite strain model

It is reminded that the first Piola-Kirchhoff stress tensor is defined by $\mathbf{P}^m = F^m S^m$, with $F^m$ the deformation gradient tensor defined by $F^m = I + \nabla_Y u^m$, $I$ the second-order identity tensor and $S^m$ the second Piola-Kirchhoff stress tensor. For instance, the isotropic Saint-Venant-Kirchhoff model is given by

$$ S^m = 2\mu^m E^m + \lambda^m \mathrm{tr}(E^m), \quad (3) $$

with $E^m = \frac{1}{2} ( (F^m)^T F^m - I )$ the Green-Lagrange strain tensor and $\mathrm{tr}(A) = \sum_{j=1}^3 A_{jj}$ the trace operator. It should be noted that Lamé parameters $\mu^m = \mu^m(Y)$ and $\lambda^m = \lambda^m(Y)$ are heterogeneous as, for instance, the stiffness of a carbon fiber is typically different from that of a polymer matrix.

For this finite strain model, it is necessary to compute the coarse scale Cauchy stress tensor $\sigma^M$ from the coarse scale first-order Piola-Kirchhoff stress tensor $\mathbf{P}^M$. This is achieved through the formula:

$$ \sigma^M = \frac{1}{\det(F^M)} \mathbf{P}^M (F^M)^T. \quad (4) $$

### 2.4. Example of small strain model

Firstly, under the small strain assumption, the coarse scale problem simplifies to

$$ \nabla_X \cdot \sigma^M (\mathbf{G}(\mathbf{X})) = \mathbf{0}, \forall \mathbf{X} \in \Omega_0^M, \quad (5) $$

where $\sigma^M$ the Cauchy stress tensor and the notation X is still used to homogenize notations with the finite strain model. Secondly, the fine
scale problem simplifies to

$$
\left\{
\begin{array}{l}
\nabla_Y \cdot \boldsymbol{\sigma}^m \left( \varepsilon \left( \nabla_Y \boldsymbol{u}^m(\boldsymbol{Y}) \right) \right) = \alpha, \forall \boldsymbol{Y} \in \Omega_0^m, \\
\boldsymbol{\sigma}^m \left( \varepsilon \left( \nabla_Y \boldsymbol{u}^m(\boldsymbol{Y}) \right) \right) \cdot \mathbf{N}^m(\boldsymbol{Y}) = \beta \cdot \mathbf{N}^m(\boldsymbol{Y}), \forall \boldsymbol{Y} \in \partial \Omega_0^m, \\
\frac{1}{|\Omega_0^m|} \int_{\Omega_0^m} \boldsymbol{u}^m(\boldsymbol{Y}) \mathrm{d}\boldsymbol{Y} = \boldsymbol{u}^M(\boldsymbol{X}), \\
\frac{1}{|\Omega_0^m|} \int_{\Omega_0^m} \nabla_Y \boldsymbol{u}^m(\boldsymbol{Y}) \mathrm{d}\boldsymbol{Y} = \mathbf{G}(\boldsymbol{X}).
\end{array}
\right.
\qquad (6)
$$

The von Mises elasto-plastic model with isotropic elasticity and com-
bined isotropic and kinematic hardening, for instance, is given by

$$
\begin{align}
\boldsymbol{\sigma}^m &= 2\mu^m \boldsymbol{\epsilon}^{m,e} + \lambda^m \operatorname{tr}(\boldsymbol{\epsilon}^{m,e}), \nonumber \\
\boldsymbol{\epsilon}^{m,e} &= \boldsymbol{\epsilon}^m - \boldsymbol{\epsilon}^{m,p} \nonumber \\
\boldsymbol{\epsilon}^m &= \varepsilon (\nabla_Y \boldsymbol{u}^m(\boldsymbol{Y})) = \frac{1}{2} (\nabla_Y \boldsymbol{u}^m(\boldsymbol{Y}) + \nabla_Y^T \boldsymbol{u}^m(\boldsymbol{Y})), \tag{7}
\end{align}
$$

where an additive decomposition into an elastic part $\boldsymbol{\epsilon}^{m,e}$ and a plastic part $\boldsymbol{\epsilon}^{m,p}$ has been assumed for the small strain tensor $\boldsymbol{\epsilon}^m$. Plastic flow is ruled by the condition $\bar{f}(\boldsymbol{\sigma}^{m,dev} - \boldsymbol{\gamma}^m) - \bar{\sigma}(\bar{\boldsymbol{\epsilon}}^{m,p}) \le 0$, with the von Mises equivalent stress $\bar{f}$ given by

$$
\begin{align}
\bar{f}(\boldsymbol{\sigma}^{m,dev} - \boldsymbol{\gamma}^m) &= \sqrt{\frac{3}{2}(\boldsymbol{\sigma}^{m,dev} - \boldsymbol{\gamma}^m) : (\boldsymbol{\sigma}^{m,dev} - \boldsymbol{\gamma}^m)}, \nonumber \\
\boldsymbol{\sigma}^{m,dev} &= \boldsymbol{\sigma}^m - \frac{1}{3} \operatorname{tr}(\boldsymbol{\sigma}^m), \tag{8}
\end{align}
$$

and the yield stress $\bar{\sigma}$ given by

$$
\begin{align}
\bar{\sigma}(\bar{\boldsymbol{\epsilon}}^{m,p}) &= \sigma_Y + K \bar{\boldsymbol{\epsilon}}^{m,p}, \nonumber \\
\bar{\boldsymbol{\epsilon}}^{m,p} &= \int_0^T \bar{\boldsymbol{\epsilon}}^{m,p}(s) ds, \nonumber \\
\bar{\boldsymbol{\epsilon}}^{m,p} &= \sqrt{\frac{2}{3} \dot{\bar{\boldsymbol{\epsilon}}}^{m,p} : \dot{\bar{\boldsymbol{\epsilon}}}^{m,p}}. \tag{9}
\end{align}
$$

Material parameter $\sigma_Y$ is the initial yield stress and $K$ is the isotropic
hardening modulus, while internal variable $\bar{\boldsymbol{\epsilon}}^{m,p}$ is the equivalent plastic
strain. The latter is defined as the integral over time of the equivalent
plastic strain rate $\dot{\bar{\boldsymbol{\epsilon}}}^{m,p}$. Under plastic flow (e.g., when $\bar{f}(\boldsymbol{\sigma}^{m,dev} - \boldsymbol{\gamma}^m) - \bar{\sigma}(\bar{\boldsymbol{\epsilon}}^{m,p}) = 0$), the plastic strain rate is given by

$$
\dot{\boldsymbol{\epsilon}}^{m,p} = \bar{\boldsymbol{\epsilon}}^{m,p} \frac{3}{2} \frac{\boldsymbol{\sigma}^{m,dev}}{\bar{f}}. \qquad (10)
$$

Kinematic hardening is modeled using a linear law where

$$
\dot{\boldsymbol{\gamma}}^m = H \dot{\boldsymbol{\epsilon}}^{m,p}, \qquad (11)
$$

$H$ being the kinematic hardening modulus.

The Next Increment Corrects Error (NICE) scheme is used to integrate
this nonlinear material law [13].

**3. Reference multiscale algorithm**

Computational homogenization with an FE² approach consists in discretizing and solving the problems in Eqs. (1) and (2) using the FE method. The first step is to discretize the coarse scale domain $\Omega_0^M$ with an appropriate FE mesh. Then, at each integration point of this coarse scale mesh, a new discrete fine scale model should be created with its own and independent FE mesh. This is the reason why the computational complexity of an FE² scheme is of $\mathcal{O}(N_q^M N_{DOF}^m)$, where $N_q^M$ is the number of coarse scale integration points, and $N_{DOF}^m$ is the number of fine scale degrees of freedom. The whole point of the PFE² method is to reduce this complexity to $\mathcal{O}(N_c^M N_{DOF}^m)$, where $N_c^M$ is the number of clusters with, ideally, $N_c^M \ll N_q^M$.

In addition to the meshes themselves and the associated coordinate
and connectivity data, the FE² scheme requires to store an unknown
displacement vector for each node of the coarse scale mesh, and for
each node of each fine scale mesh, as well as an internal variables vec-
tor for each integration point of each fine scale mesh. This is the strict
minimum of information needed in theory, but in practice the numeri-
cal implementation usually requires more data such as basis functions,
mappings, strains, stresses, etc.

3.1. Coarse scale solver

In this work, although all considered material models are rate-
independent, loading is applied over several time increments $t^n = n\Delta t \in$
$[0, T]$, with $\Delta t$ the time step. At each time increment, the coarse scale
solver given in Algorithm 1 is used.

**Algorithm 1:** Reference coarse scale solver

Run any pre-processing operation such as the management of
internal variables for each fine scale problem;

Initialize $\mathbf{u}^{M,n+1,0} = \mathbf{u}^{M,n}$;

for i = 0, ..., N<sub>NR</sub> do

Interpolate the displacement and its gradient at each
integration point of the coarse scale mesh from the current
guess $\mathbf{u}^{M,n+1,j}$;

Solve each fine scale problem as per Eq. (2) to obtain the fine
scale displacement, the coarse scale stress tensor $\mathbf{P}^M(\mathbf{G}^{n+1,j})$
and the first derivative tensor $\frac{\partial \mathbf{P}^M}{\partial \mathbf{G}}(\mathbf{G}^{n+1,j})$ as per Eq. (12);

Assemble the linearized coarse scale problem;

Test for convergence and stop if the criterion is met under
some tolerance;

Solve the linearized coarse scale problem to obtain the new
guess $\mathbf{u}^{M,n+1,i+1}$;

end

Run any post-processing operation and write output files;

The convergence criterion checks whether the residual is small
enough under some absolute tolerance, here $10^{-6}$.

On the one hand, if convergence is not met within the prescribed
number of Newton-Raphson iterations (here $N_{NR}$), the time step $\Delta t$ is
divided by two and then the algorithm is restarted. This also happens if
the linearized coarse scale problem could not be solved (ill-conditioned
matrix), if the coarse scale residual increased too much (more than thou-
sand times the initial residual), or if a fine scale problem could not be
solved. Note that a similar Newton-Raphson algorithm with restarts is
used for fine scale problems, due to the same issues, to which is added
a possible failure of the return mapping algorithm used for plasticity.
On the other hand, if convergence is met without any restart, then the
time step $\Delta t$ is multiplied by $\frac{3}{2}$ up to an upper bound which is the user-
defined initial time step $\Delta t_0$. Note that this is a pure Newton-Raphson
algorithm. It does not, in particular, rely on any line search procedure
neither at the coarse scale nor the fine scale. Restarting with time step
decrease can be seen as a substitute for line searching.

On the other hand, if convergence is met without any restart, then the time step $\Delta t$ is multiplied by $\frac{3}{2}$ up to an upper bound which is the user-defined initial time step $\Delta t_0$. Note that this is a pure Newton-Raphson algorithm. It does not, in particular, rely on any line search procedure neither at the coarse scale nor the fine scale. Restarting with time step decrease can be seen as a substitute for line searching.

3.2. Computation of tangent moduli

Details on the assembly and solution of the fine scale problem
in Eq. (2) can be found in Refs. [22,23], especially regarding the
role of the Lagrange multipliers. Each fine scale problem, in addition,
should be solved several times in order to compute the tangent modu-
lus $\frac{\partial \mathbf{P}^M}{\partial \mathbf{G}} = \frac{\partial \mathbf{P}^M}{\partial \mathbf{G}}(\mathbf{G}^{n+1,j})$ through numerical differentiation using the
forward formula:

$$
\frac{\partial \mathbf{P}_{ab}^M}{\partial \mathbf{G}_{kl}^M} (\mathbf{G}^{n+1,i}) \approx \frac{\mathbf{P}_{ab}^M (\mathbf{G}^{n+1,i} + \epsilon \mathbf{1}^{kl}) - \mathbf{P}_{ab}^M (\mathbf{G}^{n+1,j})}{\epsilon}, \quad (12)
$$
where n is the time increment number, i is the coarse scale solver iteration number, $\epsilon = 10^{-4}$ and $\mathbf{1}^{kl}$ is a second order tensor, with

$$ \mathbf{1}_{ab}^{kl} = \delta_{a,k} \delta_{b,l}, a, b = 1, \dots, 3, $$

and the Kronecker delta:

$$ \delta_{a,k} = \begin{cases} 1, & a = k, \\ 0, & a \neq k. \end{cases} $$

## 4. Partitioned multiscale algorithm

In this section, the novel k-medoids-based PFE² approach is presented. The idea of the PFE² approach as proposed in Refs. [2,4] is to group coarse scale integration points into a user-defined number of clusters. In the present work, it is proposed to do so using the k-medoids clustering algorithm. This clustering is performed at each iteration of the coarse scale solver for the first $N_{clustering}$ iterations, and then clustering is frozen for the remaining iterations. At each iteration, fine scale problems are only solved for the cluster medoid, and an approximation is used to compute stresses for all non-medoid integration points.

### 4.1. Clustering

The input variable C for the k-medoids clustering algorithm is simply the concatenation of the current coarse scale displacement gradient guess at each coarse scale integration point, and the averaged fine scale internal variables over each fine scale domain. For the finite strain model described in Section 2.3, therefore, the input variable is:

$$ C(X_p) = \left\{ \left( \mathbf{G}_{kl}^{n+1,j}(X_p) \right)_{k,l=1,...,3} \right\}, \quad (13) $$

where n is the time increment number, i is the coarse scale solver iteration number and p is the integration point number. For the small strain model in Section 2.4, the input variable is changed to:

$$ C(X_p) = \left\{ \left( \mathbf{G}_{kl}^{n+1,j}(X_p) \right)_{k,l=1,...,3} \right\}, \quad (14) $$

where $\bar{\epsilon}^{M,p,n+1,i} = \int_{\Omega_0^m} \bar{\epsilon}^{m,p,n+1,j} d\Omega_0^m, \bar{\epsilon}_{kl}^{M,p,n+1,i} = \int_{\Omega_0^m} \bar{\epsilon}_{kl}^{m,p,n+1,j} d\Omega_0^m, \Omega_0^m = \Omega_0^m(X_p)$ are computed only to incorporate history information into clustering and do not have a physical meaning. Since $\bar{\epsilon}^{M,p,n+1,i}$ symmetric, it is unnecessary to include non-diagonal components twice, but this optimization is not implemented to simplify the code.

In the k-medoids clustering algorithm, each cluster is identified by its medoid, which is an integration point of the coarse scale FE mesh. The algorithm requires the definition of a distance between cluster medoids. Here the following Euclidean distance is used:

$$ d(X_p, X_q) = \sum_k (C_k(X_p) - C_k(X_q))^2, \quad (15) $$

where $C_k$ is the component k of C, with the first nine components being occupied by the displacement gradient. The sum over k includes all components: nine for the finite strain model described in Section 2.3, and 19 for the small strain model described in Section 2.4. All distances are pre-computed and stored into a symmetric distance matrix. As presented in Ref. [7], the k-medoids clustering algorithm given in Algorithm 2 is then used.

As detailed in Ref. [7], the clustering convergence criterion consists of checking for changes in cluster medoids over consecutive iterations and also over small cycles of a few iterations to ensure that the algorithm is not cycling over the same periodic solution.

It should be noted that the input variable C might be constant or lead to a distance matrix mostly filled by zeroes, especially at the first nonlinear iteration of the first time increment. Due to the random initialization of clusters, the clustering algorithm still respects the user-defined number of clusters in such cases, even though cluster sizes might be uneven.

### Algorithm 2: K-medoids clustering algorithm.

Recover the last cluster assignment, if any, or randomly assign integration points to clusters;

Compute the medoid of each cluster as the integration point with the smallest sum of distances to the other integration points in the cluster;

repeat

Compare the distance between each integration point and all current cluster medoids;

Reassign each integration point to the closest cluster (e.g., the one whose current medoid is the closest);

Recompute cluster medoids;

until Convergence criterion is met under some tolerance;

### 4.2. Computation of stresses and tangent moduli

As proposed in Ref. [4], stresses of all integration points within a cluster are computed through the approximation:

$$ \mathbf{P}^M (\mathbf{G}(X_p)) = \mathbf{P}^M (\mathbf{G}(X_c)) + \frac{\partial \mathbf{P}^M}{\partial \mathbf{G}} (\mathbf{G}(X_c)) : (\mathbf{G}(X_p) - \mathbf{G}(X_c)), \quad (16) $$

or, introducing $\mathbf{G}_p = \mathbf{G}(X_p)$ and $\mathbf{G}_c = \mathbf{G}(X_c)$:

$$ \mathbf{P}^M (\mathbf{G}_p) = \mathbf{P}^M (\mathbf{G}_c) + \frac{\partial \mathbf{P}^M}{\partial \mathbf{G}} (\mathbf{G}_c) : (\mathbf{G}_p - \mathbf{G}_c). \quad (17) $$

It is clear from this expression that

$$ \frac{\partial \mathbf{P}^M}{\partial \mathbf{G}_p} (\mathbf{G}_p) = \frac{\partial \mathbf{P}^M}{\partial \mathbf{G}} (\mathbf{G}_c). \quad (18) $$

The authors of Ref. [4], however, did not notice that

$$ \begin{aligned} \frac{\partial \mathbf{P}^M}{\partial \mathbf{G}_c}(\mathbf{G}_p) &= \frac{\partial \mathbf{P}^M}{\partial \mathbf{G}}(\mathbf{G}_c) - \frac{\partial \mathbf{P}^M}{\partial \mathbf{G}}(\mathbf{G}_c) + \frac{\partial^2 \mathbf{P}^M}{\partial \mathbf{G}^2}(\mathbf{G}_c) : (\mathbf{G}_p - \mathbf{G}_c) \\ &= \frac{\partial^2 \mathbf{P}^M}{\partial \mathbf{G}^2}(\mathbf{G}_c) : (\mathbf{G}_p - \mathbf{G}_c). \end{aligned} \quad (19) $$

During the linearization, therefore, cross contributions arise. This is specific to the clustering approach which relies on approximation (16), and they would hence not arise in a conventional FE² method. In Ref. [4], the authors used k-means clustering and defined $\mathbf{G}_c$ as the weighted average of all displacement gradients throughout each cluster. Cross contributions would have paired all integration points of a cluster with each other, thus reducing drastically the sparsity of the linearized coarse scale problem.

In the present work, because $\mathbf{G}_c$ is the displacement gradient at the cluster's medoid, cross contributions only pair each integration point with its cluster's medoid. It is reminded that, at each nonlinear iteration of the coarse scale solver, the weak form of Eq. (1) features the term $\int_{\Omega^M} \mathbf{P}^M(\mathbf{G}_p^{n+1,j+1}) : \nabla_X \mathbf{v}^M(X_p) dX_p$. As a consequence of Eqs. (18) and (19), it is linearized as:

$$ \begin{aligned} \int_{\Omega^M} \mathbf{P}^M(\mathbf{G}_p^{n+1,j+1}) : \nabla_X \mathbf{v}^M(X_p) dX_p & \approx \int_{\Omega^M} \mathbf{P}^M(\mathbf{G}_p^{n+1,j}) : \nabla_X \mathbf{v}^M(X_p) dX_p \\ & + \int_{\Omega^M} \left( \frac{\partial \mathbf{P}^M}{\partial \mathbf{G}_p}(\mathbf{G}_p^{n+1,j}) : \Delta \mathbf{G}_p^{n+1,j+1} \right) : \nabla_X \mathbf{v}^M(X_p) dX_p \\ & + \int_{\Omega^M} \left( \frac{\partial \mathbf{P}^M}{\partial \mathbf{G}_c}(\mathbf{G}_p^{n+1,j}) : \Delta \mathbf{G}_c^{n+1,j+1} \right) : \nabla_X \mathbf{v}^M(X_p) dX_p, \end{aligned} \quad (20) $$

where $\Delta \mathbf{G}_p^{n+1,j+1} = \mathbf{G}_p^{n+1,j+1} - \mathbf{G}_p^{n+1,j}$, $\Delta \mathbf{G}_c^{n+1} = \mathbf{G}_c^{n+1} - \mathbf{G}_c^i$, and $\mathbf{v}^M$ is the test displacement field. If the last term of Eq. (20) is not neglected, then lines of the linearized coarse scale problem associated to the nodes of the element containing integration point $X_p$ receive non-zero contributions for the columns associated to the nodes of the element containing cluster
medoid $X_c$. This should reduce slightly the sparsity of the linearized coarse scale problem.

The numerical differentiation formula to compute the first derivative tensor $\frac{\partial P^M}{\partial G}(G_c^{n+1,j})$ is the same as the one presented in Section 3.2. The computation of the second derivative tensor $\frac{\partial^2 P^M}{\partial G^2}(G_c^{n+1,j})$, moreover, is quite elaborated and requires a great number of fine scale solves:

$$ \frac{\partial^2 \mathbf{P}_{ab}^M}{\partial \mathbf{G}_{kl} \partial \mathbf{G}_{rs}}(G_c^{n+1,j}) \approx \frac{\frac{\partial \mathbf{P}_{ab}^M}{\partial \mathbf{G}_{kl}}(G_c^{n+1,j} + \frac{\sqrt{\epsilon}}{2} \mathbf{1}^{rs}) - \frac{\partial \mathbf{P}_{ab}^M}{\partial \mathbf{G}_{kl}}(G_c^{n+1,j} - \frac{\sqrt{\epsilon}}{2} \mathbf{1}^{rs})}{\sqrt{\epsilon}}, $$

$$ \frac{\partial \mathbf{P}_{ab}^M}{\partial \mathbf{G}_{kl}}(G_c^{n+1,j} + \frac{\sqrt{\epsilon}}{2} \mathbf{1}^{rs}) $$

$$ \approx \frac{\mathbf{P}_{ab}^M(G_c^{n+1,j} + \frac{\sqrt{\epsilon}}{2} \mathbf{1}^{rs} + \frac{\sqrt{\epsilon}}{2} \mathbf{1}^{kl}) - \mathbf{P}_{ab}^M(G_c^{n+1,j} + \frac{\sqrt{\epsilon}}{2} \mathbf{1}^{rs} - \frac{\sqrt{\epsilon}}{2} \mathbf{1}^{kl})}{\sqrt{\epsilon}}, $$

$$ \frac{\partial \mathbf{P}_{ab}^M}{\partial \mathbf{G}_{kl}}(G_c^{n+1,j} - \frac{\sqrt{\epsilon}}{2} \mathbf{1}^{rs}) $$

$$ \approx \frac{\mathbf{P}_{ab}^M(G_c^{n+1,j} - \frac{\sqrt{\epsilon}}{2} \mathbf{1}^{rs} + \frac{\sqrt{\epsilon}}{2} \mathbf{1}^{kl}) - \mathbf{P}_{ab}^M(G_c^{n+1,j} - \frac{\sqrt{\epsilon}}{2} \mathbf{1}^{rs} - \frac{\sqrt{\epsilon}}{2} \mathbf{1}^{kl})}{\sqrt{\epsilon}}. \quad (21) $$

While Eq. (12) requires $3 \times 3 + 1 = 10$ fine scale solves, Eq. (21) requires $3 \times 3 \times 2 \times 3 \times 3 \times 2 = 324$ solves. Note that it is possible to save more than half these solves by carefully identifying and eliminating redundant ones within Eq. (21) and between Eqs. (12) and (21). It is possible to save even more solves by smartly choosing either forward or centered differentiation formulas depending on components. No optimization is done in the present work to keep the implementation as simple as possible. Alternatives to numerical differentiation through analytical derivation could also be considered to optimize the computation of both first and second derivatives [3].

### 4.3. Update of internal variables

An elaborate mapping algorithm was proposed in Ref. [4] to map internal variables by identifying origin and destination clusters when an integration point was assigned to a new cluster. In the present work, it is chosen to systematically solve all fine scale problems after convergence of the partitioned coarse scale solver. It is crucial to keep all integration points up to date, indeed, since after each run of the k-medoids clustering algorithm, any integration point may become a cluster medoid. With this scheme, therefore, partitioning does not completely eliminate the need to solve all fine scale problems. It nevertheless saves many fine scale solves during iterations of the coarse scale solver as compared to the reference FE² algorithm in Algorithm 1, including those required to compute tangent moduli.

### 4.4. Partitioned coarse scale solver

The partitioned coarse scale solver used at each time increment is summarized in Algorithm 3.

As in the non-partitioned coarse scale solver, the time step is automatically decreased (with restarts) if any issue is encountered during the prescribed number of Newton-Raphson iterations (here $N_{NR}$). If convergence is not met at the last iteration, however, it is chosen to fall back to the full FE² method. This is estimated to be preferable and a better solution than the retry procedure proposed in Ref. [4]. The latter consists in retrying nonlinear solves with slightly changed cluster assignments with the hope that convergence might be reached. In the present work, it is preferred to fall back to the full FE² method in order to get a solution in a reasonable number of expensive nonlinear iterations, instead of accumulating a great number of cheap nonlinear iterations. The occurrence of fall backs and its impact on the computational cost of the new k-medoids-based PFE² method is analyzed in

### Algorithm 3: Partitioned coarse scale solver.

Run any pre-processing operation such as the management of internal variables for each fine scale problem;

Initialize $\mathbf{u}^{M,n+1,0} = \mathbf{u}^{M,n}$;

for $i = 0, \dots, N_{NR}$ do

Interpolate the displacement and its gradient at each integration point of the coarse scale mesh from the current guess $\mathbf{u}^{M,n+1,j}$;

if $i < N_{clustering}$ then

Construct the input variable $C$ from the coarse scale displacement gradient (Eq. 13) and fine scale internal variables averages (Eq. 14);

Pre-compute the distance matrix (Eq. 15);

Run the k-medoids clustering algorithm (Algorithm 2);

end

Solve each fine scale problem associated to a cluster medoid as per Eq. (2) to obtain the fine scale displacement, the coarse scale stress tensor $\mathbf{P}^M(G_c^{n+1,j})$, the first derivative tensor $\frac{\partial P^M}{\partial G}(G_c^{n+1,j})$ as per Eq. (12) and, if cross contributions are to be assembled, also the second

derivative tensor $\frac{\partial^2 P^M}{\partial G^2}(G_c^{n+1,j})$ as per Eq. (21);

Compute homogenized stresses as per Eq. (16) for all fine scale problems which are not cluster medoids;

Assemble the linearized coarse scale problem;

Test for convergence and stop if the criterion is met under some tolerance;

Solve the linearized coarse scale problem to obtain the new guess $\mathbf{u}^{M,n+1,i+1}$;

end

Solve all fine scale problems as per Eq. (2) to update internal variables;

Run any post-processing operation and write output files;

Section 5. Cluster freezing is also implemented to improve convergence, by choosing $N_{clustering} < N_{NR}$.

It should be noted that if cross contributions are omitted, the coarse scale solver is actually of quasi-Newton type, but this subtlety is voluntarily ignored throughout the paper.

## 5. Results

Both the FE² and the new k-medoids-based PFE² schemes have been implemented within the open source FEMS software [20,21]. To keep computation times reasonable for the reference FE² simulations, two-dimensional problems are considered under the plane strains assumption. Linear triangular elements are chosen for all simulations because a single integration point can be used per element, which simplifies the visualization of cluster assignments. In the following, therefore, the words integration point may be substituted by element and vice-versa. It is reminded that all the considered material models are rate-independent.

Unless otherwise mentioned, second derivative tensors and the associated cross contributions are neither computed nor assembled. Numerical simulations are conducted to evaluate the accuracy of the k-medoids-based PFE² method as compared to reference FE² results, and its interest in the point of view of computation time. The number of iterations of the partitioned coarse scale solver is set to $N_{NR} = 32$, and the number of iterations before freezing $N_{clustering} = 7$. It is reminded that if the convergence criterion is not satisfied within $N_{NR} = 32$ iterations of the partitioned coarse scale solver, the algorithm falls back to the full FE² algorithm.

Accuracy is evaluated based on the following relative $L^2$ norm error on the displacement vector field $\mathbf{u}^M$ for some k-medoids-based PFE² simulation as compared to a reference FE² displacement vector field $\mathbf{u}^{M,ref}$:
$$ Error(\mathbf{u}^M) = \sqrt{\frac{\int_{\Omega_0^M} \sum_{j=1}^2 (\mathbf{u}_j^M - \mathbf{u}_j^{M,ref})^2 d\Omega_0^M}{\int_{\Omega_0^M} \sum_{j=1}^2 (\mathbf{u}_j^{M,ref})^2 d\Omega_0^M}}, \qquad (22) $$

and also a similar error indicator for the Cauchy stress tensor field:

$$ Error(\sigma^M) = \sqrt{\frac{\int_{\Omega_0^M} \sum_{k,l=1}^2 (\sigma_{kl}^M - \sigma_{kl}^{M,ref})^2 d\Omega_0^M}{\int_{\Omega_0^M} \sum_{k,l=1}^2 (\sigma_{kl}^{M,ref})^2 d\Omega_0^M}}. \qquad (23) $$

The equivalent stress field computed with the formula $\sqrt{\sum_{k,l=1}^2 (\sigma_{kl}^M)^2}$ is shown for some results. Computation time is evaluated using a high-performance computing workstation with an Intel(R) Xeon(R) W-2175 multi-processor with 14 CPUs and 64GB of RAM. The code is sequential. Multi-threading or distributed parallel programming could be interesting in the future to take advantage of all the CPUs and reduce the computation time even more. Linearized problems are solved at both scales using the direct solver UMFPACK [6], which operates an LU decomposition. Computation times are always mentioned in terms of speed-up ratio relatively to the corresponding FE² simulation, with one meaning there was no speed-up, two meaning the PFE² simulation finished in half the time required for the FE² simulation, etc.

## 5.1. Hyperelastic beam

The first setup is similar to the one presented in Ref. [4]. Geometries, meshes and boundary conditions are illustrated in Fig. 2. Under the assumption of plane strain, a beam of length 400 mm and height 100 mm is fixed on its left side and lifted at its top right corner with a vertical displacement of 125 mm. The fine scale domain is a square of size $1 \times 1$ mm² composed of a soft matrix reinforced with an inclusion of radius 0.2 mm at its center. The Saint-Venant-Kirchhoff hyperelastic model presented in Section 2.3 is used with a Young's modulus of 210 GPa and a Poisson's coefficient of 0.3 for the inclusion, and a Young's modulus of 2 GPa and a Poisson's coefficient of 0.25 for the matrix.

The loading is applied over five increments, and the mesh is composed of 736 linear triangular elements for the coarse scale, and 916 linear triangular elements for the fine scale. This first setup does not involve any internal variable.

### 5.1.1. Analysis of medoid change

A first simulation is conducted with only 1 cluster over the whole beam. Because clustering cannot change with a single cluster, the objective is to analyze medoid changes during the simulation. Every time clustering is triggered (in the first nonlinear iterations of each time increment), the medoid position is registered. This sequence of medoid coordinates is analyzed in terms of mean and standard deviation. This

analysis is indirect, since the clustering input variable in Eq. (13) directly depends on the displacement gradient and not coordinates, but the displacement gradient itself is a spatial field.

The mean position is (164 mm, 47.1 mm) in the initial configuration. It is located at the center of the beam with a slight shift towards the fixed side. The standard deviation is (21.5 mm, 6.43 mm), which is small with respect to beam dimensions, and also to average mesh size (e.g., 11.2 mm).

A large standard deviation would have revealed large fluctuations and numerical errors due to the approximation in Eq. (16), and an inability of the clustering algorithm to converge on a stable medoid for a monotonic loading. Since the standard deviation is small, the strategy proposed in the present work is shown to be relevant.

### 5.1.2. Convergence with an increasing number of clusters

Simulations are conducted with an increasing number of clusters in order to evaluate improvements in terms of accuracy. Equivalent stress fields are presented in Fig. 3(a)–(c) for some of those simulations. The reference FE² result is shown in Fig. 3(d). Although there is a slight stress overestimation at the top left side of the beam, it is corrected by using more clusters.

This can be explained by looking at cluster assignments in Fig. 4. The number of clusters is first increased along the beam, which helps to capture stress evolution in the horizontal direction. Then, a specific cluster is created for the top right corner which is a stress concentration spot, and some clusters become split in the vertical direction. These assignments are overall similar to those found in Ref. [4], especially when using a low number of clusters.

Errors are reported in Table 1. The relative displacement error is already quite low (below 10 %) using 1 cluster. For the stress, however, 23 clusters are required to reach the same bound. Regarding the convergence rate, as shown in Fig. 5, the displacement error is divided by around two when multiplying the number of clusters by two from 1 to 2, 11 to 23, 46 to 92 and 184 to 368. The same observation applies to the stress error from 1 to 2 clusters, 5 to 11 and 11 to 23. Overall, nevertheless, convergence is much slower.

### 5.1.3. Computation time

Speed-up ratios are also reported in Table 1. For simulations with less than 11 clusters, the most expensive operation is solving all fine scale problems after convergence of the partitioned coarse scale solver: 50–92 % of the total computation time. It should be noted that this operation could be omitted for this setup because there are no internal variables. With more clusters, computation time is dominated by the partitioned coarse scale solver and, except for the case with 184 clusters, it increases linearly with the number of clusters. Except for the simulation with 184 clusters, there were no convergence issues and no restarts, including for the reference FE² simulation.

<figure><img src="image_4.png" /><figcaption>Fig. 2. Geometries, meshes and boundary conditions for hyperelastic beam simulations: (a) coarse scale, (b) fine scale with colors distinguishing matrix and inclusion.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 3. Equivalent stress (in MPa) for the hyperelastic beam using: (a) 5 clusters, (b) 11 clusters, (c) 23 clusters, (d) the FE² scheme.</figcaption></figure>

The reference FE² simulation required 13,470 s to finish. The simulation with 23 clusters, which leads to an error below 10% for both displacement and stress fields, is interesting since it reduces the computation time by a factor of 16.7. This is similar to the speed-up ratio found in Ref. [4]. Even for the simulation with 184 clusters, where a convergence issue led to a fall back and eight (among 82) coarse scale solver iterations without clustering, there is still a computation time reduction factor of 1.66.

Computation time spent clustering (including pre-computing distances), in addition, is negligible as it represented less than 1 % of the total computation time for all simulations. This is explained by the fact that the last assignment is used as initial guess for each run of the k-medoids clustering algorithm, and that, as analyzed above, assignments and medoids do not change much during these simulations with monotonic loading.

### 5.1.4. Importance of cross contributions

For all results reported in Table 1, second derivative tensors and cross contributions were neither computed nor assembled. Since no convergence issue was encountered for most of these simulations, it can be concluded that assembling cross contributions is not a necessary condition for convergence of the coarse scale solver.

In order to assess the importance of these contributions, two analysis steps are proposed. First, a focus is made on the first nonlinear iteration of the simulation with a single cluster. A reference Jacobian matrix $J^{ref,n+1,i}$ is computed using the PETSc library's finite difference Jacobian approximation tool [1] which conducts numerical differentiation with a formula similar to Eq. (12) but over the whole nonlinear system using

the formula:

$$J_{ak}^{ref,n+1,i} \approx \frac{R_a^M(U^{n+1,i} + \epsilon \mathbf{1}^k) - R_a^M(U^{n+1,i})}{\epsilon}, \quad a, k = 1, \dots, N_{DOF}^M \quad (24)$$

where $\mathbf{1}^k$ is a vector, with

$$\mathbf{1}_a^k = \delta_{a,k}, \quad a = 1, \dots, N_{DOF}^M.$$

Moreover, $R^M$ is the global right hand side vector and $U$ the global displacement vector. They both have as dimension the number of coarse scale degrees of freedom $N_{DOF}^M$, which is equal to twice (2D) the number of nodes. The computation time for this procedure is nearly prohibitive but it ensures that influences of displacements across the mesh are captured (e.g., a nodal displacement affecting the force of a node with which it does not share an element). Jacobian matrices $J^{n+1,i}$ are then computed and assembled with Eq. (12) by enabling or disabling cross contributions, and compared to the reference one $J^{ref,n+1,i}$ through the error measure:

$$\text{Error}(J^{n+1,i}) = \frac{\|J^{n+1,i} - J^{ref,n+1,i}\|_2}{\|J^{ref,n+1,i}\|_2}, \quad (25)$$

An error of $1.04 \times 10^{-1}$ is obtained without cross contribution, and $5.64 \times 10^{-5}$ with cross contributions. This clearly shows that the Jacobian matrix is not complete when cross contributions are neglected.

The second analysis step is to evaluate the influence of cross contributions on the number of iterations of the partitioned coarse scale solver, and the occurrence of convergence issues. The simulations with 92 and 184 clusters are reran with cross contributions. The number of
<figure><img src="image_6.png" /><figcaption>Fig. 4. Final cluster assignments for the hyperelastic beam using: (a) 2 clusters, (b) 5 clusters, (c) 11 clusters, (d) 23 clusters.</figcaption></figure>

<caption>Table 1</caption>

Errors and speed-up ratios for hyperelastic beam simulations using an increasing number of clusters.

<table><thead><tr><th># clusters</th><th>Error(u<sup>M</sup>)</th><th>Error(σ<sup>M</sup>)</th><th>Speed-up</th></tr></thead><tbody><tr><td>1</td><td>5.56 × 10<sup>-2</sup></td><td>2.53 × 10<sup>-1</sup></td><td>36.5</td></tr><tr><td>2</td><td>2.53 × 10<sup>-2</sup></td><td>1.35 × 10<sup>-1</sup></td><td>36.3</td></tr><tr><td>5</td><td>1.83 × 10<sup>-2</sup></td><td>1.20 × 10<sup>-1</sup></td><td>32.3</td></tr><tr><td>11</td><td>1.05 × 10<sup>-2</sup></td><td>6.23 × 10<sup>-2</sup></td><td>24.8</td></tr><tr><td>23</td><td>4.87 × 10<sup>-3</sup></td><td>3.74 × 10<sup>-2</sup></td><td>16.7</td></tr><tr><td>46</td><td>3.69 × 10<sup>-3</sup></td><td>3.35 × 10<sup>-2</sup></td><td>11.3</td></tr><tr><td>92</td><td>1.80 × 10<sup>-3</sup></td><td>2.88 × 10<sup>-2</sup></td><td>5.34</td></tr><tr><td>184</td><td>1.21 × 10<sup>-3</sup></td><td>2.82 × 10<sup>-2</sup></td><td>1.66</td></tr><tr><td>368</td><td>6.76 × 10<sup>-4</sup></td><td>2.78 × 10<sup>-2</sup></td><td>1.64</td></tr></tbody></table>

partitioned coarse scale solver iterations with 92 clusters are reduced from 57 to 41 by assembling cross contributions. The result is more impressive with 184 clusters, as the reduction is from 82 to 43. More importantly, no convergence issue and hence no fall back occurred in the simulation with cross contributions, as opposed to the one without.

One could think avoiding fall backs to the full non-partitioned coarse scale solver would be quite beneficial in terms of computation time. Computation of second derivatives is, nevertheless, very expensive, as the simulation with 184 clusters and cross contributions took seven times longer to complete. This is even worse for the simulation with 92 clusters, as it took 11 times longer.

As a conclusion, although they are not always necessary to reach convergence, cross contributions could be helpful to avoid some issues. The computational cost associated to second derivatives, is however, too prohibitive.

<figure><img src="image_7.png" /><figcaption>Fig. 5. Errors for hyperelastic beam simulations using an increasing number of clusters.</figcaption></figure>

## 5.2. Elasto-plastic bracket

The second setup is also inspired from Ref. [4]. Under the assumption of plane strains, the bracket shown in Fig. 6(a) is subjected to a non-proportional loading path where holes A and B are fixed, while hole C is fixed in the horizontal direction and subjected to the time-dependent displacement shown in Fig. 6(c) in the vertical direction. The fine scale domain is the same as in the previous setup, but the material
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Geometries, meshes and loading for elasto-plastic bracket simulations: (a) coarse scale, (b) fine scale with colors distinguishing matrix and void, (c) vertical displacement applied at C. Radii are 70 mm for holes A and B, and 120 mm for hole C. Hole centers for A and B are positioned 200 mm away from the left side and 150 mm away from the top and bottom side, respectively. Hole C's center is positioned 200 mm away from the right side, and D = 80 mm.</figcaption></figure>

model is replaced by the small strain elasto-plastic model presented in Section 2.4, and the inclusion is replaced by a void. The latter is modeled by using properties such as a Young's modulus of 2 KPa, a Poisson's coefficient of 0.3 and an initial yield stress large enough to prevent plasticity. Matrix properties are a Young's modulus of 2 GPa, a Poisson's coefficient of 0.3, an initial yield stress of 24 GPa, and a kinematic hardening modulus of 80 GPa. Isotropic hardening is neglected.

Loading is applied over 61 increments with a time step of 1 s, except for the first increment which uses a very small step to get an elastic response and initialize the NICE scheme. The mesh is composed of 8130 linear triangular elements for the coarse scale, and 916 linear triangular elements for the fine scale. This second setup involves an internal variable to track matrix plasticity at the fine scale, and a change in the loading direction at the coarse scale.

### 5.2.1. Convergence with an increasing number of clusters

Equivalent stress fields at the loading minimum are presented in Fig. 7. They are all qualitatively similar, with the only difference being some slight noise when using less clusters.

To compare results to Ref. [4], the local shear stress $\sigma_{12}^M$ history is analyzed for the point E defined in Fig. 6(a). Shear stress–strain curves at

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Equivalent stress (in MPa) at $t = 45$ s for the elasto-plastic bracket using: (a) 81 clusters, (b) 162 clusters, (c) 325 clusters, (d), the FE$^2$ scheme.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Stress-strain curves at point E defined in Fig. 6(a): components $\epsilon_{12}^M$ and $\sigma_{12}^M$.</figcaption></figure>

<caption>**Table 2**</caption>

Errors and speed-up ratios for elasto-plastic bracket simulations using an increasing number of clusters.

<table><thead><tr><th># clusters</th><th>Error(u<sup>M</sup>)<br>t = 45 s</th><th>Error(σ<sup>M</sup>)<br>t = 45 s</th><th>Error(u<sup>M</sup>)<br>t = 60 s</th><th>Error(σ<sup>M</sup>)<br>t = 60 s</th><th>Speed-up</th></tr></thead><tbody><tr><td>40</td><td>6.86 × 10<sup>-3</sup></td><td>6.10 × 10<sup>-2</sup></td><td>6.81 × 10<sup>-1</sup></td><td>5.55 × 10<sup>-1</sup></td><td>7.48</td></tr><tr><td>81</td><td>6.28 × 10<sup>-3</sup></td><td>6.75 × 10<sup>-2</sup></td><td>5.48 × 10<sup>-1</sup></td><td>5.39 × 10<sup>-1</sup></td><td>4.07</td></tr><tr><td>162</td><td>4.51 × 10<sup>-4</sup></td><td>5.12 × 10<sup>-3</sup></td><td>2.78 × 10<sup>-1</sup></td><td>4.79 × 10<sup>-1</sup></td><td>3.48</td></tr><tr><td>325</td><td>2.07 × 10<sup>-3</sup></td><td>2.42 × 10<sup>-2</sup></td><td>2.41 × 10<sup>-1</sup></td><td>3.29 × 10<sup>-1</sup></td><td>5.81</td></tr><tr><td>677</td><td>1.00 × 10<sup>-3</sup></td><td>1.58 × 10<sup>-2</sup></td><td>9.72 × 10<sup>-2</sup></td><td>2.21 × 10<sup>-1</sup></td><td>3.93</td></tr><tr><td>1355</td><td>6.58 × 10<sup>-4</sup></td><td>1.33 × 10<sup>-2</sup></td><td>1.63 × 10<sup>-1</sup></td><td>3.51 × 10<sup>-1</sup></td><td>2.32</td></tr></tbody></table>

this point are drawn in Fig. 8. The original method proposed in Ref. [4] could compute converged results with around 100 clusters. Clearly, the k-medoids-based PFE² method requires more clusters. The main issue is path direction change, which occurs twice in these simulations, and leads to severe inaccuracies with 40, 81 and 162 clusters. This issue was also observed in Ref. [4] when using a small number of clusters. As shown in Fig. 8, with more than 162 clusters, results seem to improve.

Errors for simulations with various numbers of clusters are reported in Table 2. Convergence rates are low for both displacement and stress errors and it is difficult to obtain accurate results for both fields at the final state. This is mainly due to the fact that both displacement and stress magnitudes are small at the final state. At $t = 45$ s, however, errors are already below 10 % for both fields using 40 clusters. This demonstrates that even though the proposed k-medoids-based PFE² method might not be appropriate to predict residual stresses, it could still be quite accurate in many situations. At $t = 45$ s, as shown in Fig. 9, convergence rates are perturbed by an unexpectedly very accurate result with 162 clusters and then an increase of errors with 325 clusters. This can be explained by fall backs to the full FE² scheme, since there were three of them for the simulation with 40 clusters, seven with 81, eight with 162, one with 1355, but none for the other simulations.

### 5.2.2. Computation time

The reference FE² simulation required 437,787 s to finish. The best compromise in terms of error and computation time (Table 2) is obtained with 325 clusters, which leads to a reduction factor of 5.81. This is far below the theoretical factor of 25. This poor performance of the k-medoids-based PFE² method is due to convergence issues. Although

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Errors at $t = 45$ s for elasto-plastic bracket simulations using an increasing number of clusters.</figcaption></figure>

the k-medoids-based PFE² simulation with 325 clusters encountered no restart nor fall back to the full FE² scheme, it required 771 coarse scale solver iterations, against 299 for the reference FE² simulation, which can explain its under-performance.

Computation times in Table 2, in addition, might seem quite unexpected. Since no simulation led to a restart, this is due to convergence issues leading to fall backs. As mentioned previously, there were three of them for the simulation with 40 clusters, seven with 81, eight with 162 and one with 1355.

The comparison with Ref. [4] is difficult. The authors mentioned computation times between 115,200 s using 80 clusters and 145,980 s using 120 clusters, which could be multiplied by eight given that the authors used parallel computing. They also mentioned a speed-up ratio of 19, but they did not actually conduct the reference FE² simulation and rather estimated its computation time to a month. As mentioned previously, the FE² simulation actually took around five days with the present implementation, and the simulation with 81 clusters took 107,358 s without parallel computing. These differences could be due to the use of different workstations and, more likely, the use of Matlab in Ref. [4] against C in the present work [21].

Another important aspect is the time spent updating internal variables by solving all fine scale problems after convergence of the partitioned coarse scale solver. With 40 clusters, the computation time spent on this task is 32 %, and with 325 clusters it is 25 %. For all other simulations, it is 15 %. As could be expected, this is highly dependent on the number of clusters and also coarse scale solver iterations. It is important, nevertheless, to observe that the majority of the total computation time is not spent on this task.

Similarly as for the hyperelastic beam setup, computation time spent clustering is negligible as it represented less than 2 % of the total computation time for all simulations. For these elasto-plastic bracket simulations which involve an internal variable, this means that a more advanced input variable would be worth considering. The averaged internal variable in Eq. (14) could potentially be replaced by the whole fine scale field.

### 5.2.3. Importance of cross contributions

Due to convergence issues mentioned above, cross contributions are assembled, with the hope that the computation time spent computing second derivative tensors will pay off. This is conducted for the setups with 40 and 81 clusters.
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Geometries, meshes and loading for double-notched test simulations: (a) coarse scale, (b) fine scale (square size of 1 mm) with colors distinguishing matrix and inclusions (radius of 0.13 mm), (c) horizontal and vertical displacements applied at A.</figcaption></figure>

Surprisingly, the number of fall backs increases with cross contributions from three to nine with 40 clusters, and from seven to eight with 81. Although fall backs are related to a stagnation of the residual when using clustering, even with frozen cluster assignments, it seems that implementing cross contributions does not always help decreasing this residual. This is confirmed by analyzing the total number of coarse scale solver iterations, which is increased from 818 to 1053 with 40 clusters and from 995 to 1071 with 81. This is detrimental to computation time, which is already affected by the computation of second derivative tensors, as it is multiplied by around four in both cases.

All in all, convergence issues of the coarse scale solver do not seem to be systematically caused by the absence of cross contributions. This is unexpected based on the promising results obtained in Section 5.1.4. The main difference here is that there is an internal variable due to the elasto-plastic behavior. More investigation is needed to understand what is the cause of convergence issues with this elasto-plastic material model. Based on the current analysis, given that second derivative tensors are very expensive to compute, it can be concluded that their implementation is not advisable.

### 5.3. Double-notched test

The third setup is the well-known double-notched test [19]. Under the assumption of plane strain, the specimen shown in Fig. 10(a) is fixed on its two bottom and right boundaries B, and subjected to the non-proportional loading path shown in Fig. 10(c) on its two top and left boundaries A. The fine scale domain is more complex than previous ones as it contains seven inclusions. The small strain elasto-plastic model presented in Section 2.4 is used. Inclusion properties are a Young's modulus of 210 GPa, a Poisson's coefficient of 0.3 and an initial yield stress large enough to prevent plasticity. Matrix properties are a Young's modulus of 2 GPa, a Poisson's coefficient of 0.3, an initial yield stress of 24 MPa, and an isotropic hardening modulus of 80 MPa. Kinematic hardening is neglected.

Loading is applied over 21 increments with a time step of 0.1 s, except for the first increment which uses a very small step to get an elastic response and initialize the NICE scheme. The mesh is composed of 849

<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Cluster assignments (colors) and medoids (highlighted elements) using 6 clusters for the double-notched test at: (a) $t = 0.5$ s, (b) $t = 0.9$ s, (c) $t = 1.5$ s, (d) $t = 2.0$ s.</figcaption></figure>

linear triangular elements for the coarse scale, and 2081 linear triangular elements for the fine scale.

#### 5.3.1. Analysis of medoid change

A first simulation is conducted with 6 clusters over the whole specimen. Cluster assignments and medoids are shown in Fig. 11. In the first stage, shown in Fig. 11(a) and (b), clusters are concentrated close to holes and in the shear band between the holes. Medoids also move towards the holes and the shear band between Fig. 11(a) and (b). Shear stress values are indeed significant in these regions, as illustrated in Fig. 12(a) and (b), which is associated with matrix plastic flow at the fine scale.

In the second stage, the stress field changes, as illustrated in Fig. 12(c) and (d), which also affects clustering assignments, as shown in Fig. 11(c) and (d). Medoids, as a result, move away from the shear band.

This double-notched test problem is clearly challenging, mainly due to matrix plastic flow at the fine scale and loading path direction change. The k-medoids clustering algorithm changes both cluster assignments and medoids during the simulation to try to represent the full complexity of the problem.

#### 5.3.2. Convergence and computation time

Errors and computation times are analyzed in two steps. First, only the first loading stage is considered. Even though there is plasticity, loading is proportional because the loading direction remains constant. The reference FE² simulation took 35,245 s. With a speed-up ratio of 5.34, the k-medoids-based PFE² method computes results with an accuracy below 10 % for both displacements and stresses using 6 clusters, as shown in Table 3. There is a fall back to the full FE² scheme in this simulation, which, by comparison to the simulation with 13 clusters, appears detrimental to computation time but beneficial to accuracy. Errors close to 10 % are also obtained with 26 clusters and no fall backs, which leads to a computation time divided by 10.3.

Second, both loading stages are considered. Loading direction change is challenging especially regarding the conservation of internal variables. The reference FE² simulation took 74,426 s. As shown in Table 4, a good accuracy is achieved for all simulations except the one with
<figure><img src="image_14.png" /><figcaption>Fig. 12. Shear stress component $\sigma_{xy}^M$ using 6 clusters for the double-notched test at: (a) $t = 0.5$ s, (b) $t = 0.9$ s, (c) $t = 1.5$ s, (d) $t = 2.0$ s.</figcaption></figure>

<caption>**Table 3**</caption>

Errors and speed-up ratios for double-notched test simulations up to $t = 1$ s using an increasing number of clusters.

<table><caption>**Table 4**</caption><thead><tr><th># clusters</th><th>Error(u<sup>M</sup>)</th><th>Error(σ<sup>M</sup>)</th><th>Speed-up</th><th># fall backs</th></tr></thead><tbody><tr><td>6</td><td>1.40 × 10<sup>-2</sup></td><td>4.51 × 10<sup>-2</sup></td><td>5.34</td><td>1</td></tr><tr><td>13</td><td>2.64 × 10<sup>-2</sup></td><td>1.37 × 10<sup>-1</sup></td><td>13.3</td><td>0</td></tr><tr><td>26</td><td>1.64 × 10<sup>-2</sup></td><td>1.08 × 10<sup>-1</sup></td><td>10.3</td><td>0</td></tr><tr><td>53</td><td>4.45 × 10<sup>-3</sup></td><td>3.03 × 10<sup>-2</sup></td><td>3.50</td><td>1</td></tr></tbody></table>

Errors and speed-up ratios for double-notched test simulations up to $t = 2$ s using an increasing number of clusters.

<table><thead><tr><th># clusters</th><th>Error(u<sup>M</sup>)</th><th>Error(σ<sup>M</sup>)</th><th>Speed-up</th><th># fall backs</th></tr></thead><tbody><tr><td>6</td><td>7.98 × 10<sup>-2</sup></td><td>3.04 × 10<sup>-1</sup></td><td>2.31</td><td>5</td></tr><tr><td>13</td><td>3.89 × 10<sup>-2</sup></td><td>9.70 × 10<sup>-2</sup></td><td>2.10</td><td>6</td></tr><tr><td>26</td><td>6.89 × 10<sup>-3</sup></td><td>2.00 × 10<sup>-2</sup></td><td>1.73</td><td>7</td></tr><tr><td>53</td><td>3.83 × 10<sup>-3</sup></td><td>1.32 × 10<sup>-2</sup></td><td>1.26</td><td>9</td></tr></tbody></table>
<figure><img src="image_15.png" /><figcaption>Fig. 13. Fine scale shear stress component $\sigma_{xy}^M$ values and absolute errors for the fine scale domain located at point E in the double-notched test: (a) FE² scheme at $t = 0.9$ s, (b) PFE² method with 13 clusters at $t = 0.9$ s, (c) absolute error between (a) and (b), (d) FE² scheme at $t = 2.0$ s, (e) PFE² method with 13 clusters at $t = 2.0$ s, (f) absolute error between (d) and (e). Note that different color scales are used for (c) and (f).</figcaption></figure>

only 6 clusters. Unfortunately, this is mainly attributed to a significant number of fall backs to the full FE² scheme after loading direction change. The speed-up ratio is hence of 2.10 when using 13 clusters.

To assess errors at the fine scale as well, the fine scale solution for the point E defined in Fig. 10(a) is shown in Fig. 13(a) and (d). Point E is the point of maximum shear stress. By comparing with the FE² reference in Fig. 13(b) and (e) and analyzing absolute errors in Fig. 13(c) and (f), it can be seen that errors are larger in inclusions as compared to the matrix. This can be explained by the fact that stress values are also larger in inclusions. All in all, errors are relatively low at the fine scale since absolute errors as shown by the color scales in Fig. 13(c) and (f) are at least seven times lower than stress values at $t = 0.9$ s, and 54 times lower at $t = 2.0$ s.

Regarding the time spent solving all fine scale problems after convergence of the partitioned coarse scale solver, it varies significantly. It is maximum for the simulations with 13 and 26 clusters in the first loading stage, where it reaches respectively 48 % and 35 % of the total computation time. It is minimum for the simulation with 53 clusters in both stages, where it represents only 3 % of the total computation time, mainly due to the number of fall backs. This shows once again that the cost of the strategy consisting of solving all fine scale problems for updating internal variables is acceptable.

### 5.3.3. Path-dependence

An important aspect of plasticity is path-dependence. In the PFE² method, this is dealt with by carefully updating internal variables at the fine scale. The latter play a key role for fine scale problems associated to newly-made medoids when clusters and/or medoids change. To evaluate how the PFE² method can capture path-dependence, simulations are conducted for the same double-notched test as in Fig. 10 but with inverted loading stages. The final total displacement at A is therefore identical, but the path is different. Final shear stress fields for the two loading paths are presented in Fig. 14(c) and (d). Comparison with reference FE² results in Fig. 14(a) and (b) reveals that with only 13 clusters, the change of the shear band's angle is predicted by the k-medoids-based PFE² method.

Cluster assignments at the end of each stage using 6 clusters are shown in Fig. 15(a) and (b). Clustering itself seems to be path-dependent. Comparing Fig. 15(a) and (b) and Fig. 11(b) and (d), clearly, there is a change in the orientation of the elongated cluster capturing the shear band.

The FE² simulation with inverted loading path finished in 85,370 s. The PFE² simulations with inverted loading path encountered less convergence issues as there were three fall backs with 6 clusters, four with 13, three with 26 and seven with 53. Speed-up ratios were consequently improved, respectively, to 4.00, 3.06, 3.19 and 1.54.
<figure><img src="image_16.png" /><figcaption>Fig. 14. Shear stress component $\sigma_{xy}^M$ for the double-notched test using: (a) reference loading path and FE$^2$ scheme, (b) inverted loading path and FE$^2$ scheme, (c) reference loading path and PFE$^2$ method with 13 clusters, (d) inverted loading path and PFE$^2$ method with 13 clusters.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 15. Cluster assignments (colors) and medoids (highlighted elements) using 6 clusters for the double-notched test with inverted loading path at: (a) $t = 0.9$ s, (b) $t = 2.0$ s.</figcaption></figure>

## 6. Conclusions

An original k-medoids-based PFE$^2$ approach has been introduced in this paper. As the computation time of FE$^2$ schemes typically scales with the number of fine scale problems to solve, this approach directly tackles the issue by effectively reducing the number of fine scale solves. At each of the coarse scale nonlinear iterations, coarse scale integration points are partitioned *a priori* based on the current coarse scale displacement gradient and fine scale internal variables using the k-medoids clustering algorithm. Fine scale solves are only triggered for cluster medoids, and obtained stresses and tangent moduli are extended to the remaining non-medoid coarse scale integration points.

Convergence issues are dealt with by combining cluster freezing, restarts with time step reduction, and fall backs to the full FE$^2$ scheme.
Cross contributions involving second derivatives and pairing each coarse scale integration point with its cluster's medoid have also been introduced in this paper. Their positive effect on convergence and the number of coarse scale nonlinear iterations, however, was not deemed sufficient to justify their employment, given the number of supplementary fine scale solves that they require. More optimal techniques to compute first and second derivatives are worth considering.

As compared to previous works, the proposed method has the advantage of solving fine scale problems associated to actual coarse scale integration points (e.g., the medoids) instead of creating virtual ones. This has the added benefit of simplifying the assembly of cross contributions. Fine scale displacement fields and internal variables, moreover, are systematically updated for all integration points. More optimal techniques to transfer internal variables when cluster assignments and medoids are changed are worth considering. One possibility could be to establish cluster assignments and medoids once and for all in a pre-processing step, for instance by considering a linear elastic analysis of the structure. Cluster assignments and medoids would hence not change during the simulation, which would drastically reduce the computational cost and memory consumption of the approach.

Numerical results with nonlinear material behavior such as hyperelasticity and elasto-plasticity with either kinematic or isotropic hardening show that the new method is accurate and is a promising candidate for reducing the computational cost of FE² simulations. Its implementation is also less intrusive and simplified as compared to previous works. Convergence issues, however, deteriorate significantly the performance of this method and should be the object of more research. Since the computation time spent clustering is negligible, more information on fine scale internal variables could also be integrated into the clustering input variable.

## CRediT authorship contribution statement

**Modesar Shakoor:** Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization.

## Funding

This work was supported by the French Agence Nationale de la Recherche (ANR) through the MISSA project (ANR-22-CE46-0002).

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

I have shared the code at https://hal.science/hal-03781711v3.

## References

[1] Balay S, Abhyankar S, Adams MF, Brown J, Brune P, Buschelman K, et al. Petsc web page. 2020. https://www.mcs.anl.gov/petsc

[2] Benaimeche MA, Yvonnet J, Bary B, He Q-C. A k-means clustering machine learning-based multiscale method for an elastic heterogeneous structures with internal variables. Int J Numer Methods Eng 2022;123(9):2012-41.

[3] Blanco PJ, Sánchez PJ, de Souza Neto EA, Feijóo RA. Variational foundations and generalized unified theory of RVE-Based multiscale models. Arch Comput Methods Eng 2016;23(2):191-253.

[4] Chaouch S, Yvonnet J. An unsupervised machine learning approach to reduce nonlinear FE² multiscale calculations using macro clustering. Finite Elem Anal Des 2024;229:104069.

[5] Chaouch S, Yvonnet J. Unsupervised machine learning classification for accelerating FE² multiscale fracture simulations. Comput Methods Appl Mech Eng 2024;432:117278.

[6] Davis TA. Algorithm 832: UMFPACK V4.3—an unsymmetric-pattern multifrontal method. ACM Trans Math Software 2004;30(2):196-99.

[7] de Hoon MJL, Imoto S, Nolan J, Miyano S. Open source clustering software. Bioinformatics 2004;20(9):1453-54.

[8] Feyel F. Multiscale FE² elastoviscoplastic analysis of composite structures. Comput Mater Sci 1999;16(1-4):344-54.

[9] Fritzen F, Kunc O. Two-stage data-driven homogenization for nonlinear solids using a reduced order model. Eur J Mech A Solids 2018;69:201-20.

[10] Gao J, Shakoor M, Domel G, Merzkirch M, Zhou G, Zeng D, et al. Predictive multiscale modeling for unidirectional carbon fiber reinforced polymers. Compos Sci Technol 2020;186:107922.

[11] Geers MGD, Kouznetsova VG, Brekelmans WAM. Multi-scale computational homogenization: trends and challenges. J Comput Appl Math 2010;234(7):2175-82.

[12] Gierden C, Kochmann J, Waimann J, Svendsen B, Reese S. A review of FE-FFT-Based two-scale Methods for computational modeling of microstructure evolution and macroscopic material behavior. Arch Comput Methods Eng 2022;29(6):4115-35.

[13] Halilović M, Vrh M, Štok B. NICE—an explicit numerical scheme for efficient integration of nonlinear constitutive equations. Math Comput Simulat 2009;80(2):294-313.

[14] Hashin Z, Shtrikman S. A variational approach to the theory of the elastic behaviour of multiphase materials. Acta Metallurgica 1963;11(2):127-40.

[15] Le BA, Yvonnet J, He Q-C. Computational homogenization of nonlinear elastic materials using neural networks. Int J Numer Methods Eng 2015;104(12):1061-84.

[16] Masi F, Stefanou I. Multiscale modeling of inelastic materials with thermodynamics-based artificial neural networks (TANN). Comput Methods Appl Mech Eng 2022;398:115190.

[17] Matouš K, Geers MGD, Kouznetsova VG, Gillman A. A review of predictive nonlinear theories for multiscale modeling of heterogeneous materials. J Comput Phys 2017;330:192-220.

[18] Mori T, Tanaka K. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metallurgica 1973;21(5):571-74.

[19] Seabra MRR, Šuštarč P, Cesar de Sa JMA, Rodić T. Damage driven crack initiation and propagation in ductile metals using XFEM. Comput Mech 2013;52(1):161-79.

[20] Shakoor M. FEMS—A mechanics-oriented finite element modeling software. Comput Phys Commun 2021;260:107729.

[21] Shakoor M. FEMS—finite element modeling software. 2022. https://hal.science/hal-03781711

[22] Shakoor M, Park CH. Computational homogenization of unsteady flows with obstacles. Int J Numer Methods Fluids 2023;95(4):499-527.

[23] Shakoor M, Park CH. A pure Stokes approach for coupling fluid flow with porous media flow. Finite Elem Anal Des 2024;231:104106.

[24] Wu L, Nguyen VD, Kilingar NG, Noels L. A recurrent neural network-accelerated multi-scale model for elasto-plastic heterogeneous materials subjected to random cyclic and non-proportional loading paths. Comput Methods Appl Mech Eng 2020;369:113234.

[25] Xu R, Yang J, Yan W, Huang Q, Giunta G, Belouettar S, et al. Data-driven multiscale finite element method: from concurrence to separation. Comput Methods Appl Mech Eng 2020;363:112893.

[26] Yvonnet J, He Q-C. The reduced model multiscale method (R3M) for the nonlinear homogenization of hyperelastic media at finite strains. J Comput Phys 2007;223(1):341-68.