A component-free Lagrangian finite element formulation for large strain elastodynamics

Standard finite element formulation and implementation in solid dynamics at large strains usually relies upon and indicial-tensor Voigt notation to factorized the weighting functions and take advantage of the symmetric structure of the algebraic objects involved. In the present work, a novel component-free approach, where no reference to a basis, axes or components is made, implied or required, is adopted for the finite element formulation. Under this approach, the factorisation of the weighting function and also of the increment of the displacement field, can be performed by means of component-free operations avoiding both the use of any index notation and the subsequent reorganisation in matrix Voigt form. This new approach leads to a straightforward implementation of the formulation where only vectors and second order tensors in R3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbb {R}}^3$$\end{document} are required. The proposed formulation is as accurate as the standard Voigt based finite element method however is more efficient, concise, transparent and easy to implement.


Introduction
It is well known that equations in the field of Continuum Mechanics such as conservation of mass, conservation of momentum, constitutive relations, etc., might be formulated in a component-free ( [13,29] where no reference to a basis, axes or components is made, implied or required. It is also well recognised that such a component-free notation is a convenient and concise tool to manipulate most of the relations in Continuum Mechanics. However, particularly in the field of Computational Mechanics, it seems to be essential to refer vectors and tensors quantities to a basis and thus developing a discrete formulation based on indicial-tensor notation. Most of the variational based numerical techniques adopt this indicial-tensor notation in order to factorized the weighting functions and then the indicial-tensor notation is usually arranged in a matrix Voigt form to take advantage of the symmetric structure of the algebraic objects involved. Aside from the matrix Voigt form, there are other organisations in indicial-tensor notation like the one proposed by Gupta and Mohraz [12]. However, as explained in [14] (p. 156), Gupta and Mohraz approach [12] can not be applied in the context of non-linear analysis. On the contrary, the matrix Voigt form can be adjusted to formulate and compute non-linear analysis.
Thus, standard textbooks on finite elements for solid dynamics at finite strains [3][4][5]28,38,43] consider Voigt's notation and related algebra [42] as the appropriate convention to represent symmetric second order tensors and fourth order tensors (with major and minor symmetries) to perform a finite element formulation and implementation. The main reason argued for this choice, aside from being eligible in the context of non-linear analysis, is the fact that under Voigt's notation symmetric tensors are represented by reducing its order and therefore 'simplifying' its algebraic manipulation. In other words, Voigt's notation takes advantage of the symmetric structure of the algebraic objects involved to reduce the computational effort. Also, it is often maintained that similar algebraic objects are obtained under Voigt's notation for different scenarios such as three dimensional, plane-strain, plane-stress and axisymmetric problems, where the symmetric gradient of the shape functions B and the matrix of material elasticities D always appear [27].
Due to the paramount influence of the finite element method in Computational Mechanics field, Voigt's notation has been adopted in most of the variational based numerical techniques, such us Discontinuous Galerkin Method [2], the Material Point Method [19,20,40] or Optimal Transportation Mesh-free [17,23], where symmetric gradient of the shape functions B and the matrix of material elasticities D also appear. Voigt based finite element formulation and implementation is also widespread considered in other branches of continuum mechanics like the theory of porous media including large strain setting [16,24,25,35,36,44].
In the present work, the authors propose a completely different approach, as compared to the standard Voigt's notation and also to the indicial-tensor notation derived by Gupta and Mohraz [12], to perform the finite element formulation and implementation for isothermal elastodynamics at large strain in the initial configuration. Our principal claim is that the component-free notation, where no reference to a basis, axes or components is made, implied or required, can be further considered in the finite element formulation. Moreover, the factorisation of the weighting function can be performed by means of component-free operations, avoiding both the use of any index notation and the subsequent reorganisation in matrix Voigt form. The absence of any index notation at this point of the finite element formulation, brings a concise and transparent nodal equilibrium equations and related linearisation that are easily implemented. The original idea of this component-free approach was set by Planas and co-workers in 2012 [31] with their B free formulation, and it is further extended and clarified for large strain elastodynamics under Lagrangian description in the present work. Even though, the original designation of the component-free approach was "B free", the authors of the present manuscript consider appropriate to adjust the given name to "component-free" formulation, to further clarify that the proposed formulation is not only free from the so-called B and D matrix, but, what is further more paramount, from any basis, axes or components, being the key ingredient of the propose approach. The "B free" finite element approach has been also adopted in [1,20,21,32,39].
The main benefits of the proposed component-free notation with respect to the classical Voigt notation are as follows: 1. The discrete equilibrium equations and related linearisation in the component-free approach resemble their continuum counterparts. Therefore, this approach has a more natural interpretation than the classical Voigt notation. 2. The component-free finite element formulation and implementation, being as accurate as the classical Voigt form, is more efficient in terms of computation time. 3. It is well known that first Piola-Kirchhoff stress P and its power-conjugateḞ can be considered to develop a finite element formulation in the Lagrangian description under Voigt notation. However, these second order tensors are non-symmetric and Voigt's notation leads to larger arrays of dimension 9 instead of dimension 6. Therefore, in this circumstance, no reduction of the computational effort can be gained in Voigt form. On the contrary, in the component-free approach the computational cost is the same irrespective of the conjugate pair considered P anḋ F or S andĖ. This is possible due to the fact that no reference to a basis, axes or components is made in the proposed formulation. 4. The component-free approach exploits the fact that fourth order tensors never appear by itself in the finite element equations, but always operating on second order tensors, avoiding both the use of any index notation and the subsequent reorganisation in matrix Voigt form. Therefore, the so called B and D matrices never need to be constructed or implemented and Voigt vectors never need to be employed. In the authors' opinion this is the main reason why the component-free finite element formulation and implementation is more efficient than the classical Voigt form. 5. As no reference to a basis, axes or components is made, implied or required, the component-free approach can naturally be applied to general classes of constitutive models such as a general isotropic hyperelastic materials or even multiplicative plasticity models. 6. As no reference to basis, axes or components is made, implied or required, the proposed finite element formulation is the same irrespective of the orthogonal coordinate system considered in practical application cases. This option is not accessible for any indicial-tensor formulation as are basis dependent.
Evidences to support most of these claims will be provided in the present article. The differences of the present manuscript with respect to [31] are as follows: 1. Planas and co-workers [31], mainly focused in the small strain setting where elastostatics, J 2 plasticity and mean dilatation approach were considered. Only a brief outline for large strain elastostatics was presented in [31]. Moreover, many steps in the derivation as well as a detailed description of the tensors involved for large strains were omitted for the sake of brevity in [31]. On the contrary, in the present work, a detailed derivation of the nodal equilibrium equations and related linearisation for large strain elastodynamics in the initial configuration is presented.
Most of the researches in the Computational Mechanics community, are familiar with the Voigt based approach and, in the authors' opinion, they might appreciate a detailed derivation of the equations and tensors appearing in the component-free formulation for large strains. 2. In the present work, a careful comparison between the component-free formulation and the Voigt based formulation for large strain elastodynamics in the initial configuration is included. This comparison was also omitted in [31] for large strain elastostatics. Again, in the authors' opinion, this comparison is adequate to illustrate the scope of the component-free formulation and its differences with the Voigt based approach. 3. The component-free finite element formulation and implementation proposed in this work is developed for a general isotropic hyperelastic material in the initial configuration, with a strain energy function per unit volume, W , expressed in terms of the principal invariants I C , I I C , I I I C of the right Cauchy-Green tensor, i.e. W = W (I C , I I C , I I I C ). Then, particularisation to a compressible Neo-Hookean material is presented [43]. However, although in [31] the spatial tangent density is outlined for a general isotropic material in the current configuration, no initial configuration counterpart is presented nor particularisation to any well known hyperelastic material is adopted. 4. In the present document, the component-free finite element formulation leads to different expressions to be implemented than those obtained by the standard Voigt based approach. Therefore, accuracy and efficiency comparison with respect to standard Voigt based approach must be performed for different initial boundary value problems. In the present work, five different benchmarks at finite strains are proposed for accuracy and efficiency assessment. However, no benchmarks were proposed in [31] to verify the performance of the B free formulation. 5. Planas and co-workers [31] only derived preliminary analysis to show that the component-free formulation is never slower than the traditional Voigt based approach, i.e. they did not show definitive results regarding compared computational costs. On the contrary, in the present work, the five different benchmarks proposed at finite strains enable the authors of the present work to con-clude that the component-free finite element formulation and implementation shows a clear quantitative improvement in terms of computation time with respect to the classical Voigt form.
The performance of the component-free approach for large strain isothermal elastodynamics is illustrated by five numerical examples, one in 1D, three in 2D and one in 3D.
All the proposed benchmarks are particularized for a compressible Neo-Hookean material [7,13,18,29,43]. The 1D application case is related with a Riemann problem at finite strains [9,11]. On the other hand, the 2D benchmarks are related with the plane strain Cook's membrane [37] under a dynamic load, a plane strain solid block with low stiffness under the action of gravity and a large deformation vibration of a cantilever beam under its own weight. The 3D benchmark is related with a thick-walled cylinder pinched by two opposite line dynamic loads. The 1D benchmark allows the verification of the performance of the proposed component-free approach regarding the propagation of discontinuities in non-linear materials at large strains. Moreover, the reproduction of the 2D and 3D benchmarks allows the verification of the performance under a shearing dominant behaviour (Cook's membrane), compression dominant behaviour (solid block), bending dominant behaviour (cantilever beam) and a shearing-compression-bending combined behaviour (pinched cylinder) . Comparison with the results obtained by the standard Voigt based finite element implementation is included. The rest of the paper is organized as follows. For the sake of completeness, in Sect. 2, the Lagrangian equations for large strain isothermal elastodynamics are presented. In Sect. 3, the variational formulation of the governing equations, time integration and linearisation of the variational form is outlined. In Sect. 4, a detailed derivation of the nodal equilibrium equations and related linearisation and factorisation within the component-free approach is provided in the initial configuration. In this section, key comparisons with the standard Voigt base finite element formulation are described. In Sect. 5 different expressions of the tangent matrix density are derived for a general isotropic hyperelastic material within the component-free formulation. In Sect. 6, the proposed elastodynamic component-free Lagrangian finite element formulation is considered to reproduce five different benchmarks including a Riemann problem, a plane strain Cook's membrane under a dynamic load, a plane strain solid block with low stiffness under the action of gravity, a large deformation vibration of a cantilever beam under its own weight and a pinched thick-walled cylinder. This section ends with an efficiency analysis in terms of computation time. Relevant conclusions, as well as future and promising research lines, are drawn in Sect. 7.

