Shape-Programming in Hyperelasticity Through Differential Growth

This paper is concerned with the growth-driven shape-programming problem, which involves determining a growth tensor that can produce a deformation on a hyperelastic body reaching a given target shape. We consider the two cases of globally compatible growth, where the growth tensor is a deformation gradient over the undeformed domain, and the incompatible one, which discards such hypothesis. We formulate the problem within the framework of optimal control theory in hyperelasticity. The Hausdorff distance is used to quantify dissimilarities between shapes; the complexity of the actuation is incorporated in the cost functional as well. Boundary conditions and external loads are allowed in the state law, thus extending previous works where the stress-free hypothesis turns out to be essential. A rigorous mathematical analysis is then carried out to prove the well-posedness of the problem. The numerical approximation is performed using gradient-based optimisation algorithms. Our main goal in this part is to show the possibility to apply inverse techniques for the numerical approximation of this problem, which allows us to address more generic situations than those covered by analytical approaches. Several numerical experiments for beam-like and shell-type geometries illustrate the performance of the proposed numerical scheme.


Introduction
Soft robotics is a biologically-inspired groundbreaking technology that aims to mimic mechanical deformations, which take place in humans, animals, or plants, through actuated soft materials: dielectric elastomers or magneto-active polymers, for instance.Several actuation mechanisms, such as fluidic, heat, electric, or magnetic, may be used to control these materials [1].The range of potential applications of this new generation of robots includes, among many others, medical assistance and ocean exploration [2,3].
In addition to the development of new manufacturing technologies, mathematical modelling, analysis, and numerical simulation are tools of paramount importance to speed up progress in this field.Being composed of soft matter, nonlinear continuum mechanics is the appropriate physical theory to model the kinematics and dynamics of these materials.However, the mathematical control theory of hyperelastic materials is scarce.Indeed, the first mathematically rigorous study for control problems in the hyperelasticity setting appears to be [4].More recently, several papers have addressed the control of soft materials from the viewpoints of mathematical analysis and numerical simulation [5][6][7][8][9].See also [10] for a recent survey.
Growth is another biological process susceptible to being mimicked by artificial soft materials.As a matter of fact, [11] reports on a class of soft pneumatic robots whose movements are driven by growth.As for the mathematical modelling of growth, A. Goriley provides, in his seminal book [12], the required ingredients.Doubtless, the topic of mathematical analysis and numerical simulation of growth control is in its infancy, insofar as the mathematical analysis of soft materials actuated by growth is missing in the literature, and only a few works have addressed the numerical simulation counterpart.In this regard, it is worth mentioning [13,14].Both papers tackle the socalled shape-morphing problem, where the goal is to find the growth tensor that can produce the deformation of a given soft continuum to a desired shape.The paper [13] copes with the complexity of the activation as well, and provides explicit solutions in the case of affine shape changes.In a complementary manner, [14] focuses on the case of shells and also finds analytical solutions under the stress-free assumption.
Fostered by [13,14], this paper sets up the shape-morphing problem within the framework of optimal control theory.Indeed, the control variable is the growth tensor.We consider both cases in which the growth tensor is globally compatible, meaning that it is a deformation gradient over the undeformed domain, and the incompatible one, where it is no longer a gradient.From the viewpoint of mathematical modelling, the former case is expressed as the composition of two mechanical deformations: one of them accounts for growth and the other one incorporates boundary conditions and other possible effects like external loads.The latter case relies on the theory of Morphoelasticity, where a local elastic tensor restores the compatibility that is lost by the growth tensor.The state variable is the deformation of the actuated soft continuum.As usual in hyperelasticity theory, that deformation is a minimiser of a polyconvex energy functional.The cost function uses the Hausdorff distance to account for dissimilarities between the desired shape and the final configuration.It also includes a term to deal with the complexity of the activation.In our work, we do not neeed the stress-free hypothesis of [14].
The outline of this paper is as follows.Section 2 contains the modelling details.Section 3 performs a rigorous mathematical analysis of the shape-programming problem in the globally compatible case, which is more involved analytically than the incompatible one.More precisely, firstly, we prove that, for a given growth tensor, there exist minimisers of the underlying energy functional.Secondly, we establish the existence of solutions for the optimal control problem.We rely on the Direct Method of Calculus of Variations in both cases.Section 4 straightforwardly extends these existence results to the incompatible case.Eventually, Sect. 5 addresses the numerical approximation of the shape-programming problem.Our purpose here is to show how inverse techniques may be used for the numerical resolution of this problem, thus addressing more generic situations than those covered by analytical approaches.We adopt a pragmatic point of view in this part as we are not concerned about the compatible or incompatible nature of the growth tensor; in practise, this amounts to accept the incompatibility and, hence, the study lies in the theory of Morphoelasticity.Likewise, we take the right Cauchy-Green deformation tensor as the main variable since, by a suitable parametrisation, it highly simplifies the numerical approximation of the problem.We also derive explicit formulae for the gradients of the functional involved, and transfer those to cutting-edge optimization algorithms that use gradients, in particular, the interior-point method, to obtain the desired solutions.Several numerical examples for beam-like and shell-type applications, as well as a problem converting a square into a circular geometry, illustrate the performance of the proposed numerical scheme.

