Sensitivity analysis based multi-scale methods of coupled path-dependent problems

In the paper, a generalized essential boundary condition sensitivity analysis based implementation of FE2 and mesh-inelement (MIEL) multi-scale methods is derived as an alternative to standard implementations of multi-scale analysis, where the calculation of Schur complement of the microscopic tangent matrix is needed for bridging the scales. The paper presents a unified approach to the development of an arbitrary MIEL or FE2 computational scheme for an arbitrary path-dependent material model. Implementation is based on efficient first and second order analytical sensitivity analysis, for which automaticdifferentiation-based formulation of essential boundary condition sensitivity analysis is derived. A fully consistently linearized two-level path-following algorithm is introduced as a solution algorithm for the multi-scale modeling. Sensitivity analysis allows each macro step to be followed by an arbitrary number of micro substeps while retaining quadratic convergence of the overall solution algorithm.

Index of the last and the current macro step n , n+1 Index of the last and the current micro step s Index of micro step at the end of the last macro step (r ) Index of micro problem φ set of variables calculated at the selected macro element and transferred to the corresponding micro problem (also sensitivity parameters) S Set of variables calculated at the selected micro problem and returned to the corresponding macro element φ Me

Introduction
Multi-scale methods are a growing trend in computational mechanics, especially with increasing capabilities of computers. Multi-scale methods originate from the demand to accurately model heterogeneous materials, like fiber reinforced composites, particle reinforced adhesives, concrete and even metal [4,8,17]. The goal of multi-scale modeling is to design a combined macroscopic-microscopic computational algorithm that is numerically more efficient than solving of the full microscopic model directly. At the same time, it gives us the information we need with the desired accuracy. For the overview of multi-scale methods reader is referred to [1,3,6,30]. The use of different kinds of multiscale methods is limited by the characteristics of the problem to be solved. Roughly, we can separate multi-scale methods in two groups: on methods that are based on homogenization techniques and on domain decomposition methods. A basic hypothesis of homogenization techniques is a complete separation of scales, where the size of heterogeneities is assumed to be infinitely smaller than the structural dimensions. Homogenized material behavior of representative volume elements (RVEs), which contains microstructure, is considered to be representative of the entire or part of the structure. Standard two-level finite element homogenization approach (FE 2 ) described in [16,29] is appropriate for the problems where scales are separated far enough and are only weakly coupled, see [2,5,23,27]. If the difference between two scales is finite, in the region of high gradients or in the case of localization [25] the FE 2 approach fails. Then some sort of domain decomposition method can be applied. One such method is substructuring method, now more commonly named mesh-in-element (MIEL) scheme, which was described in [9].
Within the standard implementation of nonlinear multiscale methods, only the macro scale is parametrized by the load factor. Consequently, each macro step is followed by exactly one step at the micro level and a path-following algorithm is applied only at the global level. The first aim of this paper was to develop a nonlinear multi-scale computational scheme with two interacting path-following methods at two levels. An algorithm will be derived for consistent parametrization of both macro and micro problems leading to a two-level path-following algorithm. For the purpose of convergence comparison, one method from each group was implemented, FE 2 and MIEL.
In the literature, a lot of attention has been paid to the computation of macroscopic tangent as an essential and numerically demanding part of any multi-scale simulation. The possibilities vary from expensive and inaccurate but general finite difference approximation of macroscopic tangent, to various ways how to derive corresponding analytical expressions (for discussion on methods see e.g. [28]). One of the alternatives is also standard sensitivity analysis of coupled path-dependent problems, as introduced in [11,21]. In the primal analysis, the response of the system is evaluated, whereas in sensitivity analysis the derivatives of the response, e.g. displacements, strains, stresses or energy, with respect to arbitrary design parameter are sought. For the automation of the multi-scale methods, sensitivity analysis with respect to parameters used to define micro level boundary conditions is needed. It will be shown that the consistent linearization of the two-level path-following algorithm requires the implementation of relative sensitivity analysis instead of a full one and that a second order sensitivity analysis is also needed.
The second aim of the paper is to present advantages of analytical essential boundary condition sensitivity analysis based implementation in comparison with the classical ways of implementing multi-scale methods based on the calculation of Schur complement of micro tangent matrix (see e.g. [16,22,28] for FE 2 method and [9] for MIEL method). This is especially important for path-dependent problems such as finite strain plasticity, where consistent linearization is of high importance. It will be shown that for the MIEL type of methods the analytical second order sensitivity analysis is numerically superior with respect to Schur complement implementation.
Another motivation was to create a computational environment, where the multi-scale program code is automatically derived and various types of multi-scale and single-scale approaches can be freely mixed while retaining quadratic convergence of the Newton-Raphson procedure. To achieve the goal, the introduced method uses an advanced feature of software tool AceGen [12]. AceGen is an automatic code generator, where automatic differentiation technique [7], automatic code optimization and code generation are combined with computer algebra system Mathematica [20]. The size of the code is reduced through control of expression swell [10]. The automatic-differentiation-based (ADB, [11]) formulation enables unification and automation of various multi-scale approaches for arbitrary nonlinear, path-dependent material models (e.g. general finite strain plasticity). A short overview of multi-scale methods based on first and second order boundary condition sensitivity analysis for the linear case was already given in [15]. Numerical simulations were performed with AceFEM numerical environment [14] that has built-in support for numerically efficient first and second order analytical sensitivity analysis.
The paper is organized as follows. Micro problem formulation based on finite strain elasto-plastic model is briefly described in Sect. 1.1 as a basis for all subsequently derived multi-scale methods. After an introduction to automatic differentiation based notation (ADB) in Sect. 1.2, a generalized two-level path-following multi-scale algorithm is derived in Sect. 2. The basis for the multi-scale formulations is primal and sensitivity analysis of the micro problem as given in Sect. 3. MIEL and FE 2 multi-scale methods are then described in detail, focusing on their implementation in Sect. 4, followed by numerical examples presented in Sect. 5.