Initial configuration large strain isothermal elastodynamics governing equations
This section presents the equations, in the initial configuration, for large strain elastodynamics, symbols and notations for subsequent sections. For further details the interested reader is referred to the specialized literature [6,7,13,18,29,41].
In what follows B 0 ⊂ R 3 is the initial configuration of a deformable body. As the body deforms, its configuration changes with time. Let t ∈ I ⊂ R denote time and for each and every t ∈ I associate a unique configuration of the body B t ⊂ R 3 , refer to as the current configuration at time t. Material points of the body are labelled by its position X in the initial configuration B 0 and by x in the current configuration B t at time t. The deformation χ : B 0 × I → B t , which is considered twice-continuously differentiable with respect to position and time, is a bijection for each and every t ∈ I that maps points X ∈ B 0 to points The deformation gradient F is a second order tensor that maps material line elements dX in the initial configuration B 0 to line elements dx in the current configuration B t , and is defined through any of the following equalities Here and henceforth, Grad and Div denote gradient and divergence operators in the initial configuration, i.e. with respect to X. It is assumed F to be non-singular and therefore J = det(F) = 0 (> 0). By introducing the displacement vector the deformation gradient F can be expressed as where I is the second order identity tensor. The velocity v and acceleration a of a material point are defined, respectively, by for a fixed X, i.e., by the material time derivative. Usually, the large strain isothermal elastodynamic equations in the initial configuration can be derived for conjugate pairs P andḞ or S andĖ, where P is the first Piola-Kirchhoff stress tensor, S is the second Piola-Kirchhoff stress tensor and E is the Green-Lagrange strain tensor. Both scenarios will be considered in the present document.
For the analysis of a non-linear initial boundary value problem for large strain isothermal elastodynamics in the initial configuration, a coupled system of partial differential equations has to be solved, which consist of the local balance of linear momentum and a constitutive equation. In the present article, a general isotropic hyperelastic material in the initial configuration [10], with a strain energy function per unit volume, W , expressed in terms of the principal invariants I C , I I C , I I I C of the right Cauchy-Green tensor C = F T F, is adopted first. Then, particularisation to the well known compressible Neo-Hookean material [43] is considered for the benchmarks. The Lagrangian formulation of the conservation of mass is simply stated byρ 0 = 0. This implies that the initial density ρ 0 of the material is constant and therefore does not need to be considered as part of the unknown. Thus, the strong form of these equations reads where g is the gravity acceleration vector.
To complete the mathematical model, boundary and initial conditions need to be described. Boundary conditions have to be prescribed on ∂B 0 . More specifically, accepting the decomposition ∂B 0 = Γ u ∪ Γ σ where Γ u ∩ Γ σ = ∅ and the overline denoting closure, boundary conditions for displacements are prescribed on Γ u while boundary conditions for tractions are considered on Γ σ . This leads to u = u on Γ u (8) P N = t on Γ σ (9) where N is the outward pointing unit normal at each point on the boundary ∂B 0 . Finally, initial conditions for the displacements u and velocities v at time t = t 0 ,where usually t 0 = 0, need to be prescribed Eqs. (3), (6) and (7) together with the boundary conditions, Eqs. (8) and (9), and initial conditions, Eqs. (10) and (11), completely define the non-linear initial boundary value problem for large strain isothermal elastodynamics in the initial configuration where the primary unknown variable is the displacement vector u defined in Eq. (2). Finally, u · v denotes the standard scalar product of two vectors u, v ∈ R 3 . A second order tensor A is a linear map that acts on a vector u giving a vector v = Au. The tensor product of two vectors u, v ∈ R 3 is denoted by u ⊗ v and is a second order tensor defined in component-free notation by (u ⊗ v) w = (v · w) u. Furthermore, composition of two second order tensors A, B is a second order tensor defined by (AB) u = A (Bu) while the double contraction A : B is defined by A : B = tr A T B where tr (u ⊗ v) = u · v and A T is the transpose of the second order tensor A and is the second order tensor defined by u · A T v = Au · v.

