Héloïse Dandin

Adrien Leygue

Laurent Stainier

Graph-based representation of history-dependent material response in the Data-Driven Computational Mechanics framework

Keywords: Data-Driven Computational Mechanics, computational plasticity, graph theory i

Introduction

Model-free Data-Driven Computational Mechanics (DDCM), first introduced by [START_REF] Kirchdoerfer | Data-driven computational mechanics[END_REF], is a method for solving boundary value problems in solid mechanics without constitutive models. In this approach the material behaviour is represented by a discrete set of points, instead of a set of explicit mathematical expressions as in classical constitutive modelling (see, e.g., [START_REF] Lubliner | Plasticity Theory[END_REF]) or machine learning methods (see, e.g., [START_REF] Ghavamian | Accelerating multiscale finite element simulations of history-dependent materials using a recurrent neural network[END_REF], [START_REF] Flaschel | Discovering plasticity models without stress data[END_REF]). In DDCM, the solution is then formulated as a distance minimisation problem, which becomes a mixed-integer quadratic problem (MIQP) exhibiting combinatorial complexity [START_REF] Kanno | Mixed-integer programming formulation of a data-driven solver in computational elasticity[END_REF].

The data sets that describe the material behaviour, also referred to as material databases, can be obtained in different ways. On the one hand, multiscale modelling has been used to generate data with constitutive models identified at the microscopic scale [START_REF] Karapiperis | Data-Driven multiscale modeling in mechanics[END_REF][START_REF] Korzeniowski | A multi-level method for data-driven finite element computations[END_REF][START_REF] Platzer | From tailored data to adaptive solvers for multiscale simulations[END_REF]. On the other hand, identification methods such as Data-Driven Identification (DDI) allow estimating stresses from experimental strain measures without constitutive models [START_REF] Leygue | Non-parametric material state field extraction from full field measurements[END_REF][START_REF] Dalémat | Measuring stress field without constitutive equation[END_REF][START_REF] Vinel | Characterization of the thermomechanical behaviour of metals for high strain-rates, using ultra-high speed imaging cameras[END_REF]. In that case, DDI can be seen as an inverse problem of DDCM for material properties identification.

These material databases are likely to evolve during an incremental simulation, for example to locally enrich a sparse database [START_REF] Platzer | From tailored data to adaptive solvers for multiscale simulations[END_REF] or due to the irreversibility of the behaviour [START_REF] Eggersmann | Model-Free Data-Driven inelasticity[END_REF]. In the latter case, the solver must, at each point of the structure and for each loading increment, monitor irreversibility and propose a set of admissible data points. [START_REF] Eggersmann | Model-Free Data-Driven inelasticity[END_REF] review typical representational paradigms for history-dependent materials in continuum mechanics: general materials with memory, differential representation, internal and history variables. All these methods, especially the first one, imply to increase significantly the dimensionality of the problem. For instance, the work by [START_REF] Valdés-Alonzo | Identification of material properties and phase distribution of heterogeneous materials through data-driven computational methods : Towards an enhaced constitutive space[END_REF] develops an extension of DDI to dissipative behaviours within the framework of materials with memory, while [START_REF] Pham | Tangent space Data Driven framework for elasto-plastic material behaviors[END_REF], following [START_REF] Ciftci | Model-free data-driven simulation of inelastic materials using structured data sets, tangent space information and transition rules[END_REF], and [START_REF] Langlois | Non-parametric stress field estimation for history-dependent materials: Application to ductile material exhibiting Piobert-Lüders localization bands[END_REF] use differential representations. [START_REF] Poelstra | A data-driven framework for evolutionary problems in solid mechanics[END_REF] and Bartel et al. [2023] adopted a different strategy based on history surrogates that also requires a consequent dimensionality augmentation.

In this work, we develop a graph-based representation of history-dependent material databases. Graph theory is a branch of discrete mathematics interested in networks and pairwise relations between objects, that has already been applied in solid mechanics, for example to granular materials [START_REF] Goddard | From Granular Matter to Generalized Continuum[END_REF] and dislocation kinematics [START_REF] Starkey | Development of mean-field continuum dislocation kinematics with junction reactions using de Rham currents and graph theory[END_REF]. It offers a wide range of tools and methods designed to handle large amounts of data and solve optimisation problems. We show that the history dependence can be encoded in a material directed graph (abbr. digraph). This is first achieved by expanding the usual strain-stress constitutive space to other thermodynamical quantities like internal variables, and second by making use of graph arcs to record dissipation levels. Standard efficient graph algorithms can then be used to extract the material database from the material digraph and perform the simulation in strain-stress space.

