A comparison of the primal and semi-dual variational formats of gradient-extended crystal inelasticity

In this paper we discuss issues related to the theoretical as well as the computational format of gradient-extended crystal viscoplasticity. The so-called primal format uses the displacements, the slip of each slip system and the dissipative stresses as the primary unknown fields. An alternative format is coined the semi-dual format, which in addition includes energetic microstresses among the primary unknown fields. We compare the primal and semi-dual variational formats in terms of advantages and disadvantages from modeling as well as numerical viewpoints. Finally, we perform a series of representative numerical tests to investigate the rate of convergence with finite element mesh refinement. In particular, it is shown that the commonly adopted microhard boundary condition poses a challenge in the special case that the slip direction is parallel to a grain boundary.


Introduction
Crystal plasticity is the accepted model framework for incorporating microstructural information in continuum theory with application to crystalline, metals where dislocations constitute the mechanism behind inelastic deformation. In order to account for the size effects, e.g. the Hall-Petch effect (cf. Hall [18] Petch [26]), due to the existence of grain boundaries in a polycrystal, it is convenient to include some sort of gradient-extension of the flow properties along the slip direc-B Kristoffer Carlsson kristoffer.carlsson@chalmers.se 1 Department of Applied Mechanics, Chalmers University of Technology, Gothenburg, Sweden tions, either in the dragstress or backstress. These could come from for example, geometrically necessary dislocations, cf. Ashby [1], from here on referred to as GND's, which are generally of two types, edge and screw dislocations. Various explicit models based on this conceptual background have been proposed, for example in a series of papers by Gurtin and coworkers [13,15,17]. An attempt to unify and generalize several models for gradient enhanced crystal plasticity (within the framework of large deformations) was presented by Svendsen and Bargmann [28], who also compared the models from a conceptual viewpoint. However, several modeling issues still await their resolution. An elegant way of unifying gradient theory for different application models, including inelasticity, damage and phase-field models, was presented by Miehe [22].
The so-called primal format uses the displacements and slip of each slip system as the primary unknown fields (see e.g. Bittencourt et al. [6], Borg [7], Okumura et al. [25]). An alternative format is coined the semi-dual format, which, in addition, includes the microstresses as primary fields and thereby defines a mixed variational problem where the explicit dependence on the gradient fields are avoided. Similar mixed formulations have been used by Bayley et al. [4], Yefimov et al. [29], Evers et al. [11], but typically in these references, the GND's are the primary variables instead of the microstresses. A comparison of different FE-approximations for use in crystal plasticity was carried out by Kuroda [20] using the GND's as primary fields. We note that the mixed method used extensively in our research group in recent years, see Bargmann et al. [2,3], Ekh et al. [8,9], bears resemblance with the semi-dual format, however, without possessing a well-defined variational structure.
In this paper we focus on issues related to the theoretical and computational format of material models with gradient variables. As a prototype model we choose crystal visco-elasticity of the Norton type with kinematic gradient hardening. More specifically, the back-stress resulting from GND's (representing kinematic hardening) is of purely gradient origin. We compare the primal and semi-dual variational formats in terms of pros and cons from modeling as well as numerical viewpoints.
The paper is outlined as follows. A class of gradientenhanced dissipative materials is considered in Sect. 2 and are presented in the time-continuous and time-discrete formats. The resulting space-continuous boundary value problem is described in Sect. 3. In Sects. 4 and 5, we establish the primal and semi-dual weak format as well as the resulting Finite Element problem. In particular, the semi-dual format is established in terms of a micro-stress energy that is obtained after a partial Legendre transformation of the free energy. We show that, for a common choice of the gradient hardening modulus, it is possible to decompose the residual equation including the vectorial micro-stress into one and two scalar equations, for a problem in two and three dimensions, respectively. In Sect. 6, we show computational results for both a single crystal and a polycrystal. The convergence rate with respect to mesh refinement is investigated for the unknown fields in different norms, and the results for the primal and semi-dual formulations are compared. In particular, we investigate the degenerate case when the slip direction of a slip system is parallel with a grain boundary. Finally, concluding remarks and an outlook to future work are given in Sect. 7.

