On the Single-Point Calculation of Stress–Strain Data under Large Deformations with Stress and Mixed Control

Stress–strain data with a given constitutive model of material can be calculated directly at a single material point. In this work, we propose a framework to perform single-point calculations under large deformations with stress and mixed control, to test and validate sophisticated constitutive models for materials. Inspired by Galerkin–FFT methods, a well-defined mask projector is used for stress and mixed control, and the derived nonlinear equations are solved in Newton iterations with Krylov solvers, simplifying implementation. One application example of the single-point calculator in developing sophisticated models for anisotropic single crystal rate-independent elastoplasticity is given, illustrating that the proposed algorithm can simulate asymmetrical deformation responses under uni-axial loading. Another example for artificial neural network models of the particle reinforced composite is also given, demonstrating that the commonly used machine learning or deep learning modeling frameworks can be directly incorporated into the proposed calculator. The central difference approximation of the tangent is validated so that derivative-free calculations for black-box constitutive models are possible. The proposed Python-coded single-point calculator is shown to be capable of quickly building, testing, and validating constitutive models with sophisticated or implicit structures, thus boosting the development of novel constitutive models for advanced solid materials.


Introduction
Constitutive models are important mathematical tools to describe material responses to complex mechanical loading conditions. With increasing demands for studying in depth and simulating deformations of different types of advanced materials [1][2][3][4], their corresponding constitutive models become increasingly complex. For example, crystal plasticity models for metal materials [5] under large deformations describe several microscopic deformation mechanisms, such as dislocation motions, twinning, and phase transformation, using multiple internal variables, corresponding equations, and even more than two loops of iteration procedures to build a sophisticated model structures [6]. Constructing an advanced constitutive model is challenging even for well-trained researchers. In addition, a developed complex constitutive model needs to be tested with experimental results or microscopic calculations, before it is regarded as a reliable model. Thus, it is not trivial to have a simple calculation framework to test advanced constitutive models, debug the program for models, and facilitate parameter calibration [7]. In addition, a potential application of the sophisticated constitutive models is to generate microscopic calculation data, for macroscopic model-free data-driven mechanics [8] or for establishing macroscopic data-driven machine learning models [9][10][11][12][13]. Thus, a calculation framework that quickly calculates and outputs data for demanded field variables (e.g., stress-strain data) is highly recommended.
To implement one constitutive model, a popular method is to code the model as one subroutine to the program of finite element methods (FEM), such as ABAQUS [14][15][16][17][18] and LS-DYNA [19][20][21]. To validate the constitutive model, the structure simulation or oneelement calculation is performed in FEM programs with the constitutive model subroutines and then compared with experimental results or microscopic calculations. For FEM simulations, nonlinear geometries, nonlinear boundary conditions, and chosen interpolation functions in elements may produce extra calculation errors, which sometimes need much more effort to overcome than a nonlinear constitutive model does. To validate constitutive models or calibrate parameters, different paths of deformation loading for one material point are needed to benchmark the models. In that case, a single-point calculator can be applied as an alternative but more efficient calculation framework to FEM.
The single-point calculator runs in a single material point, calculating the stress response with the prescribed strain directly with the constitutive model. For example, when a simple shear is simulated, a prescribed deformation gradient can be given in a form (γ = 0) as follows: indicating no change of volume, as its determinant equals one. The stress can then be calculated directly in the single-point calculator, with a chosen strain measurement calculated from this prescribed deformation gradient. Calculating stress directly from the strains is called the problem of strain control. However, the inverse problem in the singlepoint calculator, calculating the strain from the stress, is not direct. It is often impossible to inverse advanced constitutive models, which use complicated mechanisms to calculate the stress from the strain. In addition, in some cases of loading, for example, uni-axial loading, mixed components of stress and strain are given in the single-point calculator to determine the other components. Especially for the anisotropic incompressible solid whose determinant of deformation gradient equals 1, the components of the deformation gradient corresponding to the other two dilatation directions are hard to prescribe accurately and directly in the case of uni-axial loading [22]. The problem to solve is to calculate the deformation gradient or the missing components of the stress/deformation gradient with given stress or mixed components, which corresponds to the problem of stress and mixed control. Therefore, to extend the applications of the single-point calculator, the stress and mixed control boundary conditions must be considered.
Inspired by the Galerkin-FFT methods [23][24][25], we propose an algorithm for the problem of stress and mixed control in the single-point calculator. Analogously to the algorithm for stress and mixed control in Galerkin-FFT methods [25,26], our calculator uses a well-defined mask projector to maintain a unified form of core equations to solve for strain control, stress control, or mixed control. As in the Galerkin-FFT methods, the core equations are to be solved with a Krylov solver, such as conjugate gradient (CG) [24,27] and minimum residual (MinRes) [28], which allows circumventing the difficulties associated with the singularity of the tangent stiffness tensor. These two characteristics (the usage of a mask projector and a Krylov solver) of our proposed algorithm make the calculator simple to implement. Details of the algorithm are given in Section 2. The proposed singlepoint calculator and the associated constitutive models are recommended to be coded in Python, to facilitate the procedure of testing new algorithms, validating sophisticated mathematical structures, and identifying parameters for constitutive models. With the built-in scientific calculation package, Numpy and Scipy [29], the speed of calculation with Python is acceptable even for industrial demands. Two examples of applications of the single-point calculator for the development of constitutive models are given in Section 3. The first example for anisotropic single crystal plasticity (Section 3.1) illustrates that the proposed calculation framework is simple to develop models with sophisticated algorithms and mathematical structures. The second example of the artificial neural network elasticity model (Section 3.2) illustrates the single-point calculator's compatibility with machine learning methods or deep learning methods [30], which can be implemented in commonly used Python-entranced frameworks, such as TensorFlow [31] or PyTorch [32]. Even though the exact analytical consistent tangent is hard to obtain for implicit networks of constitutive models, the central difference approximation of the tangent [33,34] can be used in the single-point calculator. The derivative-free calculations for black-box constitutive models are validated in the second example.
Notation. In this work, scalars, second-order tensors, and fourth-order tensors are denoted as s, T, and K, respectively. The double dot product between two second-order tensors is denoted as P : F and calculated as R = P ij F ij (i, j = 1, 2, 3), where the Einstein summation convention is used. The double dot product between one fourth-order tensor and one second-order tensor is denoted as C : e and calculated as R ij = C ijmn e mn . The operator ( * ) T denotes the transpose, while ( * ) −1 means the inverse of a tensor. The operator det( * ) represents the determinant of the matrix corresponding to the tensor.