Micro problem definition
At the micro level we will consider finite element formulation of rate-independent nonlinear problems in solid and structural mechanics, such as an arbitrary finite strain rateindependent elasto-plastic material model. Here, only essential equations are presented; for more details see e.g. [26]. Finite strain, isotropic, elasto-plastic model is defined by its elastic strain energy function, plastic evolution equations and the method for time integration of evolution equations. Some of the possible variants are presented in [13]. The actual material model used to run numerical examples is summarized in Box 1. However, all the methods presented are general, independent of the specific material model.

Box 1. Micro problem material model
Formulation is based on multiplicative split of micro deformation gradient F m , the components of an inverse right Cauchy plastic strain tensor C −1 p as plastic state variables, elastic left Cauchy strain tensor b e (1), Neo-Hookean strain energy function W (b e ) (2) and Mises yield function with exponential hardening law (4). Backward Euler is combined with the exponential map for a stable, volume conserving integration of evolution equations [26]. Discretized evolution equations (5), together with yield condition φ = 0, form a set of algebraic equations Q mg (6) for a set h mg (7) of unknown components of plastic strain tensor C −1 p and plastic multiplier γ . C −1 p n and γ n are values of plastic strain tensor and plastic multiplier at the end of the last load step. The dependency of Eq. (6) on the values of variables at the end of the last load step (h mg n ) makes the whole problem path-dependent.
Standard weak form of equilibrium equations is then written as where first Piola-Kirchhof stress tensor P m can be obtained from the elastic strain energy W by P m = ∂ W /∂F m . After finite element discretization of deformation gradient F m = F m p me , where p me is a set of nodal degrees of freedom of eth micro element at the current load step. The variation δF m = ∂F m /∂p me δp me leads from (9) together with the standard Gauss integration of weak form and standard procedure of assembly of element contributions (denoted here with A operator) to a set of algebraic equilibrium equations (10). Equations are at each Gauss point coupled with an additional set of equations Q mg (6). The result is the following integration point coupled system of nonlinear algebraic equations where R me is contribution of eth element to global residual R m and R ext m is a vector of external forces. p m denotes a set of micro level nodal unknowns, h m = n tg g h mg is a set of unknowns of all Gauss point problems. G e is a set of Gauss points of eth element and w gp is Gauss point weight. The Gauss point contribution R mg to the element residual R me leads from (9) to

Automatic differentiation based (ADB) notation
The automatic differentiation (AD) technique can be used for the evaluation of the exact derivatives of any arbitrary complex function defined by an algorithm via chain rule and represents an alternative solution to the numerical differentiation and symbolic differentiation. The result of AD procedure is called "computational derivative" and is written here aŝ δa has a dual purpose to indicate the mathematical operation of differentiation as well as to indicate that the AD technique is used to obtain the required quantity. If, for example, alternative or additional dependencies for a set of intermediate variables b have to be considered for differentiation, then the AD exception is indicated by the following formalism which indicates that during the AD procedure, the total derivatives of variables b with respect to variables a are set to be equal to matrix M. b in (13) my or may not be algorithmically a function of a. When b is algorithmically function of a then (13) defines that the true derivatives ∂b ∂a are neglected and replaced by a matrix M. When b is not algorithmically a function of a then (13) introduces from the algorithmic point of view an artificial dependency between a and b. The automatic differentiation exceptions are the basis for the ADB formulation of computational problem.
For example, the Gauss point residual R mg is defined by Eq. (12). However, form (12) is not numerically efficient from the automatic differentiation point of view. Numerically efficient ADB form of (12) is derived as As a side effect of the iterative solution of Gauss-point equations (6), there exist an implicit (algorithmic) dependency of h mg on F m . The AD exception Dh mg DF m = 0 in (14) hides this dependency from automatic differentiation procedure and ensures correct evaluation of the weak form equations. In (14) we start with the scalar and make only one call to AD procedure, which is optimal for the backward mode implementation of automatic differentiation as shown in [11].
The introduced ADB notation is abstract, thus any sufficiently sophisticated software for automatic differentiation can be used for the actual implementation. The ADB notation can be directly translated to the AceGen input and is part of automatic generation of numerically efficient program codes. Details of the method and of the corresponding software Ace-Gen together with numerous examples of AceGen inputs can be found in [10][11][12]14]. The actual AceGen and AceFEM inputs are for the complex multi-scale problems addressed in the paper too lengthy to be included in the paper. However, they are freely available at http://symech.fgg.uni-lj.si/ Examples/MultiScale.pdf, in a form of Mathematica notebook at http://symech.fgg.uni-lj.si/Examples/MultiScale.nb or as a part of software documentation available at http:// symech.fgg.uni-lj.si/Download.htm.