Variational formulation of the governing equations
Following standard variational principles [14,18], we define for each t ∈ I ⊂ R the space of trial solutions and the space of weighting functions where H 1 is the first-order vector value Sobolev space. Defining the residual R (u) = Div P + ρ 0 (g −v), multiplying this residual by a weighting function η ∈ V , integrating the residual over the initial configuration B 0 , applying integration by parts, making use of the divergence theorem and introducing the traction boundary condition Eq. (9), the variational form of the linear momentum reads as Then, the variational form of the non-linear initial boundary value problem for large strain isothermal elastodynamics in the initial configuration reads as follows while the initial conditions Eqs. (10) and (11) must be fulfilled.

Time integration of the variational form
Time integration of Eq. (14) is performed to obtain a discrete evolution of the displacement vector u [33]. To this end, let I = [0, T ] be the time interval where integration of the equations take place and let t n = t 0 + nΔt with n = 0, . . . , N be a sequence of discretization nodes. Let u n , v n , a n be the approximation, at node t n , of the displacement vector u (X, t n ), the velocity vector v (X, t n ) and the acceleration vector a (X, t n ), respectively. In the present work the well known Newmark method [26,43] is adopted. This method is based on the following approximation of the velocities and accelerations at time t n+1 upon the displacements where the constants α i for i = 1, . . . , 6 are defined by [22,43] Δt (18) After evaluation of the variational expression Eq. (14) at time t n+1 , and keeping in mind Eqs. (16) and (17), the time discretization of the variational form of the non-linear initial boundary value problem reads as follows Knowing u n , v n and a n at time t = t n find u n+1 ∈ S n+1 at time t n+1 such that subjected to the initial conditions u 0 = u 0 and v 0 = v 0 . The choice of the well known implicit Newmark method [26] is justified for the sake of clarity in the description of the component-free methodology applied to the large strain isothermal elastodynamics problem in the initial configuration, avoiding elaborated time integration schemes that could obscure the process.

Linearisation of the variational form
In order to solve the non-linear initial boundary value problem for large strain isothermal elastodynamics in the initial configuration, linearisation of the variational form must be performed. The linear part of the variational form at u = u n+1 with respect to the initial configuration is given by (20) where G u n+1 , η is just the variational expression Eq. (14) evaluated at time t n+1 while is the directional derivative of G (u, η) at u = u n+1 in the direction Δu n+1 . Assuming the tractions t and gravity are independent of the displacement vector u (dead loads) and keeping in mind that the initial density ρ 0 of the material remains constant during the deformation process, the directional derivative DG u n+1 , η · Δu n+1 consists of the linearised internal virtual work and the linearised inertial term, thus where P u n+1 + εΔu n+1 means evaluation of the first Piola-Kirchhoff stress tensor P at u n+1 + εΔu n+1 , while a u n+1 + εΔu n+1 means evaluation of the acceleration vector a at u n+1 + εΔu n+1 . As the linearisation process is performed with respect to the initial configuration while P = P (F) and applying the chain rule, the linearisation of the internal virtual work reads where A n+1 is the fourth order elasticity tensor defined by If the relation P = FS between the first and the second Piola-Kirchhoff stress tensors is now considered, keeping in mind that S = 2 ∂ W ∂C and by further use of the chain rule, the linearisation of the internal virtual work then reads where C n+1 is the fourth order elasticity tensor defined by with E the Green-Lagrange strain tensor. Minor symmetries of the fourth order tensor C have been used in Eq. (24). Finally, as the initial density ρ 0 of the material remains constant during the deformation process and considering the approximation Eq. (16) of the accelerations at time t n+1 upon the displacements, the linearisation of the inertial term reads as follows where I is the second order identity tensor. Then if Eqs. (23) and (25) are considered, the directional derivative of G (u, η) at u = u n+1 in the direction Δu n+1 reads as follows On the contrary, if Eqs. (24) and (25) are taken in to account, the directional derivative of G (u, η) at u = u n+1 in the direction Δu n+1 reads as follows