Modelling Differential Growth in Nonlinear Elasticity
Let 0 ⊂ R N , N = 2, 3, be an open, bounded and connected domain which represents the reference (or undeformed) configuration of an elastic and soft body.If 0 experiences a growth effect (as happens in plants or in human tissues, for instance), then it changes its size or shape.
There are two ways to understand and model this phenomenon.In the first, we postulate that there is an underlying deformation that produces the growth.Let us then denote by g : 0 → R N the deformation mapping induced by this phenomenon, and by g := g ( 0 ) the deformed body once growth has taken place.It is assumed that g is a Sobolev map.Although driven by growth, g is still assumed to be a mechanical deformation, hence it satisfies the properties required to any such deformation; in particular, it preserves the orientation and does not interpenetrate [15].Let us denote by G = G(X) the deformation gradient tensor associated with g , i.e., G := ∇ g , where ∇ is the material gradient operator with respect to X ∈ 0 .The orientationpreserving condition is modeled with the constraint det G(X) > 0 for almost every X ∈ 0 , ( while the non-interpenetration is modelled by imposing that g is injective almost everywhere (hereafter abbreviated to a.e.), so that the restriction of g to the complement of a set of measure zero is injective [16].
For the second possibility, and according to the general modelling of growth and morphoelasticity in [12], the postulate that there is an underlying deformation g responsible for growth is discarded, so the tensor G is not assumed to come from any deformation, though still (2.1) is retained.
Since the first alternative is more involved analytically, we will keep our general discussion (this section and Sect.3) in that context, and defer some comments on the second one (Sect.4), once the main analysis has been performed.Even so, numerical experiments in Sect. 5 are explored in the morphoelasticity scenario.
Since g is an elastic and soft material, it has an internal elastic energy, which is able to induce a new deformation e on the body g .For this initial exposition of the problem, we can think that e is Lipschitz, but this assumption is not necessary in the analysis.Although, in principle, the elastic energy might depend on the configuration g and the growth deformation g , this is not the case in the current context, since g represents a growth that does not change the elastic properties of the material.Thus, we assume that the constitutive parameters of the body occupying 0 and g are the same.This assumption requires additionally that the body is homogeneous, i.e., its mechanical properties are the same at each point.This is modeled through an energy function that does not depend explicitly on material points and is the same for both configurations 0 and g , regardless of the growth deformation g .This stored energy function is denoted by W 0 : R N ×N + → R, where R N ×N + designates the set of square N × N matrices with strictly positive determinant.The precise assumptions on W 0 will be listed in Sect.2.3.
As is well known, equilibrium configurations e are minimisers of the functional (to which one may add external forces) over a suitable class of admissible deformations to be specified later.Here F e is the deformation gradient of the elastic deformation e .The variables in g have been denoted by Y , while the variables in 0 by X, so that Y = g (X).
Taking into consideration both growth g and elastic e deformation, the total deformation of the body 0 is expressed as the composition of both mappings, i.e., Fig. 1 The mapping between reference 0 and deformed configurations is expressed as the composition of the growth deformation g and the elastic deformation e = e • g .Accordingly, the deformation gradient tensor F associated with is given by F(X) = F e ( g (X))G(X), (2.3) where F e is the deformation gradient of e .The three maps involved, , e and g , are assumed to be orientation-preserving and injective a.e.(see Fig. 1 for a graphical representation).
By (2.1) and the fact that g is injective a.e., the change of variables Y = g (X) allows rewriting (2.2) in the undeformed configuration 0 as 0 where Note that the tensor G breaks the symmetries of W 0 .Indeed, if W 0 has a symmetry group (for example, it is isotropic), then W G does not, in general.Boundary conditions will be imposed in , but not in g or e independently.We will assume that the boundary 0 of 0 is Lipschitz and is decomposed into two disjoint parts: 0 D and 0 N , with 0 D of positive (N − 1)-dimensional area.On the Dirichlet part 0 D , it is imposed = ¯ for a given deformation ¯ : 0 → R N , while on 0 N we prescribe the Piola-Kirchhoff stress vector s 0 : 0 N → R N .The latter is not explicitly stated in the admissible set but it is automatically satisfied for minimisers when the surface energy term is added to the total energy.The cases 0 N = ∅ or s 0 = 0 are not excluded.In fact, volume forces can also be added, whose simplest form is linear: In view of (2.4), the total energy is Other boundary conditions are also possible (see, e.g., [17,Ch. 5]), as well as more general external forces.Finally, we fix an exponent p > 1 related to the growth at infinity of the function W 0 (see Sect. 2.3 for details) and define the class U of admissible deformations in (2.8) as where W 1, p is the notation for the Sobolev space.Naturally, we suppose that U is not empty and that G is not identically infinity in U, which amounts to assuming that ¯ ∈ U and G ( ¯ ) < ∞.

Setting of the Shape-Programming Problem
Having in mind potential applications in soft robotics, the so-called shape-programming problem [14] amounts to finding the growth tensor field G in the initial configuration 0 such that the final configuration is as close as possible to a desired target configuration target .Besides reaching this goal, and in order to facilitate its implementation in a possible soft robot, the computed growth tensor field should be as simple as possible.
Inspired by [13], a general form of a complexity functional should have a regularisation term (typically, a squared gradient of G) plus a term penalising the difference between the actual G and the target G target ; such target may well be the identity.These two ingredients should give rise to a simple growth tensor; indeed, the regularising term avoids oscillations, while the penalising term makes G similar to G target , which is chosen to be simple, too.In Sect.2.4 we will describe some possibilities of complexity functions, but for the moment we can think of the functional for α 1 , α 2 > 0, and develop the mathematical theory for general functionals of the form for a certain appropriate density φ.Here above, the norm of a second-order tensor A is defined by Similarly is defined the norm of vectors and third-order tensors.
We will see in Sect. 3 that, in order to prove existence of minimisers ∈ U of (2.8), one needs that G be in L ∞ and that det G is bounded away from zero.These facts are in agreement with the requirement that G is easily reachable.The most general way of expressing these assumptions is to fix a compact set K ⊂ R N ×N + and impose that G(X) ∈ K for a.e.X ∈ 0 .
Two relevant examples of the set K are for some M, m > 0, and Concerning the goal that the final configuration is as close as possible to a desired one, there are several options as for the distance between shapes.The simplest, but least realistic, is to consider an L 2 distance between the actual and a target deformation; the disadvantage of this choice is that, in general, a distance between deformations is only vaguely related to a distance between shapes.Introduced by [18] and analysed in [8], we consider instead the Hausdorff distance between the image domain and the target set as an adequate way of measuring distances between those sets.As in those works, we use in fact the following smooth approximation of the Hausdorff distance.Let (G) = ( 0 ) be the image domain and target the target domain.We fix three exponents α, β, γ > 0 and a continuous and strictly decreasing function ϕ : [0, ∞) → (0, ∞).We define (where |A| is the volume of the set A), so that an approximation of the Hausdorff distance is Putting all things together, the formulation of the shape-programming problem is: subject to: G ∈ H 1 ( 0 ; K ) is the gradient of an a.e.injective map, (G) = ( 0 ) , with ∈ U a minimiser of (2.8). (2.10)

Choice of the Energy Density
The classical assumptions in nonlinear elasticity for the energy function are polyconvexity and coercivity [17,19,20].To be precise, we will assume the following conditions for W 0 : (W1) W 0 is polyconvex, i.e., there exists a convex function e : If N = 2, the dependence on cof F can be dispensed with.(W2) There exist exponents p ≥ N − 1 and q ≥ N N −1 with p > 1, and a constant c > 0 such that Of course, cof denotes the cofactor matrix.Conditions (W1)-(W3) are standard [19,21].Condition (W2), in fact, implies that any ∈ U with G ( ) < ∞ satisfies cof ∇ ∈ L q ( 0 , R N ×N ) and, thanks to a well-known inequality [21, Eq. (1.4)], det ∇ ∈ L Nq N −1 ( 0 ).Condition (W4) is not standard, but a similar assumption has been used, e.g., in [22].In the following lemma we show sufficient conditions for the fulfillment of (W4).

Lemma 2.1
The following statements hold: Assume that there exists C > 0 such that for all t 1 , t 2 > 0,

.12)
If h 1 , h 2 are monotone increasing, and m ∈ R, then the function and, analogously, Therefore, so W 0 satisfies the assumptions of part (a).Part (c).We define for i = 1, 2, Thus, g 1 , g 2 , g 3 and m satisfy the assumptions of part (b).