Generalized two-level path-following multi-scale algorithm
For highly nonlinear problems in general the solution cannot be achieved in one step. More efficient procedures can be derived when the resulting system of algebraic equations can be naturally parametrized. Various path-following algorithms, such as constant load-stepping, adaptive load stepping or arc-length methods, can then be applied to solve the nonlinear problem. Within the standard implementation of multi-scale methods only the macro scale is parametrized. Consequently, each macro step is followed by exactly one step at the micro level and a path-following algorithm is applied only at the global level. Here, an algorithm is derived for consistent parametrization of both macro and micro problems leading to two-level path-following algorithm. For the sake of simplicity, the two-level constant load stepping algorithm is derived. However, it can be easily extended to other path-following approaches. Let k be the index of the last calculated macro load step and k + 1 the current macro load step, as shown in Fig. 1. Furthermore, let n be the index of the last converged micro load step, n + 1 the current micro load step, s the index of the micro load step at the end of the last converged macro load step, n m the number of micro level steps within the macro level step and s + n m the index of the micro load step at the end of the macro load step as presented in Fig. 1.
As an example, problems in solid mechanics and nonlinear structural mechanics subjected to quasi-static proportional load are frequently parametrized by introducing loading parameter λ M . λ M will be used to parametrize macro problem. The final value of parameter λ M is usually determined by the problem at hand, e.g. as total given load factorλ M . In this case the total load is split into n M macro steps. The finite element discretization of macro level then leads to a set of nonlinear equations R M at the current load level λ M = λ M k+1 Fig. 1 Generalized two-level path-following, multi-scale algorithm where R Me denotes the contribution of internal forces of eth macro element to the nodal force vector and R ref M is the reference load vector associated with the pattern of the applied nodal forces. p M represents a set of nodal unknowns of the problem at macro level and S Me is a set of variables transferred from the micro level problems to eth macro element. S Me is composed of contributions of one or several micro problems. Thus, S Me = r ∈M e S (r ) , where S (r ) is the contribution of the r th micro problem and M e is a subset of micro problems that contribute to the eth macro element. For a general scheme it is irrelevant what the data represents physically.
Within various multi-scale methods the coupling of the scales can be done in several ways. The paper addresses methods where micro-macro coupling is achieved by expressing the essential boundary conditions of micro level problem as a function of data calculated at macro level. Letp m be a set of micro problem nodal unknowns with imposed homogeneous essential (Dirichlet) boundary conditions, φ a set of variables calculated at macro level for the current macro load level λ M on which a selected micro problem depends andp m (φ) a function such that at the and of the macro stepp m =p m (φ). φ is composed of components of macro deformation gradient in the case of FE 2 method and of components of nodal displacements of macro element in the case of MIEL method. The actual form ofp m (φ) depends on the multi-scale scheme and is given in the following sections.
Let λ m be a current value of a micro level parameter. At the end of the last macro step, we additionally define λ m s as a value of a micro level parameter andp m s as a value of nodal unknowns with imposed essential boundary conditions. Linear interpolation ofp m within the macro step then leads to the following parametrization of micro level problem The total increment p m of the micro essential boundary conditions within the macro load step is defined by The micro level parameter introduced with (16) ensures continuous parametrization of micro problem and has the following properties for the kth micro step: (λ m −λ m s ) ∈ [0, 1], λ m s = k and λ m = k + 1 at the end of macro step. With the introduction of parameter λ m , the solution of micro problem within the kth macro step is achieved in n m micro steps with associated solution vectors. The finite element discretization at micro level leads from (10) to the following integration point coupled system of nonlinear algebraic equations for the chosen micro problem (19) where equilibrium equations R m are coupled with discretized evolution equations at the Gauss point level Q mg . A standard two level Newton--Raphson method can be used to solve the resulting Gauss point coupled system of algebraic equations for the unknown p m and h m , as described in [11]. In order to have quadratically convergent multi-scale solution algorithm, we also need consistently linearized macro equilibrium equations (15). The linearization of (15) leads to where φ Me = r ∈M e φ (r ) is composed of variables calculated at the eth macro element and transferred to a subset of micro problems M e that effects the eth macro element. Partial derivatives in (20) are explicit and can be easily derived analytically for a specific multi-scale scheme.

Derivation of total derivative DS
where first order derivatives Dp m /Dφ, Dh m /Dφ and second order derivatives D 2 p m /Dφ 2 , D 2 h m /Dφ 2 are implicit and require differentiation of complete path-following algorithm for the solution of selected micro problem. This can be done using analytical sensitivity analysis procedures, such as described in [11]. φ represents input data for the selected micro level simulation and is used to calculate essential boundary conditions (16). Thus, for the evaluation of implicit derivatives, a boundary condition sensitivity analysis is needed with components of φ as sensitivity parameters. The solution of a path-dependent micro problem, in general depends on all variables transferred from the macro level to the micro level in all macro steps. Consequently, a complete set of sensitivity parameters of the selected micro problem would be composed of all variables transferred from the selected macro element to the selected micro problem. Sensitivity analysis for a large number of parameters requires significant amount of memory as well as computation time. However, it is not actually needed. The variables transferred in kth step (φ = φ k+1 ) affect the selected micro problem only from the micro step at the beginning of the macro step (micro step with the index sth) and implicit derivatives with respect to φ are not needed any more after the completion of the macro step. Consequently, and implicit derivatives with respect to φ do not appear in the macro problem after the completion of the macro step. Thus, at any given time only a set of sensitivity parameters φ that belongs to the current macro step has to be considered, provided that (22) holds. Since Q mg depends only on h mg n , it is sufficient for (22) to hold to set at the start of each micro problem increment (at the sth micro step).