Strain Control in the Single-Point Calculator
We consider two conjugated stress and strain tensors, the first Piola-Kirchhoff stress P and the deformation gradient F. For the strain control problems, the single-point calculator is a numerical tool to calculate the corresponding values of P from F with a given constitutive model. If values of F corresponding to finite points in the loading path are given, the single-point calculator can simply calculate the stress responses P one by one. When the prescribed tensor is divided into N segments, the n-th stress response is calculated with the n-th deformation gradient as follows: Equation (1) can be used for the elasticity-like solids. As for the materials depending on the loading history, such as visco-elastic or plastic solids, the rate of the deformation gradient is also needed. We can use the deformation gradient in the former step to calculate the current step: where ∆t is the time increment and α represents internal variables. Here, the rate of the deformation gradient can be approximated as follows: However, the rate of deformation gradient is often replaced by the velocity gradient, which can be approximated as follows: Compared with reliable experimental results of the same loading path, the calculated responses can be used to validate the developed constitutive model (Equations (1) and (2)) or to identify model parameters.

Stress and Mixed Control in the Single-Point Calculator
When the single-point calculator is controlled by the strain, i.e., calculating the stress with the prescribed deformation gradient, the calculation is direct and simple. However, an iterative algorithm is needed to solve the inverse responses, in the case of stress and mixed control. For example, the mixed prescribed components of stress and the deformation gradient are given: (11,12,13,21,22,23,31,32,33), where A denotes the set of indices of components of the prescribed deformation gradient, while B denotes those of the prescribed stress. The corresponding components of the stress and the deformation gradient, F I J and P ij , are unknowns to be solved. It is noted that nine components in different positions are needed in order to obtain unique responses to the problem. In fact, only the unknown components of the deformation gradient F I J are to be solved in the iteration algorithm. Once all of the components of the deformation gradient are known, the components of stress can be easily calculated with Equation (1) or Equation (2).
Analogously to the algorithm for stress and mixed control in the Galerkin-FFT methods [25,26], we define a mask operator N: where δ is the Kronecker delta. Thus, it can project an arbitrary tensor into its masked counterpart, such that only the components in the position corresponding to the stress control are not zeros: Suppose that the components of the n-th deformation gradient are prescribed as follows: The n-th deformation gradient in the loading path can be separated into two parts: where F is an arbitrary deformation gradient increment. If the n-th deformation gradient F (n) and the stress P (n) are known, the unknown components of the deformation gradient and the stress in Step n + 1, where an iterative algorithm is applied, can be determined. The iteration equation is established by linearizing the expression of the stress with Equations (1) and (8): where the fourth-order tensor K denotes the consistent tangent, K = ∂P/∂F. The admissible condition for the stress is that the components should be equal to those prescribed values: Substituting Equation (9) into Equation (10), we arrive at the system of linear equations to be solved in the procedure of iterations: where the increment of an arbitrary deformation gradient δ F (k) is the only unknown.