Part (d).
Define where A ∈ R is chosen so that h 3 ≥ 0.Then, ) .

An analogous bound holds for h 2 .
As for h 3 , it is easy to check that there exist We argue by cases so as to show inequality (2.12) for h We have thus shown that h 1 , h 2 , h 3 and m satisfy the assumptions of part (c).
The numerical simulations of Sect. 5 will use the Mooney-Rivlin material in (2.13) when N = 3.In this section we have shown that this W 0 satisfies conditions (W1) and (W4).Condition (W3), on the other hand, is clear, while condition (W2) is easily seen to hold for the exponents p = q = 2.

Choices of Complexity Functionals
The work [13] introduces some examples of complexity functionals in the context of different active materials.In the presence of isotropy, their functionals are based on the right Cauchy-Green deformation tensor C g = G T G associated with growth.However, we have found several advantages to treat G, as opposed to C g as the main variable.Indeed, dealing with C g involves the use of the so-called intrinsic elasticity [23,Sect. 4.2] and needs to incorporate the constraint that C g is a metric tensor, which is difficult to handle.
One of the examples presented in [13] is with C target a given target and α 1 , α 2 > 0 weighting parameters.Although not explicitly mentioned in [13], similar in spirit is the functional with det C target = 1, in which the penalizing term only accounts for the dissimilarity of C g from C target in shape but not in volume.Since in our context, we have decided to work with G as the main variable, the counterparts of C 1 and C 2 are and respectively, for a given target G target , which in C2 satisfies det G target = 1.Here SO(N ) denotes the set of (proper) rotations, and dist the distance between a matrix a set of matrices, i.e., the minimum distance between the matrix and any element of the set of matrices.Note that in C1 we wrote dist(G, SO(N )G target ) instead of G − G target to guarantee frame-indifference and isotropy.Analogously for C2 .A final comment refers to the regularising term with integrand ∇G 2 that is to be used in any complexity functional involved.As a matter of fact, from a practical point of view, it may be advantageous to substitute it by a standard regularising Helmholtz filter Ĝ of the form Ĝ − Here l > 0 acts as a length-scale parameter controlling the amplitude of the regularisation, denotes the Laplacian operator, and N stands for the outer unit normal vector to ∂ 0 .In this case, we replace the term 0 ∇G 2 d X by the L 2 -norm of Ĝ.Note that the operation G → G enjoys much better analytical properties than G → ∇G: the former is a compact operation with nice properties even from the approximation perspective, while the latter is not even continuous.In addition, as just remarked, parameter l can be directly associated with the length-scale of the regularization, a feature that is very convenient form the practical viewpoint.At any rate, this filter has performed quite well in the simulations below.In fact, a different version of the Helmholtz filter more suitable for the implementation will be finally adopted in the numerical simulations.We will explain later how to adapt the proof of existence to these cases.