Two-level path-following algorithm
The algorithm that summarizes the above considerations is presented in Box 2. First, the micro level equations (18) (19) are solved for unknown p m and h m at fixed p M with the use of Newton method, which is also applied to solve the macro equilibrium equation (15) in an outer loop leading to a nested iteration-subiteration solution scheme for unknown p M , p m and h m . For the sake of simplicity, the algorithm is written for the constant time stepping with n M macro steps and n m micro steps per macro step. However, it can be easily extended to an arbitrary adaptive time stepping scheme. It is assumed here that the micro problem is path-dependent, thus the state of all micro problems has to be stored somewhere at the end of the solution of each macro step and restored at the beginning of each macro iteration. The basic idea of this paper is that any FE code that supports first and second order sensitivity analysis can be turned into a fully consistent, numerically efficient, quadratically convergent nonlinear multi-scale code with minimal or even without any additional coding at the level of micro finite elements. Of course, one can use finite difference approximation to evaluate macro tangent modulus K M . However, the resulting code is numerically efficient and inexact tangent can affect the rate of convergence of iterative procedure. In any case, one has to write additional code for data management, solution of the macro problem and for parallelization of the multi-scale algorithm. Let us assume that the code supports primal, first and second order sensitivity analysis. It is then needed for the implementation of the particular multi-scale scheme to define the following quantities and expressions: • micro problem sensitivity parameters as a function of macro element unknowns (φ(p Me )), • boundary conditions of micro problem as a function of sensitivity parameters (p m (φ, λ m )), • derivatives of boundary conditions with respect to sensitivity parameters ( ), • micro level variables that are passed to macro level (S), • total derivative of micro level variables with respect to sensitivity parameters (DS/Dφ), p m s ←p m ; // update and store macro and micro data Box 2. Two-level path-following multi-scale algorithm

Primal and sensitivity analysis of micro problem
The general procedures for primal and first order sensitivity analysis are for an arbitrary problem presented in detail in [11]. Here, we focus on essential boundary condition sensitivity analysis that is needed for the implementation of multi-scale methods. For this purpose, first order essential boundary condition sensitivity analysis is extended to second order essential boundary condition sensitivity analysis.

Solution and automation of primal problem
The primal problem at micro level (18), (19) represents a system of Gauss point coupled nonlinear algebraic equations that can be solved using standard nested iterationsubiteration Newton-Raphson iterative procedure (see e.g. [11]). First, the Gauss point equation (19) is solved for h mg at fixed p me with the use of Newton method, which is also applied to solve the equilibrium equation (18) in an outer loop for unknown p m . The tangent operator for the inner loop K Q is given by and consistent micro tangent matrix K m is written as Evaluation of the micro tangent matrix (25) requires proper consideration of the implicit dependency h mg (p me ) introduced by the local iterative procedure. The missing implicit derivative ∂h mg ∂p me can be easily obtained from (19) and introduced as an AD exception in (25) (for more details see [11]).

Solution and automation of sensitivity problem
For the essential boundary condition sensitivity analysis we define the residuals and the vectors of unknowns in (18) and (19) as a function of sensitivity parameters φ by wherep m (φ, λ m ) is a set of nodal DOF with prescribed essential boundary conditions defined by (16).
At the level of individual finite element, the set of nodal unknowns p me =p me ∪p me includes degrees of freedom with prescribed essential boundary conditionp me and true degrees of freedomp me , because they are at the element level indistinguishable. The first order sensitivity problem can be obtained from the primal problem by differentiating equations (26) and (27) Tangent matrix K m is already evaluated and factorized from the primal problem. Therefore, only vector IR m on the righthand side of Eq. (30) has to be calculated in order to obtain the resulting system of linear equations. This vector is called independent first-order sensitivity pseudo-load vector. After obtaining where I Z g is an additional auxiliary variable introduced to increase numerical efficiency. It can be evaluated during the evaluation of IR mg , stored in memory and used later for the evaluation of Dh mg /Dφ I . Functionp m (φ, λ m ) can be arbitrary complex and, in general, cannot be input data of the finite element analysis. However, it is not the relationp m (φ, λ m ) itself that is needed within the sensitivity analysis, but its first and second derivatives. Let φ I and φ J be an arbitrary essential boundary condition sensitivity parameters. The rate of change of essential boundary conditions are called first and second order essential boundary condition velocity fields. The values of first and second order essential boundary condition velocity fields at the nodes of eth element are defined by where n p is the total number of element nodal DOFs. The velocity field is zero for the true degrees of freedom. Thus, the proper definition of element velocity fields is sufficient to make the difference between the degrees of freedom with prescribed essential boundary condition and true degrees of freedom at the finite element level. The actual analytical expressions for (32), (31) and (33) are rather lengthy, but they can be obtained automatically using the automatic differentiation. For this purpose, an automatic differentiation based notation or ADB notation of the terms is needed. A general ADB notation of first order terms follows from (32), (31) and (33) where all implicit derivatives are replaced by appropriate AD exceptions and leads to where again only vector I JR m on the right-hand side of Eq. (39) has to be calculated. This vector is called independent second-order sensitivity pseudo-load vector. After Dφ I Dφ J is given here by Dh mg n Dφ J = J H n g ,