To illustrate the concepts, we present numerical results obtained for one-dimensional elements, in which case all quantities of interest are scalar. We first put aside the combinatorial issues arising from the alternating minimisation solver by modelling the behaviour of a single element to demonstrate the ability of the method to generate a material database from the material digraph allowing us to solve a simple problem. For larger structures, we show the accuracy issues and challenges coming both from the combinatorial complexity of DDCM and the richer material behaviour. An enhanced solver that takes advantage of the thermodynamical information encoded in the material digraph is proposed. We show that it is possible to retrieve a good response from the system both at global and local scales.

Data-Driven Computational Mechanics framework for inelastic material response

This section recalls the mathematical framework for model-free data-driven inelasticity introduced by [START_REF] Kirchdoerfer | Data-driven computational mechanics[END_REF] and [START_REF] Eggersmann | Model-Free Data-Driven inelasticity[END_REF].

We consider the discrete representation of an inelastic body composed of N nodes and M integration points, which undergoes displacements u = {u a } N a=1 and loads f = {f a } N a=1 . In a time-discrete setting, we seek to approximate solutions at times {t 0 , . . . , t k , t k+1 , . . . }.

The internal state of the system can be described at each time step by strain and stress pairs {(ε e,k+1 , σ e,k+1 )} M e=1 , with ε e,k+1 , σ e,k+1 ∈ R me and m e the dimension of the tensors at integration point e for linearised kinematics. In a three-dimensional setting, m e = 6 in Voigt notation, while m e = 1 for bar elements. A local state is thus defined as z e,k+1 = (ε e,k+1 , σ e,k+1 ) ∈ Z e,k+1 with Z e,k+1 = R me×me the local phase space1 , while z k+1 = {z e } M e=1 is a point in the global phase space Z k+1 = M e=1 Z e,k+1 , with meaning the tensor product of spaces.

The mechanical admissibility of a state is material-independent and is given by the compatibility and equilibrium constraints at time step t k+1 :

ε e,k+1 = B e u k+1 , ∀e = 1 . . . M, ( 1a
) M e=1 w e B T e σ e,k+1 = f k+1 , ( 1b
)
with w e the weights of integration points and B e the discrete kinematic operator matrix related to the integration point e. These relations define the set of mechanical constraints,

E k+1 = {z ∈ Z | Eqs. (1a) and (1b)} ⊂ Z k+1 .
(2)

In the model-free data-driven approach the material response is described by a collection of states obtained through experiments or numerical simulations at a finer material scale. These local material data sets denoted

D e,k+1 = {y i = (ε i , σ i ) ∈ Z e,k+1 , i = 1 . . . N * | past local history} ⊂ Z e,k+1 , (3)
for inelastic behaviours, with N * the number of material strain-stress pairs, also extend to the global material database at time t k+1 as

D k+1 = D 1,k+1 × • • • × D M,k+1 ⊂ Z k+1 .
The following norm can be defined in the local phase space:

∥z e,k+1 ∥ C = 1 2 (C e,k+1 : ε e,k+1 ) : ε e,k+1 + 1 2 (C -1 e,k+1 : σ e,k+1 ) : σ e,k+1 1/2 , ( 4
)
with C a positive-definite 4th-order tensor that is not related to any material property and can be updated during the simulation. Hence, the (squared) distance between two states z e,k+1 , y e,k+1 ∈ Z e,k+1 is:

d 2 C (z e,k+1
, y e,k+1 ) = ∥z e,k+1 -y e,k+1 ∥ 2 C .

(5)

Finally, the extension to the global phase space Z k+1 gives the distance

d 2 C (z k+1 , y k+1 ) = M e=1 w e d 2 C (z e,k+1
, y e,k+1 ) .

(6)

The data-driven problem is then reformulated into a double minimisation problem, whose solution at step t k+1 is given by: S = arg min

z k+1 ∈E k+1 min y k+1 ∈D k+1 d 2 C (z k+1 , y k+1 ) , (7)
i.e. the pair (z k+1 , y k+1 ) of states, respectively mechanically admissible and from the material database constrained by history, which are closest to each other according to distance d C . The solution is typically obtained by alternating minimisation over continuous (z) and discrete (y) variables, as illustrated in Fig. 1, that consists in the fixed point iteration:

z (i+1) k+1 = P E k+1 P D k+1 z (i) k+1 , ( 8
)
with i the iteration number, y

(i) k+1 = P D k+1 z (i)
k+1 the closest point projection onto D k+1 and P E k+1 y (i) k+1 the projection of a material state onto E k+1 . The discrete nature of the material database leads to high combinatorial complexity as M and N * increase. The proportion of local minima also increases, implying a greater probability that the solver will converge on local minima that are far from the global minima [START_REF] Kanno | Mixed-integer programming formulation of a data-driven solver in computational elasticity[END_REF].

The main challenge thus lies in the selection of the local material databases D e,k+1 subject to the history of strain and stress as defined in Eq. 3.

history-dependent materials