Notation
Throughout the paper we use direct notation. Contractions are defined by where a, b are vectors (rank 1 tensors), A, B are rank 2 tensors and E is a rank 4 tensor and repeated indicies are summed over (Einstein summation convention). Outer products are defined by the relations The vector product between two vectors is defined by The symmetric part of a rank-2 tensor is expressed as where E e is the elasticity tensor, = (u ⊗ ∇) sym is the strain and p is the plastic strain. The plastic deformation is modeled by the slip γ α on the slip system α, where s α and m α are the slip direction and the normal to the slip plane, respectively. These slip directions and slip normals are assumed to stay fixed (unaffected by the deformation). We treat γ = (γ 1 , γ 2 , . . . , γ M ) as a set of internal variables. 1 We shall investigate the inherent properties of the kinematics associated with crystallographic slip in some further detail. To this end, we introduce the Cartesian axes (defined by three orthogonal unit vectors) as (m α , s α , k α ) where k α def = m α × s α , for any slip system. We may then decompose the gradient of an arbitrary function f as where we have defined the tangential gradient operator ∇ α tan , operating on a function f as Next, we introduce Nye's dislocation density tensor G (cf. [24]) as follows: We may derive an alternative representation of G α as where ρ GND are the edge dislocation and screw dislocation densities, respectively (Fleck et al. [12]). Here, b is the magnitude of the Burgers vector. The decomposition in (7) thus represents continuous distributions of the edge and screw dislocations. In the literature, these densities are called geometrically necessary dislocations (GNDs). We note that G K def = −G T is Kröner's dislocation tensor ( [19]).

Crystal visco-elasticity with kinematic gradient hardening
Associated with each slip system, we introduce kinematic hardening that is of purely gradient origin and is represented by the gradient variables g α def = ∇γ α . The free energy density ψ is then proposed as the additive decomposition ψ , γ , g = ψ e , γ + ψ g g (9) where ψ e is the contribution from elastic (stored) energy, whereas ψ g are the contributions from gradient hardening in the slip systems. They are defined as where we introduced the gradient hardening tensors H gra α . We note that the gradient contributions vanish when the length parameters l α tend to zero.
Remark A more complex expression of ψ g was proposed in, e.g, Bargmann et al. [3]: Omitting gradient cross-hardening i.e. restricting the sum in (11) to α = β, we obtain This expression also represents a special case of (10b), c.f. Sect. 2.3.
The (equilibrium) stress is given by A dual dissipation function φ * α (τ di α ), that is associated with each slip system, is chosen as where τ di α are (scalar) "dissipative microstresses" that are energy-conjugated toγ α , F is an overstress function and t * is a material parameter corresponding to a relaxation time. The evolution rules for γ α thus reaḋ where the plastic multiplier λ α is given as Remark For the numerical examples later, we choose where C is a reference stress and m is the Norton exponent. We note that increasing the value of m makes the response tend towards perfect plasticity with a yield stress C.
In order to "close" the problem formulation, we need to establish the Biot equations, see e.g. Biot and Romain [5], Nguyen [23] and Ziegler and Wehrli [30], which in this case read where the "energetic microstresses" are defined as Equation (18) is sometimes denoted microforce balance; Gurtin [14]. We note that τ α = σ : [s α ⊗ m α ] is the classical Schmid (resolved) stress. Henceforth, we adopt the more concise notation ξ α = ξ en α corresponding to the absence of a dissipative counterpart. From Biot's equations (18), we obtain the relation (20) where −χ α has the form of a back-stress. In view of (14), it is clear that the gradient effect embedded in τ di α represents the sole source of "back-stress". In fact, we may rewrite (14) more explicitly as

Crystal visco-elasticity with kinematic gradient hardening: special case
A simple choice of the gradient hardening tensors H gra α (see e.g. Gurtin et al. [17], Gurtin [16], Erturk et al. [10]) is by which the gradient part of the free energy, given in (10b), is replaced by where we note the identities g ⊥,α = −bρ GND ⊥,α and g ,α = bρ GND ,α . Clearly, this simple expression does not allow for cross-hardening. We may now express ξ α as where we introduced the "edge" and "screw" components of ξ α as follows: Furthermore, we may decompose χ α as

Time-continuous strong format
We consider a body occupying the domain and, for simplicity, restrict the problem formulation to quasistatic conditions. The strong format of the coupled problem of finding the fields u(x, t) : The boundary conditions on the body with boundary are defined as follows: Remark Each of the given boundary conditions in (28b) is of either the Dirichlet (essential) or Neumann (natural) type. The appropriate classification depends on the actual choice of variational format, as discussed further below.

Time-discrete strong format
Adopting the standard Backward Euler rule for integrating the system in (27) over the time interval (t n , t n+1 = t n + t), we obtain the solution for the time-updated fields u(x), (29c)