MIEL method
MIEL (mesh-in-element) method is a multi-scale finite element method that can be classified as a domain decomposition method. This method is appropriate for cases where the difference between two scales is finite and the scales remain coupled, or when in the region of high gradients the FE 2 multi-scale approach fails. The MIEL scheme was described in [18,19,24]. Next, we developed an automatized sensitivity analysis based version of the MIEL method. At the macro level, we have compatible interpolation of unknown fields at the boundary of macro elements, whereas material characteristics, inhomogeneities, inner structure, such as openings, incisions of different materials, are defined only at micro scale. In Fig. 2 the MIEL procedure is presented. Let assume the standard interpolation of displacements u M on the boundary of the macro element where N i ( ) are finite element shape functions, = (ξ, η, ζ ) reference coordinates and u Me i are displacements in ith macro element node. To ensure compatibility of displacements at macro and micro level, we impose the essential boundary conditions at the complete boundary of the micro mesh bȳ whereū m s ( ) are displacements at the boundary at the end of the last macro step. The derivatives of (44) with respect to components of macro element nodal displacements are given by A set of macro element unknowns is p Me = j,k u Me j k and p m is composed of the micro mesh nodal displacements. Thus, (44) defines the dependency between the degrees of freedom with prescribed essential boundary condition at the micro levelp m =p m (p Me , λ m ) and macro element unknowns p Me .
The macro element residual R Me is in the case of MIEL obtained by the integration of the internal forces, part of weak form (9), over the micro mesh, where the micro deformation where again K Mg is a contribution to the macro element tangent evaluated at micro mesh Gauss points. The residual and tangent matrix are for each macro element obtained directly from the micro scale problem and each macro element is associated with exactly one micro problem. Macro element performs only proper transfer of components of the macro element residual vector and tangent matrix from micro scale to macro scale finite element assembly procedure.  implementation of (48) and (50), we also need ADB form of (48) and (50). From (14) ADB form of (48) and (50) leads to

Sensitivity analysis is required for the evaluation of implicit
where Y φ = Dp me Dp Me and Y φφ = D 2 p me Dp 2 Me are first end second order sensitivities calculated and stored during the analysis.

Schur complement based implementation of MIEL
Let us consider formulations where the solution is within one macro step path-independent, such as hyper-elastic problems solved with an arbitrary number of micro steps or elastoplastic problems solved at the micro level in one load step. In this case, an alternative formulation of MIEL based on the calculation of Schur complement is possible, as originally presented in [18]. Let us form, at the converged state of the micro problem, a full set of equations that include unconstrained p m and constrainedp m unknowns by Schur complement of (53) leads to reduced set of equations K cc p m = R c , where K cc , and R c are condensed tangent matrix and residual of micro problem, respectively. Since the relationp m =p m (p Me , λ m ) is linear [see (44)], we can writē where T is a transformation matrix (for details see [9]). The macro element residual and tangent matrix are then expressed by With R Me and K Me known, one can apply the algorithm presented in Box 2, with sensitivity analysis related parts omitted. The size of K cc is equal to the number of constrained DOFs at the boundary of the mesh and grows with the micro mesh density. For densely meshed microstructure the calculation of the Schur complement inflicts high memory allocation and is time consuming. Contrary, the number of sensitivity parameters is the same as the number of nodal unknowns of the macro element, thus independent of micro mesh density. Schematic comparison can be seen for 2D case discretized with 4 nodded elements in Fig. 3. For Schur complement implementation, condensation is done with respect to DOFs of 20 border nodes. The dimension of the resulting matrix K cc is 40 × 40. To get macro element tangent matrix K Me with dimension 8 × 8, additional transformations (55), (56) need to be performed. With the growth of mesh density, also the number of micro-structure border nodes grows and with that the dimension of the matrix to be calculated. In the case of sensitivity based implementation, the second order sensitivity analysis is needed with respect to 8 DOFs in macro element corner nodes and summation of (51), (52) over the micro mesh integration points. The comparison of the computational cost of the two implementations is done for the 3D case, which is more computationally demanding than the 2D case. In Fig. 4, the calculation time for the Schur complement and for the second order sensitivity analysis is presented in relation to the density of micro mesh. The example is composed of one 3D, hexahedral macro element. The macro element is uniformly subdivided into n × n × n micro mesh. Two micro material models are considered, finite strain elasto-plastic as defined in Sect. 1.1 and hyper-elastic based on hyper-elastic strain energy (2). The Schur complement's computational time grows polynomially, whereas sensitivity calculation retains approximate linearity with the number of equations at the micro level. The timing of the sensitivity analysis increases with the complexity of the material model and the number of DOFs of the macro element. However, overall behavior remains the same.

FE 2 method
Standard two-level finite element homogenization approach FE 2 is appropriate for the problems where scales are separated far enough and are only weakly coupled, see [16]. The FE 2 method was already implemented using sensitivity analysis in [22,27], but without two interacting path-following schemes. Within the FE 2 approach we have one micro FE model, also called a representative volume element (RVE), at each macro mesh integration point as shown in Fig. 5. All information about micro-structure is obtained from computations at the micro level by averaging the material response characterized by an appropriate stress measure and constitutive tangent matrix over RVE. With the Gauss point contribution to the macro level weak form (P M · δF M ) and macro level discretization of deformation gradient δF M = ∂F M ∂p Me δp Me , the macro element residual leads to where the macro level first Piola-Kirchoff stress tensor P M is obtained by averaging the micro level first Piola-Kirchoff stress tensor P M = {P m }. The operation of averaging is here denoted by {·}. Several types of boundary conditions can be imposed on the RVE: e.g., fully prescribed displacements and fully prescribed traction, which are based on the uniform strain and stress assumptions and periodic boundary conditions that enforce a displacement constraint, which is suited for periodic media. Here, periodic boundary conditions are achieved (see e.g. [16]) by applying first the prescribed displacements in the corners of RVE bȳ where F M s is macro deformation gradient at the end of the last macro step. The derivatives of (59) with respect to components of F M are given by Thus, (59) defines the dependencyp m =p m (F M , λ m ) between the set of micro nodal unknowns with prescribed essential boundary conditionp m and the macro deformation gradient F M . For the unconstrained boundary nodes, the periodicity of boundary conditions is adopted with the use of Lagrange multipliers (for details see [27]). Note that the introduction of Lagrange constraints only extends the vector of micro level unknowns p m with Lagrange multipliers and micro level residual R m with constraint equations and it does not change the primal and sensitivity analysis procedures described in Sect. 3 .
Differentiation of (48) then leads to the macro element tangent matrix where DP M DF M = { ∂P m ∂p me Dp me DF M } is macroscopic constitutive matrix obtained by averaging the microscopic constitutive matrices. where M e is a set of micro problems that corresponds to G e . For the numerically efficient implementation of (58) and (62), we also need the ADB form of (58) and (62). From P M · δF M = S · δF M the ADB form of (58) and (62) leads to S=const. (64) and where Y φ = Dp me Dp Me are already calculated and stored first order sensitivities.