In this section, we generalise the material database to the material digraph concept. We also show how to construct this digraph and extract the local database D e,k+1 from it.

Representation of the discrete material behaviour with a material digraph

The definitions and concepts of graph theory provided here can be found in [START_REF] Harary | Structural models: an introduction to the theory of directed graphs[END_REF], Bondy and Murty [2008] and [START_REF] Rigo | Advanced Graph Theory and Combinatorics[END_REF].

A directed graph G (abbr. digraph) consists of a pair (V (G), A(G)), or shortly, (V, A), with V a set of vertices and A a set of arcs disjoint from V . An incidence function ψ G links each arc to an ordered pair of vertices in V :

ψ G : A → (V, V ) a → (uv), (9)
with u the tail and v the head of a. A cost c(a) is assigned to each arc a ∈ A. G is connected if, for every partition of its vertex set into two non-empty subsets X and Y , there is an arc with one end in X and the other in Y ; otherwise, it is disconnected. We note G = (V, A) the digraph encoding the material behaviour, illustrated on Fig. 2. Its vertices are the material states and its arcs are thermomechanically consistent transitions between these states such that G contains no isolated vertex. Given an arc (y i y j ) ∈ A between vertices y i , y j ∈ V , the cost c D (y i y j ) is the dissipative cost of the transition from states y i to y j :

c D (y i y j ) = D i→j , ( 10
)
with D i→j the dissipation level of the transition. In this way, a non-dissipative or reversible transition is such that c D (y i y j ) = 0 and is encoded with two symmetric arcs in the digraph, i.e.

(y i y j ) and (y j y i ), while an irreversible transition is represented with one arc directed such that c D (y i y j ) > 0.

Figure 2: A digraph G for an elastoplastic material response with loading and unloading.

We define a non-dissipative directed subgraph E = (V (E), A(E)) ⊆ G such that A(E) is the set of non-dissipative arcs of A(G) and V (E) the set of vertices joined by arcs in A(E). E might be disconnected and all the vertices belonging to the same connected component can be linked together via zero-cost arcs only (see Fig. 3 (a)). As a result, each connected component is associated with an elastic domain, which fundamentally corresponds to a clique, i.e. a digraph in which every pair of vertices is joined by exactly two arcs, one in each direction. Furthermore, we define the dissipative directed subgraph P = (V (P ), A(P )) ⊆ G such that A(P ) = A(G) \ A(E) and V (P ) has no isolated vertex. A(P ) thus contains all thermodynamically irreversible transitions: all arcs in A(P ) have strictly positive dissipative cost, as in Fig. 3 (b).

For a database representing a purely elastic material behaviour, the digraph G is a clique and is equal to the non-dissipative subgraph E. Therefore, the local database D e,k+1 contains all material states.

From this point on, we will use the terms of computational mechanics and graph theory interchangeably to refer to a state of the material database or the corresponding vertex in the digraph and a thermomechanical transition or the equivalent arc. We also call material digraph and write G the digraph associated with the material data set.

Distinction between identical states in the constitutive space

Let us now consider two identical states in the constitutive space y k = (ε, σ) and y κ = (ε, σ) but with different histories of strain and stress, denoted {y l } l⩽k and {y λ } λ⩽κ respectively.

A path is a special type of graph made of an ordered sequence of arcs such that the tail of each arc of the sequence coincides with the head of the preceding and where all the visited vertices are pairwise distinct. A path with the same first and last vertices is called a cycle.

Histories {y l } l⩽k and {y λ } λ⩽κ exist in the digraph as paths whose vertices are the states of histories. Yet according to its definition, P is acyclic. Thus, if either {y l } l⩽k or {y λ } λ⩽κ or both contain any arcs a k ∈ A(P ) and a κ ∈ A(P ) that belong to different connected components of P , then y k and y κ are different vertices of G (see Fig. 4). Otherwise, y k and y κ belong to the same elastic domain and are strictly equal (same characteristics in the constitutive space and same loading histories). In conclusion, the digraph framework makes it possible to distinguish between states with the same localisation in the constitutive space but different past histories and thus to represent complete and repeated loading cycles with various loading directions without any modification in the formalism. In that case, simulating the material response is only achievable in the sampled areas of the constitutive space with adequate dissipation levels. For example, to predict the mechanical behaviour for a loading-unloading-reloading cycle, the material digraph must encode the material response for at least one cycle and for the same strain range as implied by the applied loading.

Local database selection

The solution of the data-driven problem at time step k + 1 and integration point e is given by the pair (z e,k+1 , y e,k+1 ) of respectively mechanical and material states computed with Eq. 7. This first requires knowing the local material database D e,k+1 , which contains all the admissible states in the data set given the history of e.

As shown in Section 3.1, the history of a state is encoded in the global material database digraph. Therefore, selecting D e,k+1 now comes down to searching for admissible futures in the material database or equivalently to searching for existing paths in the digraph.