Component-free finite element discretization
In this section the nodal equilibrium equations and related linearisation are derived applying the component-free Bubnov-Galerkin finite element discretization process [27,31] to the variational form and its linearisation.
To obtain the finite element discretization of the variational problem Eq. (19), the domain B 0 is subdivided into n e elements Ω e connecting n nd nodes. Each node has as many degrees of freedom as the number of spatial dimensions. The space of trial solutions S n+1 at time t n+1 is restricted to a finite dimensional subspace S h n+1 spanned by isoparametric finite element shape functions N α : Ω → R. Standard Bubnov-Galerkin procedure is considered, in which the space of weighting functions V is restricted to a finite dimensional subspace V h spanned by the same shape functions N α used to approximate the primary variable u n+1 at time t n+1 . Spaces S h n+1 and V h might be express as where ξ is the local orthogonal coordinate system at the reference element Ω . The finite element discretization of the variational problem Eq. (19) then reads as follows Knowing u n , v n and a n at time t = t n Substituting the approximation η = N α (ξ ) η α in the variational form of the linear momentum Eq. (14) at time t n+1 yields where Einstein summation convention of repeated indices is considered while superscript n + 1 means evaluation at time t n+1 . In the case of P n+1 and a n+1 this evaluation is done through its dependency upon u n+1 = N β (ξ ) u n+1 β . In the case of P n+1 by means of Eqs. (3) and (7) and in the case of a n+1 by means of Eq. (16). For t n+1 only direct evaluation at time t n+1 is required.
As the evaluation at time t n+1 of the different terms in expression Eq. (30) have been clarified, superindex n + 1 are omitted for brevity from now on.
Using standard component-free tensor algebra, the internal virtual work density P : Grad N α η α becomes and keeping in mind that Eq. (30) must hold ∀η ∈ V h the following nodal equilibrium equation for each node α is obtained At this point a comparison between the component-free approach and the standard Voigt based finite element procedure for large strain elastodynamics in the initial configuration might be beneficial. Before the comparison, the process to compute the nodal internal virtual work in the standard finite element method, in terms of the Voigt's notation, is briefly outline: 1. Starting with the natural internal virtual work density P : Gradη, the following symmetrization is performed [43] P : Gradη = FS : Gradη = S : F T Gradη = . . .

S :
where the last equality is a definition for the variation of the Cauchy-Lagrange strain tensor δE, while S is the second Piola-Kirchhoff stress tensor and F is the deformation gradient. 2. The weighting functions are approximated by η = N α η α and applying standard tensor algebra the following expression of the gradient in the initial configuration is achieved Substituting back into the symmetrized internal virtual work, the following expression is obtained Then, going back to index notation and letting S = where factorisation of the nodal weighting function η αa is obtained by the commutative property of scalar multiplication.
where the symmetry of S and δE is considered. 6. The variation of the Cauchy-Lagrange strain tensor δE can be further express by Bη α where B is the operator defined by [27,43] 7. Finally, taking into account the arbitrariness of the weighting functions η, the nodal internal virtual work density is expressed by B T S in the standard Voigt based finite element method.
After comparison of the nodal internal virtual work density P GradN α obtained with the component-free approach and the one given by the finite element method in Voigt form B T S, the following important remarks follow. As we are dealing with a variational approach, a key ingredient in the formulation of the nodal internal virtual work density is how and when the factorisation of the nodal weighting functions η α is performed. In the case of a Voigt notation approach, this factorisation is obtained by means of the commutative property of scalar multiplication, after the index notation is considered (step 4) while the remainder terms are reorganized in the matrix Voigt form B T S, where a dense 6 × 3 matrix B is computed and then a 3 × 6 times a 6 × 1 vector multiplication is finally performed.
On the contrary, in the component-free approach the factorisation of the nodal weighting functions η α is performed by means of the component-free operation Eq. (31), avoiding the use of any index notation for this factorisation and the subsequent reorganisation in matrix Voigt form. The absence of any index notation in the factorisation of the nodal weighting functions η α leads to a concise and transparent nodal internal virtual work density expression P GradN α . In contrast to the Voigt's notation, the nodal internal virtual work term P GradN α in Eq. (32) is computed directly as a 3 × 3 matrix 3 × 1 vector multiplication, once an orthogonal basis E A and e A are defined in the initial and in the current configuration, respectively. Thus if P = P a A e a ⊗ E A and A crucial point is that this basis plays no role in the factorisation of the nodal weighting functions η α in the component-free approach. A direct consequence of this novel procedure to performed the factorisation of the nodal weighting function is that less number of operations per element is involved to compute the nodal internal virtual work density in the component-free approach than in the classical Voigt based finite element method.
The differences between the component-free and the Voigt based approaches become even more clear when dealing with the discretisation of the linearised variational form of the linear momentum. This will be the topic for the rest of the present Sect. 4.
The starting point to derive the component-free discretisation of the linearised variational form of the linear momentum is either Eqs. (26) or (27) (superindex n + 1 omitted for brevity), depending on the stress tensor adopted in the formulation. If the first Piola-Kirchhoff stress tensor P is considered, then the discretisation of the linearised variational form starts from Eq. (26). On the contrary, if the second Piola-Kirchhoff stress tensors S is used, then the discretisation starts from Eq. (27).
where A GradN α , GradN β and C GradN α , GradN β are second order tensors that can be defined in component-free notation as follows: For any vectors a, b, c ∈ R 3 and fourth order tensor T, the second order tensor T {c, b} is defined by Substituting Eq. (38) in Eqs. (36) and (39) in Eq. (37) while keeping in mind that Eqs. (36) and (37) must hold for every η α , the following expressions are achieved As Δu β are nodal directions at time t n+1 , they can be taken off the integral. Then, from Eq. (41), the following nodal directional derivative of the equilibrium equation, for the pair of nodes (α, β), is obtained where K mat αβ is the tangent matrix related with the constitutive relation defined by while K iner αβ is the tangent matrix related with the inertial term defined by Similarly, from Eq. (42), the following nodal directional derivative of the equilibrium equation, for the pair of nodes (α, β), is obtained whereK geo αβ is the geometric tangent matrix defined bỹ whileK mat αβ is the tangent matrix related with the constitutive relation defined bỹ and again, K iner αβ is the tangent matrix related with the inertial term, already defined in Eq. (45) The Eq. (45) is the standard mass matrix under Newmark time integration and requires no further clarification. However, this is not the case for Eqs. (44) and (48) that still require an explicit expression for the second order tensors A GradN α , GradN β and C GradN α , GradN β for a particular constitutive relation. This is done in Sect. 5 for a general isotropic hyperelastic material in the initial configuration and particularised afterwards to the well known compressible Neo-Hookean material.
At this point, a comparison between the component-free and Voigt based approaches regarding the process to formulate the tangent matrix is now performed. For the sake of completeness, the process to derive the tangent matrix in the standard Voigt form [43] is now outline: [DP (u) · Δu] : Gradη (49) Considering the relation P = FS between the first and the second Piola-Kirchhoff stress tensors, keeping in mind that S = S (C) with C = F T F, that the following relations hold where C = 2 ∂S ∂C and E = 1 2 (C − I) and applying the product rule for the directional derivative, then Eq.
in the initial stress, the following expression is obtained 3. The second term δE : C : ΔE in Eq. (50) depends upon the constitutive relation. After substitution in this second term the approximations Eq. (51) and Eq. (52), the standard finite element method again goes back to index notation.
This time the index notation allows the factorisation of both, the nodal weighting function η αa and the increment of the displacement field Δu βa , by the commutative property of scalar multiplication. 4. Voigt's transformation enters the scene and the fourth order constitutive tensor C is transformed into a 6 × 6 symmetric matrix D being where the components C ABC D = 2 ∂ S AB ∂C C D depend upon the constitutive relation considered. 5. Also by the Voigt's transformation, the expression also in Eq. (54) becomes BΔu β where B is given by Eq. (33). Therefore, Eq. (54) becomes η α B T DBΔu β 6. Taking into account the arbitrariness of the weighting functions η, and that Δu β are nodal directions and can be taken off the integrals, in the standard Voigt based finite element method the tangent density is finally expressed by including the initial stress term and the material term.
After comparison of the tangent densities obtained with the component-free approach and the Voigt notation 1. Once more, as we are dealing with a variational approach, a key ingredient in the computation of the tangent matrix, for a pair of nodes (α, β), is how and when the factorisation of both, the nodal weighting functions η α and the increment of the displacement field Δu β , is performed.
In the case of a Voigt notation approach, this factorisation is obtained, again, by means of the commutative property of scalar multiplication, after the index notation is considered (step 3) while the remainder terms are reorganized in the matrix Voigt form B T DB, to take advantage of the symmetries of the second order tensor E and the fourth order tensor C, by means of the dense 6 × 3 matrix B and the symmetric 6 × 6 matrix D, respectively. 2. On the contrary, in the component-free approach the factorisation of the nodal weighting functions η α and the increment of the displacement field Δu β is performed by means of the component-free operations Eqs. (38) and (39), avoiding both, the use of any index notation and the subsequent reorganisation in matrix Voigt form. The absence of any index notation in this double factorisation leads to the second order tensor A GradN α , GradN β or FC GradN α , GradN β F T , depending on the stress P or S adopted in the formulation, respectively. The specific expression for both second order tensors, A GradN α , GradN β and C GradN α , GradN β , will be clarify for a general for a general isotropic hyperelastic material in the initial configuration and then particularised to the well known compressible Neo-Hookean material in Sect. 5.