Schur complement based FE 2 implementation
As in the case of the MIEL method, the Schur complement of constrained nodal DOF at the micro level can be used to calculate macro element residual and tangent matrix. The method leads to the traditional implementation of the FE 2 method, as introduced in [16], and it will not be repeated here. The number of RVE corner nodes is constant, which makes the cost of calculating the Schur complement independent of the density of the micro mesh thus, the advantages of using the sensitivity analysis are less pronounced than for MIEL. Note that the standard method is only consistently linearized for the problems that are path independent within a single macro step.

Numerical examples
Numerical examples were calculated using program packages AceGen and AceFEM [12]. Finite element user subroutines for primal and analytical first and second order sensitivity analyses were automatically derived, optimized and written in C language with the use of AceGen automatic code generator. The MIEL and FE 2 methods based on sensitivity analysis as well as the one based on the Schur complement were implemented within AceFEM environment according to algorithm defined in Box 2. Intel MKL sparse linear algebra numerical library was used for the lin-ear algebra operations (calculation of the Schur complement and the solution of linear systems of equations).
For 2D examples, nine nodded, isoparametric, quadrilateral, plane strain elements are used integrated with 3 × 3 Gauss integration, and for 3D examples, eight nodded, isoparametric, hexahedral elements are used, integrated with 2 × 2 × 2 Gauss integration. For all examples, a finite strain elasto-plastic material model as described in Sect. 1.1 is used at the micro level.
The abbreviations used to indicate specific combination of methods solution procedures are structured as follows where method can be MIEL or FE 2 , n M is the number of macro steps or "Adaptive" for adaptive macro time stepping, n m is the number of micro steps for each macro step or "Adaptive" for adaptive micro time stepping, and implementation is "Sens." for sensitivity analysis based implementation, "SchurMMA" for the Schur complement based formulation implemented in Mathematica and "SchurMKL" for the Schur complement based formulation implemented with Intel MKL Library. The Schur complement based implementation is computationally identical to the sensitivity analysis based implementation for n m = 1. Although Mathematica and MKL both calculate theoretically the same Schur complement, the algorithm implemented in MKL performs perturbation of the zeros at the main diagonal resulting in slightly imprecise tangent matrix as shown and explained on examples in Sect. 5.2.

Validation of implementation of multi-scale algorithm
The first numerical example is a three-dimensional cantilever with clamped right and left end as shown in Fig. 6. Uniform pressure p = 10 in the vertical z direction was imposed at the top surfaces of the middle part of the cantilever. The dimensions of the cantilever are 12 × 2.4 × 2.4. 3D, hexahedral elements are used at both levels. Material properties are E = 21,000, ν = 0.3, σ y0 = 24 and K h = 100. A homogeneous mesh is used at both levels, thus for the purpose of validation the micro level is uniform and no microstructure is present. The simulations were performed with adaptive time stepping at both levels. The displacements in the z direction of nodes on line AB are presented for all simulations in Fig. 7. The extent of the plastic zone at the end of the simulation is shown in Fig. 6, where red color indicates the plastic region.
Multi-scale results are compared with the single-scale results.
The same finite elements are used for the single scale mesh as for the micro level mesh of the multi-scale simulation. First the FE 2 method is verified by comparing the results obtained by single scale analysis with mesh 20 × 4 × 4 with Fig. 6 Clamped cantilever with macro and micro mesh and enforced natural and essential boundary conditions multi-scale analysis for macro mesh density 20 × 4 × 4 and micro mesh densities 2 × 2 × 2 and 10 × 10 × 10. The results must be independent of micro mesh density due to the homogeneous micro mesh and exact enforcement of periodicity of micro mesh boundary conditions. Multi-scale results must also be exactly the same as single-scale results, which is shown in Fig. 7 (curves 1, 2 and 3). This verifies the FE 2 implementation.
For the MIEL method, the results of single scale simulation and multi-scale simulation can be exactly the same only for micro mesh density 1 × 1 × 1. This is shown in Fig. 7  (curves 1, 4). This verifies also the MIEL implementation. With the change of micro mesh density to 2 × 2 × 2, 5 × 5 × 5 and finally to 10 × 10 × 10 (curves 5, 6, 7), the MIEL results get close to single scale FEM solution obtained with the mesh 80 × 16 × 16 (curve 8). This is the consequence of a better description of the deformation field over the domain of macro element, which partially eliminates the locking behavior of the isoparametric hexahedral element. The effect is similar to that of enhanced strain finite elements, where additional degrees of freedom are added inside the elements.