Mathematical Analysis
This section aims at providing a rigorous mathematical analysis of the shape programming problem (2.10).We shall proceed in two steps.We will first prove that for a given growth tensor G, there exist minimisers of (2.8).Then, existence of solutions for (2.10) is established.
The following lemma is an easy consequence of formula and The following lower semicontinuity result will help in the final steps of the main proof.Recall that for each measurable G : Proof Lemma 3.1 yields convegences (3.1) and (3.2).By a standard fact on the product of two sequences, one factor converging weakly in L 1 and the other one a.e. with and L ∞ bound (see, e.g., [24, Prop.2.61]), we obtain thanks to Lemma 3.1 that To sum up, we have the convergences as well as det G j → det G a.e. with and L ∞ bound, which allows us to apply a standard lower semicontinuity result for polyconvex functions (see, e.g., [25,Th. 5.4] or [24,Cor. 7.9]) and conclude that 0 This proves the result.

Existence of 8 Given G
Before presenting the existence theorems, we recall a property stating that the limit of injective a.e.functions is injective a.e.
The following fundamental existence theorem in nonlinear elasticity will be used throughout.Its proof is the sum of deep and fundamental results in Analysis that are indicated below.
Assume that 0 D is an (N − 1)-rectifiable subset of ∂ of positive (N − 1)dimensional measure and that ¯ : Let the functional 123 be defined in U. Assume that U = ∅ and that I is not identically infinity in U. Then there exists a minimiser of I in U.
Proof The treatment of the linear terms (2.6) and (2.7) is standard (e.g., [19,Sect. 7] or [17,Ch. 5]), so we can assume f = 0 and s 0 = 0. Since U = ∅ and I is not identically infinity in U, there exists a minimising sequence { j } j∈N of I in U. Thus, {I ( j )} j∈N is bounded and, by condition (c), we have that {∇ j } j∈N is bounded in L p and {cof ∇ j } j∈N is bounded in L q .Poincaré's inequality shows that { j } j∈N is bounded in W 1, p .Since p > 1, there exist ∈ L p and a subsequence (not relabelled) such that j in W 1, p .The continuity of traces shows that satisfies the boundary condition.By [21,Lemma 4 for any compact F ⊂ 0 , as j → ∞.By the lower semicontinuity of polyconvex functionals (see, e.g., [25,Th. 5 Since this is true for all compact F ⊂ 0 , by monotone convergence, we obtain Now we show that det ∇ > 0 a.e.Since det ∇ j det ∇ in L 1 and det ∇ j > 0 a.e. for all j ∈ N, we have that det ∇ ≥ 0 a.e.Let A be the set of X ∈ such that det ∇ (X) = 0. We have that det ∇ j → 0 a.e. in A. If |A| > 0, by Fatou's lemma and (d), which is a contradiction.Therefore, |A| = 0 and det ∇ > 0 a.e.By Proposition 3.3, is injective a.e.Therefore, ∈ U, and it is a minimiser of I in U.
In the following result we show how the properties of W 0 are transferred to W G .
for a.e.X ∈ 0 and all Proof We start by proving (a).As G is measurable, there exists a Borel function and so is the function (X, F) → W 0 (F Ḡ(X) −1 ) det Ḡ(X).Therefore, the function Now we show (b).By definition of polyconvexity, there exists a convex function Fix X ∈ 0 .Since e is convex, so is the function ē : as a composition of a linear map with a convex function.Therefore, the function As for (c), by Lemma 3.1 and using elementary properties of the algebra of square matrices, for a.e.X ∈ 0 , Therefore, there exists c > 0 such that Property (d) is immediate.
The existence of minimisers of (2.8) for each given, feasible G is now a straightforward consequence of Theorem 3.4 and Lemma 3.5.
Assume that U = ∅ and that G is not identically infinity in U. Then there exists a minimiser of G in U.
An important issue, which we overlook here, is the potential non-uniqueness of minimiser for given G.A much more delicate analysis would be required to deal with potential bifurcation problems as the tensor G moves in the iterative, approximation procedure implemented in Sect. 5 seeking an optimal G.However, if one sticks to a selected continuous branch of solutions, one would end up with an optimal tensor G.We have to report no difficulties here in the numerical approximations performed.

