# Data-driven non-linear elasticity: constitutive manifold construction and problem discretization

Ruben Ibañez, Domenico Borzacchiello, Jose Vicente Aguado, Emmanuelle Abisset-Chavanne, Elías Cueto, Pierre Ladevèze, Francisco Chinesta

## ► To cite this version:

Ruben Ibañez, Domenico Borzacchiello, Jose Vicente Aguado, Emmanuelle Abisset-Chavanne, Elías Cueto, et al.. Data-driven non-linear elasticity: constitutive manifold construction and problem discretization. Computational Mechanics, 2017, 60 (5), pp.813-826. 10.1007/s00466-017-1440-1. hal-01647641

HAL Id: hal-01647641

https://hal.science/hal-01647641

Submitted on 8 Dec 2017

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
# Data-driven non-linear elasticity: constitutive manifold construction and problem discretization

Ruben Ibañez<fnref n="1" /> · Domenico Borzacchiello<fnref n="1" /> · Jose Vicente Aguado<fnref n="1" /> ·

Emmanuelle Abisset-Chavanne<fnref n="1" /> · Elias Cueto<fnref n="2" /> · Pierre Ladeveze<fnref n="3" /> ·

Francisco Chinesta<fnref n="1" />

**Abstract** The use of constitutive equations calibrated from data has been implemented into standard numerical solvers for successfully addressing a variety of problems encountered in simulation-based engineering sciences (SBES). However, the complexity remains constantly increasing due to the need of increasingly detailed models as well as the use of engineered materials. Data-Driven simulation constitutes a potential change of paradigm in SBES. Standard simulation in computational mechanics is based on the use of two very different types of equations. The first one, of axiomatic character, is related to balance laws (momentum, mass, energy,...), whereas the second one consists of models that scientists have extracted from collected, either natural

or synthetic, data. Data-driven (or data-intensive) simulation consists of directly linking experimental data to computers in order to perform numerical simulations. These simulations will employ laws, universally recognized as epistemic, while minimizing the need of explicit, often phenomenological, models. The main drawback of such an approach is the large amount of required data, some of them inaccessible from the nowadays testing facilities. Such difficulty can be circumvented in many cases, and in any case alleviated, by considering complex tests, collecting as many data as possible and then using a data-driven inverse approach in order to generate the whole constitutive manifold from few complex experimental tests, as discussed in the present work.

**Keywords** Data-driven computational mechanics · Data-intensive simulation · Inverse problems · Constitutive manifold

## 1 Introduction

Machine and manifold learning techniques, and more specifically nonlinear dimensionality reduction, as for example locally linear embedding (LLE), kernel-PCA (the nonlinear counterpart of principal component analysis—PCA), referred as *k*-PCA, local-PCA, among many other choices, allows us to remove correlations in data [10, 17, 19–21]. Such data, free of correlation, constitute the real information, often very limited when compared with the big data from which it was extracted.

This information is then translated into knowledge, and from it to decision making. For the human being knowledge is primordial: we are interested in understanding the intimate and subtle mechanisms about the nature of things. However, when dealing with machines, these intellectual needs are not

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

Francisco Chinesta
Francisco.Chinesta@ec-nantes.fr

Ruben Ibañez
Ruben.Ibanez-Pinillo@ec-nantes.fr

Domenico Borzacchiello
Domenico.Borzacchiello@ec-nantes.fr

Jose Vicente Aguado
jose.aguado-lopez@ec-nantes.fr

Emmanuelle Abisset-Chavanne
Emmanuelle.Abisset-Chavanne@ec-nantes.fr

Elias Cueto
ecueto@unizar.es

Pierre Ladeveze
ladeveze@lmt.ens-cachan.fr

<footnote n="1">¹ ESI GROUP Chair and High Performance Computing Institute, Ecole Centrale de Nantes, 1 rue de la Noe, BP 92101, 44321 Nantes Cedex 3, France</footnote>

<footnote n="2">² Aragon Institute of Engineering Research, Universidad de Zaragoza, Maria de Luna s/n, 50018 Zaragoza, Spain</footnote>

<footnote n="3">³ LMT, ENS Paris-Saclay, 61 Avenue du President Wilson, 94230 Cachan, France</footnote>
inherent to their nature, and decisions can be made from a new kind of (artificial) intelligence that, more than based on mathematical expressions, are based on data via data mining and data analytics.

In many models, the extraction of uncorrelated parameters remains a tricky issue. It is the case of parameters describing microstructures or shapes for example, often referred to as *latent parameters*. As soon as the uncorrelated parameters are extracted, two main options have been considered to date: (1) when a new case, not included in the data, must be analyzed, its solution is simply interpolated on the manifold (constructed from the training data) from its closest neighbors [12] so that decisions can be taken in real time; and (2) an explicit parametric solution could be constructed by using the just extracted uncorrelated parameters so that it could be particularized in real-time [6, 7].

Data-Driven simulation constitutes another appealing opportunity. Furthermore, in our humble opinion, it constitutes a real change of paradigm in simulation-based engineering sciences (SBES), with plenty of potential [1, 11, 13–16, 18].

Standard simulation in classical mechanics is based on the use of two very different types of equations. The first one, of axiomatic or epistemic character, is related to balance (conservation) laws (momentum, mass, energy...). The second one consists of models extracted from collected data.

Data-driven simulation consists of directly employing data in order to perform numerical simulations. These simulations will employ universal laws while minimizing the need of explicit, often phenomenological, models. They are based on manifold learning methodologies able to extract the uncorrelated behavior of constitutive relations from a huge amount of collected data [8, 9].

This approach is especially interesting when considering complex engineered materials (meta-materials), for which constitutive relations become hard to write, because there are (too) many possible designs, and the intimate nature of most of them remains inaccessible and/or confidential.

The main drawback of such an approach is the huge amount of required data, some of them inaccessible from the nowadays testing facilities. Such difficulty can be circumvented in many cases, and in all cases alleviated, as proved in the present work, by considering complex tests, collecting as many data as possible and then using a data-driven inverse approach in order to generate the whole constitutive manifold [8], in a subtle alliance of testing machines, devices for collecting data and powerful computers for treating these huge amount of data in a variety of ways (machine and deep learning).

To better understand the data-driven rationale here addressed, let us consider, for the sake of clarity, a very simple problem: linear elasticity. In that case, the balance of (linear and angular) momenta leads to the existence of