Convergence rate of two-level path-following iterative procedure
The convergence rate of the two-level path-following iterative procedure defined by an algorithm in Box 2 is investigated on an example from the previous section. The simulation is performed in 10 steps with a constant load increment λ M = 0.1. Homogeneous micro mesh 5 × 5 × 5 is used for all cases. Each macro step is followed by one micro step (denoted by -10/1-) or 5 micro steps (denoted by -10/5-). Results of the convergence rate of the two-level pathfollowing iterative procedure for the last macro load step, where most of integration points are already in the plastic regime, are shown in Tables 1 and 2. The effect of the number of micro steps and the type of implementation (the Schur complement or sensitivity analysis based) is investigated.    Table 1 shows that in the case when 1 macro load step is followed by one micro step (MIEL-10/1-SchurMKL and MIEL-10/1-Sens.), convergence is quadratic and the results are exactly the same regardless of implementation. Sensitivity based implementation retains quadratic convergence also for n m = 5, while the SchurMKL based implementation converges very slowly. The column denoted by MIEL-10/5-Sens.end contains a special case, where the sensitivity equations given in Sect. 3.2 are not resolved after each micro step, but only at the end of micro solution. This is again equivalent to implementation MIEL-10/5-SchurMKL. It shows that only a fully consistent sensitivity analysis ensures quadratic convergence of the overall MIEL algorithm.
Secondly, the convergence rate of FE 2 scheme was investigated in Table 2. The same conclusions can be drawn as for MIEL. Only fully consistent sensitivity analysis ensures quadratic convergence of the overall FE 2 algorithm. The last two columns contain the results of the Schur complement based formulation implemented directly in Mathematica. This is not numerically efficient, but it is necessary to show that the FE 2 -10/1-SchurMMA implementation is numerically identical to the FE 2 -10/1-Sens. implementation. The imposition of periodic boundary conditions using Lagrange constrains results in the loss of positive definiteness of the tangent matrix as well as produces zeros at the main diagonal. Some algorithms for the evaluation of the Schur complement, such as the one implemented in the Intel MKL library, perform perturbation of the zeros on the main diagonal resulting in imprecise Schur complement. This imprecision is sufficient to alter, although not significantly, the convergence behavior. This case is shown in fourth column in Table 2, designated as FE 2 -10/1-SchurMKL.

Numerical efficiency of two-level path-following iterative procedure
The numerical efficiency of the two-level path-following iterative procedure is investigated on an example from Sect. 5.1. All simulations were performed on PC with Intel i9 2.8GHz,16 Core processor and 128 GB RAM. Micro problems were solved in parallel on 14 cores. Mathematica was used only as a steering application for parallelization and the control of the iterative procedure, while all computationally intensive operations were performed with compiled C codes. The material model used is a 3D finite strain elasto-plastic model based on an exact exponential map (see e.g. [13]), which is by itself computationally intensive. Consequently, the administrative cost turns out to be negligible when compared to the actual computational cost. The effect of implementation on computational time is for FE 2 formulation presented in Table 3. An example introduced in Sect. 5.1 is solved with the macro mesh density 20×4×4 in 10 load steps with a constant load increment λ M = 0.1. The computational time normalized with respect to FE 2 -10/1-Sens. formulation is presented along with the total number of Newton--Raphson iterations for all load steps and the total number of micro problems solved during the complete simulation. The simulation using FE 2 -10/1-Sens. formulation took 1968.5 s of real time. The results are presented for the number of micro steps n m = 1 and n m = 5 and the density of the micro mesh 5×5×5 and 10×10×10. The first order sensitivity analysis based formulation is in all cases faster than the corresponding Schur complement based formulation. The loss of quadratic convergence of the FE 2 -10/5-SchurMKL formulation results in more iterations per load step, which is the most influencing factor. The density of the micro mesh influences the total computational time; however, the relation between the sensitivity based formulation and the Schur complement based formulation remains the same.
The effect of implementation, micro mesh density and material model on computational time is for the MIEL formulation presented in Table 4. The example introduced in Sect. 5.1 is in this case solved with the macro mesh density 10 × 2 × 2 in 5 load steps with a constant load increment λ M = 0.2. The computational time is normalized with respect to MIEL-5/1-Sens. formulation, which, for micro mesh 5 × 5 × 5, took 37 s of real time. Two micro material models are considered: finite strain elasto-plastic model as defined in Sect. 1.1 and hyper-elastic model based on hyper-elastic strain energy (2). For a sparse micro mesh (5×5×5) the Schur complement based formulation is faster than the second order sensitivity analysis based formulation. The advantages of the sensitivity based implementation become apparent with denser micro meshes (30 × 30 × 30). As already shown in Fig. 4, the cost of the calculation of the Schur complement grows much faster with the density of the micro mesh than the cost of the second order sensitivity analysis. While the cost of the Schur complement does not depend on the material model used, the cost of sensitivity analysis does. Consequently, the difference between the numerical efficiency of the Schur and sensitivity based formulations is greater for the hyper-elastic material model than for the elasto-plastic material model.

Effect of non-linearity of micro-structure
The example demonstrates how the use of a two-level path-following procedure improves the numerical efficiency of multi-scale simulation in the case of highly nonlinear

(b) (a)
microstructure response and relatively monotonic response of the macro structure. To demonstrate that, a 2D, plane strain, uni-axial test is simulated using FE 2 multi-scale method based on a fully consistent sensitivity analysis. The macro geometry and mesh together with the RVE geometry and mesh are shown in Fig. 8. The macro domain is discretized with 4 × 2 macro elements and displacement u max = 0.6 is prescribed at the end. The RVE of periodic micro-structure is composed of hyper-elastic rim with material properties E = 21,000, ν = 0.3 and a narrow elasto-plastic incision with properties E = 21,000, ν = 0.3, σ y0 = 24. The incision has an additional small imperfection in the middle. At RVE level Q2, nine nodded, isoparametric, quadrilateral, plane strain elements were used, to avoid the locking effect. At a certain load level, a strongly nonlinear process of necking of the incision starts and requires very small load steps. If n m = 1 (FE 2 -Adaptive/1-Sens. approach), the maximum micro level load increment limits the macro load increment resulting in small macro steps. Due to the elastic rim, the global response remains relatively unaffected. If an adaptive path-following procedure is applied also at the micro level (FE 2 -Adaptive/Adaptive-Sens. approach), significantly larger load steps can be performed at the macro level. The FE 2 -Adaptive/Adaptive-Sens. approach needs 13 macro load steps, whereas for the FE 2 -Adaptive/1-Sens. approach 33 macro load steps are done. Figure 9a shows the macro reaction force F and Fig. 9b the absolute contraction at point B at the micro level with respect to global load factor λ M for both cases. The response curves are almost the same for both cases. Thus, an efficient solution of strongly nonlinear problems requires two level adaptive time stepping procedures where the maximum size of load increments at the micro level determines the overall efficiency of the simulation.