To do that, we use the concept of tree as defined in graph theory. A tree is a connected graph that contains no cycles. As for digraphs, trees can be written in terms of coupled sets of arcs and vertices as T = (V (T ), A(T )). Starting from a root r ∈ V (G), a tree T (r) contains all vertices of G that can be reached from r, and the arcs leading to them. Hence a vertex v ∈ V (G) belongs to V (T ) if there is a path in G leading from r to v.

Knowing that digraph G encodes the thermodynamically admissible transitions between states of material database D and given a root y e,k , we build a rooted tree from this vertex and denote it T (y e,k ). Finally, the local material database at k + 1, illustrated in Fig. 5

Local database reduction

The alternating minimisation solver used for elastic problems is highly sensitive to local minima [START_REF] Kirchdoerfer | Data Driven Computing with noisy material data sets[END_REF], the number of which increases in inelasticity. In the incremental problem developed above, we usually consider small loading increments. We therefore assume small increments of the material solution at a material point e both in terms of distance in the constitutive space and dissipative cost of the transition.

To this end, we introduce the concept of path cost defined as the sum of every arc's cost:

c(p) = a∈A(p) c(a), (12)
with p a specific path and A(p) the set of its arcs. A minimum cost path is the shortest path with respect to cost c(•) and is denoted [u, v] c with u and v the tail and head vertices respectively. As part of an elastoplastic local database, u is always a tree's root and v is a vertex of T (u), which ensures that there is at least one path from u to v. Finally, three conditions on the states belonging to D e,k+1 as defined in Eq. 11 are applied to reduce local material databases:

1. squared data-driven distance:

d 2 C (y k , y i ) ⩽ TOL1, (13)
2. cumulated squared data-driven distance along the path:

c δ (p = [y k , y i ] δ ) = a∈A(p) c δ (a) ⩽ TOL2, (14)
with c δ (a) = d 2 C (tail(a), head(a)), 3. path dissipative cost:

c D (p = [y k , y i ] D ) = a∈A(p) c D (a) ⩽ TOL3, (15)
with c D (a) = D head(a)→tail(a) . It should be noted that the first two conditions involve algorithmic criteria while the latter is a material criterion. Furthermore, for any path from y k to y i , the data-driven distance is always smaller or equal to the distance along the path used in condition 2. The effect of these criteria and their combinations will be discussed in more detail in Section 4.3.

Fig. 6 details the new data-driven procedure for the incremental inelastic problem. As stated by [START_REF] Eggersmann | Model-Free Data-Driven inelasticity[END_REF], the difference with the elastic case lies in the definition of local databases while the following algorithmic steps remain the same.

Figure 6: DDCM algorithm for rate-independent elastoplastic material response at loading step k + 1.

Numerical representation of the material digraph