Primal variational format 4.1 Time-discrete primal variational format
In order to establish the proper variational format, we introduce the solution and test spaces 3 : i.e. superindex n+1 is dropped for brevity of notation. 3 Superscript 0 denotes that functions are homogeneous on the Dirichlet boundary, e.g.
The primal variational format corresponding to the strong format in (29) becomes: Find u ∈ U, γ α ∈ G α , and τ di α ∈ T that solve where we introduced the data in terms of the linear functionals As to the boundary conditions, those given on the boundary part In both these special cases the boundary integral on the RHS in (31b) vanishes. Next, we introduce the volume-specific incremental potential π as follows: The global incremental potential = int + ext becomes The stationarity conditions of (û,γ ,τ di ) are those in (31).

Nested iteration strategy: primal variational format
Firstly, we choose to satisfy Eq. (31c) in the strong sense, i.e. we introduce the local residual in each spatial point x ∈ and we note that r (τ ) α (γ α , τ di α ) = 0 can be solved for τ di α {γ α } 4 for any given γ α . The (remaining) global residuals corresponding to Eqs. (31a), (31b) are given as In order to obtain (36), we used the functional dependence of σ , τ en α and ξ α as follows: where we introduced the new operators Global equations: Local equations: The problem in (31) is thus reformulated as that of finding the fields u ∈ U, γ α ∈ G α , and τ di α ∈ T that solve together with the local implicit relation τ di This is a non-linear system of equations, that in general requires some sort of iteration strategy to solve. The most straightforward approach would be to solve for all the fields (u, γ , τ di ) in each iteration step in a "monolithic" fashion as is done for gradient plasticity by Liebe et al. [21]. However, due to the multiple slip systems in crystal plasticity models this approach would lead to a large number of global degrees of freedoms and, consequently, the strategy would be computationally heavy. Hence, it is convenient to resort to so-called nested iterations, whereby a two-level procedure of Newton iterations is designed as follows: • Global iterations Iterations are carried out on the global fields z = (u, γ ), while assuming that the constitutive equations that govern the fields τ di α {γ α } are satisfied when carrying out linearization of (39) in order to obtain the pertinent Algorithmic Tangent Operators (see below).
• Local iterations For given value of in the global iteration k that is thus known in each spatial point x, iterations are carried out on the constitutive Eq. (40) in order to obtain τ di α for given value of γ α , We remark that the particular model adopted in this paper, (17), is simple enough to allow a closed form solution for τ di α (γ α ); hence, no local iterations are needed.
A schematic figure of the nested iteration strategy and the interaction of information between the global and local levels is shown in Fig. 1.

Linearization and algorithmic tangent tensors
The system (39) is linearized while it is assumed that the local Eq. (40) are satisfied exactly. These latter conditions are used in order to compute the pertinent sensitivities of each τ di α for a perturbation of γ β . These sensitivities are defined as the directional derivatives dτ di We note that the sensitivities are uncoupled in the sense that (dτ di α ) γ β = 0 only when β = α. The sensitivities are computed from a linearized form of the residual equation in (40), as shown in Appendix 1.

Global Newton iterations
For given value z (k) = (u (k) , γ (k) ) in the global iteration step k, we compute updated values The pertinent linearized forms are given as The FE-discretized versions of (36) and (44) are given in Appendix 1.

Local iterations
For given values of the fields z = (u, γ ), or only γ (in fact), in any given spatial point x, we compute the corresponding values of τ di . For the model in (14), the pertinent constitutive problem becomes: For given γ α , compute τ di α from The solution is obtained conveniently in two steps: Firstly, |τ di α | is solved from the equation Remark For the specific power law (Norton type) model in (17), we obtain the explicit solution Secondly, the sign of τ di α is determined as follows:

Time-discrete semi-dual variational format
In the (semi-)dual variational format we exploit a partial Legendre transformation of the free energy density ψ w.r.t. to the gradient variables g α . Our aim is to obtain the compliance format of the relation in (19b), which we repeat here, In this section, we consider the situation that H gra α is positive definite, such that it is possible to invert (48) in straightforward fashion. It is then possible to define the Legendre transformation and define the semi-dual free energy density as The stationarity condition of (49) gives the constitutive equations for g α expressed as where H * gra α is the inverse of H gra α in the classical sense, i.e. it satisfies the condition H * gra α · H gra α = I. Upon inserting (51) into (49), we obtain Next, we introduce the appropriate solution and test spaces The dual variational format can thus be stated as follows: find u ∈ U, ξ α ∈ X α , γ α ∈ T, and τ di α ∈ T that solve where we introduced the data in terms of the linear functionals Remark When l → 0, the gradient part of the free energy vanishes. Hence, the expression in (49) will render ψ * g to enforce ξ α = 0. Consequently, the entire Eq. (54b) and the divergence term in (54c) are removed, and the system of equations degenerates to exactly those arising in standard non-gradient crystal plasticity where displacements are solved for in the global iterations and the slips are solved for in the local iterations. In the primal formulation, however, the slip fields would remain global.
As to the boundary conditions, those on boundary parts (γ ) t are now of the Dirichlet type, whereas those on (γ ) u are of the Neumann type. The difference from the primary format is that the type of the boundary conditions on the different parts of have switched their roles. We note that in a finite element implementation of the dual format, enforcing microfree conditions ξ α · n = 0 would require a local coordinate transformation on the form ξ α = ξ α,n n+ξ α,t 1 t 1 +ξ α,t 2 t 2 , where t 1 and t 2 are two orthogonal vectors in the plane with normal n. The condition ξ α · n = ξ α,n can then be enforced strongly.
The corresponding volume-specific incremental potential π is introduced as follows: where we recall the notation The stationarity conditions of (û,ξ ,γ ,τ di ) are those in (54). Comparing (56) with (33), we note that ψ has been replaced by ϕ and that an additional term involving χ α is present in (56).

Time-discrete semi-dual variational format: special case of singular gradient hardening
The special case that H gra α is of the form given in (22), and thus singular, requires special care. When the gradient part of the free energy ψ g α takes the form in (23), it is possible to replace the Legendre transformation in (49) with the problem of reduced dimensionality: We define the semi-dual free energy density as The stationarity condition of (58) gives the constitutive equations Upon inserting (60) into (58), we obtain In this case the appropriate solution and test spaces are defined as The dual variational format thus becomes: Find u ∈ U, ξ ⊥,α ∈ X ⊥,α , ξ ,α ∈ X ,α , γ α ∈ T, and τ di α ∈ T that solve where we introduced the data in terms of the linear functionals The corresponding volume-specific incremental potential π is introduced as follows: where we recall the notation χ ⊥,α = ∇ξ ⊥,α · s α and χ ,α = ∇ξ ,α · k α . The global potential = int + ext becomes ,α (ξ ,α ) .
Remark It is interesting to note that the formulation in (63) becomes identical to the "mixed" formulation, originally proposed by Svedberg and Runesson [27] and used by Ekh et al. [8] in the context of crystal plasticity, when screw dislocations are ignored (or the analysis is restricted to 2D, in which case screw dislocations never appear). In such a case (63c) becomes obsolete. To show the identity, we recall the relations 5 Furthermore, we assume that c α is homogeneous. As a result, we may rephrase (63b), (63d) as where we introduced the space of test functions These two equations are precisely those used by Ekh et al. [8].

Nested iteration strategy: dual variational format
Let us now return to the more general situation where H gra α is a constant, positive definite, tensor. We first choose to satisfy (54c), (54d) in the strong sense, i.e. we introduce the local residuals in each spatial point: 5 The abbreviated notation g α def = g ⊥,α , ξ α def = ξ ⊥,α is used here.
[u], χ [ξ] γ{u, ξ} The remaining residuals in (54) are treated as global: The problem in (54) is thus reformulated as that of finding u ∈ U, ξ α ∈ X, γ α ∈ T, and τ di α ∈ T that solve together with the implicit relations γ {u, ξ } and τ di α {u, ξ } satisfying Like in the primal format, the system in (72), (73) is conveniently solved via a nested iterations strategy: • Global iterations Iterations are carried out on the global fields z = (u, ξ ), while assuming that the constitutive equations that govern the fields γ , τ di are identically satisfied when carrying out linearization of (72) in order to obtain the pertinent Algorithmic Tangent Operators. • Local iterations For given value of z (k) in the global iteration k, iterations are carried out on (73) in order to obtain γ , τ di for given value of u, ξ , i.e. γ α = γ α {u, ξ }, A schematic figure of the nested iteration strategy and the interaction between the global and local level is shown in Fig. 2. We note that, for the dual formulation, the fields τ di α are purely local.