Effect of path dependency of microstructure
The example demonstrates how the use of two-level pathfollowing procedure improves the numerical efficiency of multi-scale simulation in the case of strongly path-dependent problems. The accuracy of the integration of evolution equations depends on the micro step size, thus limiting the size of micro load steps. This, for the n m = 1 case, limits also the macro load step size similarly as in the previous example. Again, the two-level adaptive path-following procedure proves to be numerically more efficient than the standard approach where each macro step is followed by one micro step.   Long clamped beam with dimensions 20 × 1 and with macro mesh division 80 × 4 has prescribed vertical displacement v max = 0.25 at the middle, as shown in Fig. 10. The beam is perforated with 320 perforations with the radius that gives 30% perforation of the beam. Perforations are evenly distributed in a way that each perforation is placed at the center of the corresponding macro element. Two cases are investigated. In the first case, a MIEL multi-scale computational scheme is employed. Due to the even distribution of perforations, all the MIEL micro meshes look the same, as shown in Fig. 10a. In the second case, infinitely small perforations were assumed with the same 30% perforation ratio as in the first case. The second case is simulated by the FE 2 scheme with RVE, as depicted in Fig. 10b. The RVE mesh is identical to the MIEL micro mesh due to the evenly distributed perforations. Material properties of the microstructure are in both cases E = 21,000, ν = 0.3, σ y0 = 24, K h = 21, R ∞ = 12 and δ = 30. The value of strain tensor component E xy in point A is for various solution strategies compared for the MIEL scheme in Fig. 11 and for the FE 2 scheme in Fig. 12.
In Fig. 11a the response curve E A xy (λ M ) is shown for the MIEL-Adaptive/1-Sens. approach with different prescribed maximal size of macro load step λ M max , adaptive time stepping at macro level and one micro step per each macro step. Converged solution is achieved for λ M max = 0.01. Secondly, Fig. 11b displays the results for fixed λ M max = 0.2 and 1, 2, 5 and 10 micro steps per each macro step. It can be seen that the evolution equation integration error is significantly reduced with the increased number of micro steps, without the need for costly additional macro steps. There is, of course, a limit to which additional micro steps can improve the overall results, as shown in Fig. 11b.
For the FE 2 scheme the results are compared in the same way as for the MIEL scheme. In Fig. 12a the response curve E A xy (λ M ) is shown for FE 2 -Adaptive/1-Sens. with respect to the prescribed maximal size of macro load step λ M max and in Fig. 12b for λ M max = 0.2 and with different number of micro steps. Adaptive time stepping at macro level is used in all cases. Conclusions are the same as for MIEL. With a two-level path-following scheme, the same accuracy is achieved with 20-times fewer macro steps. With additional micro steps, the method was able to capture also fine details of response curve near λ M = 0.2.
Point A is in the corner, close to the boundary where deformation gradients are high. Consequently, the converged curve E A xy (λ M ) is different for the MIEL and the FE 2 scheme.

Conclusions
In the paper a generalized essential boundary condition sensitivity analysis based implementation of FE 2 and MIEL multi-scale methods was derived and described in detail, as an alternative to more traditional implementations of multi-scale analysis, where the calculation of the Schur complement of the microscopic tangent matrix is needed for bridging the scales. The paper shows that the derivation of algorithmic macroscopic tangent requires for the FE 2 multi-scale method the first order essential boundary condition sensitivity analysis and for the MIEL multi-scale method the second order essential boundary condition sensitivity analysis. Thus, a general automatic differentiation based formulation (ADB) is introduced for the first and second order essential boundary condition sensitivity analysis that can be applied on an arbitrary coupled, path-dependent micro material model or element formulation. It has been shown in the paper that ADB brings several advantages. The first advantage is the ability to naturally introduce a fully consistently linearized two-level path-following algorithm as a solution algorithm for the multi-scale modeling. Sensitivity analysis allows that each macro step can be followed by an arbitrary number of micro substeps while retaining quadratic convergence of the overall solution algorithm. Using examples, it has been shown that this cannot be achieved with the standard Schur complement based methods. Additionally, the method completely avoids the evaluation of the Schur complement of the micro tangent matrix as numerically demanding mathematical operation which, especially for the MIEL multi-scale methods, results in large dense matrices.
The advantages of the sensitivity analysis based implementation in comparison with the traditional one were recognized and verified on numerical examples. The convergence of results with respect to the size of the macro load step and the number of substeps at the micro level was investigated. With additional micro steps, the accuracy of the global response can be improved without costly additional macro steps. This is especially evident in the case of a strongly nonlinear micro response, which, for some reason, does not reflect the global response, but it still limits the size of the macro load steps. Similarly, a strongly path-dependent micro material model requires small micro load steps that limit within the standard implementation also the size of the macro load step. This restriction is again relaxed with the introduction of the two-level path-following algorithm.