In practice, material data can be generated in two different ways. On the one hand, it can come from numerical simulations, for example in the context of multiscale modelling [START_REF] Karapiperis | Data-Driven multiscale modeling in mechanics[END_REF][START_REF] Gorgogianni | Adaptive goal-oriented data sampling in Data-Driven Computational Mechanics[END_REF]. In this case, a constitutive model provides all the information needed to build the material digraph, i.e. strains, stresses and dissipation. On the other hand, experimental tests can also be used to generate this data, using a full-field strain measurement technique coupled with DDI for strains and stresses. Dissipation can be inferred from calorimetric measurements with assumptions on potential heat exchanges [Chrysochoos et al., 2010[START_REF] Seghir | An improved lagrangian thermography procedure for the quantification of the temperature fields within polycrystals[END_REF][START_REF] Vinel | Characterization of the thermomechanical behaviour of metals for high strain-rates, using ultra-high speed imaging cameras[END_REF].

A graph is usually encoded by an adjacency matrix, i.e. a square matrice A that represents adjacent vertices of the graph. A component A uv usually equals the number of arcs from vertex u to v.

In the present work, the numerical representation of the material digraph G is an N * × N * adjacency matrix A ij = exp(c(y i y j )) where N * is the number of points in the data set and c(y i y j ) is the cost of arc (y i y j ). The exponential function allows to distinguish between the absence of arc from a zero-cost arc. We use two matrices to encode the costs defined above:

• data-driven distance matrix:

A δ ij = exp(c δ (y i y j )), • dissipative cost matrix: A D ij = exp(c D (y i y j )
). Both matrices represent the same material digraph (same vertices and same arcs) but with different arcs' costs. Therefore they have the same sparsity pattern that might be stored only once.

It is to be noted that the material digraph should be designed in such a way that it contains a sufficient number of arcs to provide enough information about the material behaviour and thus improve the solver's precision, and yet as few arcs as possible to ensure adjacency matrices' sparsity and speed up computations (e.g. graph search). In particular, cliques should not be fully encoded as such as long as at least one path exists between all vertices of the corresponding elastic domain. For instance, an enhanced representation of cliques could be developed to speed up computations.

Finally, local databases can be obtained thanks to generic graph search algorithms [Bondy and Murty, 2008] or shortest-path algorithms like Dijkstra's [START_REF] Dijkstra | A note on two problems in connexion with graphs[END_REF]. We implement the reduction criteria introduced in Section 3.4 within the chosen routine as a boundary or limit of path cost.

Numerical implementation and investigation of a single element problem

Material digraph construction

We generate an artificial database thanks to an elastoplastic material model with linear kinematic hardening. The yield function is given by:

f (σ, ε p ) = |σ -Hε p | -σ y (16)
with hardening modulus H, yield limit σ y and plastic strain ε p . We build a regular data set representing the material response under elastoplastic loading, elastic unloading and plastic reloading. The strain increment in the non-hardened elastic domain and the first dissipative part is 0.01 %. Unloading and reloading paths are generated with a plastic strain increment of 0.2 % and an elastic strain increment of 0.04 %.

We define the arcs of the material digraph as follows:

• elastic domains, or sets of states that could be represented by cliques, are encoded as minimum spanning subgraphs with respect to the data-driven distance, i.e. such that all vertices are connected in both directions and with arcs that minimise the total (data-driven) cost of each connected component, • elastic domains are linked together with dissipative arcs that encode irreversibility, i.e. such that

t 0 |ε p |dt > 0.
We rely on the order of the generated sequence of material points and on the dissipation level to determine the arcs' direction and whether they should be encoded. The final digraph is similar to the one presented in Fig. 2. According to this definition, superimposed states in the (ε, σ) configuration (see Fig. A local database contains potential states for a given loading increment. The criteria defined in Section 3.4 tighten eligibility conditions and therefore downsize the database. The data-driven distance criterion is implemented with sklearn's Nearest Neighbor algorithm and a maximum radius equal to TOL1. The path distance and cost conditions are computed with Dijkstra's algorithm implementation from Python library scipy.sparse.csgraph and graph adjacency matrices A δ and A D encoding respectively data-driven distance and transition dissipation as arcs' cost. Again, TOL2 and TOL3 represent the maximum admissible path cost for the last node of the sequence to be in the local database. Note that TOL represent relative tolerances taken with respect to the chosen metric C. The second criterion limits the data-driven distance along the shortest path from the root to a point in the local database (see Fig. 8 (c)). Note that the notion of path included in this cost excludes points corresponding plastic unloading, contrarily to previous criterion. It is equivalent to a radius in the graph representation and seems more suitable to our study as it better reflects the effective loading path leading to a given transition.

Finally, the dissipative criterion in Fig. 8 (d) adds all vertices in elastic domains obtained with dissipative increments lower than the chosen tolerance to the database. It is used in the following simulations in association with the second criterion. Other criteria could also be developed, for instance, to explicitly forbid non-monotonic increments (in dissipation). This assumption is for example used in the classical return-mapping algorithm for elastoplastic constitutive modelling.

Evaluation of results quality

DDCM solution strongly depends on database sparsity, especially when DDCM metric C is not optimised [Kirchdoerfer andOrtiz, 2017, Eggersmann et al., 2019]. This is the case here as we choose a fixed value for this parameter throughout the simulation despite slope discontinuities in the material data. This work aims to demonstrate the ability of the present graph-based approach to represent the mechanical solution independently of the convergence study, which should be improved by an enhanced solver, e.g. through game theory as suggested by [START_REF] Weinberg | Data-Driven Games in Computational Mechanics[END_REF]. To decouple the influence of these parameters, we seek to limit the data set's density effect by filtering the FE reference solution with the data set. We therefore introduce a new indicator: the data-driven projection of the FE reference solution onto the material database. This value is computed in two steps:

1. projecting the FE solution z ref onto the material database (reference material states),

y ref = P D z ref ,
2. re-projecting the reference material states onto equilibrium (projected reference mechanical states),

z proj = P E k+1 y ref .
The projected states represent the best solution achievable by a data-driven solver with given data and chosen metric. This is particularly important as the material data set might contain a non-uniform discretisation of the constitutive space with possibly very different sampling rates in the elastic and dissipative domains. Even for a regular data set as in Fig. 9, gaps appear in the projected data when the database is too sparse in the neighbourhood of the FE solution.

Hence the distance between the DDCM solution and the projected reference solution provides a neutral indicator regarding data sparsity and is used in the following to evaluate the results accuracy.

It is noteworthy that the distance between DDCM solution and projected reference states is not totally independent of the database density, since the convergence of the fixed point algorithm depends on this density. However, this indicator is more neutral than the comparison with the FE reference solution usually used.

Resolution of a data-driven problem

We illustrate the approach with the one-dimensional spring-bar element system on Fig. 10 subject to a strain ε up to 15 % followed by elastic unloading and plastic reloading in 11 (a) shows the mechanical and material states obtained at all time steps. The DDCM solver is able to retrieve a good approximation of the path shape even during the unloading phase, where the material database does not contain points with the exact reference dissipation level. The projected reference solution is well recovered except from the onset of unloading (time step 188) until reloading (time step 225) where the distance between DDCM mechanical states and projected reference mechanical states increases (see Fig. 11 (b) and(c)). As stated by [START_REF] Kanno | Mixed-integer programming formulation of a data-driven solver in computational elasticity[END_REF], this phenomenon is characteristic of the solver and occurs even for simple problems. Moreover, absolute strain and stress errors, computed as

err ε = |ε k -ε proj k | , err σ = |σ k -σ proj k | . ( 17
)
are given on Fig. 11 (c) between mechanical and projected reference mechanical states. The evolution of these values during loading is coherent with the chosen metric which affects a lower weight on stresses compared to strains and thus tends to approximate the latter better. 5 Numerical experiment on a truss problem

Truss structure and related challenges

The simulation of a truss structure is a combinatorial problem that involves coupled minimisations of mechanical and material states at all integration points. The DDCM alternating minimisation solver ensures global convergence but is likely to fall into local minima, which could lead to a completely unreliable solution both at local and global scales. As an example, Fig. 12 (a) shows a 2D truss composed of 102 nodes and 252 elements, subject to displacement-driven boundary conditions. The top and left boundaries are fixed while right and bottom nodes are subject to the same displacement up to 7 mm and back to 6 mm with 135 increments, as illustrated in Fig. 12 (b). The material data set is artificially generated as described on Section 4.1 with constitutive model parameters E = 217.5 GPa, H = 1 GPa, σ y = 250 MPa, a strain increment of 0.001 % and the metric is set to C = E. These parameters are representative of steel and close (in adimensional form) to those used in the previous section. The local database reduction tolerances are set to TOL2 = 5 × 10 3 × C and TOL3 = 10 -5 × C. DDCM mechanical and material states obtained at all loading steps for two plastically deformed elements of the mesh are illustrated on Fig. 12 (c). Although prediction is satisfactory for the yellow element, the solution computed for the red element does not match the expected elastoplastic response. Strain "sliding" is observed for strains from 0 Pa 1/2 to 6 × 10 3 Pa 1/2 and DDCM underestimates maximum strain, even though the overall path shape is consistent. These effects are related to the constitutive space discretisation with finer strain increments along with a metric value that gives strains a higher weight in the distance calculation. Besides, Fig. 13 (a) highlights the median of relative strain and stress errors, calculated as:

err ε = |ε e,k -ε proj e,k |/ε proj e,k , err σ = |σ e,k -σ proj e,k |/σ proj e,k . (18)
The difference in magnitude between the strain and stress medians is consistent with the strain and stress discretisation: with the data-driven distance, which transforms into the euclidean distance in the (ε √ E, σ/ √ E) space, the stress range is much lower than that of strains, and stress discretisation is denser.

Fig. 13 (b) shows how relative errors are distributed in the mesh with elements sorted by increasing strain. The plotted time step corresponds to the loading phase, which implies that both strain and stress should grow together. For some elements labelled in range 220-252 nonetheless, stress levels are lower than expected, which corresponds to the stress underestimation observed on Fig. 12 (c) for the red element. At this specific loading step, 19 % of the elements have strain relative errors greater than 100 %, half of which reaching strain levels smaller than 500 Pa 1/2 ≈ 0.1 %. Additionally, combinatorial resolution optimises the global objective function to the detriment of local behaviours. The macroscopic response, illustrated on Fig. 13 (c) by the evolution of the resulting force along the truss right boundary as a function of the applied displacement, also exhibits the consequence of local sliding for lower strains. The DDCM maximum force however is close to reference as local over-and underestimations of strains offset each other.

Even though the laws of mechanics and thermodynamics are fulfilled, respectively with the mechanical fields satisfying kinematic compatibility and equilibrium, and the material fields satisfying the condition of non-negativity of dissipation, the lack of a unique solution within local material databases combined with the sensitivity of the alternated minimisation solver to local minima leads to defective results at local and global scales. We therefore discuss a possible improvement of the solver for combinatorial elastoplastic problems.

Predictor-corrector algorithm

We suggest to implement a two-step non-dissipative predictor-dissipative corrector algorithm inspired by resolution methods for constitutive models. This first requires to define two types of local databases:

• a non-dissipative prediction local database D pred e,k+1 ⊂ D e,k+1 that only contains states such that the dissipative cost of any path from subgraph root y e,k to a state of the local database is null, or equivalently, local databases only contain states belonging to the same elastic domain as the roots, • a dissipative correction local database D corr e,k+1 = D e,k+1 obtained as defined in Section 3.3 by taking the prediction's material state y pred e,k+1 as root. Local databases D pred e,k+1 and D corr e,k+1 might also respect the reduction criteria based on cumulated data-driven distance along path and path dissipative cost, introduced in Section 3.4. In practice, D pred e,k+1 is calculated by setting the tolerance on dissipative path cost TOL3 to zero and the tolerance on path data-driven distance TOL2 for D corr e,k+1 must be chosen such that D pred e,k+1 ⊆ D corr e,k+1

to allow permanent deformation of elements. A data-driven iteration is thus performed in two steps as illustrated on Fig. 14. First, a non-dissipative step, called predictor, where the resolution is performed using D pred e,k+1 with root y e,k until convergence. Then a dissipative computation or correction is performed using D corr e,k+1 , equal to the local database defined in Section 3.3 and Section 3.4, obtained by taking the predicted material state y pred e,k+1 as root vertex. In addition, the material states in the correction stage are initialised with the predicted material states. The corrected states are also the final solution of the loading increment. This modification of the DDCM algorithm does not involve any additional hypothesis or parameter but is solely based on physics.

The truss problem introduced above is now solved with the enhanced algorithm. The DDCM solution for the yellow element depicted in Fig. 12 solution improved significantly compared to Fig. 12. The material response is well predicted in the reversible part, although the solver still tends to underestimate strains in the dissipative domain, leading to early unloading.

Fig. 16 confirms that relative errors also improved compared to the previous simulation. The median values were reduced down to 20 % of those in the previous simulation. Additionally, at loading step 33, high relative errors mostly restrict to slightly deformed elements: only 11 % of the elements exhibit strain relative errors greater than 100 %, half of which with strain levels smaller than 500 Pa 1/2 ≈ 0.1 %. Finally, the macroscopic response of the mesh, evaluated by the resulting force against displacement curve on Fig. 16 (c), is recovered as well.

Conclusion

This work extends the concept of material database defined in the original DDCM approach to that of material digraph, suitable for representing history-dependent behaviours. The former is merely a sampling of the strain-stress response of the material, and therefore sufficient to represent elastic behaviours. The latter also encodes thermomechanically admissible transitions between states through the arcs of the digraph. In an incremental scheme, a local database representing the local material behaviour has to be selected at each loading step. We show that this selection operation can be efficiently performed with graph tools by computing a rooted tree from the last local material state. The obtained local database represents every possible future states that are compatible with the thermomechanical history. It therefore contains states that might be distant regarding different criteria. A reduction of the local database is proposed by filtering states according to three of those criteria. Additional conditions could also be introduced to filter states obtained from dissipation-wise non-monotonic paths. This restriction, similar to the hypothesis used by the classical return-mapping algorithm in constitutive modelling, would ensure that the local database only contains states belonging either to the same elastic domain as the root or to a linked dissipative part of the material database.

The presented numerical results were obtained for material databases computed with an elastoplastic constitutive model with linear kinematic hardening. However, neither the material digraph nor the local databases used by the DDCM solver rely on any internal variable. We thus believe that the method can handle any hardening type and possibly other inelastic behaviours, although the latter may require adapting the constitutive space. Cyclic loading can be predicted as long as the constitutive space discretisation (i.e. the database) covers the range of strain and dissipation that needs to be studied. For trusses, the standard DDCM solver is strongly affected by the problem combinatorial complexity and can reach a spurious solution. Model-based computational mechanics usually circumvents such solutions by limiting the material behaviour to a tangent space and controlling loading increments. Within the DDCM framework, we take advantage of the material digraph to compute two local databases: the thermomechanically admissible database, as defined above, and a smaller iso-dissipation database. The first one is called dissipative database, while the latter is non-dissipative. An enhanced two-stage solver, composed of a prediction step with the non-dissipative local database and a correction step with the dissipative database, has been developed that exhibits promising results.

Fundamentally, there is a priori no theoretical obstacle to the extension of the method to full 3D elastoplasticity. However, the extension implies an increase of the constitutive space dimensionality and a larger number of points is required to accurately represent the material behaviour. This problem is not specific to inelastic behaviours since it even arises in elasticity, and has been addressed, e.g. with adaptive approaches [START_REF] Gorgogianni | Adaptive goal-oriented data sampling in Data-Driven Computational Mechanics[END_REF]. Additionally, the amount of possible dissipative paths will increase even faster and the material digraph will necessarily have to sample the possible transitions. An efficient description of the material digraph then becomes necessary, and graph theory provides a wide range of tools and algorithms that have not yet been explored in detail. In particular, a hierarchical representation of the digraph based on cliques, that play an important role in the approach, could be of use.

Nevertheless, once this operation is achieved, graph theory provides optimised algorithms for efficient graph search and database construction. For instance regarding local databases selection, an operation that has to be repeated at every integration point, a worst case scenario gives Dijkstra's algorithm (with Fibonacci heap could be used to explore the vertices from one depth to the next, starting from the root (as opposed to Depth-First Search algorithm, which explores the complete branch before moving to the next). These numerical optimisations will be the subject of future work.

Figure 1 :3

 1 Figure 1: DDCM algorithm for elastic material response at loading step k + 1

Figure 3 :

 3 Figure 3: Subgraphs of G (Fig. 2): (a) non-dissipative subgraph E, (b) dissipative subgraph P .

Figure 4 :

 4 Figure 4: Identical states y k and y κ in constitutive space (ε, σ) but with different histories.

  Figure 5: (a) Tree T (y e,k ) built from G (Fig 2). (b) Material database, root y e,k and local database D e,k+1 .

  7 (a)) might have different histories and dissipation levels (see Fig. 7 (b)) and are different vertices of the material digraph.

