Crossing Scales: Data-Driven Determination of the Micro-scale Behavior of Polymers From Non-homogeneous Tests at the Continuum-Scale

We propose an efficient method to determine the micro-structural entropic behavior of polymer chains directly from a sufficiently rich non-homogeneous experiment at the continuum scale. The procedure is developed in 2 stages: First, a Macro-Micro-Macro approach; second, a finite element method. Thus, we no longer require the typical stress-strain curves from standard homogeneous tests, but we use instead the applied/reaction forces and the displacement field obtained, for example, from Digital Image Correlation. The approach is based on the P-spline local approximation of the constituents behavior at the micro-scale (a priori unknown). The sought spline vertices determining the polymer behavior are first pushed up from the micro-scale to the integration point of the finite element, and then from the integration point to the element forces. The polymer chain behavior is then obtained immediately by solving a linear system of equations which results from a least squares minimization error, resulting in an inverse problem which crosses material scales. The result is physically interpretable and directly linked to the micro-structure of the material, and the resulting polymer behavior may be employed in any other finite element simulation. We give some demonstrative examples (academic and from actual polymers) in which we demonstrate that we are capable of recovering “unknown” analytical models and spline-based constitutive behavior previously obtained from homogeneous tests.


INTRODUCTION
Modern applications and the easiness of 3D printing of polymers even at the microscale (e.g., via dual-photon polymerization), have renewed the interest in large deformation modeling of these entropic materials. Polymeric materials can now be found in a wide range of biomedical applications (stents, sutures, spinal cages, soft tissue implants, and tissue engineering scaffolds, … ), see Bergström and Hayman (2016). Even most human soft biological tissues, which are made of a matrix (elastine, proteoglycans) plus fibers (collagen), withstand large reversible deformations within the physiological range, and therefore, use hyperelasticity as ground for more complex aspects; Chagnon et al. (2015), Chagnon et al. (2017). As the simplest procedure to guarantee true elasticity (reversible, non-dissipative processes) the cornerstone in hyperelasticity is the free energy function, the state function, from which the stresses are uniquely derived from the strains (or vice-versa), regardless of path. Since the 3D strain energy function cannot be measured directly, the classical approach in constitutive modeling establishes a predefined form for the free energy. This function typically contains some parameters that are adjusted according to the experimentally (stress-strain) observed material behavior. Although it is relatively simple to tune model parameters to predict (up to a desirable precision) a single experimental curve, determining the parameters that produce accurate results for different modes of deformation is not trivial, as it is apparent from the unaccountable number of hyperelastic models available; Volokh (2016). Theoretically, if the proposed model is correct, this set of parameters should exist and although determined for specific tests they should predict well other modes of deformation. In practice, when parameters are obtained from a single experimental curve, they fail to generalize to other deformation states; this is the reason why in practice multiple tests are recommended to determine the parameters of the free energy function (Marckmann and Verron, 2006, 3, p.12). Using multiple tests to calibrate the parameters alleviates the deviations of the model for other modes of deformation (at the cost of accuracy for a given test), but at the same time it raises the question of whether the proposed form for the free energy really captures the physical phenomena behind experimental data or this assumed form is just a complex interpolation scheme that adapts its parameters to fit the curves used during calibration. We remark that if the physics behind were accurately represented, a single curve should be sufficient to capture the general multiaxial behavior of isotropic, incompressible polymers under reversible deformations.
This kind of problems has encouraged many researchers to pursue different approaches. One of them is the modelfree data-driven computing paradigm. In this approach, basic conservation laws and essential constraints are satisfied but the constitutive laws are eliminated in the benefit of data; Kirchdoerfer and Ortiz (2016), Kirchdoerfer and Ortiz (2018), Eggersmann et al. (2019), Ibañez et al. (2017), Ibañez et al. (2018). Regarding the leading role of data for some of these references, works that address the efficient handling of data have also been published, see Zheng et al. (2020), Korzeniowski and Weinberg (2021). On the other hand, other approaches attempt to surrogate the constitutive law with an input-output relation through Artificial neural networks (ANNs), Nguyen-Thanh et al. (2020), Liu et al. (2020). Both approaches (modelfree data-driven and surrogate-like ones) show promising results, however, the predictive capability of models that just rely on data is strongly dependent on the amount and quality of data being employed. In addition, since there is no expression for the free energy function most of those models are very difficult to interpret from a physical standpoint. It seems clear that an approach solely based on data might not be the best option for this kind of problems (the more the model needs to learn, the more data is required). This need for introducing physics information in full data-driven models has led to other works based on Physics-informed neural networks (PINNs) Liu et al. (2020) and on thermodynamically consistent datadriven approaches; González et al. (2018). Still, the increased generality of those approaches increase the amount of data required when compared to classical constitutive modeling techniques and their interpretability is much less direct. To summarize, an optimal approach for the constitutive modeling of polymers should: 1) include information about the physical equations without assuming a fixed given form for the free energy function; 2) use data to complement what we know about the physical phenomena and fill in the gaps in our knowledge, but without using more data than actually needed; 3) interpretability is also very important because understanding the solution and being able to identify its physical meaning avoids many pitfalls allowing us to search for the answer within a smaller solution space and identify spurious solutions. Interpretability also facilitates the imposition of desired (physics-based) requirements to the sought solution, for example, smoothness, monotonic increase or decrease, isotropy, etc.
With all those requirements in mind we developed the WYPiWYG (What-You-Prescribe is What-You -Get) approach to constitutive modelling, Latorre and Montáns (2013), based on some seminal ideas from the Sussman-Bathe model for isotropic, incompressible materials, Sussman and Bathe (2009). The WYPiWYG approach determines the free energy function or its contributions, but in contrast to classical phenomenological models, which presume a form for the energy function and fit the model parameters to the experimental data, our approach starts with some basic fundamental assumptions about the material behavior (isotropy/anisotropy, Valanis-Landel decomposition, invariant-based contributions to the energy function) and then obtains numerically the constitutive equation from equilibrium using a local approximation scheme based on splines. This local approximation philosophy is similar to the way shape functions in finite elements interpolate the displacement field, instead of computing coefficients of predefined analytical functions as in the Navier and Rayleigh methods. The generality of this approach is demonstrated on the models elaborated for anisotropic materials, Latorre and Montáns (2014), auxetic materials Crespo and Montáns (2018) and models for the active and passive response of skeletal muscle, Moreno et al. (2020). So far phenomenological WYPiWYG hyperelasticity circumvents the need to prescribe the shape of the energy function while maintaining the same model interpretability that the phenomenological models have. However, the amount of data required to characterize the behavior of polymeric materials is similar to the amount of data required by other parametric phenomenological models like the Ogden Model, Ogden (1972).
Other alternative to phenomenological models are those based on the micro-structure which employ additional information about the structure of polymers to get better predictive capabilities with fewer data. Most micro-structural models assume that the polymer is fully entropic and thus all the work employed in its deformation directly translates into a variation of its entropy. This physical insight has been exploited by researches, leading to well known models as the Neo-Hookean model, Flory and Rehner (1943); Gent (1989); Treloar (1975), and the 8chain model, Arruda and Boyce (1993). The expressions of those models depend just on the first invariant of the Green-Cauchy tensor, I C 1 and some material parameters, but in contrast to the phenomenological ones, the material parameters are linked to the micro-structure resulting in some additional physical insight. Although they were conceived to be characterized from a simple extension test, the results on other modes of deformation are not satisfactory, even with the additional information about the microstructural behavior. This fall from expectations added up to the conclusions of Mooney and its Mooney plots Mooney (1940), which showed that I C 1 was not the right (or at least the unique) variable to describe the polymer behavior due to the controversial slope C 2 that consistently appeared on simple extension experimental data.
The extension of the WYPiWYG approach to microstructural modelling with the aim of overcoming those dificulties resulted in a Macro-Micro-Macro (MMM) approach to obtain the polymer constituents behavior directly from experimental data with no assumptions about its analytical form or parameters to calibrate; Amores et al. (2020). Just some basic assumptions were made: 1) homogenization of the chains free energy to obtain the free energy of the continuum Ψ (λ 1 , λ 2 , λ 3 ) = ∫ S ψ ch (λ ch ) dS/S, and 2) the computation of the micro-stretch variable is obtained from the continuum stretch tensor, λ ch = r ⋅ U ⋅ r (a non-affine measure of deformation in agreement with the lack of relevant contribution of chains orientation change in the entropy reduction). The MMM approach predicts well any general deformation mode requiring just a single experimental curve to characterize the chain behavior, see (Amores et al., 2020, Figure 3). Since a similar Data-Driven MMM framework using the affine deformation measure (λ C ch ) 2 = r ⋅ C ⋅ r was not able to offer the same results, in Amores et al. (2021) we questioned the affine micro stretch assumption from theoretical grounds, which seems to be the most popular in micro-structural models; Treloar (1975), Arruda and Boyce (1993), Alastrué et al. (2009), Sáez et al. (2011, Khiêm and Itskov (2016). In that same work the non-affine measure of deformation λ ch = r ⋅ U ⋅ r did show to be in accordance with the "controversial" C 2 slope observed by Mooney that up to the date could not be successfully explained from the classical statistical theory.
The framework presented in Amores et al. (2020) seems to be in accordance with both the chain statistical theory and experimental results, but needs homogeneous tests to characterize the chain behavior. Hence, our work here is to pursue a more general approach by employing arbitrary continuum nonhomogeneous tests and using Digital Image Correlation (DIC), crossing scales from the continuum to the polymer constituent macromolecules.
The procedure consists of linking two stages. One is the previously introduced MMM method, and the other one is to link that method to a finite element analysis of non-homogeneous continuum problems continuum problems, see (Cite to Figure  1) outline. We assume in the latter that the non-homogeneous field of displacements (via DIC), plus the test loads (via load cell) are known, the input data could be either 1D, 2D or 3D depending on the case, but it is important to note that in 2D and 1D, it should be possible to employ reasonable assumptions to determine the principal stretches and stresses in the eliminated directions (incompressibility plus plane stress allow to determine both the stretch and the stress out of the plane just from the plane information). Then, the polymer chain behavior is modeled by P-splines, which vertices are to be determined-P-splines are penalyzed interpolating B-splines to guarantee smoothness; see Eilers and Marx (2021). That structure (the unknown vertices) are transferred to the continuum scale via integration in all the material directions and the result attached to the finite element integration point (the continuum constitutive behavior). Hence, the nodal forces of the finite element are set as a direct, explicit function of the unknown P-spline vertices of the polymer chains and the prescribed deformation gradient. By a least squares formulation, a linear system of equations is established, which allows for the immediate determination (i.e., simply solving a linear system of equations) of the P-spline vertices of the chain behavior from the macroscopic loads in the specimen and the macroscopic field of deformation. The physics equations present on the procedure include at the FEM level the compatibility equations (computing the strain quantities from the displacement field) and the equilibrium equations (null force residual), incompressible hyperelasticity with volumetricdeviatoric decoupling at the integration point level, properly including the micro-macro connections (energy homogenization and affine micro-stretch) at the chain level.
In the following sections we introduce the procedure, first using a continuum hyperelastic formulation and then the micromechanical one. We also demonstrate the applicability through an academic example and an example using the wellknown Treloar's rubber.

Detailed Procedure Description
One way to obtain the displacement-based finite element formulation is through the principle of virtual work, which for the quasi-static case reads δW = δW int − δW ext = 0. Considering a conventional FEM discretization and interpolation of the displacement field, both virtual works (internal and external) could also be expressed in terms of the internal and external nodal forces, To conclude, a weighted integration of the equilibrium equation leads to an alternative expression for the virtual work, see (Eq. 1), that if compared with δW = δu ⋅ (f int − f ext ) = 0 provides an expression for the internal and external nodal forces. The expression for the weighted integration of the equilibrium equation in the reference configuration can be expressed as: where δu and u are the virtual displacement field, and the displacement field, respectively, P, the 1st PK (First Piola-Kirchhoff) stress tensor, b the volumetric forces and t the surface forces. The symbol ⊗ represents the dyadic production, so [δu ⊗ ∇] ∶ P (u) = P ∶ ∇δu is the internal virtual work density.

Internal Force Term
As it has already been mentioned, a FEM discretization for the internal force term, ∫ Ω P ∶ ∇δudΩ, and the interpolation of the displacement field in the reference unit element, δu = ∑ n n a=1 h a ( 0 ξ) δu a leads to the expression for the nodal internal forces, see (Eq. 13). Note that h a ( 0 ξ) are the shape functions, δu a the virtual displacement vector for the node a, and that in δu a i , a = 1, …, n n and i indicates the spatial dimension (i = 1, 2 in 2D). The expressions for the internal forces are described with depth underneath, (Eq. 2) and (Eq. 3): where all the variables with subscript j are computed in the integration point j. The previous equations can be rewritten doing the sum over the DOF of the element (n dofs = n n × 2 in the 2D solid elements) instead of doing the sum over the nodes: In Eq. 3, i is the index that runs through the local degrees of freedom in the element, for the cases studied here (2D plane stress problems) i = 1, …, 2n n , δu i is the virtual displacement at the local degree of freedom i and ∇h i j is a second order tensor that projects the contribution of P j to the ith local degree of freedom. Regarding ∇h i j , for i = 2, ∇h 2 j is the projector for the second local degree of freedom of the element, this DOF corresponds to the first node of the element a = 1 and the second dimension 2, therefore, ∇h 2 j = e 2 ⊗ ∇h 1 ( 0 ξ j ). If i = 3 instead, ∇h 3 j is the projector for the third local degree of freedom of the element, which corresponds to the second node of the element a = 2 and the first dimension 1, therefore, ∇h 3 j = e 1 ⊗ ∇h 2 ( 0 ξ j ). Now Eq. 3 is rewritten in matrix form to identify the components of the internal force vector: where Since the PK1 tensor is a two-leg tensor placed in 2 configurations at the same time (material in the right and spatial in the left), it might be more suitable to rewrite the term P j ∶ ∇h i j in terms of the PK2 (Second Piola-Kirchhoff) stress tensor, S, which lies completely in material configuration, P = XS, being X = ∂ t x/∂ 0 x the deformation gradient: In order to simplify the notation we define l i j ≔ sym (X T j ∇h i j ), on the other hand the double contraction will be computed using the Voigt notation. Note that S 3 = 0 since we are considering the case of plane stress and that the components in the Voigt notation are with respect to the principal directions of deformation/stress. Then, we write in principal directions (denoted as the X ppal system of representation) To compute the principal stress components of a polymeric and quasi-incompressible material, we assume that the volumetric and deviatoric contributions of the energy can be separated (a typical assumption in quasi-incompressibility), the volume ratio, and A = 1/2 (X T X − I) the Green-Lagrange strain tensor. A typical choice for the penalization could be U (J) = κ/2 (J − 1) 2 with the bulk modulus, κ, selected such that the stress produced by a volumetric deformation grows rapidly. If an estimation for the shear modulus of the material, μ, is known, κ/μ∼10 4 , might suffice to ensure the satisfaction of the quasiincompressibility condition J = λ 1 ⋅ λ 2 ⋅ λ 3 ≈ 1. The decoupling of the energy in a volumetric and deviatoric energy results in a similar decoupling for the PK2 stress tensor, S = S v + S d . We note that S v and S d are not the volumetric and deviatoric part, respectively, of a second order tensor S (as obtained from the respective mathematical operators), but S v is the contribution to S that comes from the volumetric energy U(J), and S d is the contribution to S that comes from the deviatoric energy, where N i are the principal referential directions of deformation, W k = ∂W/∂λ d k and λ d i = λ i J −1/3 , the isochoric stretches. As we have already mentioned the test cases are plane stress and therefore, S 3 = 0 = S ∶ (N 3 ⊗ N 3 ): From the previous equation, the pressure term could be isolated and introduced in the equations for S 1 and S 2 : Once the principal stretches are known for a particular integration point, the principal PK2 stress components in that same point could be determined through evaluation of the deviatoric contribution derivatives. When the stress tensor has been computed in all the integration points of the element, the internal nodal forces for that element are computed using (Eq. 4). On the other hand, if the functions for the deviatoric contribution derivatives are unknown, a cubic P-Spline local approximation can be employed, see Amores et al. (2019), Eilers and Marx (2021). The B-splines (or its penalized version, Psplines) are one of the approaches for expressing any general function y(x) as the product of a set of known basis functions (cubic in this case), B i (x) and it's a priori unknown corresponding weights (also called vertices),v i , if the number of vertices of the B-Spline is n vert , the expression for the unknown function is a scalar product, y ( In contrast to the basis of features proposed in Flaschel et al. (2021), where a basis of preassumed functions is used, the local B-Spline basis is general with no assumptions about the possible function expressions (the only assumption would be that locally the function is at most a cubic function when cubic B-Splines are used). Additionally, the local P-Spline approximation in this work is typically performed directly on the derivatives so there is no need to compute derivatives of the approximated function. With the mentioned P-Spline for the derivative of the energy function, the principal stresses could then be represented as row is a vector defined with 1D P-Splines basis vector (B (λ d 1 ), B (λ d 2 )), see Section 3.1 and Section 3.2 for specific expressions of this vector. Therefore, depending on the approach that is used to compute the continuum principal stresses, the unknown vertices,v , could represent the derivative of the energy either on the macro-structure or in the micro-structure: and for the element local degree of freedom i which means that every component of the elemental nodal force vector is obtained as the product of a known row multiplied by the unknown vertices, or expressed in a different manner: where F e int is a known matrix, given in square brackets in Eq. 13.

External Forces Term
In case the body forces and the tractions are not zero, the external force vector for each element, f e ext has also to be computed and assembled into f ext , looking at the right hand side of (Eq. 1).
On the previous equations, the reaction forces on the boundary where displacements are imposed are accounted in the traction term, for the sake of simplicity we are going to consider that those nodal reaction forces are known, while this is not typically the case in a experimental setting, instead, the total reaction forces are known rather than the nodal forces. To deal with this fact, 2 different approaches can be followed: 1) take another artificial boundary far enough from the original one in which we can suppose that the reaction force is evenly distributed according to the Saint-Venant's principle or 2) consider two independent set of equations one for the free dofs and other for the DOF with fixed displacements, in the fixed displacement DOF the resultant of the internal forces equals the reaction force at each of the boundaries.