Direct comparison of the component-free tangent density
for S and the tangent density in Voigt notation give rise to the following equality while the Voigt expression B T DB is a 3 × 6 times 6 × 6 times 6 × 3 matrix multiplication, the component-free expression FC GradN α , GradN β F T is a 3 × 3 times 3 × 3 times 3 × 3 matrix multiplication, once an orthogonal basis E A and e A are defined in the initial and in the current configuration, respectively. Therefore, less number of operations per element is involved to compute the nodal tangent density in the component-free approach than in the classical Voigt based finite element method. 4. Again, a crucial point is that the basis considered to compute the component-free tangent density plays no role in the factorisation of the nodal weighting functions η α neither in the factorisation of the increment of the displacement field Δu β . This novel procedure to performed this double factorisation, leads to a concise and transparent nodal tangent density in terms of the second order tensor A GradN α , GradN β or C GradN α , GradN β .
Finally, the following general remarks are stated: • As isoparametric finite element shape functions N α : Ω → R are considered, its gradient GradN α at an arbi-trary node α is computed by where the gradientF is defined bŷ whileĜradN β means gradient of the isoparametric shape function N β : Ω → R with respect to the local orthogonal coordinate system ξ at the reference element Ω and X β is the position of the node β in the initial configuration B 0 . • In order to compute integrals in Eqs. (32,44,45,47,48) standard assembly process is considered. • The final non-linear algebraic system of equations in Eq. (29) can be solved by means of an iterative procedure, like the Newton-Raphson method, to determine the unknown displacements u n+1 at time t n+1 . Combining Eqs. (32) and (43) or (32) and (46), the Newton-Raphson method might be expressed as where subindex (i) is the iteration index.

General isotropic hyperelastic material
In this subsection, explicit expression of the second order tensor A GradN α , GradN β that appears in Eq. (44) and the second order tensor C GradN α , GradN β that appears in Eq. (48), are derived for a general isotropic hyperelastic material, with a strain energy function per unit volume, W , expressed in terms of the principal invariants I C , I I C , I I I C of the right Cauchy-Green tensor C, i.e. W = W (I C , I I C , I I I C ).
The steps followed to obtain the material tangent matrix density for a general isotropic hyperelastic material are: 1. The expression of the fourth order elasticity tensor C related with the second Piola-Kirchhoff stress S, i.e. C = 2 ∂S ∂C = 4 ∂ W ∂C∂C , is derived by means of the chain rule. 2. The expression of the second order tensor C {v, b} for any vectors v, b ∈ R 3 is derived by means of a component-free manipulation through the definition Eq. (40). At the end of this step, an explicit expression for C GradN α , GradN β will be derived.
3. Based on the expression for C GradN α , GradN β , an explicit expression for A GradN α , GradN β is finally obtained.
It is well known that the second Piola-Kirchhoff stress tensor S for a hyperelastic material with a strain energy function per unit volume W that fulfils objectivity requirement can be obtained [13,29] by means of S = 2 ∂ W ∂C where C is the right Cauchy-Green tensor.
If the hyperelastic material is isotropic, then the strain energy function, W , may be regarded as function of the principal invariants I C , I I C , I I I C . Therefore, by the chain rule, the second Piola-Kirchhoff stress tensor S reads [43] by further use of the chain rule and making use of the component-free product rule ∂(αS) ∂C = S ⊗ ∂α ∂C + α ∂S ∂C for a scalar valued α and second order tensor valued S tensor functions, the following expression for the fourth order elasticity tensor C is obtained [10] where the coefficients Γ i are scalar functions of the principal invariants I C , I I C , I I I C .
The second order tensor C {v, b}, for any vectors v, b ∈ R 3 , is defined by Eq. (40). This definition is a componentfree definition of a second order tensor, i.e., the second order tensor C {v, b} is understood as a linear map that acts over a vector a, i.e. [C {v, b}] a, giving as a result the vector As the space of fourth order tensors is a vector space itself, then Taking into account the following three component-free relations for any vectors a, b ∈ R 3 , expression Eq. (60) reads Applying the second order tensor obtained in Eq. (61) to an arbitrary vector v ∈ R 3 , i.e.[C : (a ⊗ b)] v, and considering the following component-free expressions the following expression follows As the space of second order tensors is a vector space itself, the vector a can be extracted as a common factor. Moreover, considering now definition Eq.  (a ⊗ b)] v, and keeping in mind that C is symmetric, the following final expression is obtained for the second order tensor C {v, b} In the authors' opinion, the formulation and implementation with Voigt notation, by means of the so-called matrix D, of the material tangent matrix density for a general isotropic hyperelastic material would be non-trivial as compared with the component-free approach shown in this Subsect. 5.1.