a symmetric second-order tensor $\boldsymbol{\sigma}$ (the so-called Cauchy's stress tensor) verifying equilibrium, expressed in the absence of body forces and inertia effects, as

$$ \int_{\Omega} \boldsymbol{\varepsilon}^* : \boldsymbol{\sigma} \, d\mathbf{x} = \int_{\Gamma_N} \mathbf{u}^* \cdot \mathbf{t} \, d\mathbf{x}, \qquad (1) $$

$\forall \mathbf{u}^*$ regular enough and vanishing on $\Gamma_D$ (portion of the domain boundary $\Gamma \equiv \partial\Omega$ where the displacement is prescribed), being the tractions $\mathbf{t}$ known in the complementary boundary region $\Gamma_N$, with $\Gamma_D \cup \Gamma_N = \Gamma$ and $\Gamma_D \cap \Gamma_N = \emptyset$.

In order to solve problem (1) some relationship linking kinematic and mechanical (static) variables is required, the so-called constitutive equation. The simplest one, giving rise to isotropic linear elasticity, is known as Hooke's law (even if, more than a law, it is simply a constitutive model), and reads

$$ \boldsymbol{\sigma} = \lambda \mathrm{Tr}(\boldsymbol{\varepsilon})\mathbf{I} + \mu \boldsymbol{\varepsilon}, \qquad (2) $$

where $\mathrm{Tr}(\bullet)$ denotes the trace operator, $\boldsymbol{\varepsilon}$ is the strain tensor, and $\lambda$ and $\mu$ are the Lamé coefficients directly related to the Young modulus $E$ and the Poisson coefficient $\nu$.

By introducing the constitutive model, Eq. (2), into the weak form of the balance of momentum, Eq. (1), the so-called virtual work principle, a problem is obtained that can be formulated entirely in terms of the displacement field $\mathbf{u}$. By discretizing it, using standard finite element approximations, for instance, and performing numerically the integrals involved in Eq. (1), we finally obtain a linear algebraic system of equations, from which the nodal displacements can be obtained.

In the case of linear elasticity there is no room for discussion: the approach is simple, efficient and has been applied successfully to many problems of practical interest. Today, there are numerous commercial codes making use of this mechanical behavior and nobody doubts about its pertinence in engineering practice. However, there are other material behaviors for whom simple models fail to describe any experimental finding. These models lack of generality (universality) and due to this reason a mechanical system is usually associated to different models that are progressively adapted and/or enriched from collected data.

The biggest challenge could then be formulated as follows: can simulation proceed directly from data by circumventing, or at least alleviating, the necessity of establishing a constitutive model? In the case of linear elasticity it is obvious that such an approach lacks of interest. However, in other branches of engineering science and technology it should be an appealing alternative to standard constitutive model-based simulations.

In [8] we proposed some methods for performing simulation employing the just developed concept of constitutive manifolds arising from data. However, in that work the issue
related to the data-collection and curation was not addressed. The present work constitutes a first tentative in that direction.

As argued in [8], the main issue related to data-based approaches lies in the huge amount of data required in order to represent the mechanical behavior. Thus, constructing the constitutive manifold by carrying out a sequence of homogeneous tests with the purpose of activating all the possible strain states, seems today too expensive, but probably not in the future where data is expected playing a major role. In the present paper we consider an alternative route. If instead of performing simple tests, we consider one involving complex and evolving loads applying on a quite complex geometry. In this way numerous mechanical states will coexist in the part, and having access for example to the strain in a region of the specimen, we could by using an inverse identification strategy, identifying a large part of the constitutive manifold. This is the main idea explored in the present paper.

However, this inverse technique can be performed in different ways. In this paper we consider two strategies. The first consists in gradually constructing the manifold from data collected during the loading. Thus, at each loading step the elastic tensor for a new strain value is identified. However, such a procedure has as main drawback the fact of using the elastic tensor as main mechanical variable as well as its complexity in the case of nonlinear behaviors, as discussed later. Another appealing possibility consists of constructing a polynomial approximation of the elastic energy, whose second derivative results in the elastic tensor, and whose identification from collected data seems to be more robust.

Both strategies will be presented in Sects. 2 and 3 respectively, and then their efficiency checked as discussed in Sect. 4.

# 2 Progressive construction of the constitutive manifold

## 2.1 Linear setting

We consider first, for simplicity, mechanical tests conducted on a perfectly linear elastic material, in a specimen exhibiting uniform stresses and strains. We will later consider issues related to data generation and curation. Thus, for M randomly applied external loads, we assume ourselves able to collect M couples $(\sigma_m, \varepsilon_m)$, $m = 1, \dots, M$. Each stress-strain couple could thus be represented as a single point $P_m$ in a phase space of dimension $D = 12$ (the six distinct components of the stress and strain tensors, respectively). In the sequel Voigt notation will be considered, i.e. stress and strain tensors will be represented as vectors and consequently the fourth-order elastic tensor reduces to a $6 \times 6$ square matrix.

Each vector $P_m$ thus defines a point in a space of dimension D and, therefore, the whole set of samples represents a set of

M points in $\mathbb{R}^D$. We conjecture that all these points belong to (or can be embedded into) a certain low-dimensional manifold embedded into the high-dimensional space $\mathbb{R}^D$ allowing for a nonlinear dimensionality reduction as discussed in [8].

As soon as the elastic manifold $C(\varepsilon)$ is determined from a locally-linear interpolation, we can proceed from the standard weak form

$$ \int_{\Omega} \varepsilon^*(\mathbf{x}) : \sigma(\mathbf{x}) \, d\mathbf{x} = \int_{\Gamma_N} \mathbf{u}^*(\mathbf{x}) \cdot \mathbf{t}(\mathbf{x}) \, d\mathbf{x}, \quad (3) $$

that using Voigt notation and the behavior derived from the constitutive manifold, becomes

$$ \int_{\Omega} \varepsilon^*(\mathbf{x}) \cdot (C(\varepsilon(\mathbf{x}))\varepsilon(\mathbf{x})) \, d\mathbf{x} = \int_{\Gamma_N} \mathbf{u}^*(\mathbf{x}) \cdot \mathbf{t}(\mathbf{x}) \, d\mathbf{x}. \quad (4) $$

By using an appropriate linearization, this last expression allows one to compute (at convergence) every mechanical field.

However, as previously argued, prior to proceed with the calculations summarized above and analyzed in detail in [8], one must accomplish the construction of the so-called constitutive manifold.

Taking as reference the strain- and stress-free reference configuration of the solid, the problem can be expressed in the incremental form (particularly interesting in the nonlinear case addressed in Sect. 2.2)

$$ \int_{\Omega} \Delta \varepsilon^*(\mathbf{x}) \cdot \Delta \sigma(\mathbf{x}) \, d\mathbf{x} = \int_{\Gamma_N} \Delta \mathbf{u}^*(\mathbf{x}) \cdot \Delta \mathbf{t} \, d\mathbf{x}, \quad (5) $$

with $\Delta \varepsilon^*(\mathbf{x})$ the virtual strain field related to the kinematically admissible increment of displacement test field $\Delta \mathbf{u}^*$. By introducing the linear behavior the previous integral equation reads

$$ \int_{\Omega} \Delta \varepsilon^*(\mathbf{x}) \cdot (C_T \Delta \varepsilon(\mathbf{x})) \, d\mathbf{x} = \int_{\Gamma_N} \Delta \mathbf{u}^*(\mathbf{x}) \cdot \Delta \mathbf{t} \, d\mathbf{x}, \quad (6) $$

where the tangent matrix $C_T$ (that coincides with the secant one C because of the assumed linearity) is unknown. However, because of the assumed linear elastic behavior, it remains constant everywhere in the domain.

Using a c-term parametrization of $6 \times 6$ matrices (the more general consisting of canonical matrices fulfilling symmetry constraints) we can write

$$ C_T = \sum_{i=1}^{c} \alpha_i M_i, \quad (7) $$

with coefficients $\alpha_i$ unknown.

By introducing this tangent matrix representation into the equilibrium weak form it results
$$ \int_{\Omega} \Delta \boldsymbol{\varepsilon}^*(\mathbf{x}) \cdot \left( \left( \sum_{i=1}^{c} \alpha_i \mathbf{M}_i \right) \Delta \boldsymbol{\varepsilon}(\mathbf{x}) \right) d\mathbf{x} = \int_{\Gamma_N} \Delta u^*(\mathbf{x}) \cdot \mathbf{t} d\mathbf{x}, \quad (8) $$

whose discrete form reads

$$ \Delta \mathbf{U}^* \cdot \left( \sum_{i=1}^{c} \alpha_i \mathbf{K}_i \right) \Delta \mathbf{U} = \Delta \mathbf{U}^* \cdot \mathbf{T}, \quad (9) $$

with $\mathbf{K}_i$ the stiffness matrices corresponding to the canonical behaviors and $\Delta \mathbf{U}$ the nodal vectors of incremental displacements.

We assume that local displacements, and consequently their associated strains, are accessible (experimentally measurable) at a certain region of the domain (in general a portion of its boundary). Their associated degrees of freedom are hereafter referred to with the superscript $\bullet^O$. Thus, making use of a partition of the displacement vector $\Delta \mathbf{U}^O$ and $\Delta \mathbf{U}^H$ referring to the observable and hidden displacements, respectively, the previous discrete system reads

$$ \left( \sum_{i=1}^{c} \alpha_i \mathbf{K}_i^{\mathcal{H}\mathcal{H}} \sum_{i=1}^{c} \alpha_i \mathbf{K}_i^{\mathcal{H}O} \quad \sum_{i=1}^{c} \alpha_i \mathbf{K}_i^{\mathcal{O}\mathcal{H}} \sum_{i=1}^{c} \alpha_i \mathbf{K}_i^{\mathcal{O}O} \right) \begin{pmatrix} \Delta \mathbf{U}^{\mathcal{H}} \\ \Delta \mathbf{U}^{\mathcal{O}} \end{pmatrix} = \begin{pmatrix} \mathbf{T}^{\mathcal{H}} \\ \mathbf{T}^{\mathcal{O}} \end{pmatrix}. \quad (10) $$

This system of equations is obviously complemented with appropriate Dirichlet boundary conditions on $\Gamma_D$. In the previous algebraic system, vectors $\mathbf{T}^O$ and $\mathbf{T}^H$ refer to the nodal traction contributions at nodes related to the observable and hidden displacements, respectively.

The algebraic system (10) has as unknowns the hidden displacements $\Delta \mathbf{U}^H$ and the constitutive coefficients $\alpha_i$, being known the observable displacements $\Delta \mathbf{U}^O$. If the number of known displacements that corresponds with the size of vector $\Delta \mathbf{U}^O$ is large enough (in all cases larger than the number of alpha-coefficients, c) it is thus possible to solve the resulting nonlinear algebraic problem to compute both the unknown displacements $\Delta \mathbf{U}^H$ and the coefficients defining the material behavior $\alpha_i$. In the opposite case it is always possible to apply some regularization to solve the undetermined resulting problem (e.g. Tikhonov regularization). In the sequel we focus in the former scenario.

System (10) can be rewritten as follows

$$ \left( \sum_{i=1}^{c} \alpha_i \mathbf{K}_i^{\mathcal{H}\mathcal{H}} \mathbf{K}_1^{\mathcal{H}O} \Delta \mathbf{U}^O \dots \mathbf{K}_c^{\mathcal{H}O} \Delta \mathbf{U}^O \quad \sum_{i=1}^{c} \alpha_i \mathbf{K}_i^{\mathcal{O}\mathcal{H}} \mathbf{K}_1^{\mathcal{O}O} \Delta \mathbf{U}^O \dots \mathbf{K}_c^{\mathcal{O}O} \Delta \mathbf{U}^O \right) \begin{pmatrix} \Delta \mathbf{U}^H \\ \alpha_1 \\ \vdots \\ \alpha_c \end{pmatrix} = \begin{pmatrix} \mathbf{T}^H \\ \mathbf{T}^O \end{pmatrix}, \quad (11) $$

or, by defining vector $\boldsymbol{\alpha}$ and matrices $\boldsymbol{\kappa}^{\mathcal{H}O}$ and $\boldsymbol{\kappa}^{\mathcal{O}O}$ as

$$ \left\{ \begin{array}{l} \boldsymbol{\alpha} = (\alpha_1, \dots, \alpha_c)^T \\ \boldsymbol{\kappa}^{\mathcal{H}O} = (\mathbf{K}_1^{\mathcal{H}O} \Delta \mathbf{U}^O, \dots, \mathbf{K}_c^{\mathcal{H}O} \Delta \mathbf{U}^O), \\ \boldsymbol{\kappa}^{\mathcal{O}O} = (\mathbf{K}_1^{\mathcal{O}O} \Delta \mathbf{U}^O, \dots, \mathbf{K}_c^{\mathcal{O}O} \Delta \mathbf{U}^O) \end{array} \right. \quad (12) $$

the previous system can be rewritten as

$$ \left( \sum_{i=1}^{c} \alpha_i \mathbf{K}_i^{\mathcal{H}\mathcal{H}} \boldsymbol{\kappa}^{\mathcal{H}O} \quad \sum_{i=1}^{c} \alpha_i \mathbf{K}_i^{\mathcal{O}\mathcal{H}} \boldsymbol{\kappa}^{\mathcal{O}O} \right) \begin{pmatrix} \Delta \mathbf{U}^{\mathcal{H}} \\ \boldsymbol{\alpha} \end{pmatrix} = \begin{pmatrix} \mathbf{T}^{\mathcal{H}} \\ \mathbf{T}^{\mathcal{O}} \end{pmatrix}, \quad (13) $$

that represents an overdetermined nonlinear algebraic system.

By premultiplying by the transpose of the matrix, a square algebraic system is obtained,

$$ \left( \sum_{i=1}^{c} \alpha_i \mathbf{K}_i^{\mathcal{H}\mathcal{H}} \boldsymbol{\kappa}^{\mathcal{H}O} \right)^T \left( \sum_{i=1}^{c} \alpha_i \mathbf{K}_i^{\mathcal{H}\mathcal{H}} \boldsymbol{\kappa}^{\mathcal{H}O} \right) \left( \sum_{i=1}^{c} \alpha_i \mathbf{K}_i^{\mathcal{O}\mathcal{H}} \boldsymbol{\kappa}^{\mathcal{O}O} \right) \begin{pmatrix} \Delta \mathbf{U}^{\mathcal{H}} \\ \boldsymbol{\alpha} \end{pmatrix} = \left( \sum_{i=1}^{c} \alpha_i \mathbf{K}_i^{\mathcal{H}\mathcal{H}} \boldsymbol{\kappa}^{\mathcal{H}O} \right)^T \begin{pmatrix} \mathbf{T}^{\mathcal{H}} \\ \mathbf{T}^{\mathcal{O}} \end{pmatrix}, \quad (14) $$

that allows us to calculate $\Delta \mathbf{U}^H$ and $\boldsymbol{\alpha}$ by using an adequate nonlinear solver (e.g. fixed point, Newton, etc.).

When considering a linear behavior the resulting displacements, strains and stresses can easily be derived from $\mathbf{U} \equiv \Delta \mathbf{U} = (\Delta \mathbf{U}^H, \Delta \mathbf{U}^{\mathcal{O}})^T$ by considering

$$ \left\{ \begin{array}{l} \mathbf{u}(\mathbf{x}) = \sum_{i=1}^{n_d} \mathbf{U}_i N_i(\mathbf{x}) \\ \boldsymbol{\varepsilon} = \nabla_s \mathbf{u} \\ \boldsymbol{\sigma} = \mathbf{C} \boldsymbol{\varepsilon} \end{array} \right. , \quad (15) $$

where $\mathbf{C} \equiv \mathbf{C}_T$, $n_d$ is the number of nodes considered to approximate the displacement field $\mathbf{u}(\mathbf{x})$ and $N_i(\mathbf{x})$ the associated shape functions. $\nabla_s(\bullet)$ denotes the symmetric component of the gradient operator and $\mathbf{C}_T$ results from $\alpha_i$

$$ \mathbf{C}_T = \sum_{i=1}^{c} \alpha_i \mathbf{M}_i. \quad (16) $$

## 2.2 Nonlinear elastic behavior

In the nonlinear case a major difficulty appears: since the behavior depends on strain, and it can be different at each physical point $\mathbf{x} \in \Omega$, the procedure just proposed and described to address the linear case must be adapted accordingly.
The external traction **t** is progressively applied, that is

$$ \mathbf{t} = \sum_{j=1}^{\ell} \Delta \mathbf{t}_j, \qquad (17) $$

with traction increments small enough to ensure the accuracy of the identified behavior. Thus, from the stress- and strain-free reference state the application of the first traction increment $\Delta \mathbf{t}_1$ results in the equilibrium weak form

$$ \int_{\Omega} \boldsymbol{\varepsilon}^*(\mathbf{x}) \cdot \boldsymbol{\sigma}_1(\mathbf{x}) \, d\mathbf{x} = \int_{\Gamma_N} \mathbf{u}^*(\mathbf{x}) \cdot \mathbf{t}_1 \, d\mathbf{x}, \qquad (18) $$

or its incremental counterpart taking into account that $\boldsymbol{\sigma}_1 = \boldsymbol{\sigma}_0 + \Delta\boldsymbol{\sigma}_1$ (with $\boldsymbol{\sigma}_0 = \mathbf{0}$) and $\mathbf{t}_1 = \mathbf{t}_0 + \Delta\mathbf{t}_1$ (with $\mathbf{t}_0 = \mathbf{0}$)

$$ \int_{\Omega} \Delta \boldsymbol{\varepsilon}^*(\mathbf{x}) \cdot \Delta \boldsymbol{\sigma}_1(\mathbf{x}) \, d\mathbf{x} = \int_{\Gamma_N} \Delta \mathbf{u}^*(\mathbf{x}) \cdot \Delta \mathbf{t}_1 \, d\mathbf{x}, \qquad (19) $$

whose linearized form writes

$$ \int_{\Omega} \Delta \boldsymbol{\varepsilon}^*(\mathbf{x}) \cdot (\mathbf{C}_{T_1} \Delta \boldsymbol{\varepsilon}_1(\mathbf{x})) \, d\mathbf{x} = \int_{\Gamma_N} \mathbf{u}^*(\mathbf{x}) \cdot \Delta \mathbf{t}_1 \, d\mathbf{x}, \qquad (20) $$

where the tangent matrix $\mathbf{C}_{T_1}$, assumed unknown, can be considered almost constant everywhere in the domain as soon as the first traction increment is taken small enough to ensure that this first loading produces a linear response everywhere in the domain $\Omega$, i.e. the first tests considered for starting the construction of the behavior manifold should avoid the appearance of stress or strain localization.

Applying the same rationale that was employed in the linear case, we consider the c-term parametrization

$$ \mathbf{C}_{T_1} = \sum_{i=1}^{c} \alpha_i^1 \mathbf{M}_i, \qquad (21) $$

that, together with

$$ \left\{ \begin{array}{l} \boldsymbol{\alpha}^1 = (\alpha_1^1, \dots, \alpha_c^1)^T \\ \boldsymbol{\kappa}_1^{\mathcal{H}\mathcal{O}} = (\mathbf{K}_1^{\mathcal{H}\mathcal{O}} \Delta \mathbf{U}_1^{\mathcal{O}}, \dots, \mathbf{K}_c^{\mathcal{H}\mathcal{O}} \Delta \mathbf{U}_1^{\mathcal{O}}), \\ \boldsymbol{\kappa}_1^{\mathcal{O}\mathcal{O}} = (\mathbf{K}_1^{\mathcal{O}\mathcal{O}} \Delta \mathbf{U}_1^{\mathcal{O}}, \dots, \mathbf{K}_c^{\mathcal{O}\mathcal{O}} \Delta \mathbf{U}_1^{\mathcal{O}}) \end{array} \right. , \qquad (22) $$

leads to

$$ \left( \sum_{i=1}^{c} \alpha_i^1 \mathbf{K}_i^{\mathcal{H}\mathcal{H}} \quad \boldsymbol{\kappa}_1^{\mathcal{H}\mathcal{O}} \right)^T \left( \sum_{i=1}^{c} \alpha_i^1 \mathbf{K}_i^{\mathcal{H}\mathcal{H}} \quad \boldsymbol{\kappa}_1^{\mathcal{H}\mathcal{O}} \right) \left( \Delta \mathbf{U}_1^{\mathcal{H}} \quad \boldsymbol{\alpha}^1 \right) \\ = \left( \sum_{i=1}^{c} \alpha_i^1 \mathbf{K}_i^{\mathcal{O}\mathcal{H}} \quad \boldsymbol{\kappa}_1^{\mathcal{O}\mathcal{O}} \right)^T \left( \Delta \mathbf{T}_1^{\mathcal{H}} \quad \Delta \mathbf{T}_1^{\mathcal{O}} \right), \qquad (23) $$

that allows us to calculate $\Delta \mathbf{U}_1^{\mathcal{H}}$ and $\boldsymbol{\alpha}^1$ and, from them, displacements, strains and stresses, according to

$$ \left\{ \begin{array}{l} \mathbf{U}_1 = \mathbf{U}_0 + \Delta \mathbf{U}_1 = \Delta \mathbf{U}_1 \\ \boldsymbol{\varepsilon}_1 = \boldsymbol{\varepsilon}_0 + \Delta \boldsymbol{\varepsilon}_1 = \Delta \boldsymbol{\varepsilon}_1 \\ \boldsymbol{\sigma}_1 = \boldsymbol{\sigma}_0 + \Delta \boldsymbol{\sigma}_1 = \boldsymbol{\sigma}_0 + \mathbf{C}_{T_1} \Delta \boldsymbol{\varepsilon}_1 = \mathbf{C}_{T_1} \Delta \boldsymbol{\varepsilon}_1 \end{array} \right. . \qquad (24) $$

This last equation makes use of the constitutive matrix $\mathbf{C}_{T_1}$

$$ \mathbf{C}_{T_1} = \sum_{i=1}^{c} \alpha_i^1 \mathbf{M}_i. \qquad (25) $$

Consider now a second loading step $\Delta \mathbf{t}_2$. The process is repeated to calculate the sequence $(\boldsymbol{\varepsilon}_2, \boldsymbol{\sigma}_2), \dots, (\boldsymbol{\varepsilon}_\ell, \boldsymbol{\sigma}_\ell)$. However, in the second and subsequent iterations the situation is a bit different with respect to the first one just described, deserving some additional comments.

After the first loading step, in which a uniform stress-strain state was assumed due to the small traction increment, now the different points in $\Omega$ will be subject to a non-uniform strain field and consequently we cannot assume that a constant $\mathbf{C}_{T_2}$ will apply in the whole domain. Because the very small magnitude of the applied loading increments only the points having the maximum deformation energy at the previous iteration are potential candidates to exhibit at the present iteration a tangent behavior different to one of the previously identified. Thus, two groups of finite elements are considered: (1) first, the ones whose stress-strain couple remains close enough to any of the ones previously identified, and (2) the ones that do not fulfill that condition. For the first group the already identified tangent behavior can be considered (up to a tolerance) still valid, whereas for the second group the strains are checked to verify if they are close enough among them. If it is the case, a common unknown tangent behavior [parametrized according to Eq. (7)] is assigned to all them. This constitutive clustering can be performed by using an appropriate classifier. In our case we employed k-means with two populations.

Assume therefore that the domain $\Omega$ can be decomposed in two parts $\Omega^\mathcal{U}$ and $\Omega^\mathcal{A}$, the former involving elements belonging to the cluster whose behavior is assumed unknown (thus far from all the behaviors already identified) and the last the one that concerns elements whose behavior, already identified, is assumed to remain valid. The linearized equilibrium at the loading increment $j$ can be written accordingly as

$$ \int_{\Omega^\mathcal{A}} \Delta \boldsymbol{\varepsilon}^*(\mathbf{x}) \cdot (\mathbf{C}_T(\mathbf{x}) \Delta \boldsymbol{\varepsilon}_j(\mathbf{x})) \, d\mathbf{x} \\ + \int_{\Omega^\mathcal{U}} \Delta \boldsymbol{\varepsilon}^*(\mathbf{x}) \cdot (\mathbf{C}_{T_j} \Delta \boldsymbol{\varepsilon}_j(\mathbf{x})) \, d\mathbf{x} \\ = \int_{\Gamma_N} \mathbf{u}^*(\mathbf{x}) \cdot \Delta \mathbf{t}_j \, d\mathbf{x}. \qquad (26) $$
The first integral defines a linear contribution, whereas the second one remains nonlinear because it involves the unknown displacement vector as well as the unknown constitutive coefficients, grouped into the tangent matrix $C_{T_j}$.

### 3 Polynomial approximation of the constitutive manifold

As proposed in [4] a simple and still appealing possibility to describe the constitutive manifold consists of approximating it in an adequate polynomial basis. The simplest alternative consists of approximating the elastic energy (as a function of the strain) whose first derivative results in the stress tensor and the second one leads to the elastic tensor.

Proceeding to identify the energy density functional seems to be a better alternative than identifying the elastic tensor, for two important reasons. The first is that it ensures thermomechanical consistency, and thus all symmetries associated to the material behavior. The second reason is based on the fact that the polynomial approximation of all the components of the elastic tensor is much more expensive computationally than the approximation of a single scalar function, the energy in the present case, and the identification procedure much more robust from a computational viewpoint.

The choice of the approximation basis deserves some comments. Imagine for a while the approximation of a one-dimensional function $f(\xi)$ in $\mathcal{I} = [\xi^{-}, \xi^{+}]$. A natural possibility consists of using piecewise continuous linear functions $N_i(\xi)$ to define its approximation, as it is usual within the finite element framework, by considering a mesh composed of $\mathfrak{q}$ nodes uniformly distributed in $\mathcal{I}$, with coordinates $\xi_i$, $i = 1, \dots, \mathfrak{q}$ ($\xi_1 = \xi^{-}$ and $\xi_{\mathfrak{q}} = \xi^{+}$), from which the approximation reads

$$f(\xi) = \sum_{i=1}^{\mathfrak{q}} f(\xi_i) N_i(\xi), \qquad (27)$$

where $N_i(\xi)$, for $1 < i < \mathfrak{q}$ writes

$$N_i(\xi) = \begin{cases} \frac{\xi - \xi_{i-1}}{\xi_i - \xi_{i-1}} & \text{if } \xi \in [\xi_{i-1}, \xi_i] \\ \frac{\xi_{i+1} - \xi}{\xi_{i+1} - \xi_i} & \text{if } \xi \in [\xi_i, \xi_{i+1}], \\ 0 & \text{elsewhere} \end{cases}, \qquad (28)$$

$$N_1(\xi) = \begin{cases} \frac{\xi_2 - \xi}{\xi_2 - \xi_1} & \text{if } \xi \in [\xi_1, \xi_2], \\ 0 & \text{elsewhere} \end{cases}, \qquad (29)$$

and

$$N_{\mathfrak{q}}(\xi) = \begin{cases} \frac{\xi - \xi_{\mathfrak{q}-1}}{\xi_{\mathfrak{q}} - \xi_{\mathfrak{q}-1}} & \text{if } \xi \in [\xi_{\mathfrak{q}-1}, \xi_{\mathfrak{q}}], \\ 0 & \text{elsewhere} \end{cases}. \qquad (30)$$

If the solution is known at different positions $\Xi_j$, $j = 1, \dots, j$, Eq. (27) will read

$$f(\Xi_j) = \sum_{i=1}^{\mathfrak{q}} f(\xi_i) N_i(\Xi_j), \quad j = 1, \dots, j, \qquad (31)$$

that results in the linear system

$$\begin{pmatrix} N_1(\Xi_1) & \dots & N_{\mathfrak{q}}(\Xi_1) \\ \vdots & \ddots & \vdots \\ N_1(\Xi_j) & \dots & N_{\mathfrak{q}}(\Xi_j) \end{pmatrix} \begin{pmatrix} f(\xi_1) \\ \vdots \\ f(\xi_{\mathfrak{q}}) \end{pmatrix} = \begin{pmatrix} f(\Xi_1) \\ \vdots \\ f(\Xi_j) \end{pmatrix}. \qquad (32)$$

At this point, different situations can be found:

- An undetermined system if $j < \mathfrak{q}$;

- A determined one, if $j = \mathfrak{q}$;

- An overdetermined one, if $j > \mathfrak{q}$. However, even when $j \ge \mathfrak{q}$ the resulting system can become undetermined if at least for one node $\xi_i$, no point $\Xi_j$, $\forall j$, falls in its support, $[\xi_{i-1}, \xi_{i+1}]$.

In these circumstances different algebraic solutions exist (e.g. pseudo-inverse, “matlab backslash”, $L^2$ or $L^1$ optimization, ...). However, in this work we decided to consider global approximation functions in $[\xi^{-}, \xi^{+}]$. To avoid the issues related to high-order Lagrange approximations, we consider approximations based on the use of orthogonal polynomials, and more precisely Chebyshev polynomials. Thus, Eq. (27) is replaced by

$$f(\xi) = \sum_{i=1}^{\mathfrak{q}} \gamma_i T_i(\xi), \qquad (33)$$

where $T_i(\xi)$ refer to Chebyshev polynomials and the weights $\gamma_i$ are computed from its associated linear system

$$f(\Xi_j) = \sum_{i=1}^{\mathfrak{q}} \gamma_i T_i(\Xi_j), \quad j = 1, \dots, j, \qquad (34)$$

where singularity issues are circumvented as soon as $j \ge \mathfrak{q}$ and there are not repeated points.

However, problems arise as soon as the approximation becomes multidimensional. This is the case when approximating the elastic energy $\psi$ as a function of the 6 components of the strain tensor $\boldsymbol{\epsilon}$, using the same degree ($\mathfrak{q}$) for each component. In this case, the approximation

$$\psi(\boldsymbol{\epsilon}) \approx \sum_{ijklmn}^{\mathfrak{q}^6} \gamma_{ijklmn} T_i(\boldsymbol{\epsilon}_{11}) T_j(\boldsymbol{\epsilon}_{12}) T_k(\boldsymbol{\epsilon}_{13}) \\ T_l(\boldsymbol{\epsilon}_{22}) T_m(\boldsymbol{\epsilon}_{23}) T_n(\boldsymbol{\epsilon}_{33}), \qquad (35)$$
contains too many coefficients $\gamma_{ijklmn}$ (in fact $\mathbb{Q}^6$), and consequently the accuracy requires the same number of data points (even if sparse sampling could be an appealing alternative). Of course the approximation could be limited to a certain degree $\mathcal{D}$ by considering in the previous sum indexes verifying $i + j + k + l + m + n \le \mathcal{D}$.

An alternative approximation makes use of a separated representation (usually considered within the proper generalized decomposition (PGD) framework [2,3]) that reads

$$ \psi(\boldsymbol{\epsilon}) \approx \sum_{i}^{\mathbb{P}} E_{i}^{11}(\boldsymbol{\epsilon}_{11}) E_{i}^{12}(\boldsymbol{\epsilon}_{12}) E_{i}^{13}(\boldsymbol{\epsilon}_{13}) \\ E_{i}^{22}(\boldsymbol{\epsilon}_{22}) E_{i}^{23}(\boldsymbol{\epsilon}_{23}) E_{i}^{33}(\boldsymbol{\epsilon}_{33}). \quad (36) $$

This separated representation is specially appropriate when $\boldsymbol{\epsilon}$ is defined in the hyper-hexahedral domain $\mathcal{E} = [\boldsymbol{\epsilon}_{11}, \boldsymbol{\epsilon}_{11}^-] \times [\boldsymbol{\epsilon}_{12}, \boldsymbol{\epsilon}_{12}^+] \times \cdots \times [\boldsymbol{\epsilon}_{33}, \boldsymbol{\epsilon}_{33}^+]$. However, admissible deformations imply non separable domains. The application of separated representation in non-separable domains was deeply addressed in [5] where the use of $R$-functions succeeded to represent complex non-separable geometries.

To avoid singularity issues, functions $E_i^{kl}$, are approximated by using global Chebyshev polynomials, according to

$$ E_i^{kl}(\boldsymbol{\epsilon}_{kl}) \approx \sum_{j=1}^{\mathbb{Q}kl} \gamma_j^{kl,i} T_j(\boldsymbol{\epsilon}_{kl}). \quad (37) $$

Starting from the weak form

$$ \int_{\Omega} \boldsymbol{\epsilon}^* \cdot \boldsymbol{\sigma} \, d\mathbf{x} = \int_{\Gamma_N} \mathbf{u}^* \cdot \mathbf{t} \, d\mathbf{x}, \quad (38) $$

we substitute the constitutive relationship

$$ \int_{\Omega} \boldsymbol{\epsilon}^* \cdot \mathbf{C} \boldsymbol{\epsilon} \, d\mathbf{x} = \int_{\Gamma_N} \mathbf{u}^* \cdot \mathbf{t} \, d\mathbf{x}, \quad (39) $$

with the elastic tensor expressed as the second order derivative of the energy.

Note that functions $E_i^{kl}$, as well as the unknown nodal displacements, should be computed from the knowledge of the measurable nodal displacements accessible in a part of the domain $\Omega$, as was the case in the procedures discussed previously.

As in the case of the PGD constructor, we consider a greedy algorithm that computes sequentially these functions [3]. Thus at iteration $n$, $n < \mathbb{P}$, we assume that the rank-$n$ approximation of the elastic energy $\psi^n$ was already computed, i.e.,

$$ \psi(\boldsymbol{\epsilon}) \approx \psi^n(\boldsymbol{\epsilon}) = \sum_{i}^{n} E_{i}^{11}(\boldsymbol{\epsilon}_{11}) \cdots E_{i}^{33}(\boldsymbol{\epsilon}_{33}). \quad (40) $$

At present iteration we look for the new functional product leading to the updated enriched rank-$n + 1$ expression of $\psi^{n+1}(\boldsymbol{\epsilon})$ from

$$ \begin{aligned} \psi^{n+1}(\boldsymbol{\epsilon}) &= \psi^n(\boldsymbol{\epsilon}) + E_{n+1}^{11}(\boldsymbol{\epsilon}_{11}) \cdots E_{n+1}^{33}(\boldsymbol{\epsilon}_{33}) \\ &= \psi^n(\boldsymbol{\epsilon}) + \Delta\psi(\epsilon), \end{aligned} \quad (41) $$

that introduced into the weak form results

$$ \int_{\Omega} \boldsymbol{\epsilon}^* \cdot (\mathbf{C}^n(\boldsymbol{\epsilon}) + \Delta\mathbf{C}(\boldsymbol{\epsilon})) \boldsymbol{\epsilon} \, d\mathbf{x} = \int_{\Gamma_N} \mathbf{u}^* \cdot \mathbf{t} \, d\mathbf{x}. \quad (42) $$

where $\Delta\mathbf{C}(\epsilon)$ results from the second derivative of the energy enrichment $\Delta\psi(\boldsymbol{\epsilon})$.

As is the case when applying the PGD solver, the solution procedure consists of using an alternated direction fixed point strategy, that proceeds as follows [3]:

1. By considering $E_{n+1}^{12(r-1)}$, ..., $E_{n+1}^{33(r-1)}$ from the previous fixed point iteration $r-1$ of the nonlinear solver (initialized at $r=1$ from the functions at the previous enrichment iteration $n$), we compute $E_{n+1}^{11(r)}$.

2. The process is repeated but now with $E_{n+1}^{11(r)}$, $E_{n+1}^{13(r-1)}$, ..., $E_{n+1}^{33(r-1)}$ known. This allows to compute the unknown nodal displacements and functions involved in $E_{n+1}^{12(r)}$. The process is repeated for all the other components until computing $E_{n+1}^{33(r)}$. Then, the fixed point convergence is checked and if it is not attained we move to the next fixed point iteration $r+1$.

3. When reaching the fixed point convergence, the enrichment convergence is evaluated and if it is not attained we move to the next elastic energy approximate $\psi^{n+2}$ from the just competed $\psi^{n+1}$. We assume that at iteration $\mathbb{P}$ the enrichment process converges and consequently we have access to the elastic tensor manifold from which simulations can be carried out as described in [8].

**Remark** It is important to note that at each solution step in the fixed point loop, the unknowns are the unknown nodal displacements as well as the nodal variables related to the approximation of functions involved in $E_i^{kl}$.

## 4 Numerical results

To illustrate the capabilities of the just described procedure, we consider the simple mechanical problem depicted in Fig. 1. It consists of a two-dimensional unit squared solid, $\mathbf{x} = (x, y) \in \Omega = (0, 1) \times (0, 1)$, equipped with a nonlinear elastic material, clamped along its basis $y=0$, free of traction on its lateral boundaries $x=0$ and $x=1$ and with a uniformly distributed traction $\mathbf{t}$ on its upper boundary $y=1$. We analyze the performance of the just presented strategies.
<figure><img src="image_3.png" /><figcaption>Fig. 1 Schema of the considered mechanical problem</figcaption></figure>

## 4.1 Progressive construction of the behavior manifold

In the present case, we consider an applied traction whose orientation, i.e. $\mathbf{t} = t\mathbf{p}$, is arbitrary: $\mathbf{p}(\theta) = (\cos\theta, \sin\theta)^T$, $\theta \in [0, 2\pi)$. As just discussed, this traction is applied incrementally in magnitude and orientation. From

$$t = \sum_{j=1}^{\ell} \Delta t_j, \qquad (43)$$

we can define an intermediate traction magnitude at step $r$, $t_r$, from

$$t_r = \sum_{j=1}^{r} \Delta t_j, \qquad (44)$$

that leads to different tractions depending on the orientation

$$\mathbf{t}_r^s = \sum_{j=1}^{r} \Delta t_j \mathbf{p}(\theta^s), \qquad (45)$$

$$\theta^s = \sum_{m=1}^{s} \Delta\theta, \qquad (46)$$

with $M$ defining the angular discretization, that is the number of discrete angles considered,

$$M \Delta\theta = 2\pi. \qquad (47)$$

Thus, for each intermediate traction magnitude, the whole orientation space is fully swept before incrementing the traction magnitude, to better explore the constitutive manifold.

To discretize the mechanical problem, the domain $\Omega$ was equipped with a uniform mesh consisting of $P \times P$ square finite elements, where a bilinear approximation of the displacement field was considered. The displacement is

therefore assumed to be experimentally measurable at each finite element node located on the upper-boundary, that is, it is assumed measurable at the $P + 1$ nodes located on its upper boundary $y = 1$.

### 4.1.1 Synthetic generation of displacement measures

In order to generate pseudo-experimental displacement measurements, we consider a nonlinear elastic behavior of the type (Voigt notation is employed here)

$$C = \frac{E}{1 - \nu^2} \begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1-\nu}{2} \end{bmatrix}, \qquad (48)$$

with the elastic coefficients given by

$$\begin{cases} E = E_0 + E_1 \mathrm{Tr}(\boldsymbol{\epsilon}) \\ \nu = \nu_0 + \nu_1 \mathrm{Tr}(\boldsymbol{\epsilon}) \end{cases}, \qquad (49)$$

with $E_0$, $E_1$, $\nu_0$ and $\nu_1$ positive constants and where $\mathrm{Tr}(\bullet)$ refers to the trace operator acting on tensor $\bullet$. Coefficients $\nu_0$ and $\nu_1$ where selected such that $\nu \in (0, 0.5)$ in the range of deformations considered.

In the numerical example discussed below the material coefficients were selected as $E_0 = 10$, $\nu_0 = 0.1$, $E_1 = 10$ and $\nu_1 = 0.1$. The applied tension was $t = 0.1$ and it was applied by considering 10 loading steps, i.e. $\ell = 10$ and 10 orientations, i.e. $M = 10$. The mesh consisted of $10 \times 10$ $Q1$ finite elements.

The standard finite element solution of the resulting non-linear model allowed the calculation of the displacement at each loading step at each of the 11 nodes located on the upper-boundary, $y = 1$.

### 4.1.2 Unveiling the constitutive manifold

The fact that the constitutive law employed to generate pseudo-experimental displacements was known is now forgotten, and the behavior is assumed unknown from now on. The main objective is therefore to determine the constitutive manifold of the material, that is, its sampling stress-strain couples, with the only information provided by the mechanical test illustrated in Fig. 1 and the recorded displacements at the 11 locations at each loading step.

For this purpose we proceed as described in Sect. 3. Figures 2, 3 and 4 compare the different identified components of the stress tensor, $\boldsymbol{\sigma}_{xx}$, $\boldsymbol{\sigma}_{yy}$ and $\boldsymbol{\sigma}_{xy}$ and the reference ones obtained from Eq. (48). The stress magnitude in those figures is represented from the color bar. These figures reveal an almost perfect stress-strain couple match with the pseudo-experimental ones, with relative errors lower than 1%.
<figure><img src="image_4.png" /><figcaption>Fig. 2 $\sigma_{xx} = \sigma_{xx}(\boldsymbol{e})$: identified from data (left) and reference values (right)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3 $\sigma_{yy} = \sigma_{yy}(\boldsymbol{e})$: identified from data (left) and reference values (right)</figcaption></figure>

It was proved, that as expected, by decreasing the loading step, that is, by increasing $\ell$ and $M$ in the loadings expressed from Eqs. (43) and (47), the error with respect to the reference one [related to the constitutive Eq. (48)] decreases proving the expected convergence of the proposed inverse identification strategy as Fig. 5 reveals. The accuracy was also checked by comparing the identified components of the tangent matrix with the analytical one, and as shown in Fig. 6 the results almost match as soon as the sampling (loading increments) becomes fine enough. Errors fewer than few percent using the norm $\|C^{\text{identified}} - C\|_2$, are easily reachable (Fig. 5).

To further explore the method, we decided to apply a non-linear dimensionality reduction technique to the stress-

strain couples just obtained. By applying on them Locally Linear Embedding nonlinear dimensionality reduction strategy, see Fig. 7, we compare the dimensionality of the resulting linear and nonlinear constitutive manifolds. In the linear case two parameters seemed to be enough for visualizing and parametrizing the constitutive data (this number corresponds with the number of lowest eigenvalues before reaching the typical plateau of LLE techniques [19], see Fig. 7). However, when considering the manifold that results from the identified stress-strain couples describing the nonlinear case, the dimensionality seems to increase to three parameters. This is natural since the nonlinear behavior implies the need for more complex descriptions.
<figure><img src="image_6.png" /><figcaption>**Fig. 4** $\sigma_{xy} = \sigma_{xy}(\boldsymbol{\epsilon})$: identified from data (**left**) and reference values (**right**)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 5** Evolution of the error with the loading steps $\ell$</figcaption></figure>

Alternatively, we employed $k$-PCA nonlinear dimensionality reduction [7], that allows to visualize low-dimensional manifolds within a higher dimensional space. When applying $k$-PCA to the identified data corresponding to the nonlinear behavior, we obtain the embedding depicted in Fig. 8. Here, in order to prove that the embedded, low-dimensional data is well distributed on the slow manifold, we assigned a color to each data point corresponding to its elastic energy. In order to prove that the reduced data define an almost perfect 2D manifold, we represent in Figs. 9 and 10 two different views of the solution shown in Fig. 8.

### 4.1.3 Data-driven simulation

The nonlinear elastic problem is now solved by employing the constitutive manifold just identified, when a traction $\mathbf{t} = t\mathbf{p}$, $\mathbf{p}^T = (\cos 3\pi/2, \sin 3\pi/2)$ applies on the top boundary.

The reference displacement field calculated with the constitutive model (48)–(49) is depicted in Fig. 11 and compared with the one obtained when solving the same problem but now with the identified constitutive model whose solution is depicted in Fig. 12. Both results are in good agreement despite the coarse descriptions considered.

## 4.2 Polynomial approximation of the constitutive manifold

The second proposed procedure consisted on the polynomial representation of the constitutive manifold described in Sect. 3. In this case, the identified behavior was in perfect agreement with the reference one considered above. Figure 13 represents the relative error in the component $C_{11}$. It can be noticed that maximum relative errors remain again lower than 1%.

When considering $\forall k, l, q_{kl} = q = 5$ and only one load applied on the top boundary $\mathbf{t} = 0.1\mathbf{p}$, with $\mathbf{p} = (-1/\sqrt{2}, -1/\sqrt{2})$, and by assuming that our pseudo-experimental technique is able to provide us with nodal displacement values in the 25% of the nodal locations in the model, a perfect agreement was obtained between the identified and the reference behavior as proved in Fig. 14.

The reduction of the number of measured displacements requires the use of additional loading test cases. We also proved that the convergence is significantly enhanced with the number of considered loading cases, the number of measured displacements and the considered polynomial degree for approximating the behavior. Moreover, the use of a separated representation allows to diminish the number of experimental measurements because it involves optimal polynomial representations of the constitutive manifold.
<figure><img src="image_8.png" /><figcaption>Fig. 6 Comparing the exact and identified $C_{11}(\epsilon)$. Most of that data-point are almost superpose</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7 Dimensionality of the linear (**left**) and nonlinear (**right**) manifolds when applying the LLE nonlinear dimensionality reduction technique on the stress-strain data</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8 Constitutive manifold represented in a 3D space</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9 View of the manifold represented in Fig. 8 as a function of the two first embedding coordinates</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. 10** Lateral view of manifold represented in Fig. 8 that allows concluding on its almost perfect two-dimensional nature</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 13** Relative error related to the component C₁₁. Black points indicate the strains that were available for the identification procedure</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 11** Data-driven simulation based on the reference constitutive manifold. Displacement field: **u**<sub>x</sub> (left) and **u**<sub>y</sub> (right)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 12** Data-driven simulation based on the identified constitutive manifold. Displacement field: **u**<sub>x</sub> (left) and **u**<sub>y</sub> (right)</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Fig. 14** Identified constitutive coefficients $C^{11}(\epsilon)$ (left) versus the reference ones (right)</figcaption></figure>

**5 Conclusions**

We proved in a previous work [8] that numerical simulations can be performed from the only knowledge of data defining the material behavior. It was claimed that the main drawback of one such approach is the necessity of unveiling the whole constitutive manifold. However, at present, testing facilities are not able to explore the whole strain-stress space in a continuous way. In this paper we considered elastic behaviors (linear and nonlinear), proving that the constitutive manifold can be extracted from a data-driven inverse procedure in an effective manner.

Two procedures have been proposed: the first one is indeed a progressive construction of the constitutive manifold, while the second involves the polynomial approximation of the whole constitutive manifold. The first scheme results to be quite simple. However, the error accumulates all along the identification process. Moreover, the fact of using clustering techniques remains also a tricky issue. The second route, however, seems to be more robust from all points of view, and its immersion in a hierarchical or multi-resolution strategy seems an appealing choice for future developments.

Even if the results only concerned some simplistic behaviors, the methodology seems to be appropriate to address more complex scenarios, such as behaviors involving large strains, as well as inelastic deformations. Other points that should be considered are the ones related to existence and propagation of noise. This constitutes our current effort of research.

**Acknowledgements** This work has been supported by the ESI GROUP Chair at Ecole Centrale of Nantes as well as by the Spanish Ministry of Economy and Competitiveness, through grants number CICYT DPI2014-51844-C2-1-R and DPI2015-72365-EXP and by the Regional Government of Aragon and the European Social Fund, research group T88. This support is gratefully acknowledged.

**References**

1. Brunton SL, Proctor JL, Kutz JN (2016) Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proc Nat Acad Sci 113(15):3932-3937

2. Chinesta F, Leygue A, Bordeu F, Aguado JV, Cueto E, Gonzalez D, Alfaro I, Ammar A, Huerta A (2013) Parametric PGD based computational vademecum for efficient design, optimization and control. Arch Comput Methods Eng 20(1):31-59

3. Chinesta F, Keunings R, Leygue A (2014) The proper generalized decomposition for advanced numerical simulations: a primer. Springer, New York

4. Crespo J, Latorre M, Montans F (2017) WYPIWYG hyperelasticity for isotropic, compressible materials. Comput Mech 59(1):73-92

5. Gonzalez D, Ammar A, Chinesta F, Cueto E (2010) Recent advances in the use of separated representations. Int J Numer Methods Eng 81(5):637-659

6. Gonzalez D, Cueto E, Chinesta F (2015) Computational patient avatars for surgery planning. Ann Biomed Eng 44(1):35-45

7. Gonzalez D, Aguado JV, Cueto E, Abisset-Chavanne E, Chinesta F (2016) kPCA-based parametric solutions within the PGD framework. Arch Comput Methods Eng. doi:10.1007/s11831-016-9173-4

8. Ibanez R, Abisset-Chavanne E, Aguado JV, Gonzalez D, Cueto E, Chinesta F (2016) A manifold learning approach to data-driven computational elasticity and inelasticity. Arch Comput Methods Eng. doi:10.1007/s11831-016-9197-9

9. Kirchdoerfer T, Ortiz M (2016) Data-driven computational mechanics. Comput Methods Appl Mech Eng 304:81-101

10. Lee JA, Verleysen M (2007) Nonlinear dimensionality reduction. Springer, Berlin

11. Liu Z, Bessa MA, Liu WK (2016) Self-consistent clustering analysis: an efficient multi-scale scheme for inelastic heterogeneous materials. Comput Methods Appl Mech Eng 306:319-341

12. Lopez E, Gonzalez D, Aguado JV, Abisset-Chavanne E, Cueto E, Binetruy C, Chinesta F (2016) A manifold learning approach for integrated computational materials engineering. Arch Comput Methods Eng. doi:10.1007/s11831-016-9172-5

13. Michopoulos J, Farhat C, Houstis E (2004) Dynamic-data-driven real-time computational mechanics environment. In: Bubak M, van Albada GD, Sloot PMA, Dongarra J (eds) Computational science-ICCS 2004: 4th international conference, Krakow, Poland, June 6-9, 2004, proceedings, Part III. Springer, Berlin, pp 693-700

14. Olson GB (2000) Designing a new material world. Science 288(5468):993-998
15. Peherstorfer B, Willcox K (2015) Dynamic data-driven reduced order models. Comput Methods Appl Mech Eng 291:21–41

16. Peherstorfer B, Willcox K (2016) Data-driven operator inference for nonintrusive projection-based model reduction. Comput Meth- ods Appl Mech Eng 306:196–215

17. Polito M, Perona P (2001) Grouping and dimensionality reduction by locally linear embedding. In: Advances in neural information processing systems, vol 14. MIT Press, pp. 1255–1262

18. Raghupathi W, Raghupathi V (2014) Big data analytics in health- care: promise and potential. Health Inf Sci Syst 2(1):1–10

19. Roweis ST, Saul LK (2000) Nonlinear dimensionality reduction by locally linear embedding. Science 290(5500):2323–2326

20. Tenenbaum JB, de Silva V, Langford JC (2000) A global framework for nonlinear dimensionality reduction. Science 290:2319–2323

21. Wang Q (2012) Kernel principal component analysis and its applications in face recognition and active shape models. CoRR. arXiv:1207.3538