System of Linear Equations for the P-Spline Vertices
The problem to solve is an overdetermined linear system of equations F int ⋅v = f ext , wherev are the unknown vertices of the P-Spline that approximates the derivative of the energy function at the macro-scale (phenomenological energy function) or the derivative of the energy function of the constituents at the micro-scale. The solution process consists of finding the vertices,v , that minimizes the mean square error, MSE, in the force residual: The solution for the previous minimization problem is analytical and result in the mentioned linear system of equations: As we have already mentioned, one of the advantages of using the P-Spline-based local approximation is that although we do not assume the form of the energy function, additional requirements can be added to the solution, a typical one is smoothness, which can be translated into a penalization on the second order finite differences of the P-spline vertices, if the solution obtained is not monotonically increasing, this property could also be imposed FIGURE 1 | On the right side the FEM part (green contour), from bottom to top there are different levels, first at the solid level, the nodal displacement vector, u, the nodal internal forces, f int , and the external nodal forces, f ext , are encountered. Second, at the element level, the elemental nodal displacement vector, u e , the elemental nodal internal forces, f e int , and the elemental external nodal forces, f e ext . Finally, the integration point (at the macro-scale) characterized by the Right Green Cauchy tensor C, and the PK2 stress tensor, S, at this point. On the left hand side the Macro-Micro-Macro approach (magenta contour), the general idea would be "zooming" (hence the magnifying glass) on the integration point and consider that for polymeric materials, the chains oriented in a certain direction (representative direction) suffer a deformation described by the continuum stretches ellipsoid (upper left). The a priori unknown chain mechanical behavior, P ch (λ ch ) (bottom left), will be evaluated for each direction and the effect in each representative direction will be integrated for all the representative directions on the reference unit micro-sphere (center left) to obtain the continuum principal stress components. Finally, the Macro-Micro transition (center of the figure) dictates how the continuum energy is obtained from the constituents free energy and how the micro-structural stretches (non-affine chain stretch) are obtained from the continuum deformation measures. From a more practical perspective: 1) The Macro-Micro-Macro approach expresses the continuum principal stresses in terms of the unknown P-Spline vertices for the constituents mechanical behavior, S (v ), and 2) A classical FEM assembly process is employed to write nodal equilibrium equations from the information at the integration points, F int ⋅v = f ext . The substeps involved are: 1.1) The continuum strain, C, tensor defines the micro-structural deformations, λ ch on the representative directions; 1.2) The micro-mechanical behavior (unknown a priori) is evaluated in all the representative directions in terms of the unknown verticesv , P ch (λ ch ) = B row (x) ⋅v ; 1.3) The constituents mechanical behavior in all directions is integrated over the unit micro-sphere to obtain the continuum principal stresses at the integration point, S (v ); 2.1) The integration point stress S (v ) is integrated over the element volume to obtain the elemental nodal internal forces, f e int (v ); 2.2) The elemental nodal internal forces are assembled to obtain the nodal internal forces of the solid, f int (v ); 2.3) The unknown vertices are determined from equilibrium, f int (v ) = F int ⋅v = f ext . by iterative penalization on the vertices that do not satisfy the condition see Amores et al. (2019), Eilers and Marx (2021). With all the penalizations, the general system of equations to solve would be:v where W is a diagonal matrix that weights the relative importance of the equations in F int , D 1 and D 2 are the first and second differences matrices respectively, Ω 2 and Ω 1 , are also diagonal weight matrices for the penalizations in the second and first differences of the vertices, note that k indicates the step of the iterative penalization for monotonic smoothing. When all the intervals are monotonically increasing the weighting matrix Ω k 1 = 0 and that penalization automatically disappears. In contrast to the overdetermined system of linear equations produced by piecewise approximation presented in this work, the methodology in Flaschel et al. (2021) employs a basis of functions which extend over the whole function domain, although this procedure seems similar to the local one, the global support typically leads to dense solutions which are less interpretable. For the sake of interpretability Flaschel et al. (2021) proposes a sparsity promotion with ℓ p regularization, but at the price of requiring a fixed point iteration procedure to obtain the solution and introducing p as an hyperparameter that has also to be determined in the solution process. On the other hand, the additional matrices added to the system matrix in our work ensure smoothness of the solution (making the solution more robust against noise) and allows penalization in monotonicity if required (this penalization has not been required for the examples presented in this article).