Linearization and algorithmic tangent tensors
The system (72) will be linearized while it is assumed that the local Eq. (73) are satisfied exactly. We thus compute the sensitivities of each γ α and τ di α for perturbations of u, ξ . These sensitives are defined as the directional derivatives The sensitivities in (74) are computed from linearized forms of the residual Eq. (73), as shown in Appendix 1. In (75), G α is a tangent modulus defined in (81), whereasẼ αβ is a symmetric matrix defined in (86). Furthermore, we recall the definition of D e β in (38a).

Global Newton iterations
For given value z (k) = (u (k) , ξ (k) ) in the global iteration step k, we compute updated values whereby the increment z = ( u, ξ ) ∈ Z 0 def = U 0 × X 0 is solved from the system The pertinent linearized forms are given as (78d) Remark The symmetry, which is inherent in the use of a variational principle, is shown explicitly in (78).
The FE-discretized versions of (72) and (78) are given in Appendix 1.

Numerical results
In this section we benchmark the primal and semi-dual formulations by comparing the respective solutions for two different problems: a single crystal and a polycrystal system. The analyses are performed in two dimensions, with a plain strain assumption, and the simplified format for the gradient hardening tensors H gra α in Eq. (22) is adopted. The pertinent variational format used for the semi-dual formulation is given in Sect. 5.2. Regarding FE-interpolation, piecewise linear interpolation of all global fields on triangular meshes are used throughout. Table 1 shows the parameter values that were used. The values l α and H ⊥,α are chosen to be equal for each slip system which means that the subscript α is obsolete, and is therefore henceforth dropped. We remark that the internal length is given in dimensionless form as l/L, where L represents the length scale of the grain structure, e.g. as shown for the single grain in Fig. 3a. However, it is important to note that the relevant measure for the numerical results is l 2 H ⊥ (and not the individual values of l and H ⊥ ).

Single crystal
The first problem is a polyhedral single crystal with Dirichlet boundary conditions on the displacement field and microhard  boundary conditions on the slip field. The shape of the crystal with an imposed mesh is shown in Fig. 3b. The purpose of the analysis is to investigate how the solutions for the global and local fields converge with mesh refinement. In the primal formulation the slip fields γ α are global fields and the gradient fields g α are local, while in the semi-dual formulation the slip fields are local and the microstress fields ξ α are global. Consequently, the slip fields and the fields including the gradient of the slip are approximated differently and are, therefore, expected to converge at different rates.
We introduce two different measures of the slip fields: the norm of the slip fields denoted ||γ ||, and the weighted norm 6 of the gradient fields, denoted ||g|| H , as follows: 6 For the singular H gra α introduced in (22), which is used for the numerical investigation, ||g|| H is actually a seminorm. However, this property has no relevance for the results.
We note that the weighted norm pertains to the gradient part of the free energy when applied to the solution. Henceforth, we study the errors ||γ − γ h || and ||g − g h || H , where γ h and g h are finite element approximations. In practice, the exact solution is replaced by an overkill solution on a sufficiently fine mesh. Since the overkill and the FE-solution generally live on different meshes it is not obvious how the subtraction between these fields should be carried out. The method used in this paper is to transfer the solution from the coarse mesh to the overkill mesh using the element interpolation functions The displacement boundary condition on the crystal is chosen asū(x) = 0.01x 2 (e 1 + e 2 ). This corresponds to vertical stretching combined with shear (on the average). The inhomogeneous Dirichlet boundary condition was applied incrementally up to a time t = 10 s. The crystal contains two slip systems with the slip directions oriented 20 • and 40 • , respectively, counterclockwise from the x 1 axis, c.f. Fig. 3a. Figure 4 shows the slip field for slip system 1 taken along a vertical line between the points (0.25, 0) and (0.25, 0.15), as shown in Fig. 3a. In the primal formulation the slip is identically zero on the boundary of the crystal, while for the dual formulations the microhard condition is not fulfilled exactly. In fact, the smaller the internal length scale the more the microhard condition is violated in the dual formulation. This is an effect of the strong vs. weak enforcement of the microhard condition in the two formulations.
Next, we consider convergence of the slip and its spatial gradient with mesh refinement. We then recall that the gradient is given as g α def = ∇γ α in the primal format, whereas g α is an independent global field in the semi-dual format. Figure 5 shows that ||g|| H converges from above in the primal format and from below in the semi-dual format. This behavior mimics linear elasticity based on potential energy and complementary energy, respectively.
The convergence of the slip and gradient fields is illustrated in Fig. 6. We denote the rate of convergence by the value of the exponent p in the expression k(h/L) p where h is a characteristic length of the elements in the mesh. The slip converges at approximately twice the rate in the primal formulation as compared to the dual formulation. On the other hand, the gradient of the slip converges at approximately twice the rate in the dual formulation as compared to the primal formulation.

