Impact Statement
While phenomenological constitutive models have been extensively employed in solving solid mechanics problems, formulation and characterization of models become difficult and often lose generality for complex material systems. Data-driven modeling that directly integrates material data with physical laws into simulations has shown promising potential in recent years. This study aims to develop a data-driven modeling approach for nonlinear anisotropic materials by introducing a two-level data search scheme that utilizes both stress–strain data and orientations of material anisotropy to discover the underlying anisotropic material manifold that represents the anisotropic material behaviors. This data search scheme is implemented with a local convexity-preserving reconstruction algorithm for enhanced robustness against noise and outliers in data-driven computing. Ultimately, the direct integration of discrete material data into simulations avoids the need for pre-defining the stress–strain relationship. The proposed data-driven modeling approach has demonstrated satisfactory accuracy in predicting nonlinear anisotropic material behaviors compared with constitutive model-based solutions, showing promising potential in modeling nonlinear anisotropic material systems.
1. Introduction
Constitutive laws that explicitly describe a stress–strain relationship are one of the key components in computational mechanics. However, the phenomenological constitutive models for material modeling are typically based on empirical observations, hypothesis, and mathematical assumptions, often lack of generality to capture complex material behaviors, such as, anisotropy and nonlinearity. For example, it remains challenging to formulate conventional phenomenological constitutive laws for musculoskeletal systems consisting of multiple highly nonlinear and anisotropic biological materials (Huijing, Reference Huijing1999; Gillies and Lieber, Reference Gillies and Lieber2011; Takaza et al., Reference Takaza, Moerman, Gindre, Lyons and Simms2013; Zhang et al., Reference Zhang, Chen, He, He, Basava, Hodgson, Sinha and Sinha2020).
With advancements in computing power and significant progresses in data mining, machine learning based data-driven approaches, for example, manifold learning and artificial neural networks, have been extensively applied to various fields owing to their strong ability of feature/pattern extraction, including constitutive material modeling (Ghaboussi et al., Reference Ghaboussi, Garrett and Wu1991; Furukawa and Yagawa, Reference Furukawa and Yagawa1998; Lefik and Schrefler, Reference Lefik and Schrefler2003; Wang and Sun, Reference Wang and Sun2018; Liu et al., Reference Liu, Wu and Koishi2019; Mozaffar et al., Reference Mozaffar, Bostanabad, Chen, Ehmann, Cao and Bessa2019; Heider et al., Reference Heider, Wang and Sun2020; He and Chen, Reference He and Chen2020; He et al., Reference He, Laurence, Lee and Chen2020b), surrogate and reduced order modeling (Lee and Chen, Reference Lee and Chen2013; Chen et al., Reference Chen, Marodon and Hu2015; He et al., Reference He, Chen and Marodon2019; Lee and Carlberg, Reference Lee and Carlberg2020; Mendizabal et al., Reference Mendizabal, Márquez-Neila and Cotin2020; Rocha et al., Reference Rocha, Kerfriden and van der Meer2020; Maulik et al., Reference Maulik, Lusch and Balaprakash2020), solving partial differential equations and system identification (Raissi et al., Reference Raissi, Perdikaris and Karniadakis2019; He et al., Reference He, Barajas-Solano, Tartakovsky and Tartakovsky2020a; Goswami et al., Reference Goswami, Anitescu, Chakraborty and Rabczuk2020), and structural design (Zhang and Ye, Reference Zhang and Ye2019; Lei et al., Reference Lei, Liu, Du, Zhang and Guo2019). Recently, a data-driven computing paradigm has been developed to bypass the material modeling step by formulating boundary value problems as a new optimization problem to search for a solution that is closest to the material dataset and constrained by equilibrium and compatibility (Kirchdoerfer and Ortiz, Reference Kirchdoerfer and Ortiz2016; Ibanez et al., Reference Ibanez, Abisset-Chavanne, Aguado, Gonzalez, Cueto and Chinesta2018; Conti et al., Reference Conti, Müller and Ortiz2018). Recent developments have extended the physics-constrained data-driven computing framework to dynamics (Kirchdoerfer and Ortiz, Reference Kirchdoerfer and Ortiz2018), nonlinear elasticity (Nguyen and Keip, Reference Nguyen and Keip2018; He et al., Reference He, Laurence, Lee and Chen2020b), inelasticity (Eggersmann et al., Reference Eggersmann, Kirchdoerfer, Reese, Stainier and Ortiz2019), constitutive manifold construction (Ibañez et al., Reference Ibañez, Borzacchiello, Aguado, Abisset-Chavanne, Cueto, Ladevèze and Chinesta2017; Stainier et al., Reference Stainier, Leygue and Ortiz2019), and multiscale modeling (Xu et al., Reference Xu, Yang, Yan, Huang, Giunta, Belouettar, Zahrouni, Zineb and Hu2020; Mora-Macías et al., Reference Mora-Macías, Ayensa-Jiménez, Reina-Romo, Doweidar, Domínguez, Doblaré and Sanz-Herrera2020).
To counteract the curse of dimensionality and instability due to noise and outliers present in material datasets, He and Chen (Reference He and Chen2020) introduced locally linear embedding manifold learning and convexity-preserving reconstruction into the data-driven solver, called local convexity data-driven (LCDD) computing. The LCDD framework has been extended to model nonlinear elastic solids and applied to model mechanical responses of biological heart valve tissues with experimental datasets by He et al. (Reference He, Laurence, Lee and Chen2020b). Recently, Eggersmann et al. (Reference Eggersmann, Stainier, Ortiz and Reese2020) introduced a tensor voting machine learning technique into the entropy based distance-minimizing data-driven (DMDD) framework to construct locally linear tangent spaces in order to utilize underlying data structures to achieve high-order convergence of data-driven solutions, referred to as a second-order data-driven scheme. The aforementioned physics-constrained data-driven computing frameworks have demonstrated promising performances in various scientific fields. Despite these advances, these data-driven computing frameworks cannot effectively handle anisotropic material systems with various anisotropic orientations, that is, orientations of material anisotropy. This is due to the fact that the data-driven solvers of these frameworks do not consider information of material orientations, which prevents the applications of these data-driven computing frameworks to complex material systems with direction-dependent material behaviors, for example, musculoskeletal systems with significant differences in the muscle fiber direction leading to material anisotropy.
The objective of this study is to develop a new data-driven solver built upon the local convexity-preserving reconstruction scheme by He and Chen (Reference He and Chen2020) in order to model anisotropic nonlinear elastic solids. The remainder of this paper is organized as follows. Section 2 reviews basic formulations of physics-constrained data-driven modeling of nonlinear elastic solids and makes a comparison of the material solvers of the DMDD by Kirchdoerfer and Ortiz (Reference Kirchdoerfer and Ortiz2016) and LCDD by He and Chen (Reference He and Chen2020). A local convexity-preserving material solver designed for anisotropic solids is introduced. Particularly, a rotated material database is constructed offline and a two-level data search is integrated into the material solver to capture direction-dependent material behaviors during online data-driven computing. In Section 3, the proposed data-driven computing framework is verified by modeling the deflection of a cantilever beam with layers containing different fiber directions and inflation of a cylinder where fiber directions vary along the circumferential direction of the cylinder. The data-driven solutions are compared with the constitutive model-based reference solutions to examine the effectiveness and robustness of the proposed methods. Concluding remarks and discussions are given in Section 4.
2. Methodologies
This section revisits the basic formulations of the physics-constrained data-driven computing framework for nonlinear elastic solids. The distinction between the material solvers of the DMDD by Kirchdoerfer and Ortiz (Reference Kirchdoerfer and Ortiz2016) and LCDD by He and Chen (Reference He and Chen2020) and He et al. (Reference He, Laurence, Lee and Chen2020b) is discussed. To model anisotropic solids, we propose to integrate a two-level data search scheme into the local convexity-preserving material solver in LCDD in order to capture the anisotropic material properties with anisotropic orientations (orientations of material anisotropy) information.
2.1. Physics-constrained data-driven modeling of nonlinear elastic solids
We start by summarizing the physics-constrained data-driven computing formulation for nonlinear elastic solids. The deformation of a nonlinear elastic solid defined in a domain 
$ {\Omega}^X $ with a Neumann boundary 
$ {\Gamma}_t^X $ and a Dirichlet boundary 
$ {\Gamma}_u^X $ are described by two basic laws:
$$ \mathrm{Equilibrium}:\left\{\begin{array}{ll} DIV\left(\mathbf{F}\left(\mathbf{u}\right)\cdot \mathbf{S}\right)+\mathbf{b}=\mathbf{0},& \mathrm{in}\hskip1em {\Omega}^X,\\ {}\left(\mathbf{F}\left(\mathbf{u}\right)\cdot \mathbf{S}\right)\cdot \mathbf{N}=\mathbf{t},& \mathrm{on}\hskip1em {\Gamma}_t^X,\end{array}\right. $$
$$ \mathrm{Compatibility}:\left\{\begin{array}{ll}\mathbf{E}=\mathbf{E}\left(\mathbf{u}\right)=\left({\mathbf{F}}^T\mathbf{F}-\mathbf{I}\right)/2,& \mathrm{in}\hskip1em {\Omega}^X,\\ {}\mathbf{u}=\mathbf{g},& \mathrm{on}\hskip1em {\Gamma}_u^X,\end{array}\right. $$where the superscript “
$ X $” denotes the reference (undeformed) configuration, 
$ \mathbf{u} $ is the displacement vector, 
$ \mathbf{E} $ is the Green-Lagrangian strain tensor, and 
$ \mathbf{F} $ is the deformation gradient related to 
$ \mathbf{u} $, defined as 
$ \mathbf{F}\left(\mathbf{u}\right)=\partial \left(\mathbf{X}+\mathbf{u}\right)/\partial \mathbf{X} $, where 
$ \mathbf{X} $ is the material coordinate, 
$ \mathbf{S} $ is the second Piola–Kirchhoff (second-PK) stress tensor, and 
$ \mathbf{b} $, 
$ \mathbf{N} $, 
$ \mathbf{t} $, and 
$ \mathbf{g} $ are the body force, the surface normal on 
$ {\Gamma}_t^X $, the traction on 
$ {\Gamma}_t^X $, and the prescribed displacement on 
$ {\Gamma}_u^X $, respectively. Let 
$ \mathcal{Z} $ denote the phase space that contains strain–stress pairs 
$ \mathbf{z}=\left(\mathbf{E},\mathbf{S}\right) $. The physical laws, Equations (1) and (2), define a material independent subset:
$$ \mathcal{C}=\left\{\mathbf{z}\hskip.3em \in \hskip.3em \mathcal{Z}:(1)\;\mathrm{and}\;(2)\right\} $$which contains all states satisfying the physical laws and they are called physical state.
Consider that the behaviors of the elastic solid are described by a discrete material dataset 
$ \unicode{x1D53C}={\left\{({\hat{\mathbf{E}}}_j,{\hat{\mathbf{S}}}_j)\right\}}_{j=1}^M $, where 
$ M $ is the number of measured material data and a hat symbol “
$ \wedge $” is used to denote material quantities. The strain–stress pair 
$ \hat{\mathbf{z}}=(\hat{\mathbf{E}},\hat{\mathbf{S}})\hskip.3em \in \hskip.3em \unicode{x1D53C} $ is called the material data. Conventionally, the material dataset 
$ \unicode{x1D53C} $ is used to characterize material constants associated with a predefined stress–strain relation, 
$ \mathbf{S}=\mathbf{f}\left(\mathbf{E}\right) $ in 
$ {\Omega}^X $. The constructed material model with properly estimated model coefficients is then combined with the physical laws (Equations 1 and 2) to solve the boundary value problem. However, as material models 
$ \mathbf{f} $ are usually based on empirical experimental observations, physical principles, and mathematical assumptions, with parameters calibrated from experimental data of material samples (Ghaboussi et al., Reference Ghaboussi, Garrett and Wu1991; Sussman and Bathe, Reference Sussman and Bathe2009), the phenomenological material modeling process inevitably lacks generality, especially for complex material systems (Latorre and Montáns, Reference Latorre and Montáns2014; Ibanez et al., Reference Ibanez, Abisset-Chavanne, Aguado, Gonzalez, Cueto and Chinesta2018).
The recently developed physics-constrained data-driven computing paradigm takes an alternative path to solve boundary value problems, which directly embeds material data into physical simulations, bypassing the need of phenomenological model construction (Kirchdoerfer and Ortiz, Reference Kirchdoerfer and Ortiz2016; Ibanez et al., Reference Ibanez, Abisset-Chavanne, Aguado, Gonzalez, Cueto and Chinesta2018; He and Chen, Reference He and Chen2020). One class of data-driven computing framework consists of minimizing the distance between the material data 
$ \hat{\mathbf{z}} $ and the strain–stress state 
$ \mathbf{z} $ satisfying the physical laws in Equations (1) and (2), which is called distance-minimizing data-driven (DMDD) computing (Kirchdoerfer and Ortiz, Reference Kirchdoerfer and Ortiz2016). A global distance functional is defined to measure the distance between the material data 
$ \hat{\mathbf{z}} $ and the physical state 
$ \mathbf{z} $
$$ d\left(\mathbf{z},\hat{\mathbf{z}}\right)={\int}_{\Omega^X}{d}_z^2\left(\mathbf{z},\hat{\mathbf{z}}\right)d\Omega $$with
$$ {d}_z^2\left(\mathbf{z},\hat{\mathbf{z}}\right)={d}_E^2(\mathbf{E}\left(\mathbf{u}\right),\hat{\mathbf{E}})+{d}_S^2(\mathbf{S},\hat{\mathbf{S}}), $$
$$ {d}_E^2(\mathbf{E}\left(\mathbf{u}\right),\hat{\mathbf{E}})=\frac{1}{2}(\mathbf{E}\left(\mathbf{u}\right)-\hat{\mathbf{E}}):\mathbf{M}:(\mathbf{E}\left(\mathbf{u}\right)-\hat{\mathbf{E}}), $$
$$ {d}_S^2(\mathbf{S},\hat{\mathbf{S}})=\frac{1}{2}(\mathbf{S}-\hat{\mathbf{S}}):{\mathbf{M}}^{-1}:(\mathbf{S}-\hat{\mathbf{S}}), $$where 
$ \mathbf{M} $ is a predefined symmetric and positive-definite weight tensor used to regularize the distances between 
$ \mathbf{z} $ and 
$ \hat{\mathbf{z}} $. For a system with multiple materials, multiple material datasets are required for data-driven computing, with each dataset describing material behaviors of one material. An ensemble of these material datasets is called a material database, 
$ {\unicode{x1D53C}}_{em}={\unicode{x1D53C}}_1\times {\unicode{x1D53C}}_2\times \dots \times {\unicode{x1D53C}}_m\hskip.3em \subset \hskip.3em \mathcal{Z} $, where 
$ m $ denotes the number of material types, 
$ {\unicode{x1D53C}}_p={\left\{({\hat{\mathbf{E}}}_j^p,{\hat{\mathbf{S}}}_j^p)\right\}}_{j=1}^{M_p} $, and 
$ {M}_p $ is the number of material data in 
$ {\unicode{x1D53C}}_p $, 
$ p=1,\dots, m $.
2.2. Data-driven two-step processes: material and physical solvers
The data-driven computing problem is formulated by stating:
$$ \underset{\hat{\mathbf{z}}\hskip.3em \in \hskip.3em {\unicode{x1D53C}}_{em}}{\min}\hskip1em \underset{\mathbf{z}\hskip.3em \in \hskip.3em \mathcal{C}}{\min}\hskip1em d\left(\mathbf{z},\hat{\mathbf{z}}\right)=\underset{\mathbf{z}\hskip.3em \in \hskip.3em \mathcal{C}}{\min}\hskip1em \underset{\hat{\mathbf{z}}\hskip.3em \in \hskip.3em {\unicode{x1D53C}}_{em}}{\min}\hskip1em d\left(\mathbf{z},\hat{\mathbf{z}}\right), $$where 
$ d\left(\mathbf{z},\hat{\mathbf{z}}\right) $ is the global distance functional defined by Equation (4). The objective is to find the material data 
$ \hat{\mathbf{z}}\hskip.3em \in \hskip.3em {\unicode{x1D53C}}_{em} $ in the material database that is closest to the constraint set 
$ \mathcal{C} $ of physical states 
$ \mathbf{z} $ or equivalently to find the physical state 
$ \mathbf{z}\hskip0.80em \in \hskip0.80em \mathcal{C} $ that is closest to the material database 
$ {\unicode{x1D53C}}_{em} $. The data-driven problem in Equation (8) can be decomposed into a two-step problem:
$$ \mathrm{Physical}\hskip0.17em \mathrm{step}:{\mathbf{z}}^{\ast }=\underset{\mathbf{z}\hskip.3em \in \hskip.3em \mathcal{C}}{\arg\;\min}\;d\left(\mathbf{z},\hat{\mathbf{z}}\right), $$
$$ \mathrm{Material}\hskip0.17em \mathrm{step}:{\hat{\mathbf{z}}}^{\ast }=\underset{\hat{\mathbf{z}}\hskip.3em \in \hskip.3em {\unicode{x1D53C}}_{em}}{\arg\;\min}\;d\left({\mathbf{z}}^{\ast },\hat{\mathbf{z}}\right), $$where 
$ {\hat{\mathbf{z}}}^{\ast } $ is the optimal material data sought from the material step that is closest to the physical state 
$ {\mathbf{z}}^{\ast } $ computed from the physical step.
2.2.1. Physical solver
A physical solver is used to solve the physical step defined in Equation (9), which searches for the physical state 
$ \mathbf{z}\hskip.3em \in \hskip.3em \mathcal{C} $ that is closest to a given material data 
$ \hat{\mathbf{z}} $. It can be reformulated as a constrained minimization problem:
$$ {\displaystyle \begin{array}{ll}\underset{\mathbf{z}\hskip.3em \in \hskip.17em \mathcal{C}}{\min}\;d\left(\mathbf{z},\hat{\mathbf{z}}\right)=\underset{\mathbf{u},\mathbf{S}}{\min }& {\int}_{\Omega^X}\left({d}_E^2(\mathbf{E}\left(\mathbf{u}\right),\hat{\mathbf{E}})+{d}_S^2(\mathbf{S},\hat{\mathbf{S}}\Big)\right)d\Omega \\ {}\mathrm{subject}\hskip0.17em \mathrm{to}: DIV\left(\mathbf{F}\left(\mathbf{u}\right)\cdot \mathbf{S}\right)+\mathbf{b}=\mathbf{0}\hskip0.24em \mathrm{in}\hskip0.24em {\Omega}^X,\\ {}\left(\mathbf{F}\left(\mathbf{u}\right)\cdot \mathbf{S}\right)\cdot \mathbf{N}=\mathbf{t}\hskip0.24em \mathrm{on}\hskip0.24em {\Gamma}_t^X.\end{array}} $$ Note that the Green–Lagrangian strain tensor 
$ \mathbf{E} $ is obtained from the displacement function of sufficient smoothness so the compatibility condition is naturally imposed in the physical solver. Enforcing the physical constraints by the Lagrange multiplier 
$ \boldsymbol{\lambda} $ gives the following functional:
$$ {\mathrm{\int}}_{\Omega^X}\left[{d}_E^2(\mathbf{E}(\mathbf{u}),\hat{\mathbf{E}})+{d}_S^2(\mathbf{S},\hat{\mathbf{S}}\Big)\right]d\Omega +\boldsymbol{\lambda} \left[{\int}_{\Omega^X}\left[ DIV\left(\mathbf{F}\left(\mathbf{u}\right)\cdot \mathbf{S}\right)+\mathbf{b}\right]d\Omega -{\int}_{\Gamma_t^X}\left[\left(\mathbf{F}\left(\mathbf{u}\right)\cdot \mathbf{S}\right)\cdot \mathbf{N}-\mathbf{t}\right]d\Gamma \right]. $$The stationary condition of Equation (12) yields the following data-driven variational equations (He and Chen, Reference He and Chen2020; He et al., Reference He, Laurence, Lee and Chen2020b; Nguyen et al., Reference Nguyen, Rambausek and Keip2020):
$$ \delta \mathbf{u}:\hskip1em {\int}_{\Omega^X}\delta \mathbf{E}\left(\mathbf{u}\right):\mathbf{M}:\left(\mathbf{E}\left(\mathbf{u}\right)-\hat{\mathbf{E}}\right)d\Omega ={\int}_{\Omega^X}\delta {\mathbf{F}}^T\left(\mathbf{u}\right)\cdot \nabla \boldsymbol{\lambda} :\mathbf{S}d\Omega, $$
$$ \delta \mathbf{S}:\hskip1em {\int}_{\Omega^X}\delta \mathbf{S}:\left({\mathbf{M}}^{-1}:\mathbf{S}-{\mathbf{F}}^T\left(\mathbf{u}\right)\cdot \nabla \boldsymbol{\lambda} \right)d\Omega ={\int}_{\Omega^X}\delta \mathbf{S}:{\mathbf{M}}^{-1}:\hat{\mathbf{S}}d\Omega, $$
$$ \delta \boldsymbol{\lambda} :\hskip1em {\int}_{\Omega^X}\delta \nabla \boldsymbol{\lambda} :\left(\mathbf{F}\left(\mathbf{u}\right)\cdot \mathbf{S}\right)d\Omega ={\int}_{\Omega^X}\delta \boldsymbol{\lambda} \cdot \mathbf{b}d\Omega +{\int}_{\Gamma_t^X}\delta \boldsymbol{\lambda} \cdot \mathbf{t}d\Gamma . $$ The variational formulation in Equations (13a–c) can be solved by numerical solvers, for example, the Finite Element Method and the Reproducing Kernel Particle Method (RKPM) (Liu et al., Reference Liu, Jun and Zhang1995; Chen et al., Reference Chen, Pan, Wu and Liu1996). In this study, the RKPM numerical solver is adopted for approximating the displacement field and the Lagrange multiplier field due to its nodal approximation of state and field variables that are particularly effective for data-driven computing, see more details of RKPM in Appendix. The stabilized conforming nodal integration (SCNI) by Chen et al. (Reference Chen, Yoon and Wu2002) scheme is adopted such that the integration points share the same set of points as the discretization nodes, which allows all material data search and variable evaluation to be performed only at the nodal points, enhancing efficiency and accuracy of data-driven computing. Consistent to the nodal integration, the stress 
$ \mathbf{S} $ is approximated by indicator functions 
$ {\mathcal{X}}_i\left(\mathbf{X}\right) $ so that its nodal values are directly associated with the material stress data without introducing additional interpolation errors, see He and Chen (Reference He and Chen2020) for details.
$$ \mathbf{S}\left(\mathbf{X}\right)\hskip.5em \approx \hskip.5em {\mathbf{S}}^h\left(\mathbf{X}\right)=\sum \limits_{i=1}^N{\mathcal{X}}_i\left(\mathbf{X}\right){\mathbf{S}}_i, $$where 
$ N $ is the number of integration points, 
$ {\mathbf{S}}_i $ is the nodal stress associated with the 
$ i $th node, and 
$ {\mathcal{X}}_i\left(\mathbf{X}\right) $ is expressed as
$$ {\mathcal{X}}_i\left(\mathbf{X}\right)=\left\{\begin{array}{cc}1,& \mathrm{if}\hskip1em \mathbf{X}\hskip0.80em \in \hskip0.80em {\Omega}_i^X,\\ {}0,& \mathrm{if}\hskip1em \mathbf{X}\hskip.5em \notin \hskip.5em {\Omega}_i^X.\end{array}\right. $$Employing nodal integration and the stress approximation in Equation (14), the discrete form of Equation (13b) becomes
$$ \sum \limits_{i=1}^N{V}_i\delta {\mathbf{S}}_i^T\left({\mathbf{M}}_i^{-1}{\mathbf{S}}_i-\left({\mathbf{F}}^T{\left(\mathbf{u}\right)}_i\cdot \nabla {\boldsymbol{\lambda}}_i\right)\right)=\sum \limits_{i=1}^N{V}_i\delta {\mathbf{S}}_i^T{\mathbf{M}}_i^{-1}{\hat{\mathbf{S}}}_i,\hskip1em \forall \hskip1em \delta {\mathbf{S}}_i,\hskip1em i=1,\dots N, $$where the subscript 
$ i $ denotes the terms evaluated at 
$ {\mathbf{X}}_i $, e.g., 
$ {\mathbf{S}}_i=\mathbf{S}\left({\mathbf{X}}_i\right) $ and 
$ {\mathbf{F}}^T{\left(\mathbf{u}\right)}_i=\mathbf{F}\left(\mathbf{u}\left({\mathbf{X}}_i\right)\right) $, and 
$ {V}_i $ is the integration weight associated with the 
$ i $th node. Equation (16) yields the stress update for all integration points:
$$ \mathbf{S}=\hat{\mathbf{S}}+\mathbf{M}:\left({\mathbf{F}}^T\left(\mathbf{u}\right)\cdot \nabla \boldsymbol{\lambda} \right). $$ Substituting Equation (17) into Equations (13a) and (13c) yields the following nonlinear system of 
$ \mathbf{u} $ and 
$ \boldsymbol{\lambda} $ (He et al., Reference He, Laurence, Lee and Chen2020b)
$$ {\displaystyle \begin{array}{l}{\int}_{\Omega^X}\delta \mathbf{E}\left(\mathbf{u}\right):\mathbf{M}:(\mathbf{E}\left(\mathbf{u}\right)-\hat{\mathbf{E}})d\Omega ={\int}_{\Omega^X}\left(\delta {\mathbf{F}}^T\left(\mathbf{u}\right)\cdot \nabla \boldsymbol{\lambda} \right):\mathbf{M}:\left({\mathbf{F}}^T\left(\mathbf{u}\right)\cdot \nabla \boldsymbol{\lambda} \right)d\Omega \\ {}\operatorname{}+{\int}_{\Omega^X}\left(\delta {\mathbf{F}}^T\left(\mathbf{u}\right)\cdot \nabla \boldsymbol{\lambda} \right):\hat{\mathbf{S}}d\Omega, \end{array}} $$
$$ {\int}_{\Omega^X}\left({\mathbf{F}}^T\left(\mathbf{u}\right)\cdot \delta \nabla \boldsymbol{\lambda} \right):\left[\mathbf{M}:\left({\mathbf{F}}^T\left(\mathbf{u}\right)\cdot \nabla \boldsymbol{\lambda} \right)+\hat{\mathbf{S}}\right]d\Omega ={\int}_{\Omega^X}\delta \nabla \boldsymbol{\lambda} \cdot \mathbf{b}d\Omega +{\int}_{\Gamma_t^X}\delta \nabla \boldsymbol{\lambda} \cdot \mathbf{t}d\Gamma, $$which can be solved by the Newton–Raphson method (Belytschko et al., Reference Belytschko, Liu and Moran2000)
2.2.2. Material solver
Following the physical step, a material step is performed to search for the optimal material data 
$ {\hat{\mathbf{z}}}^{\ast } $ that is closest to the physical state 
$ {\mathbf{z}}^{\ast } $ obtained from the physical step by minimizing the distance functional. This is done by the following material solver, which decomposes the global minimization in Equation (10) into 
$ N $ local minimization problems (Kirchdoerfer and Ortiz, Reference Kirchdoerfer and Ortiz2016).
$$ {\hat{\mathbf{z}}}_i^{\ast }=\underset{{\hat{\mathbf{z}}}_i\hskip.3em \in \hskip.3em {\unicode{x1D53C}}_{em}}{\arg\;\min}\;{d}_z^2\left({\mathbf{z}}_i^{\ast },{\hat{\mathbf{z}}}_i\right),\hskip1em i=1,\dots, N, $$where 
$ i $ denotes the indices of integration points and 
$ N $ is the total number of integration points.
The data-driven solutions are obtained through fixed-point iterations of the physical step (Equations 17 and 18) and the material step (Equation 19) until the variation in material-step solutions between two consecutive iterations is within a certain tolerance, see Kirchdoerfer and Ortiz (Reference Kirchdoerfer and Ortiz2016) and Conti et al. (Reference Conti, Müller and Ortiz2018) for discussion of convergence properties of this two-step fixed-point iteration solver. The DMDD computing process is depicted in Figure 1, where 
$ \mathrm{\mathcal{E}} $ is the material admissible set, 
$ (v) $ is the iteration index, and one iteration contains one physical step and one material step. In DMDD computing, given a physical state 
$ {\mathbf{z}}^{\ast (v)} $ of a material point, the material solver searches for the closest material data 
$ {\hat{\mathbf{z}}}^{\ast (v)} $ directly from the material database 
$ {\unicode{x1D53C}}_{em} $. However, it can result in unsatisfactory accuracy of data-driven solutions when noise and outliers presented in the material data (Kirchdoerfer and Ortiz, Reference Kirchdoerfer and Ortiz2017; He and Chen, Reference He and Chen2020).
Figure 1. Geometric schematic of the distance-minimizing data-driven (DMDD) solver (Kirchdoerfer and Ortiz, Reference Kirchdoerfer and Ortiz2016).
2.3. Local convexity-preserving material solver
The solutions from the material step are critical to data-driven computing, as they represent the material behaviors of the anisotropic solids. As discussed in Section 2.2, the conventional DMDD material solver (Equation 19) that directly searches for the closest material data can lead to inaccurate data-driven solutions when the material data contains noise and outliers. To enhance the robustness of data-driven computing against noise and outliers, the LCDD computing (He and Chen, Reference He and Chen2020) was developed by introducing the underlying manifold structure of material data to the material solver, which solves the below material step:
$$ {\hat{\mathbf{z}}}_i^{\ast }=\underset{{\hat{\mathbf{z}}}_i\hskip.3em \in \hskip.3em \mathrm{\mathcal{E}}\left({\mathbf{z}}_i^{\ast}\right)}{\arg\;\min}\;{d}_z^2\left({\mathbf{z}}_i^{\ast },{\hat{\mathbf{z}}}_i\right),\hskip1em i=1,\dots, N, $$where 
$ \mathrm{\mathcal{E}}\left({\mathbf{z}}_i^{\ast}\right) $ is a local convex space formed by 
$ k $ material data points that are closest to a given physical state 
$ {\mathbf{z}}_i^{\ast } $ of material point 
$ {\mathbf{X}}_i $, providing a smooth and bounded solution space for optimal material data search, and preserving the convexity of the constructed local material manifold for enhanced robustness and convergence stability. The material step defined in Equation (20) involves two substeps. First, given a physical state 
$ {\mathbf{z}}_i^{\ast } $, 
$ k $ nearest neighbors (material data points), 
$ {\left\{{\hat{\mathbf{z}}}_{\alpha}\right\}}_{\alpha \hskip.3em \in \hskip.3em {\mathcal{N}}_k\left({\mathbf{z}}_i^{\ast}\right)}\subset \hskip.5em {\unicode{x1D53C}}_{em} $, are identified based on the distance measured by 
$ {d}_z^2\left({\mathbf{z}}_i^{\ast },\hat{\mathbf{z}}\right) $, where the indices of the nearest neighbors of 
$ {\mathbf{z}}_i^{\ast } $ are stored in 
$ {\mathcal{N}}_k\left({\mathbf{z}}_i^{\ast}\right) $. Then, a local convex space is constructed based on the collected nearest neighbors 
$ {\left\{{\hat{\mathbf{z}}}_{\alpha}\right\}}_{\alpha \hskip.3em \in \hskip.3em {\mathcal{N}}_k\left({\mathbf{z}}_i^{\ast}\right)} $ of 
$ {\mathbf{z}}_i^{\ast } $, defined as
$$ \mathrm{\mathcal{E}}\left({\mathbf{z}}_i^{\ast}\right)=\mathrm{Conv}\left({\left\{{\hat{\mathbf{z}}}_{\alpha}\right\}}_{\alpha \hskip.3em \in \hskip.3em {\mathcal{N}}_k\left({\mathbf{z}}_i^{\ast}\right)}\right)=\left\{\left.\sum \limits_{\alpha \hskip.3em \in \hskip.3em {\mathcal{N}}_k\left({\mathbf{z}}_i^{\ast}\right)}{w}_{\alpha }{\hat{\mathbf{z}}}_{\alpha}\;\right|\sum \limits_{\alpha \hskip.3em \in \hskip.3em {\mathcal{N}}_k\left({\mathbf{z}}_i^{\ast}\right)}{w}_{\alpha }=1\hskip1em \mathrm{and}\hskip1em {w}_{\alpha}\ge 0,\hskip1em \forall \alpha \hskip.3em \in \hskip.3em {\mathcal{N}}_k\left({\mathbf{z}}_i^{\ast}\right)\right\}, $$ The optimal coefficients 
$ {\mathbf{w}}_i={\left\{{w}_{\alpha}\right\}}_{\alpha \hskip.3em \in \hskip.3em {\mathcal{N}}_k\left({\mathbf{z}}_i^{\ast}\right)} $ are obtained by solving the following minimization problem:
$$ {\displaystyle \begin{array}{ll}{\mathbf{w}}_i^{\ast }=\underset{{\mathbf{w}}_i}{\arg\;\min }& {d}_z^2\left({\mathbf{z}}_i^{\ast },\sum \limits_{\alpha \hskip.3em \in \hskip.3em {\mathcal{N}}_k\left({\mathbf{z}}_i^{\ast}\right)}{w}_{\alpha }{\hat{\mathbf{z}}}_{\alpha}\right)\\ {}\hskip0.84em \mathrm{subject}\hskip0.17em \mathrm{to}:& \sum \limits_{\alpha \hskip.3em \in \hskip.3em {\mathcal{N}}_k\left({\mathbf{z}}_i^{\ast}\right)}{w}_{\alpha }=1,\\ {}& {w}_{\alpha}\ge 0,\hskip1em \forall \alpha \hskip.3em \in \hskip.3em {\mathcal{N}}_k\left({\mathbf{z}}_i^{\ast}\right),\end{array}} $$where 
$ {\mathbf{w}}_i^{\ast }={\left\{{w}_{\alpha}^{\ast}\right\}}_{\alpha \hskip.3em \in \hskip.3em {\mathcal{N}}_k\left({\mathbf{z}}_i^{\ast}\right)} $ denotes the optimal coefficients. Equation (22) is solved by means of a non-negative least-square algorithm with penalty relaxation, see details in He and Chen (Reference He and Chen2020). The optimal material data 
$ {\hat{\mathbf{z}}}_i^{\ast } $ can then be obtained by the following local convex construction:
$$ {\hat{\mathbf{z}}}_i^{\ast }=\sum \limits_{\alpha \hskip.3em \in \hskip.3em {\mathcal{N}}_k\left({\mathbf{z}}_i^{\ast}\right)}{w}_{\alpha}^{\ast }{\hat{\mathbf{z}}}_{\alpha }, $$which ensures that the optimal material data 
$ {\hat{\mathbf{z}}}_i^{\ast } $ always lie within the local convex space 
$ \mathrm{\mathcal{E}}\left({\mathbf{z}}_i^{\ast}\right) $. The LCDD computing process is depicted in Figure 2, where the material solver finds the optimal material data within the local convex space 
$ \mathrm{\mathcal{E}}\left({\mathbf{z}}_i^{\ast}\right) $ (denoted by enclosed black dash lines) formed by the 
$ k $ selected data points closest to the given physical state, which expands the feasible solution space for robustness in data-driven iterations against noise and outliers in the material data (He and Chen, Reference He and Chen2020).
Figure 2. Geometric schematic of the local convexity data-driven (LCDD) solver, where the physical state is projected onto the local convex hulls spanned by the nearest material data points located inside the polygons (He and Chen, Reference He and Chen2020).
2.4. Local convexity-preserving material solver for anisotropic solids
To capture direction-dependent material properties of anisotropic solids, a two-level local data search scheme is introduced in the local convexity-preserving material solver. Let us consider an anisotropic material dataset 
$ \unicode{x1D53C}={\left\{({\hat{\mathbf{E}}}_j,{\hat{\mathbf{S}}}_j)\right\}}_{j=1}^M $ with stress-strain data measured in a global reference frame. Every material point in a physical system is associated with its anisotropic orientations (orientations of material anisotropy), which are represented by the vector of Euler angles 
$ {\boldsymbol{\theta}}_i $ between the local fiber frame and the global reference frame, where 
$ i=1,\dots, N $ is the index of material (integration) points. For simplicity, the methodologies of the proposed material solver for anisotropic solids is illustrated by using the examples with in-plane (two-dimensional) anisotropy, which can be easily extended to three-dimensional problems.
Figure 3a shows an in-plane anisotropic material sample under testing in a global reference frame (
$ {\mathbf{e}}_1^g,{\mathbf{e}}_2^g $), where the dash lines indicate the material anisotropic orientation in 
$ {\mathbf{e}}_1^g $ direction. A bar under uniaxial stretching, as shown in Figure 3b, contains material points 
$ {\mathbf{X}}_i $ associated with certain angles (anisotropic orientations) 
$ {\theta}_i $ between the local fiber frame (
$ {\mathbf{e}}_1^l,{\mathbf{e}}_2^l $) of material point 
$ i $ and the reference frame. The corresponding rotation tensor is defined as
$$ {\mathbf{R}}_i=\left[\begin{array}{ccc}\mathit{\cos}\left({\theta}_i\right)& -\mathit{\sin}\left({\theta}_i\right)& 0\\ {}\mathit{\sin}\left({\theta}_i\right)& \mathit{\cos}\left({\theta}_i\right)& 0\\ {}0& 0& 1\end{array}\right]. $$
Figure 3. (a) Material sample under testing in a reference frame where the dash lines indicate the material anisotropic orientation; (b) uniaxial stretching of a bar. The material behaviors of the material point marked in blue are characterized by the material data from the sample shown in (a).
Applying the rotation 
$ {\mathbf{R}}_i $ to the strain–stress data 
$ \unicode{x1D53C}={\left\{({\hat{\mathbf{E}}}_j,{\hat{\mathbf{S}}}_j)\right\}}_{j=1}^M $ yields the rotated strain–stress data 
$ {\unicode{x1D53C}}_i^{\theta }={\left\{({\hat{\mathbf{E}}}_j^{\theta },{\hat{\mathbf{S}}}_j^{\theta })\right\}}_{j=1}^M $ representing the material behaviors with the anisotropic orientation 
$ {\boldsymbol{\theta}}_i $ under the reference frame:
$$ {\hat{\mathbf{E}}}_j^{\theta }={\mathbf{R}}_i\cdot {\hat{\mathbf{E}}}_j\cdot {\mathbf{R}}_i^T,\hskip1em j=1,\dots, M $$
$$ {\hat{\mathbf{S}}}_j^{\theta }={\mathbf{R}}_i\cdot {\hat{\mathbf{S}}}_j\cdot {\mathbf{R}}_i^T,\hskip1em j=1,\dots, M. $$The rotated material datasets are normally obtained offline with various angles associated with material points in the physical system. Rotated material datasets are then used to reconstruct the optimal material data in the material solver during online data-driven computing.
However, it is impractical to prepare rotated material datasets when a system contains a large number of varying anisotropic orientations associated with material points as a collection of all possible rotated material datasets requires prohibitive memory. A typical example is muscle tissues in musculoskeletal systems that involve randomly oriented fibers. To effectively model anisotropic behaviors in complex material systems, we introduce a two-level local data search scheme into the material solver. To this end, the anisotropic orientation 
$ \boldsymbol{\theta} $ is encoded as an additional feature in material data and physical states of material points. Consequently, the distance between the material data and physical state is not only measured by the distance between their strain–stress values, called state distance, but also the distance between their anisotropic orientations, called anisotropic distance, expressed as
$$ {d}_{zf}^2\left((\mathbf{z},\boldsymbol{\theta} )(\mathbf{z},\hat{\boldsymbol{\theta}},)\right)={d}_z^2\left(\mathbf{z},\hat{\mathbf{z}}\right)+{d}_f^2(\boldsymbol{\theta} \hat{\boldsymbol{\theta}}). $$ The state distance 
$ {d}_z\left(\mathbf{z},\hat{\mathbf{z}}\right) $ is computed by Equations (6) and (7) and the anisotropic distance 
$ {d}_f\left(\boldsymbol{\theta}, \hat{\boldsymbol{\theta}}\right) $ is defined as
$$ {d}_f(\boldsymbol{\theta}, \hat{\boldsymbol{\theta}})=\left\Vert \boldsymbol{\theta} -\hat{\boldsymbol{\theta}}\right\Vert, $$where 
$ \hat{\boldsymbol{\theta}} $ and 
$ \boldsymbol{\theta} $ denote the anisotropic orientations of material data and the material point, respectively.
2.4.1. Level 1 data search
Let us consider a system constituted by an anisotropic material with various anisotropic orientations, for example, a musculoskeletal system consisting muscle fibers with various fiber orientations. The anisotropic material with various anisotropic orientations in the global reference frame exhibits the same material behaviors under the local fiber frame. We first prepare a rotated material database that contains 
$ \hat{m} $ rotated material datasets obtained by rotating the original dataset with 
$ \hat{m} $ different angles (anisotropic orientations). Each rotated material datasets contains strain–stress data 
$ \hat{\mathbf{z}} $ and an anisotropic orientation 
$ \hat{\boldsymbol{\theta}} $. The range of variations in anisotropic orientations of the rotated material database is sufficiently large to cover all material points in the system. Given a physical state of a material point, 
$ \left({\mathbf{z}}_i,{\boldsymbol{\theta}}_i\right) $, where the subscript 
$ i $ denotes the index of a material (integration) point 
$ {\mathbf{X}}_i $, we first compute the anisotropic distance (Equation 27) between the material point (
$ {\boldsymbol{\theta}}_i $) and all rotated material datasets (
$ {\hat{\boldsymbol{\theta}}}_p $, 
$ p=1,2,\dots, \hat{m} $). Two rotated material datasets, 
$ {\unicode{x1D53C}}_p^{\theta } $ and 
$ {\unicode{x1D53C}}_q^{\theta } $, will then be selected so that 
$ {\hat{\boldsymbol{\theta}}}_p\hskip.5em \le \hskip.5em {\boldsymbol{\theta}}_i\hskip.5em \le \hskip.5em {\hat{\boldsymbol{\theta}}}_q $, as illustrated by the blue dash block in Figure 4. The selected material datasets have the anisotropic orientations closest to that of the material point and therefore are considered to best represent the anisotropic material behaviors of the material point.
Figure 4. Illustration of two-level local data search in the proposed material solver.
2.4.2. Level 2 data search
Given the two selected rotated material datasets by the Level 1 search scheme, we search for 
$ k $ nearest neighbors (
$ KNN $) within each dataset based on the state distance between the physical state of the material point and the selected material datasets. The lists of 
$ k $ material data points closest tco the physical state of the material point from the first and second selected material dataset are denoted as 
$ {KNN}_1 $ and 
$ {KNN}_2 $, respectively, as shown in Figure 4. To properly consider the effect of the two 
$ KNN $ datasets, we propose to use a linear weighting scheme based on the anisotropic orientation information as follows:
$$ {\hat{w}}_2=\frac{d_f({\boldsymbol{\theta}}_i,{\hat{\boldsymbol{\theta}}}_p)}{d_f({\hat{\boldsymbol{\theta}}}_q,{\hat{\boldsymbol{\theta}}}_p)}=\frac{\left\Vert {\boldsymbol{\theta}}_i-{\hat{\boldsymbol{\theta}}}_p\right\Vert }{\left\Vert {\hat{\boldsymbol{\theta}}}_q-{\hat{\boldsymbol{\theta}}}_p\right\Vert }, $$
$$ {\hat{w}}_1=1-{\hat{w}}_2. $$ The final list of 
$ k $ nearest neighbors, which is formed by int(
$ k\cdot {\hat{w}}_1 $) (rounded to the nearest integer) nearest neighbors from 
$ {KNN}_1 $ and int(
$ k\cdot {\hat{w}}_2 $) nearest neighbors from 
$ {KNN}_2 $, is used for local convexity-preserving reconstruction of the optimal material data that is closest to the physical state of the material point by
$$ ({\hat{\mathbf{z}}}_i^{\ast },{\hat{\boldsymbol{\theta}}}_i^{\ast })=\underset{({\hat{\mathbf{z}}}_i,{\hat{\boldsymbol{\theta}}}_i)\in \tilde{\mathrm{\mathcal{E}}}({\mathbf{z}}_i^{\ast },{\boldsymbol{\theta}}_i)}{\arg\;\min}\;{d}_z^2({\mathbf{z}}_i^{\ast },{\hat{\mathbf{z}}}_i)+{d}_f^2({\boldsymbol{\theta}}_i,{\hat{\boldsymbol{\theta}}}_i),\hskip1em i=1,\dots, N, $$where 
$ {\mathbf{z}}_i^{\ast } $ is the optimal strain–stress state obtained from the physical step (Equation 9) for the material point 
$ {\mathbf{X}}_i $, 
$ {\boldsymbol{\theta}}_i $ is the anisotropic orientation of the material point 
$ {\mathbf{X}}_i $, 
$ \tilde{\mathrm{\mathcal{E}}}\left({\mathbf{z}}_i^{\ast },{\boldsymbol{\theta}}_i\right) $ is a local convex space formed by the selected 
$ k $ nearest neighbors 
$ {\left\{({\hat{\mathbf{z}}}_{\alpha },{\hat{\boldsymbol{\theta}}}_{\alpha })\right\}}_{\alpha \hskip.3em \in \hskip.3em {\mathcal{N}}_k\left({\mathbf{z}}_i^{\ast },{\boldsymbol{\theta}}_i\right)} $ of the physical state 
$ \left({\mathbf{z}}_i^{\ast },{\boldsymbol{\theta}}_i\right) $ based on the distance measured by 
$ {d}_z^2\left({\mathbf{z}}_i^{\ast },\hat{\mathbf{z}}\right)+{d}_f^2({\boldsymbol{\theta}}_i,\hat{\boldsymbol{\theta}}) $. The anisotropic oritentations 
$ {\boldsymbol{\theta}}_i $ and 
$ {\hat{\boldsymbol{\theta}}}_i $ are normalized by 
$ {\hat{\boldsymbol{\theta}}}_q $ before the calculation of anisotropic distance.
Note that the construction of convexity-preserving local manifold takes the anisotropic orientations of the selected nearest neighbors into account in forming the local convex space 
$ \tilde{\mathrm{\mathcal{E}}}\left({\mathbf{z}}_i^{\ast },{\boldsymbol{\theta}}_i\right) $ so that the anisotropic orientations of nearest neighbors play an important role in reconstructing the optimal material data 
$ ({\hat{\mathbf{z}}}_i^{\ast },{\hat{\boldsymbol{\theta}}}_i^{\ast }) $ closest to the physical state 
$ \left({\mathbf{z}}_i^{\ast },{\boldsymbol{\theta}}_i\right) $. Let 
$ {\mathbf{a}}_i^{\ast } $ and 
$ {\hat{\mathbf{a}}}_i $ denote 
$ \left({\mathbf{z}}_i^{\ast },{\boldsymbol{\theta}}_i\right) $ and 
$ ({\hat{\mathbf{z}}}_i,{\hat{\boldsymbol{\theta}}}_i) $, respectively. The total distance between 
$ {\mathbf{a}}_i^{\ast } $ and 
$ {\hat{\mathbf{a}}}_i $ is denoted by 
$ {d}_{zf}^2\left({\mathbf{a}}_i^{\ast },{\hat{\mathbf{a}}}_i\right) $. The local convex space 
$ \tilde{\mathrm{\mathcal{E}}}\left({\mathbf{a}}_i^{\ast}\right) $ is constructed based on the collected nearest neighbors 
$ {\left\{{\hat{\mathbf{a}}}_{\alpha}\right\}}_{\alpha \hskip.3em \in \hskip.3em {\mathcal{N}}_k\left({\mathbf{a}}_i^{\ast}\right)} $, defined as
$$ \tilde{\mathrm{\mathcal{E}}}\left({\mathbf{a}}_i^{\ast}\right)=\mathrm{Conv}\left({\left\{{\hat{\mathbf{a}}}_{\alpha}\right\}}_{\alpha \hskip.3em \in \hskip.3em {\mathcal{N}}_k\left({\mathbf{a}}_i^{\ast}\right)}\right)=\left\{\left.\sum \limits_{\alpha \hskip.3em \in \hskip.3em {\mathcal{N}}_k\left({\mathbf{a}}_i^{\ast}\right)}{\tilde{w}}_{\alpha }{\hat{\mathbf{a}}}_{\alpha}\;\right|\sum \limits_{\alpha \hskip.3em \in \hskip.3em {\mathcal{N}}_k\left({\mathbf{a}}_i^{\ast}\right)}{\tilde{w}}_{\alpha }=1\hskip1em \mathrm{and}\hskip1em {\tilde{w}}_{\alpha}\ge 0,\hskip1em \forall \alpha \hskip.3em \in \hskip.3em {\mathcal{N}}_k\left({\mathbf{a}}_i^{\ast}\right)\right\}, $$ The optimal coefficients 
$ {\tilde{\mathbf{w}}}_i={\left\{{\tilde{w}}_{\alpha}\right\}}_{\alpha \hskip.3em \in \hskip.3em {\mathcal{N}}_k\left({\mathbf{a}}_i^{\ast}\right)} $ are obtained by solving the following minimization problem:
$$ {\displaystyle \begin{array}{ll}{\tilde{\mathbf{w}}}_i^{\ast }=\underset{{\tilde{\mathbf{w}}}_i}{\arg\;\min }& {d}_{zf}^2\left({\mathbf{a}}_i^{\ast },\sum \limits_{\alpha \hskip.3em \in \hskip.3em {\mathcal{N}}_k\left({\mathbf{a}}_i^{\ast}\right)}{\tilde{w}}_{\alpha }{\hat{\mathbf{a}}}_{\alpha}\right)\\ {}\hskip0.84em \mathrm{subject}\hskip0.17em \mathrm{to}:& \sum \limits_{\alpha \hskip.3em \in \hskip.3em {\mathcal{N}}_k\left({\mathbf{a}}_i^{\ast}\right)}{\tilde{w}}_{\alpha }=1,\\ {}& {\tilde{w}}_{\alpha}\hskip.5em \ge \hskip.5em 0,\hskip1em \forall \alpha \hskip0.80em \in \hskip0.80em {\mathcal{N}}_k\left({\mathbf{a}}_i^{\ast}\right),\end{array}} $$where 
$ {\tilde{\mathbf{w}}}_i^{\ast }={\left\{{\tilde{w}}_{\alpha}^{\ast}\right\}}_{\alpha \hskip.3em \in \hskip.3em {\mathcal{N}}_k\left({\mathbf{a}}_i^{\ast}\right)} $ denote the optimal coefficients. Equation (31) is solved by the non-negative least-square algorithm with penalty relaxation used to solve for the optimal coefficients of the local convex space in LCDD’s material solver (Equation 22), see more details in He and Chen (Reference He and Chen2020). The optimal material data 
$ {\hat{\mathbf{a}}}_i^{\ast }=({\hat{\mathbf{z}}}_i^{\ast },{\hat{\boldsymbol{\theta}}}_i^{\ast }) $ can then be obtained by the following local convex construction:
$$ {\hat{\mathbf{a}}}_i^{\ast }=\sum \limits_{\alpha \hskip.3em \in \hskip.3em {\mathcal{N}}_k\left({\mathbf{a}}_i^{\ast}\right)}{\tilde{w}}_{\alpha}^{\ast }{\hat{\mathbf{a}}}_{\alpha }, $$ The optimal material data sought from the material solver is considered to be “closest” to the physical state in terms of both state distance and anisotropic distance and thus best represent the anisotropic material properties of the material point. The optimal material data will then be input to the physical step (Equations 17 and 18) to solve for the closest physical state in the next data-driven iteration. Note that anisotropic properties are embedded in material database, which are independent of physical laws. Hence, the physical solver only requires the optimal strain–stress material data, 
$ {\hat{\mathbf{z}}}_i^{\ast } $, reconstructed from the material solver.
2.4.3. Schematic illustration
Figure 5 shows schematic examples of the material step in a two-dimensional space, where 
$ {\hat{\boldsymbol{\theta}}}_p={20}^{\circ } $, 
$ {\hat{\boldsymbol{\theta}}}_q={40}^{\circ } $, and 
$ k=6 $ are considered. Three different anisotropic orientations 
$ {\boldsymbol{\theta}}_i $ of the material point 
$ i $ are considered, that is, 
$ {36}^{\circ } $, 
$ {30}^{\circ } $, and 
$ {24}^{\circ } $. The linear weights for selecting the nearest neighbors (
$ KNN $ weights) from each material dataset are computed based on the anisotropic orientation information by Equation (28). The number of nearest neighbors from each material dataset are computed by multiplying their 
$ KNN $ weights with the total number of nearest neighbors 
$ k $, as listed in Table 1, which are rounded to the nearest integer. This approach follows the idea of instance-based learning (Mitchell et al., Reference Mitchell, Agle and Wood1997) to learn the underlying relation from the neighboring manifold of observation data, but additionally, it assigns different number of votes to the selected datasets based on their anisotropic distance 
$ {d}_f $ (Equation 27). For example, in the case that the anisotropic orientation of the material point 
$ i $ is 
$ {\boldsymbol{\theta}}_i={36}^{\circ } $, which is closer to that of Dataset 2 (
$ {\hat{\boldsymbol{\theta}}}_q={40}^{\circ } $) and thus the computed 
$ KNN $ weight for Dataset 2 is 
$ {\hat{w}}_2=0.8 $, larger than 
$ {\hat{w}}_1=0.2 $ of Dataset 1 (
$ {\hat{\boldsymbol{\theta}}}_p={20}^{\circ } $). Hence, there are 1 and 5 nearest neighbors selected from Dataset 1 and Dataset 2, respectively, which are used to form a local convex space capturing the underlying anisotropic data structure, as shown in Figure 5a. The optimal material data (red circle in Figure 5a) is then reconstructed from the local convex space. As the anisotropic orientation of the material point is closer to that of Dataset 2, more nearest neighbors are selected from Dataset 2 for local convex reconstruction and therefore the optimal material data that is closest to the given physical state will have anisotropic properties closer to that of Dataset 2. In the case that the anisotropic orientation is 
$ {\boldsymbol{\theta}}_i={30}^{\circ } $, both datasets have the same 
$ KNN $ weights and contribute the same number of nearest neighbors to the construction of local convex space, as shown in Figure 5b, indicating that they have the same effects on reconstructing the optimal material data.
Figure 5. Schematic illustration of the proposed material solver for anisotropic solids. Dataset 1 has an anisotropic orientation of 
$ {\hat{\boldsymbol{\theta}}}_p={20}^{\circ } $. Dataset 2 has an anisotropic orientation of 
$ {\hat{\boldsymbol{\theta}}}_q={40}^{\circ } $. The total number of nearest neighbors 
$ k $ is 6. The material step of a material point with different anisotropic orientations are compared: (a) 
$ {\boldsymbol{\theta}}_i={36}^{\circ } $; (b) 
$ {\boldsymbol{\theta}}_i={30}^{\circ } $; (c) 
$ {\boldsymbol{\theta}}_i={24}^{\circ } $.
Table 1. The number of nearest neighbors from each dataset in the examples shown in Figure 5 with 
$ {\hat{\boldsymbol{\theta}}}_p={20}^{\circ } $, 
$ {\hat{\boldsymbol{\theta}}}_q={40}^{\circ } $, and 
$ k=6 $.
The weights are computed by Equation (28). 
$ k\cdot {\hat{w}}_1 $ and 
$ k\cdot {\hat{w}}_2 $ are rounded to the nearest integer.
Note that the relation between the anisotropic properties and anisotropic orientations is nonlinear. More sophisticated anisotropic distance metrics and local data reconstruction schemes are required in order to achieve higher accuracy in anisotropic material data reconstruction if the selected material datasets have a large anisotropic distance, which will be investigated in our future study. In the cases that the anisotropic distance of the selected rotated material datasets is small, the 
$ {L}_2 $-based metric for anisotropic distance (Equation 27) and the linear local convexity-preserving reconstruction scheme (Equation 32) adopted in this study can yield desirable accuracy in the data-driven modeling of anisotropic materials, which will be demonstrated in the numerical examples in the following sections.
2.5. Local convexity-preserving data-driven solver for anisotropic solids
Given an anisotropic material dataset 
$ \unicode{x1D53C} $ and anisotropic orientations of material points in the system, the proposed data-driven solver for modeling anisotropic solids is summarized as follows.
Offline stage:
Step 1. Define a range for variations in anisotropic orientations for material data rotation so that it covers that of the material points in the system. Depending on the distribution of the anisotropic orientations in the system, anisotropic orientations for data rotation can be evenly distributed in the defined range or follow certain statistical distributions in order to better capture anisotropic properties of the system.
Step 2. Construct a rotated material database 
$ {\unicode{x1D53C}}_{em}={\unicode{x1D53C}}_1^{\theta}\times \cdots \times {\unicode{x1D53C}}_m^{\theta } $ by rotating the original anisotropic material dataset 
$ \unicode{x1D53C} $ with the anisotropic orientations defined in Step 1. The rotation of strain and stress data is performed by Equation (25).
Online stage:
Step 1. Randomly initialize 
$ {\hat{\mathbf{z}}}_i^{(0)} $, 
$ i=1,\dots, N $, where 
$ i $ denote the indices of material points and 
$ N $ is the number of material points in the system, and set the data-driven iteration index 
$ v=0 $.
Step 2. Solve the physical step (Equations 13–19) for physical states of all material points, 
$ {\mathbf{z}}_i^{\ast (v)} $, 
$ i=1,\dots, N $.
Step 3. For each physical state at each integration point, 
$ ({\mathbf{z}}_i^{\ast (v)},{\boldsymbol{\theta}}_i) $, perform two-level local data search (material step):
3.1. Search for two rotated material datasets with anisotropic orientations 
$ {\hat{\boldsymbol{\theta}}}_p $ and 
$ {\hat{\boldsymbol{\theta}}}_q $ so that 
$ {\hat{\boldsymbol{\theta}}}_p\le {\boldsymbol{\theta}}_i\le {\hat{\boldsymbol{\theta}}}_q $.
3.2. From each selected material dataset, search for 
$ k $ material data points that are closest to the physical state based on the state distance 
$ {d}_z^2({\mathbf{z}}_i^{\ast (v)},\hat{\mathbf{z}}) $ and obtain the final 
$ KNN $ using the weights computed by Equation (28).
3.3. Construct local convex space by solving Equation (31) and obtain the optimal material data 
$ {\hat{\mathbf{z}}}_i^{\ast (v)} $ within the local convex space by Equation (32).
Step 4. Update 
$ v=v+1 $. If 
$ \underset{i=1,\dots, N}{\mathit{\max}}{d}_z^2({\hat{\mathbf{z}}}_i^{\ast (v)},{\hat{\mathbf{z}}}_i^{\ast \left(v-1\right)})> $ tolerance, repeat Step 2–4.
Step 5. Data-driven solution: 
$ {\mathbf{z}}_i^{\ast }={\hat{\mathbf{z}}}_i^{\ast (v)} $, 
$ i=1,\dots, N $.
Remark: Modeling anisotropic materials by classical computational mechanics requires online rotation of element stiffness matrices from local fiber frames to the global reference frame before global assembly, which could lead to high computational cost especially for large systems. In the proposed data-driven modeling of anisotropic materials, online rotation is replaced with an offline rotation, where anisotropic material datasets are rotated by various anisotropic orientations under the global reference frame to form a rotated material database for efficient online data-driven computing. In the general three-dimensional orthotropic case, each anisotropic orientation is associated with three Euler angles and thus for every anisotropic orientation considered, offline rotation operations associated with three Euler angles are required for every stress–strain data, which could be CPU intensive but can be performed in parallel especially when dealing with a large amount of data. Once the rotated material database is constructed offline, it can be efficiently applied to different material points in data-driven modeling and to different systems with the same anisotropic materials. Note that when the size of anisotropic material datasets is large and the range of variations in anisotropic orientations of the system is very broad, the offline rotated material database may require a large amount of memory.
If memory resources are not available, one can also adopt online rotation in the data-driven computing for modeling anisotropic materials. In this case, at each data-driven iteration, the physical step remains unchanged. In the material step, the physical states are first rotated from the global reference frame to local fiber frames of material points and then the LCDD material solver ( Equations 20–23) can be applied to find the optimal material data from the unrotated anisotropic material dataset. The optimal material data is then rotated back to the global reference frame for the next data-driven iteration. Compared to the proposed data-driven approach with offline rotation, this online rotation approach requires higher CPU but less memory.
3. Results and Discussion
3.1. Preparation of material datasets
For the following numerical demonstration, the two-dimensional Saint Venant–Kirchhoff phenomenological model with isotropic and orthotropic elastic tensors are respectively considered as the reference models and used to generate synthetic data for data-driven computing. The plane-stress isotropic elastic stress–strain relation in the Voigt notation is expressed as
$$ \left[\begin{array}{c}{S}_{11}\\ {}{S}_{22}\\ {}{S}_{12}\end{array}\right]=\frac{E}{\left(1+\nu \right)\left(1-2\nu \right)}\left[\begin{array}{ccc}1-\nu & \nu & 0\\ {}\nu & 1-\nu & 0\\ {}0& 0& \frac{1-2\nu }{2}\end{array}\right]\left[\begin{array}{c}{E}_{11}\\ {}{E}_{22}\\ {}2{E}_{12}\end{array}\right], $$where 
$ E $ denotes the Young’s modulus and 
$ \nu $ is the Poisson’s ratio. As another reference material model, the plane-stress orthotropic elastic stress–strain relation in the Voigt notation is expressed as
$$ \left[\begin{array}{c}{S}_{11}\\ {}{S}_{22}\\ {}{S}_{12}\end{array}\right]=\frac{1}{1-{\nu}_{12}{\nu}_{21}}\left[\begin{array}{ccc}{E}_1& {\nu}_{21}{E}_1& 0\\ {}{\nu}_{12}{E}_2& {E}_2& 0\\ {}0& 0& {G}_{12}\left(1-{\nu}_{12}{\nu}_{21}\right)\end{array}\right]\left[\begin{array}{c}{E}_{11}\\ {}{E}_{22}\\ {}2{E}_{12}\end{array}\right], $$where 
$ {E}_1 $ and 
$ {E}_2 $ are Young’s moduli in the 
$ {\mathbf{e}}_1^g $ and 
$ {\mathbf{e}}_2^g $ directions of the reference frame, respectively, see Figure 3a. 
$ {\nu}_{12} $ and 
$ {\nu}_{21} $ are Poisson’s radios. 
$ {G}_{12} $ is the shear modulus.
To demonstrate the robustness of the proposed data-driven computing framework against noise presented in given material datasets, noisy material datasets are generated and the procedure is described below. First, a noiseless material dataset, 
$ {\unicode{x1D53C}}^0={\left\{({\hat{\mathbf{E}}}_i^0,{\hat{\mathbf{S}}}_i^0)\right\}}_{i=1}^M $ with 
$ M={20}^3 $, is generated, where each strain component is uniformly distributed within a certain range, for example, [−0.2, 0.2], and the stress components are obtained by using the orthotropic material model in Equation (34). 
$ M={20}^3 $ strain and stress data points with 
$ {E}_1={10}^4,{E}_2=2.5\times {10}^3,{\nu}_{21}=0.1,{\nu}_{12}=0.4 $, and 
$ {G}_{12}=4.8\times {10}^3 $ are shown in Figure 6a,c, respectively, serving as a noiseless base dataset without any noise and rotation. The anisotropic orientation of the noiseless base dataset is along the horizontal direction of the reference frame, as shown in Figure 3a. Then, each component of the noiseless base dataset is perturbed by Gaussian noise with a scaling factor, 0.4
$ {\overline{\mathbf{z}}}_{max}/\sqrt[3]{M} $, where 
$ {\overline{\mathbf{z}}}_{max} $ is a vector of the maximum values for each component among the noiseless dataset. Figure 6b,d show the strain and stress of the noisy base dataset, 
$ \unicode{x1D53C}={\left\{({\hat{\mathbf{E}}}_i,{\hat{\mathbf{S}}}_i)\right\}}_{i=1}^M $ with 
$ M={20}^3 $, generated from the noiseless base dataset 
$ {\unicode{x1D53C}}^0 $ shown in Figure 6a,c, respectively.
Figure 6. Material datasets with 8,000 data points, 
$ {E}_1={10}^4,{E}_2=2.5\times {10}^3,{\nu}_{12}=0.1,{\nu}_{21}=0.4 $, and 
$ {G}_{12}=4.8\times {10}^3 $: (a) noiseless strain data; (b) noisy strain data; (c) noiseless stress data; and (d) noisy stress data.
Given anisotropic orientations of material points in a physical system, for example, 
$ {\theta}_i $ of the material point 
$ i $ in Figure 3b, the strain and stress data of the base dataset can be rotated using Equation (25). Figure 7 shows the strain and stress data rotated from the base dataset by three different angles 
$ {30}^{\circ } $, 
$ {60}^{\circ } $, and 
$ {90}^{\circ } $. In the following numerical examples, rotation of material datasets is performed offline by various angles (anisotropic orientations) to cover the range of variations in anisotropic orientations of systems. A collection of rotated material datasets form a rotated material database, which is then used for online data-driven computing with the material solver equipped with two-level data search designed for capturing material’s directional dependence.
Figure 7. Rotated material datasets with 8,000 data points, 
$ {E}_1={10}^4,{E}_2=2.5\times {10}^3,{\nu}_{12}=0.1,{\nu}_{21}=0.4 $, and 
$ {G}_{12}=4.8\times {10}^3 $: (a) strain data rotated by 
$ {30}^{\circ } $; (b) strain data rotated by 
$ {60}^{\circ } $; (c) strain data rotated by 
$ {90}^{\circ } $; (d) stress data rotated by 
$ {30}^{\circ } $; (e) stress data rotated by 
$ {60}^{\circ } $; and (f) stress data rotated by 
$ {90}^{\circ } $.
3.2. Multi-layer anisotropic cantilever beam subjected to a tip shear load
We first investigate a multi-layer anisotropic cantilever beam benchmark problem to verify the proposed data-driven modeling framework for anisotropic nonlinear solids, where three layers of orthotropic materials are considered. Two cases are examined and compared. In Case 1, all material points in the beam have an identical anisotropic orientation, which is along the horizontal direction, as shown in Figure 8a. In Case 2, the beam is made of three layers of anisotropic materials with different anisotropic orientations. The material points within each layer have the same anisotropic orientation. Figure 8b shows that the bottom, the middle, and the top layers of the beam have anisotropic orientations of 
$ -{45}^{\circ } $, 
$ {0}^{\circ } $, and 
$ {45}^{\circ } $, respectively. A synthetic noiseless orthotropic material dataset with 8,000 strain–stress data points is first generated from the phenomenological orthotropic elastic model with 
$ {E}_1={10}^4,{E}_2=2.5\times {10}^3,{\nu}_{21}=0.1,{\nu}_{12}=0.4 $, 
$ {G}_{12}=4.8\times {10}^3 $, and a range of [−0.2, 0.2] for each strain component, as shown in Figure 6a,c. Then, a noisy orthotropic material dataset is generated from the noiseless dataset using the procedure described in Section 3.1, as shown in Figure 6b,d.
Figure 8. Schematic of cantilever beam subjected to a tip shear load: (a) Case 1: material points have only one anisotropic orientation 
$ {0}^{\circ } $; (b) Case 2: material points located in different layers of the beam have different anisotropic orientations. The anisotropic orientations of the bottom, the middle, and the top layers are 
$ -{45}^{\circ } $, 
$ {0}^{\circ } $, and 
$ {45}^{\circ } $, respectively. 
$ P=10{E}_1I/{L}^2 $, and 
$ I={H}^3/12 $.
For Case 1, no rotation of the material dataset is required as the anisotropic orientations of materials points in the beam are identical with that of the synthetic dataset, which is along the horizontal direction. For Case 2, a rotated material database is constructed by rotating the generated synthetic noiseless/noisy dataset from 
$ -{60}^{\circ } $ to 
$ {60}^{\circ } $ with a 
$ {10}^{\circ } $ interval, that is, 
$ -{60}^{\circ } $, 
$ -{50}^{\circ } $, 
$ -{40}^{\circ } $, 
$ -{30}^{\circ } $, 
$ -{20}^{\circ } $, 
$ -{10}^{\circ } $, 
$ {0}^{\circ } $, 
$ {10}^{\circ } $, 
$ {20}^{\circ } $, 
$ {30}^{\circ } $, 
$ {40}^{\circ } $, 
$ {50}^{\circ } $, and 
$ {60}^{\circ } $. The rotated material database is then used for online data-driven computing. The numerical studies of both cases are performed with noiseless or noisy material datasets. A tip shear load 
$ P=10{E}_1I/{L}^2 $ with 
$ I={H}^3/12 $ is applied and the data-driven analysis is performed with 20 loading steps.
The normalized tip deflection-loading and stress distribution of data-driven solutions are compared with those of the constitutive model-based reference solutions obtained from an in-house Finite Element code, as shown in Figures 9 and 10. The data-driven solutions of both cases show a satisfactory agreement with the model-based reference solutions, which demonstrates the effectiveness of the proposed data-driven computing framework for modeling anisotropic nonlinear elastic materials. The data-driven solutions obtained by using the noisy material datasets agree well with the model-based reference solutions, showing the robustness of the proposed framework against noise in the datasets, which is achieved by the local convexity-preserving reconstruction in the material solver.
Figure 9. Comparison of data-driven solutions with constitutive model-based reference solutions. Normalized tip deflection-loading, where 
$ I={H}^3/12 $: (a) Case 1 and (b) Case 2.
Figure 10. Comparison of data-driven solutions with constitutive model-based reference solutions. Distribution of 
$ {S}_{xx} $: (a) Case 1: reference solution; (b) Case 1: data-driven solution with noiseless data; (c) Case 1: data-driven solution with noisy data; (d) Case 2: reference solution; (e) Case 2: data-driven solution with noiseless data; and (f) Case 2: data-driven solution with noisy data.
3.3. Anisotropic cylinder subjected to internal pressure
To further evaluate the effectiveness of the proposed data-driven computing framework in handling physical systems with a large variation in anisotropic orientations, an anisotropic cylinder subjected to internal pressure is analyzed, as shown in Figure 11a, which is composed of an orthotropic material with anisotropic orientations along the circumferential direction of the cylinder. Considering axisymmetry of the geometry and the applied load, only a quarter model shown in Figure 11b is modeled. The inner and outer radius are 1 and 2, respectively. In Figure 11c, the anisotropic orientations of the 
$ 10\times 20 $ discretization points are denoted by red arrows. Two cases are examined and compared.
Figure 11. (a) Schematic of a cylinder subjected to internal pressure and a quarter model to be simulated and (b) red arrows denote nodal anisotropic orientations of material points in a discretization with 
$ 10\times 20 $ nodes.
In Case 1, an isotropic elastic material is applied, with Young’s modulus 
$ E=9\times {10}^3 $ and Poisson’s ratio 
$ \nu =0.2 $. A noiseless synthetic isotropic material dataset with 8,000 strain–stress data points is generated by the phenomenological isotropic elastic material model with a range of [−0.4, 0.4] for each strain component. In Case 2, an orthotropic elastic material is applied, with 
$ {E}_1=4\times {10}^4,{E}_2=9\times {10}^3,{\nu}_{21}=0.045,{\nu}_{12}=0.2 $, 
$ {G}_{12}=2\times {10}^4 $, and anisotropic orientations along the circumferential direction of the cylinder. A noiseless synthetic orthotropic material dataset with 8,000 strain–stress data points is generated by the phenomenological orthotropic elastic material model with a range of [−0.4, 0.4] for each strain component. Considering the uniform distribution of anisotropic orientations of material points in the cylinder, a rotated material database is constructed by rotating the generated synthetic dataset from 
$ {90}^{\circ } $ to 
$ {180}^{\circ } $ with a 
$ {5}^{\circ } $ interval. Noisy material datasets are generated from the noiseless material datasets by the procedure described in Section 3.1. Internal pressure p = 1,000 is applied and the data-driven analysis is performed with 20 loading steps.
The cross-sectional (
$ y=0 $) radial displacement (
$ {U}_r $) as well as the radial and the circumferential stress distributions of the data-driven solutions are compared with those of the constitutive model-based reference solutions, as shown in Figures 12–14. The data-driven solutions obtained from using noisy material datasets are very close to those of reference solutions, which again shows the robustness of the proposed data-driven framework to deal with noisy datasets. The results of both cases show a good agreement with the reference solutions, which demonstrates the effectiveness of the proposed data-driven framework in dealing with systems with a large variation in anisotropic orientations.
Figure 12. Comparison of data-driven solutions with constitutive model-based reference solutions. Cross-sectional radial displacement 
$ {U}_r $: (a) Case 1 and (b) Case 2.
Figure 13. Comparison of data-driven solutions with constitutive model-based reference solutions of Case 1. Distribution of 
$ {S}_{rr} $ (stress in the radial direction) and 
$ {S}_{\theta \theta} $ (stress in the circumferential direction): (a) Case 1: reference 
$ {S}_{rr} $; (b) Case 1: data-driven solution with noiseless data of 
$ {S}_{rr} $; (c) Case 1: data-driven solution with noisy data of 
$ {S}_{rr} $; (d) Case 1: reference 
$ {S}_{\theta \theta} $; (e) Case 1: data-driven solution with noiseless data of 
$ {S}_{\theta \theta} $; and (f) Case 1: data-driven solution with noisy data of 
$ {S}_{\theta \theta} $.
Figure 14. Comparison of data-driven solutions with constitutive model-based reference solutions of Case 1. Distribution of 
$ {S}_{rr} $ (stress in the radial direction) and 
$ {S}_{\theta \theta} $ (stress in the circumferential direction): (a) Case 2: reference 
$ {S}_{rr} $; (b) Case 2: data-driven solution with noiseless data of 
$ {S}_{rr} $; (c) Case 2: data-driven solution with noisy data of 
$ {S}_{rr} $; (d) Case 2: reference 
$ {S}_{\theta \theta} $; (e) Case 2: data-driven solution with noiseless data of 
$ {S}_{\theta \theta} $; and (f) Case 2: data-driven solution with noisy data of 
$ {S}_{\theta \theta} $.
4. Conclusion
In this study, we develop a new data-driven material solver built upon the local convexity-preserving reconstruction scheme by He and Chen (Reference He and Chen2020) to capture anisotropic material behaviors and enable data-driven modeling of anisotropic nonlinear elastic solids. The proposed data-driven approach assumes that the material data of anisotropic materials with a specific anisotropic orientation in a reference frame is accessible. The information of anisotropic orientations, for example, the rotation angles between local fiber frames and the reference frame of the material data are utilized to construct an offline material database, which contains rotated material datasets representing anisotropic material properties with various anisotropic orientations. The offline rotated material database can be efficiently constructed and applied to data-driven simulations of anisotropic materials.
During online data-driven computing, a two-level local data search is integrated into the local convexity-preserving material solver. In the Level 1 data search, two rotated material datasets with minimum anisotropic distance (Equation 27) to the material anisotropic orientation are identified as the local datasets. The selected rotated material datasets are considered to contain anisotropic material properties closest to that of the material point. In the Level 2 data search, 
$ k $ data points closest to the given physical state based on strain–stress state distance (Equation 5) are obtained separately from the two rotated material datasets selected from the Level 1 data search. A linear weighting scheme based on anisotropic distance is adopted to determine the number nearest data points from each of the selected material datasets. The final 
$ k $ nearest data points that contain information of strain, stress, and anisotropic orientations are used to construct a local convex space capturing the underlying anisotropic data structure, within which the optimal material data is reconstructed by the material solver. The optimally reconstructed material data is closest to the physical state of the material point in terms of both anisotropic distance and state distance and thus is considered to best represent the anisotropic properties of the material point.
The performance of the proposed data-driven computing framework is examined by two numerical examples, including deflection of a multi-layer anisotropic cantilever beam made of materials with different anisotropic orientations and inflation of an anisotropic cylinder where anisotropic orientations are along the circumferential direction of the cylinder. Synthetic noiseless and noisy material data are generated from the phenomenological material models and employed for the data-driven analysis for accuracy assessment of the proposed method. The data-driven solutions show a good agreement with the constitutive model-based reference solutions, which demonstrates its effectiveness and robustness of the proposed data-driven framework against noise present in the material data.
The proposed two-level data search can achieve high computational efficiency if the rotated material database is constructed offline such that online computation does not involve any frame transformation of states or data. For applications with a large amount of anisotropic material data and strong variations in anisotropic orientations, constructing a rotated material database with small anisotropic distance between rotated datasets requires a large amount of memory. In this case, more sophisticated metrics for anisotropic distance and reconstruction schemes will be investigated in order to achieve high accuracy in reconstructing anisotropic material properties from given material datasets that have a large anisotropic distance. In future studies, the proposed data-driven anisotropic modeling framework will be further examined with real-world data on three-dimensional material systems with aniostropic material behaviors, e.g., musculoskeletal systems consisting of muscle fibers with varying anisotropic orientations.
Funding Statement
The support of this work by the National Science Foundation under Award Number CCF-1564302 to the first, second, and the third authors, and the National Institute of Health under grant number 5R01AG056999-03 to the first, third, forth, and fifth authors is very much appreciated.
Competing Interests
The authors declare no competing interests.
Data Availability Statement
All data involved in this study are synthetically generated and data generation procedures are described in Section 3.1.
Author contributions
Conceptualization, J.-S.C., Q.H., X.H.; Methodology, X.H., Q.H., J.-S.C.; Formal Analysis, X.H.; Data Curation, X.H.; Writing–Original Draft, X.H.; Writing–Review and Editing, X.H., Q.H., J.-S.C., U.S., S.S.; Supervision, J.-S.C.; Funding Acquisition, J.-S.C., U.S., S.S.; All authors approved the final submitted draft.
Appendix Reproducing Kernel Particle Method
The reproducing kernel (RK) approximations (Liu et al., Reference Liu, Jun and Zhang1995; Chen et al., Reference Chen, Pan, Wu and Liu1996) is adopted in the physical solver of the proposed data-driven framework due to its nodal approximation of state and field variables that are particularly effective for data-driven computing. The RK approximation functions can be constructed to possess desired completeness and continuity, which are determined by basis functions and kernel functions, respectively. Figure 15a shows a domain 
$ \Omega $ discretized by a set of nodes.
Figure 15. (a) A domain 
$ \Omega $ discretized by the a set of RK nodes; (b) a cubic B-spline function widely used as a kernel function in RK approximation; and (c) an example of RK approximation function centered at 
$ X=5 $ with a support size 
$ a=1.5 $ 
$ \times $ (nodal spacing).
The RK approximation of a field variable 
$ \mathbf{u}\left(\mathbf{X}\right) $ is given by
$$ \mathbf{u}\left(\mathbf{X}\right)\approx {\mathbf{u}}^h\left(\mathbf{X}\right)=\sum \limits_{i=1}^{NP}{\Psi}_i\left(\mathbf{X}\right){\mathbf{d}}_i, $$where 
$ NP $ is the number of discretization nodes, 
$ {\mathbf{d}}_i $ is the nodal coefficient associated with the 
$ i $th node, and 
$ {\Psi}_i\left(\mathbf{X}\right) $ is the RK approximation function expressed as
$$ {\Psi}_i\left(\mathbf{X}\right)={\mathbf{H}}^T\left(\mathbf{X}-{\mathbf{X}}_i\right)\mathbf{b}\left(\mathbf{X}\right){\phi}_a\left(\mathbf{X}-{\mathbf{X}}_i\right), $$where 
$ {\mathbf{H}}^T\left(\mathbf{X}-{\mathbf{X}}_i\right)=\left[1,{X}_1-{X}_{1i},{X}_2-{X}_{2i},{X}_3-{X}_{3i},\dots, {\left({X}_3-{X}_{3i}\right)}^n\right] $ is a vector of monomial basis functions up to 
$ n $th order, and 
$ {\phi}_a\left(\mathbf{X}-{\mathbf{X}}_i\right) $ is a kernel function with a local support size “
$ a $,” which controls the smoothness of the RK approximation function. For example, the cubic B-spline kernel function is widely used as a kernel function in RK approximation, see Figure 15b,
$$ {\phi}_a(y)=\left\{\begin{array}{cc}\frac{2}{3}-4{y}^2+4{y}^3,& 0\le y<\frac{1}{2}\\ {}\frac{4}{3}-4y+4{y}^2-\frac{4}{3}{y}^3,& \frac{1}{2}\le y<1\\ {}0,& y\ge 1\end{array}\right.\hskip1em \mathrm{with}\hskip1em y=\frac{\left\Vert \mathbf{X}-{\mathbf{X}}_i\right\Vert }{a}. $$ The vector 
$ \mathbf{b}\left(\mathbf{X}\right) $ in Equation (36) is a parameter vector determined by imposing the 
$ n $-th order reproducing conditions (Liu et al., Reference Liu, Jun and Zhang1995; Chen et al., Reference Chen, Pan, Wu and Liu1996),
$$ \sum \limits_{i=1}^{NP}{\Psi}_i\left(\mathbf{X}\right){X}_{1i}^i{X}_{2i}^j{X}_{3i}^k={X}_1^i{X}_2^j{X}_3^k,\hskip1em \mid i+j+k\mid =0,1,\dots, n. $$ Substituting Equation (36) into Equation (38) yields 
$ \mathbf{b}\left(\mathbf{X}\right)={\mathbf{M}}^{-1}\left(\mathbf{X}\right)\mathbf{H}\left(\mathbf{0}\right) $, where 
$ \mathbf{M}\left(\mathbf{X}\right) $ is a moment matrix given by
$$ \mathbf{M}\left(\mathbf{X}\right)=\sum \limits_{i=1}^{NP}\mathbf{H}\left(\mathbf{X}-{\mathbf{X}}_i\right){\mathbf{H}}^T\left(\mathbf{X}-{\mathbf{X}}_i\right){\phi}_a\left(\mathbf{X}-{\mathbf{X}}_i\right). $$The RK approximation function is expressed as, see Figure 15c,
$$ {\Psi}_i\left(\mathbf{X}\right)={\mathbf{H}}^T\left(\mathbf{0}\right){\mathbf{M}}^{-1}\left(\mathbf{X}\right)\mathbf{H}\left(\mathbf{X}-{\mathbf{X}}_i\right){\phi}_a\left(\mathbf{X}-{\mathbf{X}}_i\right). $$
Open access
Comments
No Comments have been published for this article.