Existence of Optimal G
The lower semicontinuity of the function ρ H , as given by (2.9), was shown in [8].Although the framework here is somewhat different, the same proof is valid.For the convenience of the reader, we state in a precise way the result inside the proof of [8, Th. 4.2] that will be used in Theorem 3.9 below.
and j is injective a.e., for each j ∈ N. Assume that there exists ∈ W 1, p ( 0 , R N ), with det ∇ > 0 a.e., such that j → a.e., det The following result is an easy consequence of (W4).
for some constant C > 0. Similarly, for some constants c 1 , c 2 > 0. The conclusion readily follows.
Our main result is concerned with the existence of an optimal G. Theorem 3.9 Let W 0 : is the gradient of an a.e.injective map, G is not identically infinity in U, (G) = ( 0 ) , with a minimiser of G in U}.
Assume that A = ∅ and that J is not identically infinity in A. Then there exists a minimiser of J in A.
Proof We will rely on the Direct Method of Calculus of Variations.Let {G j } j∈N be a minimising sequence of J in A. The coercivity of J with respect to ∇G and the fact that Thus, we can extract a subsequence (not relabelled) such that G j G in H 1 and G j → G in L 2 and a.e., for some G ∈ H 1 .Since K is closed, we see that G(X) ∈ K for a.e.X ∈ 0 .Thanks to (b), for a.e.X ∈ 0 , and so, by Fatou's lemma, 123 By the weak convergence in In addition, we ought to check that G ∈ A.
Let us check that G is the gradient of an a.e.injective map.To this aim, we use that for each j ∈ N we have G j = ∇ g j for some Sobolev map g j : 0 → R N that is injective a.e.Without loss of generality, we can assume that 0 g j = 0.
As G j ∈ H 1 , we have that g j ∈ H 2 .By the Poincaré-Wirtinger inequality, Therefore, the sequence { g j } j∈N is bounded in H 2 , so we can extract a subsequence weakly convergent in H 2 to some g ∈ H 2 .Moreover, we can assume that the convergence g j → g also holds a.e.As G j G in H 1 , we have that G = ∇ g .Let us see that g is injective.For this, we can apply Proposition 3.3, according to which it is enough to show that g j → g a.e., det ∇ g j > 0, det ∇ g j det ∇ g in L 1 , with det ∇ g > 0 a.e.Those conditions are satisfied because of the convergence g j g in H 2 and the Sobolev embeddings.Indeed, g j ∈ W 1,N −1 because the embedding so sup j∈N cof ∇ g j L 1 < ∞.Convergence g j → g a.e. was shown earlier.Now, det ∇ g j ≥ m for some m > 0, since G j ∈ K a.e.On the other hand, for N ≤ 3 the compact embedding H 2 ⊂ W 1,r holds for 1 ≤ r < 6, so ∇ g j → ∇ g in L r and, hence, det ∇ g j → det ∇ g in L s for all s < 2. This implies the last condition since det ∇ g ≥ m.
Another main step should focus on the first contribution to the cost given in terms of the Hausdorff distance ρ H (G) , target , as well as the minimising relationship between and G in (2.10) and (2.8).To treat this step, it is mandatory to work with the minimiser j ∈ U of (2.8) corresponding to G j , for each j ∈ N. (3.4) for some C 1 > 0. On the other hand, using (3.3), we find that for some C 2 > 0. This inequality, together with (3.4) and (3.5) shows that {∇ j } j∈N is bounded in L p and {cof ∇ j } j∈N is bounded in L q .As in the proof of Theorem 3.4, we obtain the existence of a ∈ U such that j in W 1, p , together with On the other hand, using dominated convergence, bound (3.6) and Lemma 3.8, we find that lim and the arbitrariness of ˜ in U implies that is a minimiser of (2.8) in U for our limit G.
The final ingredient is provided by Proposition 3.7.Indeed, its assumptions have already been checked, so Altogether, we see that 123 and the proof is finished.

Remark 3.10
The same conclusion of Theorem 3.9 holds if we replace the term ∇G(X) with the term Ĝ(X) ; see (2.16) and note that l > 0 is given.In this case the new functional is We explain the only steps of the proof that differ from that of Theorem 3.9.Let {G j } j∈N be a minimising sequence of J in A. Then {G j } j∈N is bounded in L 2 ( ), so there exists G ∈ L 2 ( 0 ) such that, for a subsequence, G j G in L 2 ( ).Let Ĝ j ∈ H 2 ( 0 ) be the solution of (2.16) with right-hand side G j , and Ĝ ∈ H 2 ( 0 ) the solution with right-hand side G.By standard elliptic regularity theory (see, e.g., [26, Prop.9.26]), Ĝ j Ĝ in H 2 ( 0 ), and, hence, Ĝ j → Ĝ in H 1 ( 0 ).From here, the rest of the proof is identical to that of Theorem 3.9.

The Theory of Morphoelasticity
Our main source in this section is the book [12].The basic principle of morphoelasticity postulates a multiplicative decomposition of the deformation gradient in the form

F(X) = A(X)G(X). (4.1)
This decomposition replaces (2.3), the main difference being that there is no intermediate mappings g and e to account for growth and elastic deformation, respectively: tensors A and G are not associated with any deformation."The growth tensor G takes the initial configuration to a virtual stress-free state that may be incompatible.Then, a local elastic tensor A restores compatibility of the body and enforces the boundary conditions and body forces so that the body is in a compatible configuration in mechanical equilibrium" ([12, p. 355]).Yet, the elastic constitutive law is formulated through an internal energy density W = W ( A) that depends only on the elastic deformation tensor A = FG −1 , i.e., (2.4) and (2.5) are still valid, with F = ∇ .The rest of Sect.2.1 is also valid word by word.The idea of a decomposition of the form (4.1) in Mechanics can be traced back to the mid of the last century and first appeared in the contexts of anelasticity, placticity, dislocations, thermoelasticity and, more recently, biomechanics and growth mechanics.A survey of the history of this decomposition can be found in [27].In fact, the recent papers [28,29] explore when "virtual, incompatible" state actually exists as a global intermediate configuration.
Since our preceding analysis does not rely in any way on the fact that growth tensor G comes from a gradient, i.e., is globally compatible, all of our previous results and discussions are correct in this new setting as well.In particular, the shape programming problem is the same as (2.10): Notice that the only difference with the gradient case is the non-occurrence of the constraint that G must be the gradient of an a.e.injective Sobolev map, so the proof of the existence result in this case would be shorter and less technical than that of Theorem 3.9.For record purposes, we state the main existence theorem in this setting.
, with a minimiser of G in U}.
Assume that A = ∅ and that J is not identically infinity in A. Then there exists a minimiser of J in A.
Remark 4.2 According to (2.16) and Remarks 3.10, the same conclusion of Theorem holds if we replace the term ∇G(X) with the term Ĝ(X) for a given l > 0. In other words, the same result holds for the functional We anticipated in Sect.2.4 that by isotropy one can work with C g = G T G, instead of G as the main variable.In Sects. 2 and 3 we opted for G so as not to deal with the constraint that C g is a metric tensor, but in the context of this section, the theory of morphoelasticity only requires that C g is a field of positive definite symmetric matrices.Let us see why isotropy allows working with C g .Recall from Sect.2.1 that the stored energy function of the material is W 0 and that, once the growth takes place, the total energy of the deformation is given by the integral (2.4), where the new stored-energy function is W G given by (2.5).Now, if a growth tensor G 1 changes to G 2 = RG 1 for some R ∈ SO(N ), then, using (2.5) and the isotropy of W 0 , we find that Thus, by polar decomposition, the dependence of W G on G is only through C g .Likewise, the cost functional C should be isotropic (as well as frame-indifferent).The required conditions for that were written in [13,Sect. 2(a)].In the particular case of functionals of the form 0 for a.e.X ∈ 0 , all symmetric definite positive C ∈ R N ×N and all R ∈ SO(N ).
Since this is the framework of the numerical experiments in the next section, we present the programming problem in the language of (4.2): g and ∈ U a minimiser of (2.8). (4.4) In this case, C C g can take the form (2.14), (2.15) or the general form given in the following theorem, which we present without proof.
, with a minimiser of G in U}.
Assume that A = ∅ and that J is not identically infinity in A. Then there exists a minimiser of J in A.
In the theorem we have taken G as the only symmetric positive definite square root of C g , but, as explained earlier, any square root of C g gives rise to the same problem.Indeed,