Compressible Neo-Hookean material
The compressible Neo-Hookean material might be expressed by the following strain energy function per unit volume [43] W where J is the determinant of the deformation gradient F, I C is the first principal invariant of the right Cauchy-Green tensor C, while Λ and μ are the Lamé constants. Moreover, the second Piola-Kirchhoff stress tensor S = 2 ∂ W ∂C can be computed by means of Eq. (58) as while the fourth order elasticity tensor C = 2 ∂S ∂C reads as follows In order to derive explicit expression for the second order tensor C GradN α , GradN β , particularized to the well known compressible Neo-Hookean material, one could reproduce the steps follow in Subsect. 5.1.
However, under the scope of the component-free approach, the second order tensor C GradN α , GradN β has been easily obtained for a general isotropic hyperelastic material. Thus, we can take advantage of it.
Following this second approach, it would be enough to consider v = GradN α and b = GradN β in Eq. (63) together with the following values for the coefficients Γ i Moreover, if the explicit expression for the second order tensor A GradN α , GradN β is sought for the compressible Neo-Hookean material, then combination of Eqs. (63) and (64) with the values of the coefficients Γ i given by Eq. (68) would be also enough.

Validation benchmarks and efficiency analysis
The performance of the component-free approach for large strain isothermal elastodynamics is illustrated by five numerical examples, one in 1D, three in 2D and one in 3D.
All the proposed benchmarks are particularized for a compressible Neo-Hookean material. The 1D application case is related with a Riemann problem at finite strains [9,11].
On the other hand, the 2D benchmarks are related with the plane strain Cook's membrane [37] under a dynamic load, a plane strain solid block with low stiffness under the action of gravity and a large deformation vibration of a cantilever beam under its own weight [34]. The 3D benchmark is related with a thick-walled cylinder pinched by two opposite dynamic line loads. The 1D benchmark allows the verification of the performance of the proposed component-free approach regarding the propagation of discontinuities in non-linear materials at large strains. Moreover, the reproduction of the 2D and 3D benchmarks allows the verification of the performance under a shearing dominant behaviour (Cook's membrane), compression dominant behaviour (solid block), bending dominant behaviour (cantilever beam) and a shearing-compression-bending combined behaviour (pinched cylinder). Comparison with the results obtained by the standard Voigt based finite element implementation is included. In all benchmarks considered, time integration is performed by means of the Newmark method with γ = 0.5 and β = 0.25, e.g., an implicit unconditionally stable second order trapezoidal rule is considered. The final non-linear algebraic system is solved by means of Newton-Raphson iterative procedure with T O L = 10 −10 . All numerical solutions with the component-free proposed formulation are obtained at a quadratic rate of asymptotic convergence.

1D Riemann problem at finite strains
In the present case of validation, a Riemann problem, consisting of a longitudinal propagation of an initial discontinuity in the velocity field over an infinite homogeneous compressible Neo-Hookean bar, is analysed. This benchmark is of paramount importance, as numerical methods might fail to capture the complex behaviour display in time-dependent non-linear hyperbolic conservation laws, where shock and rarefaction waves might appear.
This problem might be represented in a Lagrangian configuration by a longitudinal wave propagating in the X −direction where x = (x, y, z) and X = (X , Y , Z ), on a homogeneous compressible Neo-Hookean continuum, subjected to a piecewise constant initial velocity field where v = u t while v L > v R , on an undeformed bar and without boundary conditions in X −direction.
The present initial problem leads to the following quasilinear second order partial differential equation for u (X , t).
The analytical solution might be obtained as a first order non-linear hyperbolic system in the independent variables ε = u x and v = u t [9,11] and consists of two different shock waves travelling in opposite directions from the initial discontinuity location. Clustering the independent variables as w = (ε, v) this analytical solution might be expressed as where s 1 and s 2 are the shock speed of the left and right travelling sock waves, respectively, while w M is the intermediate state connected to the left state w L = (0, v L ) through the left travelling sock wave and to the right state w R = (0, v R ) through the right travelling sock wave. For further details in this analytical solution, the interested reader is referred to the specialized literature [9,11]. The numerical set up of this benchmark is as follows. A bar of 100m length and 0.1m wide is analysed under a plane strain representation, Fig. 1. In order to obtain a 1D longitudinal Fig. 1 Initially undeformed bar of 100m × 0.1m. Velocity field discontinuity at X = 50m with v L = 10m/s and v R = 0m/s  Fig. 2 Zero initial deformation u x along the bar at time t = 0s (upward); initial discontinuity in the velocity field at time t = 0s (downward) wave propagation, zero vertical displacement is imposed in the upper and lower boundaries of the bar Fig. 1. The spatial discretization has been performed by a unique row of 1000 four-node linear isoparametric quadrilateral elements (Q4) with four Gauss points per element. The material parameters considered for the compressible Neo-Hookean material, including the initial density, are displayed in Table 1.
The Young modulus E, density ρ 0 , time step Δt = 5 · 10 −5 s and spatial discretization considered, give a Courant number of 0.98. In order to represent a bar of infinite length, Neumann homogeneous boundary conditions were placed at X = 0m and at X = 100m. The initial discontinuity in the velocity field is set at X = 50m with v L = 10m/s and v R = 0m/s on an initially undeformed bar Fig. 1.
Before the comparison between the component-free and the Voigt based methods is performed, it would be beneficial to visualised the longitudinal propagation of the initial discontinuity in the velocity field over the bar. To this end, in Figs. 2, 3, 4 and 5, snapshots of the propagation at times t = 0s, 0.001s, 0.002s, 0.003s can be observed In Fig. 6, comparison of the results obtained with the component-free finite element formulation and the standard Voigt based finite element method are shown for the deformation field u x = F 11 − 1 and the velocity field u t , along the bar, at time t = 0.003s. The initial discontinuity in the velocity field can be also appreciated in Fig. 2 as well as the analytical solution outline in the present Subsect. 6.1.
As per Fig. 6 it is clear that both approaches, componentfree and Voigt base, are able to capture quite well the sock strength, i.e. the amplitude of the jump, and the sock speed, i.e. the propagation speed of the discontinuity, both in the deformation and the velocity field. The spurious numerical oscillations observed in Fig. 2 are due to the time integration considered. As no algorithmic damping is included in the Newmark method (γ = 0.5), no dissipation of highfrequency modes is introduced.
Although the nodal internal virtual work and the tangent density matrix of the proposed component-free Lagrangian