Toy Example: Recovering an Analytical Neo-Hookean With a P-Spline for the Continuum Free Energy
In this section we are going to demonstrate just the FEM part of the methodology using a P-Spline based approximation for the functions that form the deviatoric continuum free energy. The experimental data is a virtual test (FEM results) of an analytical quasi-incompressible Neo-Hookean model with parameters μ = 3.5MPa and κ = 1 × 10 5 MPa. The FEM model employed for the virtual tests is a 2D plate with a hole on plane stress conditions. The solicitation is an imposed displacement in the upper border with u = 4. The dimensions of the plate and the mesh employed in the simulation is shown in Figure 2, the elements are quadratic quadrilaterals with nine integration points. The formulation employed for plane stress in large deformations is detailed in Supplementary Appendix SA, the reader can also find the details about the Neo-Hookean material model in Supplementary Appendix SB.
Regarding the software employed, Julia programming language, Bezanson et al. (2017) has been used, in particular the package FerriteFem.jl, see Carlsson and Ekre (2021), for the FEM simulations and Amores (2022) for the P-Splines functionality.
A typical assumption made on the phenomenological approach for isotropic incompressible solids is the Valanis Landel decomposition, . With that assumption, the expressions for the principal PK2 stresses can be obtained in terms of an unknown function ω ′ (x): Once ω ′ (x) is obtained, the principal stresses are determined for any deformation state under the plane stress hypothesis or the contribution to the principal stresses coming from the deviatoric energy for any general case (not in plane stress), see (Eq. 8), this is all that can be determined from a constitutive point of view because for an incompressible solid, the pressure does not come from the constitutive equation but from equilibrium conditions. Using P-splines, the unknown function can be written in terms of some known basis functions and some unknown vertices in a similar manner that it is done in FEM formulation with the displacement field:  With the previous expression, the principal stresses could be written as: [ Referring again to Eq. 13, internal nodal forces of the element can be written as the product of a known matrix by a vector of unknown vertices: Doing the assembly of the internal nodal forces of the elements in the mesh: From the overdetermined linear system of equations F int ⋅ω ′ = f ext we can solve forω ′ : The number of vertices has to be enough to capture the complexity of the curve, typically n vert = 14 suffice, but additional vertices can be added. If more and more vertices are added, the number of equations required to determine them increases and the problem can become ill-conditioned, this is solved by the smoothing term which adds the information of smooth transition between vertices and links the vertices to its neighbours. In the homogeneous case, just a single equation is obtained for every load/displacement step, therefore, in order to obtain information for the function on the considered domain (from the minimum principal stretch to the maximum principal stretch on the simulation), it would be necessary to sweep a whole range of load/displacement steps. On the other hand, for the non-homogeneous case, in principle, it would be possible to employ just a single load/displacement step if the step under consideration is rich enough (in this case, just the last step, u = 4 was used). In case that additional information is needed to determine the function on the considered range, it is also possible to add more steps between u = 0 and u = 4. Since the initial solution was directly monotonically increasing, Ω 1 = 0 and it was not necessary to follow an iterative process, ω ′ was directly obtained from the simple initial system of equations: With the vertices obtained, the P-Spline could be reconstructed and compared to the original one, see Figure 3.
Since the function ω ′ (x) has been reverse-engineered exactly, the values for the free energy partial derivatives for any state of deformation can be determined, and will be exactly equal to those obtained with the original Neo-Hookean model. If desired, the FEM simulations can be run again with the reversed-engineered derivatives of the energy function and the σ VM plot could be compared with the one obtained for the original analytical Neo-Hookean, see Figure 4. As it can be seen in Figure 4 the plots are exactly the same as it might be expected from having obtained the exact same ω′ as the one from the original Neo-Hookean.
In this particular example, the matrix of the system A sys ⋅ω ′ = b sys , presents a null eigenvalue with its corresponding associated eigenvector, i.e., there is a certain subspace of possible solutionsω ′ 0 (containing just the direction of the eigenvector associated to the null eigenvalue) such that added to the solution,ω ′ , produce null effect on the independent term, b sys and thereforeω ′ + αω ′ 0 , for an arbitrary scalar α, is also a solution of the system. Looking at Eq. 20, to not produce any effect, such function has to satisfy the condition Regarding the functions that satisfy this condition, there is a specific form of function that complies (which at the same time must be smooth due to the regularization) The only way two functions of different independent variables can be equal ∀ λ d Since the shape of the eigenvector associated to the null eigenvalue exactly matches the shape of the function, ω ′ (λ) = β 0 /λ, we would like to understand physically the origin of this null contribution to the system. Looking at Eq. 8, we can see that the pressure term is of the form p (λ 1 , λ 2 , λ 3 ) /λ 2 i , the same form of the term resulting from introducing an additional β 0 /λ d i to the function ω ′ (λ d i ). From a more general perspective we can see that this indetermination comes from the split of the energy into a volumetric and a deviatoric contribution: (29) where C = X T X is the right Cauchy-Green tensor and ℙ is the fourth order projector tensor, see Supplementary Appendix SB.
Since ℙ ∶ (αC −1 ) = 0, from (Eq. 29) it can be seen that there are multiple S |d = ∂W/∂C d that produce the same S d and the undetermined part has the form of a volumetric-like PK2 stress tensor. Any contribution of this form has no effect on the solution since it is going to be finally eliminated by the projector ℙ.

Micro-Mechanical Approach Based on the Non-Affine Deformation Chain Model
In this section the complete methodology described in Section 2 (Macro-Micro-Macro first, FEM second) is demonstrated. In contrast to the procedure described in Section 3.1, here the P-Spline will approximate the derivative of the chain energy function P ch (λ ch ) = dψ ch /dλ ch (mechanical behavior of the material constituents at the micro-scale) rather than the continuum contributions to the energy. The experimental data is obtained from a virtual (FEM simulation) using the analytical structure-based material model in Amores et al. (2020). To be more precise, the P ch (λ ch ) that will be employed to perform the virtual test is the one obtained for the Treloar test data for unfilled rubber, Treloar (1944), that function is displayed in (Amores et al., 2020, Figure B.1). The FEM model employed for the virtual tests is again the 2D plate with a hole on plane stress conditions with an imposed displacement in the upper border (u = 8). All the required dimensions and the mesh employed in the simulation is shown in Figure 2. The reader can find further details about the non-affine deformation chain material model in Supplementary Appendix SC and Amores et al. (2020), a discussion about the suitability of the non-affine stretch employed in the model (free-fluctuating network assumption) is presented in Amores et al. (2021). In Amores et al. (2020) a way was established to compute the strain energy function of the continuum with an homogenization of the micro-structural chain free energy function in polymeric like materials under the assumption of incompressibility (λ 1 ⋅ λ 2 ⋅ λ 3 = 1): For FEM simulations even pure incompressible solids are simulated with the quasi-incompressible material framework: in which U (J) is a penalization function that leads to incompressibility for κ → ∞ like U (J) = 1 2 k (J − 1) 2 and where W (λ d 1 , λ d 2 , λ d 3 ) is the same as the one used in pure incompressibility but with the deviatoric stretches instead of the total ones: The unknown function that will be approximated using P-Splines is P ch (λ d ch ), as described above, the P-Splines representation expresses the unknown function in terms of some known basis vector multiplied by a vector of unknown vertices: Introducing the expression for P ch (x) in (Eq. 32) and then all the derivatives of the form W k in (Eq. 11) a expression for the principal stresses in terms of the vertices that define the Polymeric-chain response function can be obtained: Note that the previous integral in the microsphere is computed by a numerical quadrature, ∫ S f (r) dS = ∑ n qS j=1 f (r) w S j an that although λ d ch depends on λ d 1 , λ d 2 and λ d 3 , since λ d 1 ⋅ λ d 2 ⋅ λ d 3 = 1, there are just 2 independent variables. Again a procedure similar to the one in the toy example, Section 3.1, is followed: Looking at the previous matrix equation and using again (Eq. 13), it is straightforward to write the internal vector force for an element as where F e int is a matrix. Now in the same way that the components that the components of f e int can be assembled to obtain the global internal forces vector, f int . The rows of F e int can be assembled to obtain a global internal force matrix F int .
At that point it is important to note that if compared to the procedure followed in Section 3.1, now the internal force term is linked with the unknown vertices that correspond to the constituents behavior at the micro-scale. With that approach we show that scales can be crossed and information in one scale can be pushed up to other scales, that of course taking into account that the macro to micro connection has already been established.
From the overdetermined linear system of equations F int ⋅P ch = f ext a solution can be obtained forP ch : The discussion about the number of vertices and steps of load/displacement required in Section 3.1 is also applicable here. In this particular case, just the last step of deformation (u = 8) was used. Again, in case that additional information is needed to determine the function on the considered range, it is also possible to add more steps between u = 0 and u = 8. Since the initial solution was directly monotonically increasing, Ω 1 = 0 and it was not necessary to follow an iterative process,P ch was directly obtained from the simple initial system of equations: With the vertices obtained, the P-Spline could be reconstructed and compared it to the original one, see Figure 5.
As it is shown in Figure 5, the P ch (λ d ch ) extracted from (Amores et al., 2020, Figure B.1) and the one obtained through the methodology described in the paper are practically identical. Looking closely, there are some small deviations close to λ d ch = 1 FIGURE 6 | Plot comparing the results of σ VM for the reference micro-mechanical model in Amores et al. (2020) that uses (Amores et al., 2020, Figure B.1) and the reverse-engineered one described in the article. (A) σ VM using the P ch (λ ch ) obtained for unfilled rubber, Treloar (1944) extracted from (Amores et al., 2020, Figure B.1). (B) σ VM using the P ch (λ ch ) reverse-engineered from the non-homogeneous virtual test using the methodology described in this article. on the compression range. As it has already been mentioned, the origin of those deviations is the eigenvector linked to an almost zero eigenvalue on the reduced matrix of the system, this matter has been carefully justified at the end of Section 3.1. With the obtained P ch (λ d ch ) the partial derivatives of the free energy can be determined for any state of deformation, with them, it is possible to re-run the FEM simulations and compare for example, the σ VM plot. Figure 6 compares the σ VM with P ch (λ d ch ) extracted from (Amores et al., 2020, Figure B.1) and the P ch (λ d ch ) with the methodology described in this paper. As we can see, although the functions are slightly different, that discrepancy has no effect on the final stresses, this is so because the additional spurious term that appears in the solution corresponds to a pressure term (in the cauchy stress sense) and therefore, it will not appear in the final computed stress due to the projector ℙ.

CONCLUSION
We have proposed a numerical method for determining both the continuum free energy and the polymer macromolecules behavior from arbitrary non-homogeneous DIC-based tests at the continuum scale. The procedure consists of a combination of our Macro-Micro-Macro approach and a Finite Element model. As a novel contribution of this approach, we show that by crossing scales transferring the microscale unknowns to the finite element formulation it is possible to determine the mechanical chain behavior from non-homogeneous experiments at the continuum scale by simply solving a linear system of equations (i.e., in an even more efficient manner than the subsequent simulations of the polymer behavior). Another key aspect is that Penalized B-splines (P-splines) preserve the general form of the energy function while retaining sufficient tools for enforcing specific desired conditions on the sought functions. The methodology at hand recovers the analytical free energies used as starting point in the virtual tests. Therefore, any finite element simulation performed with the reversed-engineered energy function will provide identical results to the original material model. From the authors perspective, this new approach opens new possibilities for data-driven characterization of the micro-structure from nonhomogeneous tests at the macro-scale. Further research has to be conducted to evaluate the generalization of the approach in Amores et al. (2020) to more complex material behaviors from which the applicability of this same methodology has to be assessed.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
VA: Conceptualization, methodology, software, and writing original draft. FM: Conceptualization, methodology, funding acquisition, and supervision. EC: Methodology and supervision. FC: Conceptualization, methodology, funding acquisition, and supervision. All authors contributed to manuscript revision, read, and approved the submitted version.

FUNDING
This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Grant Agreement No. 101007815.