Implementation of the Single-Point Calculator
The procedure of the single-point calculator is described as follows: 1. Divide the loading into N segments. Thus, the prescribed deformation gradient and stress increase in N steps accordingly.

2.
Calculate the deformation gradient and the stress with the known tensors in the last step and the prescribed values in this step. The initial deformation gradient is set as F (0) = I, and the initial stress can be set as P (0) = 0.

3.
Perform the Newton iteration within one step, as shown in Algorithm 1. In the (k + 1)-th iteration, the increment of the deformation gradient can be solved from Equation (11). The deformation gradient is then updated: The stress is updated: P (k+1) = P(F (k+1) ). The iteration goes on until the acceptable deformation gradient and stress are found.
The algorithm of the single-point calculator is implemented in Python. The coded Python program only uses scientific calculation packages Numpy and Scipy. The codes of the single-point calculator are available online (https://gitee.com/withmc/spc, accessed on 19 septembre 2022). The system of linear Equations (Equation (11)) of the singlepoint calculator is solved with Krylov solvers, such as CG or MinRes from Scipy, which are also the core solvers for the Galerkin-FFT methods. Even if the consistent tangent in Equation (11) is singular, the MinRes solver is able to find relevant solutions [28]. Before and during the Newton iteration of finding increments of the deformation gradient, Step 4 in Algorithm 1, the studied constitutive model is called to calculate the stress P and the tangent stiffness K. It is noted that the iteration is not necessary when the problem is only strain control, as the components of the mask operator N are zeros.
Algorithm 1 Iteration in the stress and mixed control single-point calculator.
(4) Calculate P and K with the constitutive model. end while (5) Accept the converged results: F (n+1) = F, P (n+1) = P end procedure Note: 1 and 2 are two positive values chosen as approximations to zero.

Applications of the Single-Point Calculator
In this section, two examples of applications of the single-point calculator are introduced. The first example illustrates the calculator's capability of dealing with a complicated constitutive model for incompressible and anisotropic single crystal plastic solids. The second example illustrates that the Python coded calculator is able to build machine learning constitutive models directly and simply.

Application in Tests of a Single Crystal Plasticity Model
As the first example, a recently proposed advanced model for single crystal plasticity [18] is tested with the single-point calculator.

Constitutive Model
A short review of the model is presented first. A single crystal plasticity model is pathdependent, so its constitutive model is in the form of Equation (2). When the constitutive model is called, the deformation gradient F in the current iteration, the deformation gradient F 0 and the internal variables α 0 in the last increment are given, in order to calculate the current stress P, the consistent tangent K, and the updated internal variables. For this single crystal plasticity model, the internal variables are the Cauchy stress σ, the total accumulated slip Γ, the rate of slip in each slip systemγ (α) , the critical resolved shear stress τ (α) c , and the orientation matrix Q. The index α takes values from 1 to the number of slip systems N s .
Given the current and the last deformation gradients, we can calculate the velocity gradient L with Equation (4). The symmetric and the skew-symmetric part of the velocity gradient are used in constructing the model. They are called the strain rate and the spin, respectively: An increment of rotation can be calculated with the spin [35] as follows: A trial orientation matrix Q trial can then be updated from the last orientation matrix Q 0 with the increment of rotation: The initial orientation matrix is determined by the orientation of the crystal before the deformation. In this single crystal plasticity model, the stress and the related internal variables are updated in the co-rotational frame or the lattice frame, so that the lattice tensors are constants during the deformation. Thus, the components of the Cauchy stress, the strain rate, and the spin need to be rotated as follows: A return mapping algorithm is then used to determine the updated stress in the lattice frame. First, a trial stress is obtained with an elastic predictor: where the fourth order tensor C denotes the constant elastic modulus with cubic symmetry in the lattice frame, characterized by only three independent material parameters C 11 , C 12 , and C 44 . If the trial stress is inside the yield surface f ( σ) ≤ 0, then the material is elastic and the stress is the same as the trial stress. Otherwise, the deformation is plastic and the return mapping should be employed to calculate the stress. The yield surface is the key function of the plasticity model, and this time the Holmedal model [36] is chosen: where n is a material parameter. In Equation (19), the tensor P (α) is the symmetric part of the Schmid factor in the α slip system, and it is constant in the lattice frame. When plasticity, it is to solve the following system of nonlinear equations to obtain the values of σ,λ, and Γ: In Equation (20b), the function ϕ depends on Γ because the critical resolved shear stress in Equation (19) is a function of Γ. The Voce equation is chosen as the hardening law: where τ ini , ∆τ sat(α) , and ∆γ sat(α) are material parameters. In Equation (20b), the rate of the accumulated slip is calculated as follows: where the rate of slip is calculated as follows: The system of nonlinear equations is solved with Newton iterations including a line-search algorithm [18,37]. After the values of σ,λ, and Γ are obtained, the consistent algorithmic modulus C alg can be calculated in the lattice frame. The stress and the consistent algorithmic modulus should then be rotated back into the global frame by the updated orientation matrix Q: The orientation matrix is updated from the last value Q 0 as follows: With an assumption of small elastic deformation, the co-rotated spin can be approximated as follows: where Ω (α) is the skew-symmetric part of the Schmid tensor in the lattice frame. At last, the first Piola-Kirchhoff stress is obtained with Equation (33). The tangent stiffness K is calculated as done by Lucarini and Segurado [38], and the formulation for the components is

Calculation and Results
This stable but complicated constitutive model was originally coded in Fortran, serving as a user-defined material (UMAT) subroutine for commercial FEM codes. We coded the constitutive model in Python to perform the single-point calculation, taking advantage of the Python's convenience in debugging and building. It is noteworthy that the Fortran subroutine can also be executed in the Python program, with the help of an interface already encapsulated in the package Numpy.
A uni-axial traction was simulated with the single-point calculator using this constitutive model. The plastic deformation preserves the volume, while the elastic deformation is not necessarily isochoric. For metals, the elastic deformation is small, so an incompressible uni-axial dilatation is often used in the strain control single-point calculation or one cubic element calculation. The prescribed values of the deformation gradient is often given as follows [26,27] where λ > 1, and its determinant J = 1. However, in the case of anisotropic solids, the component F 22 is not necessarily identical to the component F 33 , and the determinant is not necessarily 1. Therefore, a mixed control of the components of the stress and the deformation [25,39] should be used, as follows: With the material parameters in Table 1, the loading was simulated with the singlepoint calculator. After calculation, the Python program of the single-point calculator output the stress components and the deformation gradient, as well as all the internal variables. The components F 22 , F 33 , and the determinant J are depicted in Figure 1. From Figure 1, we can see that the determinant of the deformation gradient stays approximately at the value of 1, implying that the volume is nearly unchanged during deformations. The values of the components F 22 and F 33 decrease as the loading is gradually completed. A great difference between these two components F 22 and F 33 should not be ignored. This difference results from the anisotropic nature of the crystalline lattice, and this cannot be observed if the uni-axial dilatation as Equation (29) is applied in the strain control. This difference depends on the lattice orientation. Thus, 1000 randomly generated groups of Euler angles were investigated. The final values of the ratio of components F 22 to F 33 as functions of Euler angles are depicted in Figure 2. In Figure 2, we can observe the distributions of the values of the ratio of components F 22 to F 33 , which is a reasonable result: as the solid described by this single crystal plasticity model is nearly isochoric, it is inferred that F 22 F 33 ≈ 1, but F 22 = F 33 under uni-axial traction for a given lattice orientation. This numerical phenomenon can be obtained easily with the proposed single-point calculator.

Performance of the Calculator
The problem to solve for the performance of the single crystal plasticity model is highly nonlinear, as shown in Figure 1. The number of increment steps needs to be sufficient to avoid the issue of convergence. The minimum number of increments needed to converge for the numerical example presented in Section 3.1.2 was found to be 605 through the trial-and-error method. The performances of the single-point calculator was studied and compared with different total numbers of increments: 605, 1000, 2000, 5000, and 10,000. Table 1. Material parameters used in the uni-axial traction simulation with the single crystal plasticity model.

Cubic Elasticity Modulus
Hardening Parameters Euler Angles Four dimensions of the performance are compared and shown in Figure 3. The total time needed to finish the simulation (as in Figure 3a) was recorded in a desktop computer with a 10-core CPU of 3.7 GHz. The total number of iterations for the MinRes solver and that of the Newtons were also recorded. The average numbers of iterations are calculated as the total number divided by the corresponding number of increments, shown in Figure 3b,c for the MinRes solver iterations and the Newton iterations, respectively. It is worth noting that the Newton iteration is used to solve nonlinear equations at one increment, while iterations for the Krylov solvers are needed to solve the system of linear equations, as shown in Equation (11), at each step in the Newton iteration. The last dimension is the relative difference, shown in Figure 3d. To calculate the differences, the deformation gradient and the first Piola-Kirchhoff stress calculated in 10,000 increments are chosen as the reference values, denoted as F REF and P REF , respectively. The relative differences for each number of increments are calculated as follows: where || · || denotes the norm. The linearization the nonlinear expression of the stress (as in Equation (9)) implies a small increment of the deformation gradient. With a larger number of increments, the increment of the deformation gradient is smaller, so the results are more reliable, as shown in Figure 3d. In addition, smaller increments lead to a lower average number of iterations at one increment, as shown in Figure 3b,c. However, a higher computation cost is needed to achieve a better result, as shown in Figure 3d. In this case, the total number of increments of 2000 is recommended, as it can achieve an accuracy of about 99.95% of that with 10,000 increments, in less than one quarter of the time. In practical uses of the single-point calculator, an automatic stepping algorithm can be used to determine the magnitude of the increments at each step, to improve the efficiency of the calculator. The choice of the Krylov solvers is important for the proposed single-point calculator. In Equation (11), the corresponding stiffness matrix of the system of linear equations, N : K : N, is possibly singular. For example, in the case of the mixed control as in Equation (29), the matrix representation of the stiffness is with a vector representation of the increment of the deformation gradient as The bottom right non-zero components can be reformed into a new 8 × 8 matrix K minor . If the matrix K minor is non-singular, Equation (11) can be solved with its inverse, which is the elimination method. There are two flaws. First, the size of the non-zero minor matrix K minor depends on the stress control (N controlled components of stress lead to an N × N matrix). A problem-oriented treatment of the stiffness matrix without the generality of the solver is needed. Second, the minor matrix is not guaranteed to be non-singular. For example, at the beginning of the increment, the 8 × 8 matrix is singular as Despite the singularity, the system of Equations (Equation (11)) can be solved by Krylov solvers, such as CG or MinRes. The performances of CG and MinRes solvers are compared, as shown in Figure 4, from which we can see the advantage of the MinRes solver over the CG solver. In addition, an unexpected behavior of the CG solver can be observed that it does not converge when the total number of increments is superior to 4000. For the CG solver, N : K : N needs to transform to a tridiagonal matrix T by the Lanczos process. Then the Cholesky decomposition to T is applied. However, it is not always successful when N : K : N is singular. In this case, the CG solver is unstable. Different from CG, the MinRes solver finds approximated solutions to Equation (11) by minimizing the residual: A QR-decomposition of T through Givens rotations is performed, so that the residual decreases monotonically without divergence and stagnation, when N : K : N is singular [40]. With numerical experiments, it is found that the MinRes solver can guarantee the convergence of the Newton iterations, as long as the number of increment steps is sufficient. Thus, the MinRes solver is chosen as the solver of the system of linear Equations (Equation (11)) in the single-point calculator, due to its • Unified form of Equation (11) for all kinds of the stress control; • Ability to deal with the possible singular stiffness tangent; • Robustness, comparing to the CG solver.

Application in Tests of an Artificial Neural Network Elasticity Model
As the second example, constitutive models developed with machine learning techniques for elasticity are tested with the single-point calculator. A black-box Artificial Neural Network (ANN) model [4] and a white-box Tensorial Sparse Symbolic Regressed (TSSR) model [13] are compared. A "black-box model" refers to the model relating the system inputs to the outputs with implicit calculations, whose arithmetic mechanism cannot be seen explicitly. On the contrary, a "white-box model" has an explicit arithmetic mechanism. These models are designed to model the macroscopic behavior of the hierarchical materials, fitting data generated from microscopic calculations.

Data-Driven Constitutive Models
First, the calculation for generating stress-strain data is introduced. A particle reinforced composite (PRC) is considered, whose structure can be seen in Figure 5. The matrix phase (95.53% in volume) and the hard particle phase (4.47% in volume) can be characterized by the Eulerian linear relation between the Cauchy stress and the Almansi-Euler strain locally at each voxel: where µ and λ are material parameters given in Table 2. The particles are distributed randomly and uniformly in the matrix. Fifty-seven different paths of deformation of the PRC were calculated with the displacement-based-FFT (DBFFT) algorithm [27], so more than 10,000 pairs of stress-strain data were generated for establishing macroscopic data-driven models.  The ANN model is used to establish an implicit mechanism to predict the components of the stress tensor with those of the strain tensor. The Cauchy stress and the left Cauchy-Green strain were chosen as output and input, respectively. As both the chosen stress and strain are symmetric, there are only six input nodes and six output nodes. Trained with the normalized data from DBFFT calculations, a multi-layer, fully connected ANN constitutive model with three hidden layers was established. The activation function is the hyperbolic tangent function. The number of nodes in each hidden layer was tuned as 30, 30, and 24 from input to output. This established an implicit relation model between the macroscopic Cauchy stress and the macroscopic left Cauchy-Green strain, as follows: The Cauchy stress σ can be related to the first Piola-Kirchhoff stress as where J is the determinant of the deformation gradient F, J = det(F). The left Cauchy-Green deformation tensor B is calculated from the deformation gradient as For comparison, a TSSR model was also obtained with the same data. Even though TSSR is a data-driven method, it established an explicit expression of the relation between the stress and strain tensors: where I 1 and I 2 are the first and second invariants of B, respectively. These two models were coded in Python and involved with the single-point calculator. For the ANN model, the offline training process was performed with the help of a Python machine learning package Scikit-learn [41]. The well established ANN model can be saved in a file after training. The constitutive model function in the single-point calculator then loads the ANN file to incorporate the ANN model, to predict the stress responses with given deformation gradients. As both the machine learning package and the calculator are Python-based, there were no additional efforts made to link the offline training and the calculating.

Calculation and Results
In the case of stress or mixed control, the consistent tangent K is needed, as shown in Equation (11). For the TSSR model, the expression is already explicit, so the analytical expression of the tangent is easy to obtain. For multi-layer perceptrons, a kind of ANN with relatively simple structures, a consistent tangent can be constructed with the learned parameters, the layer structure information, and the derivatives of the activation functions [42,43]. However, in the general scope, it is hard to calculate the exact consistent tangent of an implicit and complicated network. A surrogate tangent can be learned with another ANN [44,45] from the data of the deformation gradient and the tangent calculated with DBFFT. One difficulty of this approach is that the tangent data is not always available. In addition, learning another ANN for the tangent is possibly very time-consuming. In our simple single-point calculator, a derivative-free approach is available with the central difference (CD) approximation of the derivatives [33,34]. The component K ijmn of the consistent tangent can be approximated as K ijmn = P ij ([F 11 , F 12 , ..., F mn + e, ..., F 33 ]) − P ij ([F 11 , F 12 , ..., F mn − e, ..., where e is a small scalar that is empirically taken as e = 0.01/N (where N is the number of loading segment). With the help of the iterative Krylov solver, the single-point calculator can converge in the problem of stress and mixed control in using the CD approximated consistent tangent. A mixed control was simulated successfully with the single-point calculator using these two data-driven constitutive models and the CD approximated consistent tangents. Components of the prescribed deformation gradient and the stress are given as follows: The calculated components P 11 and F 22 are depicted in Figure 6. To evaluate the accuracy of the two surrogate models, the results calculated with DBFFT are also shown in Figure 6 as the reference values. Figure 6 shows that the trained ANN model predicts the increments of the stress components similarly to the TSSR model. Differences still exist between the stress results calculated with these two models, especially when the deformation is large. As for the deformation gradient, the differences from the reference are more significant for the ANN model than the TSSR model, from which we can argue that the established TSSR has a greater generalization ability than the ANN model. Nevertheless, calculations with these two models were easily made with the single-point calculator, facilitating fast tests, comparisons, and evaluations of these data-driven models.

Conclusions
In this work, we formulated a framework of single-point calculation for stress and mixed control to build, validate, and calculate with constitutive models. A projection operator was defined to mask the components of the deformation gradient and the stress, unifying the expressions of the equations for stress and mixed control boundary conditions. The Krylov solvers were chosen to solve the problem within the Newton iteration, as in the Galerkin-FFT methods. Python is easy to code and debug, and it is the most preferred language for well-known machine learning (e.g., Scikit-learn) or deep learning (e.g., TensorFlow or PyTorch) frameworks. Thus, the proposed single-point calculator and its demonstration constitutive models were coded in Python. Two examples of applications of the single-point calculator were then considered: a constitutive model with sophisticated mathematical structures for anisotropic single crystal plasticity and an artificial neural network constitutive model for elasticity. For implicit models whose analytical consistent tangents are hard to calculate, an alternative derivative-free method is available in the single-point calculator, with the central difference approximation method. It is illustrated that the proposed framework is suitable for developing and calculating with novel constitutive models, with structures either of a complicated explicit mathematical formula or of implicit multi-layer neural networks, for solid materials under large deformations, thus improving understanding and promoting the development of advanced materials.