Numerical Simulation
This section presents the numerical simulations of the shape-programmig problem analyzed in the previous sections, all done in dimension N = 3.As noticed in [13] and explained in Sect.4, under the presence of isotropy it is more convenient to implement numerically the growth-driven actuation by means of the right Cauchy-Green strain tensor C g = G T G. Indeed, this choice of the control variable introduces less nonlinearity into the formulation of the problem, which facilitates its numerical approximation.Moreover, in this section we use the theory of morphoelasticity, so it is not relevant in practise whether or not the growth tensor is a deformation gradient.All in all, this section addresses the numerical resolution of the shape-programming problem (4.4).
The layout of this section is as follows.In Sect.5.1, after parametrising C g in terms of its eigenvalues and eigenvectors, we find an equivalent formulation of (4.2), which is more amenable to computing the gradients that are required in gradientbased optimisation algorithms.We also present the numerical scheme.In Sect.5.2 we perform several numerical experiments.A set of experiments deals with an initial geometry resembling a beam and another set resembling a shell.For these experiments it is enough to control the eigenvalues of the tensor C g , while keeping the eigenvectors fixed.In the final example we show that, when the initial configuration is a cube and the final configuration is a cylinder, it is necessary to consider both eigenvalues and eigenvectors as design variables to achieve a satisfactory match between the final and the target configurations.
To accomodate the notation to that widely used in Computational Mechanics, from now on in this section, we denote by H = cof F and by J = det F.

Parametrisation of the Growth Tensor
Consider the following version of the Mooney-Rivlin density energy presented in (2.13): This energy is isotropic, so it is valid to work with C g instead of G.The actuated energy density ψ(F, C g ), equivalent to W G in (2.5), adopts the expression (5.1)The eigenvalue decomposition of C g is given by where the ortonormal eigenvectors are encapsulated in the columns of V , i.e., V = v 1 v 2 v 3 , whilst encodes the eigenvalues {λ 1 , λ 2 , λ 3 } of C g .Thus, C g is rewritten as Positive definiteness of C g entails positivity of the eigenvalues λ 1 , λ 2 , λ 3 .Moreover, in some applications one may wish to impose incompressibility in C g , which is modelled by condition det C g = 1 and is equivalent to a restriction only on , namely, det = 1.This can be accomplished, for instance, by parametrising as although we have not included any experiment in this context.
It remains to define the eigenvectors {v 1 , v 2 , v 3 }.A possibility for that is to define the matrix V by using the Rodrigues formula, according to which V is parametrised in terms of a unitary vector k ∈ R 3 , and a rotation angle where E is the third-order alternating tensor (or Levi-Civita tensor), and k is defined through a spherical parametrisation as (5.5)

An Equivalent Formulation of the Optimisation Problem
This new set of design variables, λ = {λ 1 , λ 2 , λ 3 } and θ = {θ 1 , θ 2 , θ 3 }, allows us to consider the new optimisation problem where, in analogy to (2.16), the fields are regularised versions of λ and θ , respectively.More precisely, for 1 ≤ i ≤ 3, the functions λi (X) and θi (X) are in H 1 ( 0 ) and are the solutions of the boundary value problems where l is a length-scale parameter and N is the normal to 0 .Note that in (5.6) we have not included the compexity functional C described in Sect.2.4, although it can easily be incorporated.In any case, problem (5.6) fits in the theory of Theorem 4.3, just by putting φ = 0, while the regularisation term based on the L 2 norm of ∇C g is replaced by the regularisations λ and θ.
In addition, to account for the complexity of the actuation, we impose lower and upper pointwise bounds on λ i , namely λ lb i ≤ λ i ≤ λ ub i , 1 ≤ i ≤ 3, rather than an L 2norm constraint.Indeed, those pointwise contraints are more effectively handled by constrained optimisation methods as they prevent from tuning the weighting parameters that appear in the complexity functional.Similarly, the angles θ i are confined in suitable predefined intervals.Due to the maximum principle, the regularised variables λi and θi satisfy the same bounds.
We present the existence result for (5.6).