Solid block under the action of gravity
In the present case of validation, a plane strain solid block with low stiffness under the action of gravity is analysed. This benchmark allows the verification of the performance of the proposed Lagrangian formulation approach under a compression dominant behaviour. The numerical set up of this benchmark is as follows. A solid block of 1m × 1m is analysed under plane strain Fig. 7. At the base of the solid block, a zero vertical displacement is imposed being free the horizontal component Fig. 7. The rest of the boundaries are considered Neumann homogeneous. The spatial discretization is performed by means of 10 × 10 eight-node isoparametric quadrilateral finite element (Q8) with nine Gauss points per element Fig. 7.
The material parameters considered for the compressible Neo-Hookean material, including the initial density, are displayed in Table 2.
The Young modulus E, density ρ 0 , time step Δt = 10 −2 s and spatial discretization considered, give a Courant number of 0.98. Initial conditions u (X, 0) = 0 and v (X, 0) = 0 are considered. The gravity g = 10s/m 2 is applied gradually by a ramp curve during the first second followed by a constant value, e.g. the gravity is applied by means of g (t) = g · min (t, 1).
In Fig. 8, the deformed mesh at time t = 1.08s along with the contour plot of the determinant J of the deformation gradient F, can be appreciated for the component-free approach. The minimum vertical displacement at node A ( Fig. 7) is attained for the first time at t = 1.08s. The initial undeformed mesh is also included in Fig. 8. No amplification factor is considered.
In Fig. 9, comparison of the results obtained with the component-free formulation and the standard Voigt based finite element method are shown for the vertical displacement with respect to time at node A (Fig. 7).
As per Fig. 9 it is clear that both approaches, componentfree and Voigt base, obtain the same evolution of the vertical displacement at the central node of the upper horizontal boundary. As no algorithmic damping is included in the Newmark method (γ = 0.5) no dissipation is introduced and the same oscillations, same phase and amplitude, are recorded.
Again, no significant differences can be appreciated between both formulations in Fig. 9. Thus, the componentfree proposed formulation for large strain elastodynamics might be understood as a novel, concise, transparent, simple to implement, component-free tensor based approach being

Cook's membrane under a dynamic loading
In the present case of validation, a plane strain Cook's membrane under a dynamic load is analysed. This benchmark allows the analysis of the performance of the componentfree Lagrangian approach under a combination of bending and shearing. The numerical set up of this benchmark is as follows. The domain consists of a tapered cantilever [8] with the lefthand side clamped and a time-dependent shear load f (t) in vertical direction applied on the right-hand side (Fig. 10). The profile of the time-dependent shear load corresponds with a ramp curve during the first half second followed by a constant value f max = 20M Pa, e.g. the shear load follows f (t) = 2 f max · min (t, 0.5). The thickness of the Cook's membrane is chosen to be 1 mm. The spatial discretization is performed by means of 10 × 20 eight-node isoparametric quadrilateral finite element (Q8) with nine Gauss points per element.The geometry, boundary conditions and applied load are illustrated in Fig. 10.
The material parameters considered for the compressible Neo-Hookean material, including the initial density and maximum applied load f max , are displayed in Table 3.
The values considered for the Young modulus E, density ρ 0 , and maximum of the applied load f max are obtained from [37]. The time step considered for the simulation is Δt = 0.02s. Initial conditions u (X, 0) = 0 and v (X, 0) = 0 are adopted.
In Fig. 11, the deformed mesh at time t = 0.76s along with the contour plot of the σ x x Cauchy stress component can be appreciated for the component-free approach.The Cauchy  J PF T . The maximum vertical displacement at node A (Fig. 11) is attained for the first time at t = 0.76s. The initial undeformed mesh is also included Fig. 11. No amplification factor is considered.
In Fig. 12, comparison of the results obtained with the component-free finite element formulation and the standard Voigt based finite element method are shown for the vertical displacement with respect to time at node A in the upperright corner. The stationary solution provided by [37] is also included. Once more, as per Fig. 12, it is clear that both approaches, component-free and Voigt base finite elements, describe the same evolution of the vertical displacement at the upperright corner. Again, no algorithmic damping is included in the Newmark method (γ = 0.5) therefore, no dissipation is introduced and the same oscillations, same phase and amplitude, are recorded for both approaches.
Again, no significant differences can be appreciated between both formulations in Fig. 12. Thus, the componentfree proposed formulation for large strain elastodynamics might be understood as a novel, concise, transparent, simple to implement, component-free tensor based approach being as accurate as the traditional Voigt based formulation under a combination of bending and shearing.