Parallel slip direction and grain boundary with microhard boundary conditions
It is of interest to investigate the degenerate case when the slip direction and the grain boundary are parallel. To this end, we consider a square-shaped polycrystal consisting of three grains with boundaries oriented at 45 • to the x 1 -axis. Each grain has a single slip system with the same slip direction for all grains.  Fig. 7. As the slip directions of the grain boundary and slip system coincide, the grain boundary becomes "invisible" in the dual formulation due to its weak enforcement of the microhard condition. For the primary formulation, the fact that the microhard boundary condition is enforced strongly means that the slip will always be identically zero along the grain boundary. It is of interest to see if, and how much, this fact influences the error in a homogenized quantity. The convergence in the standard L 2 -norm of the deviatoric plastic strain with mesh refinement is displayed in Fig. 8 . For the primal formulation the error is quite insensitive to the slip direction, while for the dual formulation it matters more. As the difference between the direction of the grain boundary and the slip direction tends towards zero, the error increases for the primal formulation while it decreases for the dual formulation. This can be explained by the following reasoning: For the primal formulation, in the case of complete alignment between the grain boundary and slip direction, even though there is no gradient coupling perpendicular to the slip direction in the equations, the entire element adjacent to the grain boundary will be affected by the hard (Dirichlet) boundary condition. To reduce the contribution of this element the mesh has to be refined such that the elements connected to the grain boundary are small. For the dual formulation, the effect of alignment is instead that the grain boundary no longer has an effect on the solution. Recalling Eq. (64b), we see that when the slip direction and grain boundary are parallel (s α · n = 0) the load l (ξ ) ⊥,α is identically equal to zero. Large elements at the boundary do, therefore, not induce a large error in the result.
While there is a large relative difference in how the two formulations behave, the absolute error is still small so this effect can likely be concluded to be insignificant from a homogenization perspective.

Conclusions and outlook
In this paper we have presented two different formulations of a gradient-enhanced crystal viscoelasticity model: The primary format uses the displacements and slips as global fields, whereas the dual format uses the displacements and slip gradients as global fields. When gradient effects are disregarded, the residual equations for the dual format degenerate to those typically used in local, i.e. non-gradient, crystal plasticity.
For a common special choice of the gradient hardening tensors H gra α , the gradient part of the energy is quadratic in terms of the GND densities. In this case, it turns out that the residual equations in the dual format for the microstress ξ α can be reduced from a vectorial equation to one or two scalar equations in two or three dimensions, respectively, thereby saving computational time.
As to the numerical examples, we have examined the convergence rate for mesh refinement and showed that the gradient fields converge faster in the dual formulation, while slip fields converge faster in the primal formulation. Finally, we looked at the degenerate case where the grain boundary and the slip direction align. It was found that the error in the deviatoric plastic strain became significantly smaller in the dual formulation than in the primal formulation. The discrepancy was explained theoretically from the difference in type of the boundary conditions for the two formulations. The absolute error was, however, small for both formulations and it is, therefore, likely that the difference in the results for the two formulations can be ignored in practice.
Future work will focus on extending the investigation to three dimensions and thereby incorporating the influence of the screw dislocations. Further investigations should also target the upscaled / homogenized properties of polycrystals on length scales of engineering interest.
Next, consider the system (84). Proceeding in the same manner as above, we obtain dτ di and The final result is

Primal format
The displacement field u(x) and slip fields γ α (x),are approximated in finite dimensional subspaces: Upon inserting (92) into (36), we obtain the discrete residuals in matrix form where we introduced B e αβ def = E e αβ + δ αβ H α + A α = B βα and where A α are the (symmetric) tangent matrices Upon introducing the abbreviated notation: we arrive at the condensed system (107)

Dual format
In the same fashion as for the primal format in Sect. 1 the fields u(x) and ξ α (x) are approximated by: Upon inserting (108) into (71), we obtain the discrete residuals in matrix form R (u) a (u) , a where we introduced the matrices Upon using Newton's method to solve for the increments a (u) , a (ξ ) α in a given iteration, we obtain the system