Computation of Continuous Gradients
As is customary in gradient-based optimisation, in order to compute a descent direction, we use the standard Lagrangian method [31].To this end, let us consider the Lagrangian L defined as L ( , p, λ, θ ) = J ( ) − 0 P ∇ , λ, θ : ∇ p d X, (5.8) which is defined for ( , p, λ, θ ) and satisfying the boundary condition in 0 D .For the strain energy in (5.1), the first Piola-Kirchhoff stress tensor P = ∂ F ψ is (5.9) where ( A × B) i I = E i jk E I J K A j J B k K , for A, B ∈ R 3×3 , and E I J K represents the components of the third-order alternating tensor.Notice that ( , p, λ, θ ) are considered as independent variables in (5.8).The stationary condition of L (5. (5.12)where C represents the fourth order elasticity tensor, defined as which takes the form and for A ∈ R 3×3 and A ∈ R 3×3×3×3 .From the linear equation in (5.12) it is therefore possible to obtain the adjoint state p.
The directional derivative of the Lagrangian L with respect to the design variables {λ 1 , λ 2 , λ 3 } and {θ 1 , θ 2 , θ 3 } yields which, for the specific expression for P in (5.9), takes the following form where the second order tensor T is defined as (5.13)

Numerical Scheme
In order to clarify how the different equations featured in the current section have been embedded into a gradient algorithm, we summarise the steps involved in the optimisation.
Starting from an initial guess (λ 0 , θ 0 ), we proceed with the loop: (i) Solve the state equation (5.10) with expression (5.11), which yields the new deformation mapping and the new deformed configuration .(ii) Based on the new deformation mapping , compute the adjoint state field p by means of (5.12).(iii) Compute the objective function J ( (λ, θ )).(iv) Compute descent directions for each of the design variables, namely ∂ λ i L and to the gradient algorithm in order to determine the step size and hence, the new value of the design variables (λ 1 , λ 2 , λ 3 ) and (θ 1 , θ 2 , θ 3 ).

Remark 5.2
Although the lower bound conditions λ i > 0 ensuring the positive definiteness of the tensor C g have not been explicitly included in the Lagrangian L in (5.8), any standard gradient-based algorithm such as the interior-point, can easily handle this type of constraint, by augmenting the Lagrangian L in (5.8) by means of the method of Lagrange multipliers.Evaluation of this additional term and of its derivatives with respect to the design variables (λ 1 , λ 2 , λ 3 ) is therefore omitted in the derivations included in this section.In fact, as mentioned in Sect.5.1.2,pointwise bounds on λ i and θ i , 1 ≤ i ≤ 3, have been included in the optimisation method.

Numerical Experiments
The objective of this section is to demonstrate the applicability of the proposed formulation in the context of shape morphing, i.e., determining the value of the optimal design variables {λ 1 , λ 2 , λ 3 } and {θ 1 , θ 2 , θ 3 } following a gradient-based approach with the aim of attaining the closest growth-driven configuration to a given target configuration.As indicated in the introduction, one of the objectives of this paper is to present an alternative formulation to other analytical approaches that make use of simplifying assumptions such as the absence of boundary conditions, which permit to obtain a closed-form solution of the optimal growth tensor [32,33].We do not intend to claim that our formulation is more advantageous than others.On the contrary, as in other areas of continuum mechanics, analytical solutions can be extremely useful.Our purpose is to illustrate the possibility to apply inverse techniques for the optimal solution of this problem, which can address more generic situations than those covered by analytical approaches.
In the first two examples (Sects.5.2.1 and 5.2.2), we will advocate for a widely accepted formulation in engineering, according to which the eigenvectors {v 1 , v 2 , v 3 } of C g (see equation (5.2)) remain fixed, while only the eigenvalues {λ 1 , λ 2 , λ 3 } serve as the unknown fields to be determined analytically [32,33].Although this approach is less flexible compared to the more comprehensive formulation discussed in Sect.5.1, it has exhibited reasonably positive outcomes in terms of achieving the target configuration.Typically, the eigenvectors of C g are considered coincident with the tangent vectors associated with the curvilinear coordinate system that describes the geometry of the initial solid configuration.
However, in the final example of Sect.5.2.3, we will illustrate a scenario where incorporating additionally the eigenvectors as design variables (specifically, the three angular fields {θ 1 , θ 2 , θ 3 } in equation (5.5)) allows for a higher degree of flexibility.This enhanced formulation enables a significantly better approximation to the target configuration.
In all the examples, the upper and lower bounds used for the eigenvalues λ i (i = {1, 2, 3}) (see (5.6)) are λ lb i = 0.01; λ ub i = 12.With respect to the upper and lower bounds used for {θ 1 , θ 2 , θ 3 } (see (5.6)), these are We have chosen these bounds since in the performed experiments, it is not expected that rotations of more than one loop take place, but of course the above bounds can be expanded if the geometry of the problem suggests so.

Beam-Like Applications
The first examples consider applications where the geometry of the undeformed domain 0 resembles that of a beam.In particular, we consider the rectangular section beam in Fig. 2a and the beam with circular cross-section in Fig. 2b.For both cases, the eigenvectors {v 1 , v 2 , v 3 } featuring in the definition of C g in (5.2) are defined as for the case in Fig. 2a and for the case in Fig. 2b.In both cases, the boundary conditions are such that the displacements in X 1 = 0 are 0 in the three directions {E 1 , E 2 , E 3 } of the configuration {X 1 , X 2 , X 3 }.Three target configurations, target = d ( 0 ), have been prescribed: (i) Shape morphing configuration 1: rectangular cross-section beam with target configuration given by d (ii) Shape morphing configuration 2: rectangular cross-section beam with target configuration given by (iii) Shape morphing configuration 3: circular cross-section beam with target configuration given by For the case of the rectangular cross-section beam in Fig. 2a, the final configurations attained at convergence are depicted in Fig. 3, corresponding with the optimal solutions that yield the closest growth-driven configurations to the target configurations denoted as shape morphing configurations 1 and 2. In addition, Fig. 4 depicts the evolution of the cost function for the case of the shape morphing configuration 1.The interior-point algorithm has been used as the optimization method.
With regard to the circular cross-section beam in Fig. 2b, with target configuration given in Eq. (5.16), the final growth-driven configuration is displayed in Fig. 5, along with the contour plot distribution of the three design variables { λ1 , λ2 , λ3 }.The tight agreement with respect to the target configuration initially prescribed in Eq. (5.16) is shown in Fig. 5d.

Shell-Type Applications
Next, we consider the two undeformed configurations given in Fig. 6a and the beam with circular cross-section in Fig. 6b.For both cases, the eigenvectors {v 1 , v 2 , v 3 } are defined as Fig. 3 Representation of the computed optimal deformation (indistinguishable from the target configuration) and the contour plot distribution of λ1 −1 (x) , for the beam with rectangular section in Fig. 2a for target configurations: a equation ( 5 In both cases, the boundary conditions are such that the displacements vanish in X 3 and r = R 1 (for Fig. 6a) and r = R (for Fig. 6b) in the three directions {E 1 , E 2 , E 3 } of the configuration {X 1 , X 2 , X 3 }.Two target configurations, target = d ( 0 ), have been prescribed: (i) Shape morphing configuration 4: initial geometry given in Fig. 6a with target configuration given by For the case of the initial geometry in Fig. 6a, the final configuration attained can be observed in Fig. 7, corresponding with the optimal solutions that yield the closest growth driven configurations to the target configuration denoted as shape morphing configuration 4. It is worth empashising how the optimal solution is capable of, starting with a flat disk geometry, inducing a deformation of the continuum into the final conical shape illustrated in this figure.From Fig. 7d, the almost perfect match between the growth-driven and target configurations can be observed.
Finally, for the case of the initial geometry in Fig. 6b, the final configuration attained can be observed in Fig. 8.In this case, attaining the target configuration entails a considerable enlargement of the initial geometry along the X 3 direction, in addition to a bending in the X 1 and X 2 directions, yielding the conical shape illustrated in this figure.Figure 8d shows the almost perfect match between the growth-driven and target configurations can be observed.

Cube to Cylinder Geometry
The objective of this example is to evidence what was anticipated in the introductory part of Sect.5.2.Specifically, although the examples shown in Sects.5.2.1 and 5.2.2 have demonstrated that including only eigenvalues as design variables whilst maintaining the eigenvectors fixed throughout the optimisation process can yield extremely good results, this might not be the case for any predefined target configuration.In order to illustrate that, we consider now the initial and target configurations shown in Fig. 9.The problem requires a transformation from a cube into a cylinder, although the pictures are represented in 2D.
We solved the problem using two formulations: Figure 10 includes the results yielded by both formulations.As expected, the deformed configuration resulting from the second formulation (including both {λ 1 , λ 2 , λ 3 } and {θ 1 , θ 2 , θ 3 } as design variables) yields a significantly better approximation to the unattainable circular target configuration.This is corroborated by the values of the objective function ρ H attained by both formulations in the last optimisation iteration, when numerical convergence was observed.Specifically, the ratio between the values yielded by both formulations were where the value in the numerator refers to the first formulation (only λ 1 , λ 2 , λ 3 as design variables).Finally, Fig. 11 illustrates the optimal solution obtained for the eigenvectors v 1 (θ 1 , θ 2 , θ 3 ) and v 2 (θ 1 , θ 2 , θ 3 ).This figure demonstrates the necessity to modify spatially these eigenvectors in order to yield the observed higher flexibility.
any singular value σ of F} for some 0 < m 1 ≤ m 2 .In addition, one may want to model the growth given by G as incompressible; in this case, relevant sets K are K = {F ∈ R N ×N : F ≤ M and det F = 1} and {F ∈ R N ×N : det F = 1 and m 1 ≤ σ (F) ≤ m 2 for any singular value σ of F}.

Theorem 3 . 4
Assume that W : 0 × R N ×N + → [0, ∞] satisfies the following conditions: (a) W is L × B-measurable, where L denotes the Lebesgue σ -algebra in 0 , and B stands for the Borel σ

Fig. 4
Fig.4 Evolution of the objective function with the number of iterations for the target configuration given in Eq.(5.14)

Fig. 5 8 ( 7 − 2 (X 3 + 1 )
Fig. 5 Representation of the computed optimal deformation and contour plot distribution of a λ1 −1 (x) , b λ2 −1 (x) and c λ3 −1 (x) for the example with initial configuration depicted in Fig. 2. The translucid geometry represents the initial configuration.d Agreement between the target configuration (grey colour) and the deformed solid subjected to the optimal growth tensor (red) (Color figure online)

Fig. 7
Fig. 7 Representation of the computed optimal deformation and contour plot distribution of a λ1 , b λ2 and c λ3 for the example with initial configuration depicted in Fig. 6a.The translucid geometry represents the initial configuration.d Agreement between the target configuration (grey meshed domain) and the deformed solid subjected to the optimal growth tensor (red) (Color figure online)

Fig. 8 Fig. 9
Fig. 8 Representation of the computed optimal deformation and contour plot distribution of a λ1 , b λ2 and c λ3 for the example with initial configuration depicted in Fig. 6b.The grey domain represents the initial configuration.d Agreement between the target configuration (grey meshed domain) and the deformed solid subjected to the optimal growth tensor (red) (Color figure online) Fig. 9 Undeformed configuration (red), representing a square of side 4. The circular geometry represents the target configuration, with diameter 12 (Color figure online)