Large deformation vibration of a cantilever beam
In the present case of validation, a large deformation plane strain vibration of a cantilever beam under its own weight [34] is analysed. This benchmark allows the analysis of the performance of the component-free Lagrangian approach under a bending dominant behaviour. The numerical set up of this benchmark is as follows: It consists of a compressible Neo-Hookean cantilever beam clamped in the left boundary and loaded by the gravity, g = 10m/s 2 . As shown in Fig. 13, the left end is fixed while the rest of the boundaries are traction free. The large deformation vibration of the beam is induced by its own weight where gravity is fully applied at t = 0s. This large deformation vibration is analysed along 3s using time steps of Δt = 10 −3 s. The spatial discretization is performed by means of 8 × 2 eight-node isoparametric quadrilateral finite element (Q8) with nine Gauss points per element.The geometry, boundary conditions, mesh and applied load are illustrated in Fig. 13.
The material parameters considered for the compressible Neo-Hookean material are displayed in Table 4.   The values considered for the Young modulus E, Poisson ν and density ρ 0 are obtained from [34]. Initial conditions u (X, 0) = 0 and v (X, 0) = 0 are adopted.
In Fig. 14, the deformed mesh at time t = 1.4s along with the contour plot of the vertical displacement can be appreciated for the component-free approach. The initial undeformed mesh is also included in Fig. 14. No amplification factor is considered.
In Fig. 15, comparison of the results obtained with the component-free finite element formulation and the standard Voigt based finite element method are shown for the vertical displacement with respect to time at node A in the lower-right corner.
Once more, as per Fig. 15, it is clear that both approaches, component-free and Voigt based finite elements, describe the same evolution of the vertical displacement at the lowerright corner. As no significant differences can be appreciated between both formulations in Fig. 15, the component-free proposed formulation for large strain elastodynamics might be understood as a novel, concise, transparent, simple to

Pinched thick-walled cylinder
In the present case of validation, an adaptation of the traditional pinched thin-walled cylinder [30] to a pinched thick-walled cylinder with solid elements (instead of shell elements) is analysed. This benchmark allows the analysis of the performance of the component-free Lagrangian approach under a combination of shearing, compression and bending in 3D.
The numerical set up of this benchmark is as follows. The domain consists of a single-layered thick-walled cylindrical tube with rigid end-diaphragms where x and y degrees of freedom are fixed while the z degree of freedom is kept free. Geometry and initial configuration can be found in Fig. 16. The thick-walled cylindrical tube is pinched by two opposite time-dependent line loads f (t) parallel to the y axis (Fig. 17). The profile of the time-dependent line load f (t) corresponds with a ramp curve during the first half second followed by a constant value f max = 3G N /m, e.g. the line load follows the equation f (t) = 2 f max · min (t, 0.5). Due to the symmetry, the computational domain is restricted only to one eight of the cylinder as can be observed in grey colour in Fig. 17. The spatial discretization is performed by means of 50 quadratic 20-node hexahedral isoparametric elements (H20) with 14 Gauss points per element [15].
The material parameters considered for the compressible Neo-Hookean material, including the initial density and maximum applied load f max , are displayed in Table 5.
The time step considered for the simulation is Δt = 0.02s. Initial conditions u (X, 0) = 0 and v (X, 0) = 0 are adopted.
The deformed mesh for two different viewpoints at time t = 0.64s along with the component y displacement contour plot can be appreciated for the component-free approach in Fig. 18. The maximum component y displacement at node    (Fig. 16) is attained for the first time at t = 0.64s. No amplification factor is considered. Comparison of the results obtained with the componentfree finite element formulation and the standard Voigt based finite element method are shown in Fig. 19 for the component y displacement with respect to time at node A (Fig. 16).
Once more, as per Fig. 19, it is clear that both approaches, component-free and Voigt base finite elements, describe the same evolution of the component y of the displacement at node A. The same oscillations, phase and amplitude are  Component y displacement at node A for the pinched thick-walled cylinder recorded for both approaches. No significant differences can be appreciated between both formulations in Fig. 19.
Thus, the component-free proposed formulation for large strain elastodynamics might be understood as a novel, concise, transparent, simple to implement, component-free tensor based approach being as accurate as the traditional Voigt based formulation under a combination of shearing, compression and bending in 3D.

Efficiency analysis
As no significant accuracy differences can be appreciated between the component-free and Voigt formulations regarding the numerical results for each and every proposed benchmark, in this section comparison is focus in computation time. Thus, an efficiency analysis is performed for four of the benchmarks considered in the manuscript: The plane strain Cook's membrane under a dynamic load, the plane strain solid block with low stiffness under the action of gravity, the large deformation vibration of a cantilever beam under its own weight and the thick-walled cylinder pinched by two opposite time-dependent line loads. For each of these cases, four simulations have been considered. Two of them with the proposed component-free finite element formulation and implementation and other two with the classical Voigt form. Furthermore, 8-node and 4-node isoparametric quadrilateral finite elements with 9 and 4 Gauss points per element, respectively, have been considered for the 2D benchmarks. For the 3D benchmark 20-node and 8-node isoparametric hexahedral finite elements with 14 and 8 Gauss points per element, respectively, have been considered. Efficiency analysis is performed in terms of computation time. Results are displayed in Tables 6 and 7.
It can be observed that the proposed component-free finite element formulation and implementation is more efficient than the classical Voigt form. For the linear elements Q4 in 2D,the proposed component-free finite element formulation and implementation shows an average improvement of 20% of the computation time required by the classical Voigt form. This average improvement is reduced to 16% for the quadratic elements Q8 in the 2D. On the contrary, the improvement obtained in the 3D benchmark is of 21% regardless of the order considered for the elements, H8 or H20.

Conclusions
In the present work, a component-free finite element Lagrangian formulation and implementation for isothermal elastodynamics at large strain is presented. The key idea of the proposed formulation is how the factorisation of the nodal weighting function for the nodal internal virtual work and the factorisation of both, the nodal weighting functions and the increment of the displacement field for tangent matrix, is performed. This factorisations are performed in the proposed approach by means of component-free operations instead of using the commutative property in index notation reorganized with the aid of Voigt notation. This point of view leads to concise and transparent nodal internal virtual work and tangent matrix densities, that can naturally be applied to general classes of constitutive models, being as accurate as the classical Voigt form but more efficient in terms of computation time.
The proposed formulation for large strain elastodynamics has been validated against five different benchmarks giving very similar results than standard Voigt based finite element approach in terms of accuracy. On the contrary, an average improvement of 18% of the computation time require by the classical Voigt form is observed for the proposed componentfree formulation. The benchmarks were a Riemann problem at finite strains, a Cook's membrane under a dynamic load, a solid block with low stiffness under the action of gravity, a large deformation vibration of a cantilever beam under its own weight and a pinched thick-walled cylinder.
In the authors' opinion, the clarity of the componentfree formulation should be favoured with respect to Voigt based approach. Following this leading idea, the authors have recently proposed these new paradigm to the novel mesh-free variational based Material Point Model [20] with very promising results. Moreover, application of the same component-free philosophy to model and compute large strain plasticity as well as multiphase materials will be considered.
the ETSI Caminos, Canales y Puertos, from the Universidad Politécnica de Madrid, as well.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.