Figure 7 :

 7 Figure 7: Material data set under loading, unloading and reloading conditions, the black solid line represents an example history: strain and stress reach a maximum before decreasing, while dissipation only increases. (a) In stress-strain space, colouring represents the dissipation state, (b) 3d representation.

Fig. 8

 8 Fig. 8 shows local databases for different reduction criteria in space (ε √ E, σ/ √ E), which has the same dimensions as the constitutive space. In this representation and because we set the DDCM parameter C equal to Young's modulus, data-driven and euclidean distance are equivalent. Criterion 1 expresses a maximum data-driven radius centred on the root and is represented by a circle in these coordinates (see Fig. 8 (b)).

Figure 8 :

 8 Figure 8: Material database and effect of reduction criteria on the local database. (a) Material database. (b) Data-driven distance, TOL1=10 -5 × C. (c) Cumulated data-driven distance along the path, TOL2=10 -2 × C. (d) Path dissipative cost, TOL3=5 • 10 -5 × C.

Figure 9 :

 9 Figure 9: Material data set and (left) FE solution, (right) projected reference mechanical states.

Figure 10 :

 10 Figure 10: Boundary value problem: spring-bar element system

Figure 11 :

 11 Figure 11: (a) Mechanical and material solutions for the spring-bar element model under loading, unloading and reloading. (b) local material database, DDCM solution (blue star and red point) and projected reference mechanical state (green point) at maximum loading. (c) Errors between DDCM and projected reference mechanical states.

Figure 12 :

 12 Figure 12: (a) Truss geometry and boundary conditions. (b) Applied displacement. (c) Evolution of DDCM mechanical, material states and projected reference mechanical states for yellow and red elements during loading.

Figure 13 :

 13 Figure 13: (a) Median relative errors between DDCM and projected reference mechanical states. (b) Relative errors at time step 33 (black dotted line on plot (a)) between DDCM and projected reference mechanical states, the yellow and red dotted lines refer to the elements highlighted on Fig. 12 (a). (c) Evolution of resulting force against displacement for the truss problem.

  Figure 14: Illustration of the predictor-corrector elastoplastic DDCM algorithm.

]Figure 15 :

 15 Figure 15: Evolution of DDCM mechanical, material states and projected reference mechanical states for the truss problem with predictor-corrector algorithm. (a) Yellow element. (b) Red element.

Figure 16

 16 Figure 16: (a) Median relative errors between DDCM and projected reference mechanical states. (b) Relative errors at time step 33 (black dotted line on plot (a)) between DDCM and projected reference mechanical states, the yellow and red dotted lines refer to the elements highlighted on Fig. 12 (a). (c) Evolution of resulting force against displacement for the truss problem with predictor-corrector algorithm.

  ) running time complexity in O(|A| + |V | × log(|V |)), with | • | the cardinal number of the sets of arcs A and vertices V . The reduction criteria yet restrict the graph search to the root's neighbourhood, decreasing time cost. Alternatively, a Breadth-First Search (BFS) algorithm with a complexity in O(|A| + |V |)

The notation Z e,k+1 is introduced here for convenience as it will be used in the following. However, note that the local phase space is independent of k.

Acknowledgments

All authors gratefully acknowledge the financial support of the Deutsche Forschungsgemeinschaft (DFG) and French Agence Nationale de la Recherche (ANR) through the project "Direct Data-Driven Computational Mechanics for Anelastic Material Behaviours" (ANR-19-CE46-0012, RE 1057/47-1, project number 431386925) within the French-German Collaboration for Joint Projects in Natural, Life and Engineering (NLE) Sciences.

Bibliography

T. Bartel, M. Harnisch, B. Schweizer, and A. Menzel