Shape deformation from metric ’s transport

International Journal of Non-Linear Mechanics 119 (2020) 103326

Contents lists available at ScienceDirect

International Journal of Non-Linear Mechanics journal homepage: www.elsevier.com/locate/nlm

Shape deformation from metric ’s transport L. Teresi a ,∗, F. Milicchio a , S. Gabriele a , P. Piras b , V. Varano a a b

Università degli Studi Roma Tre, Italy Università degli Studi di Roma ‘‘La Sapienza", Italy

ARTICLE Keywords: Elastic metric Distortions Non linear elasticity Metric tensor

INFO

ABSTRACT We exploit the possibility of deforming a solid volume, or a surface as well, by assigning a target metric. As well known, an elastic solid may change its shape under two different kinds of actions: one are the loadings, the other one are the distortions, aka, the pre-strains. Actually, a distortion prescribes the target metric sought by the solid, and the metric effectively realized is the one that minimize the distance, measured through an elastic energy, between the target metric and the actual one. The challenge of our work is to use the metric information taken from a pair of configurations 𝑋, 𝑋 of a given body 𝑥 , to deform a different body 𝑦 , so that its deformation be as similar as possible to the one suffered by 𝑥 . In particular, we assess the metric change between the two configurations 𝑋, 𝑋 of the template 𝑥 ; then, we use this information to build a target metric for 𝑦 , exploiting the notion of distortions: given a source configuration 𝑌 , we find a target configuration 𝑌 by minimizing an appropriate elastic energy. The proposed method is very effective in deforming solids as well as surfaces, and in reproducing similar deformations even when applied to different bodies.

1. Introduction The challenge of our work is the following: gauge the metric change ̄ 𝑋 of a given material body B𝑥 , between a pair of configurations 𝑋, and use that information to morph a different body B𝑦 , so that its deformation be as similar as possible to the one suffered by B𝑥 . The problem of the shape reconstruction from intrinsic geometrical properties, such as the metric or the curvature, is common to different fields, among them Shape Analysis [1,2] and Solid Mechanics [3,4]. Among the main applications of metric’s transport in Computer Graphics, we cite Motion Capture, Style/Pose Transfer, and Frame interpolation. The first allows to capture facial expressions or body-poses from human actors and then use them to animate digital characters. The second allows to reuse many times a deformation hand-sculpted by a digital artist to deform other shapes [5]. Finally, building a suitable interpolation between shapes can be used to create an animation starting from two frames (the so-called key-frames). In Statistical Shape Analysis (for example, in Medical Image Analysis) the definition of a trajectory interpolating different shapes allows introducing a notion of ‘‘distance’’ between shapes, making it possible to quantitatively measure the differences between a group of shapes [6]. Here, we are interested in linking these two fields: our goal is to develop a reliable mathematical tool to probe the deformation between a pair of shapes of a same object, and then use this information to deform likewise a third different object.

In the last years, the concepts of target metric, non-euclidean plates [7] and non-linear distortions [8], have been widely used to describe shape formation processes in both natural (e.g. growth and remodeling of biological structures), and artificial contexts (e.g. design of actuators). In particular the local strain, in the deformation of a solid body (or a shell), can be characterized by the change of metric and curvature [9]. Our approach rests on two pillars: one is morphing from a target metric, the other one is the transport of metric between body manifolds. We morph a material body basing on the theory of finite elasticity with large distortions: a distortion induces a target metric, and the configuration which is effectively realized is the one that minimizes the distance, measured through the elastic energy, between the target metric and the actual one. We note that this technique may be used to the full hierarchy of geometric manifolds, that is, for 1D curves, 2D surfaces, and 3D solids, provided the appropriate metric information and elastic energy is used. As example, for the curve, metric information are given by arc elongation, curvature and torsion, while the energy is the so-called beam energy; for the surfaces, we have the first and second fundamental forms, and the shell energy; for the solids, the standard metric tensor and the hyper-elastic energy. The metric transport is necessary to use the metric information gauged on a pair of configurations to morph a different third body. In particular, to morph 𝑌̄ , we compute the metric tensor associated to

∗ Correspondence to: Dipartimento di Matematica e Fisica, Università Roma Tre, via della Vasca Navale 84, I-00146 Roma, Italy. E-mail addresses: [email protected] (L. Teresi), [email protected] (F. Milicchio), [email protected] (S. Gabriele), [email protected] (P. Piras), [email protected] (V. Varano).

https://doi.org/10.1016/j.ijnonlinmec.2019.103326 Received 29 March 2019; Accepted 15 October 2019 Available online 24 October 2019 0020-7462/© 2019 Elsevier Ltd. All rights reserved.

L. Teresi, F. Milicchio, S. Gabriele et al.

International Journal of Non-Linear Mechanics 119 (2020) 103326

̄ 𝑋, and then transport this tensor onto 𝑌̄ to the two configurations 𝑋, define a distortion field. The task of transporting deformations is well known [6,10–13]. In [14] a morphing based on the nonlinear elasticity with distortions has been introduced, with applications in Medical Imaging. A drawback of the cited method is the trivial deformation transport which is used. Here, we propose a better formalization and a generalization of the method in [14], by introducing a common reference configuration 𝑍 ̄ 𝑋, 𝑌̄ and 𝑌 . This common for the four different shapes under exam 𝑋, reference allows the use of a unique atlas for all the manifolds, and yields a point-wise correspondence between their points, a feature which is essential to extend the method from solid bodies to shell bodies. Fig. 1. The essence of a distortion. The volume element 𝑑𝑣𝑦̄ , attached to a point 𝑦̄ ∈ 𝑌̄ of the source configuration (top, right) is mapped by 𝐅𝑜 onto a distorted element 𝑑𝑣𝑜 , attached to the same point 𝑦̄ (top, left). To get a compatible volume element 𝑑𝑣𝑦 (bottom, left), we must add a further strain to 𝐅𝑜 , called the elastic strain 𝐅𝑒 . The embedding 𝐅 is thus given by 𝐅 = 𝐅𝑒 𝐅𝑜 .

2. The elastic metric We briefly introduce the notion of distortion and the associated elastic metric; the theoretical framework of finite elasticity with distortions is detailed in [15], while a thoroughly discussion on the elastic metric can be found in [8]. Let E be the 3D Euclidean ambient space; we denote with VE the corresponding translation space, and with Lin = VE ⊗ VE = Sym ⊕ Skw the space of double tensors on VE (linear maps of VE into itself), decomposed in the sum of the symmetric tensors Sym, and the skewsymmetric ones Skw. Finally, Rot ⊂ Lin is the space of the orthogonal transformations of VE . Given the configuration 𝑌̄ of a body manifold B𝑦 , we define a new configuration 𝑌 of the body via the deformation 𝑓 : 𝑓 ∶ 𝑌̄ → E

embedded in the 3D Euclidean space, and any actual configuration will have 𝐂𝑒 ≠ 𝐈. A necessary condition for the existence of a deformation 𝑓 such that 𝐅𝑇 𝐅 = 𝐂𝑜 , that is, which realizes a given target metric 𝐂𝑜 , is that the Riemannian curvature associated to 𝐂𝑜 be null; is such a case, the target metric, as well the associated distortion, are called compatible. 2.1. The minimization problem

(1)

Let 𝑑𝑣𝑦̄ and 𝑑𝑣𝑦 be the volume elements of 𝑌̄ and 𝑌 , respectively, and let 𝐽𝑦̄ = det 𝐅𝑦̄ , 𝐽𝑦 = det 𝐅𝑦 , 𝐽𝑒 = det 𝐅𝑒 = 𝐽 ∕𝐽𝑜 be the Jacobians of the tangent maps; moreover, let 𝑑𝑣𝑜 be the distorted volume element of 𝑌̄ . Volume elements are related by the following relations:

(2)

𝑦̄ ↦ 𝑦 = 𝑓 (𝑦) ̄ ,

The deformation 𝑓 is a smooth embedding of 𝑌̄ onto E , and its gradient is denoted by 𝐅 = ∇𝑓 ; the gradient 𝐅 is the tangent map from the tangent bundle 𝑇 𝑌̄ of 𝑌̄ , to the tangent bundle 𝑇 𝑌 of 𝑌 . A body element attached at 𝑦̄ ∈ 𝑌̄ is the tangent space 𝑇 𝑌̄𝑦̄ to 𝑌̄ at 𝑦. ̄ Distortions are described by a smooth tensor-valued field 𝐅𝑜 ∶ 𝑌̄ → Lin(𝑇 𝑌̄ , 𝑇 𝑌̄ ) ,

𝑑𝑣𝑜 = 𝐽𝑜 𝑑𝑣𝑦̄ ,

(3)

𝐅𝑒 = 𝐅 𝐅−1 𝑜 ,

Thus, the elastic-energy density 𝜑, per unit of initial volume 𝑑𝑣𝑦̄ , is given by multiplying 𝜑𝑜 for the Jacobian 𝐽𝑜 :

target metric, due to 𝐅𝑜 ,

𝐂 = 𝐅𝑇 𝐅 ,

actual metric, given by 𝐅 ,

−1 𝐂𝑒 = 𝐅⊤𝑒 𝐅𝑒 = 𝐅−⊤ 𝑜 𝐂 𝐅𝑜

elastic metric .

where 𝜑 = 𝜑(𝐅𝑜 , 𝐅) has been considered as a function of the two terms which defines 𝐂𝑒 . Given an initial configuration 𝑌̄ and a distortion field 𝐅𝑜 , the actual configuration 𝑌 can be found by minimizing with respect to 𝑓 the elastic energy (9). 𝑌 = 𝑓 (𝑌̄ ) = min 𝜑(𝐅𝑜 , ∇𝑓 ) .

(10)

𝑓

3. Morphing from metric ̄ 𝑋 of a 3D body We now consider two different configurations 𝑋, B𝑥 ⊂ E , and one configuration 𝑌̄ of a different body B𝑦 . Our challenge is to map the source 𝑌̄ onto a target 𝑌 such that the deformation between 𝑌̄ and 𝑌 be as close as possible to the one between 𝑋̄ and 𝑋. We tackle this problem by measuring the metric change from 𝑋̄ to 𝑋, and then by using this info to deform 𝑌̄ onto 𝑌 . A key issue is to establish a point-wise correspondence between the points of all ̄ 𝑋, 𝑌̄ , and 𝑌 ; thus, we introduce a common our configurations 𝑋, reference domain 𝑍, denote with 𝑧 ∈ 𝑍 its points, and consider our configurations as images of 𝑍 under appropriate deformation maps:

(5)

The name target metric, which is also known as natural or intrinsic metric [18], stems form the fact that, if the actual metric 𝐂 is equal to 𝐂𝑜 , then the elastic metric is trivial: 𝐂 = 𝐂𝑜 .

(9)

𝜑 = 𝜑(𝐅𝑜 , 𝐅) = 𝐽𝑜 𝜑𝑜 (𝐂𝑒 ) ,

(4)

𝐂𝑜 = 𝐅𝑇𝑜 𝐅𝑜 ,

(8)

𝜑𝑜 = 𝜑𝑜 (𝐂𝑒 ) .

Thus, 𝐅𝑒 is the difference between the distortion 𝐅𝑜 and the tangent map 𝐅, in the sense of the multiplicative composition (see Fig. 1). Basing on the three tangent maps 𝐅, 𝐅𝑜 , and 𝐅𝑒 , we may define the following metric tensors, having the role of right Cauchy–Green strain-measures:

⇔

(7)

Any elastic-energy density 𝜑𝑜 per unit distorted volume 𝑑𝑣𝑜 has to be a function of the elastic metric 𝐂𝑒 :

with Lin(𝑇 𝑌̄ , 𝑇 𝑌̄ ) an endomorphism of 𝑇 𝑌̄ , having positive Jacobian determinant 𝐽𝑜 ∶= det 𝐅𝑜 > 0. The essence of a distortion can be briefly summarized as follows: being 𝐅 the gradient of a deformation, it embeds body elements onto E ; conversely, 𝐅𝑜 is not required to be the gradient of any deformation from 𝑌̄ , and in general is not compatible, that is, it does not exist an embedding for the distorted body elements. Thus, an extra strain must be added to 𝐅𝑜 , called elastic strain, whose introduction dates back to [16,17]

𝐂𝑒 = 𝐈

𝑑𝑣𝑦 = 𝐽 𝑑𝑣𝑦̄ = 𝐽𝑒 𝐽𝑜 𝑑𝑣𝑦̄ .

(6)

In general, 𝐂𝑜 is not Euclidean, that is, flat Riemannian; as a consequence, the body manifold 𝑌̄ , endowed with such a metric, cannot be

𝑥̄ = 𝑓𝑥̄ (𝑧) , 2

𝑥 = 𝑓𝑥 (𝑧) ,

𝑦̄ = 𝑓𝑦̄ (𝑧) ,

𝑦 = 𝑓𝑦 (𝑧) .

(11)

L. Teresi, F. Milicchio, S. Gabriele et al.

International Journal of Non-Linear Mechanics 119 (2020) 103326

Fig. 2. Diagram of the deformation maps from the reference domain 𝑍 to the four configurations under exam. Given 𝑓𝑥 , 𝑓𝑥̄ and 𝑓𝑦̄ , we look for a map 𝑓 such that the metric change between 𝑌̄ and 𝑌 be as close as possible to the one between the pair ̄ 𝑋. 𝑋,

The first three maps are given, or can be inferred by the given configurations, while the last one 𝑓𝑦 is the unknown of the problem; by composition, we may define the two maps ̄ 𝑓𝑥𝑥 ̄ ∶𝑋 →𝑋, 𝑓 ∶ 𝑌̄ → 𝑌 ,

Fig. 3. Diagram of the tangent maps; the endomorphisms of a same bundle are denoted with an arc (bar 𝐅𝑜 ). It is worth noting that 𝐅𝑜 is an endomorphism of 𝑇 𝑌̄ : it distorts the volume elements of 𝑌̄ . In general, 𝐅𝑜 is not compatible: 𝐅𝑒𝑦̄ is the extra strain to ̄ and of 𝑌 add in order to deform 𝑌̄ into 𝑌 . The relative stretch of 𝑋 with respect to 𝑋, with respect to 𝑌̄ , can be measured by 𝐔𝑥𝑥 ̄ and 𝐔𝑥𝑥 ̄ , respectively, both endomorphisms of 𝑇 𝑍. The target metric 𝐂𝑜 is an endomorphism of 𝑇 𝑌̄ : it can be transport to 𝑇 𝑍, yielding 𝐂𝑜𝑦 .

such that 𝑥 = 𝑓𝑥 (𝑧) = 𝑓𝑥𝑥 ̄ ◦𝑓𝑥̄ (𝑧) , such that 𝑦 = 𝑓𝑦 (𝑧) = 𝑓◦𝑓𝑦̄ (𝑧) .

(12)

Note that for the map 𝑓 , the unknown of our problem, we use a simplified notation without subscripts (see Fig. 2). Given a deformation 𝑓𝑘 , with 𝑘 = 𝑥, ̄ 𝑥, 𝑦, ̄ 𝑦, we denote with 𝐅𝑘 = ∇𝑓𝑘 its gradient, that is, the linear map between the corresponding tangent spaces: 𝐅𝑥 ∶ 𝑇 𝑍 → 𝑇 𝑋 ;

𝐅𝑥̄ ∶ 𝑇 𝑍 → 𝑇 𝑋̄ ,

(13)

𝐅𝑦 ∶ 𝑇 𝑍 → 𝑇 𝑌 ,

𝐅𝑦̄ ∶ 𝑇 𝑍 → 𝑇 𝑌̄ ,

(14)

Moreover, both 𝐂𝑥 and 𝐂𝑥𝑥 ̄ correspond to the standard right Cauchy– Green tensor used in finite elasticity; the latter measures the metric of the embedding 𝑓𝑥 with respect to the metric of 𝑓𝑥̄ . Let us consider the polar decomposition of both 𝐅𝑥 and 𝐅𝑥̄ : 𝐅𝑥 = 𝐑𝑥 𝐔𝑥 ,

We use a same notation to denote the volume elements 𝑑𝑣𝑘 of the four configurations and the Jacobian 𝐽𝑘 = det 𝐅𝑘 of the tangent map 𝐅𝑘 ; moreover, we denote with 𝑑𝑉 the reference volume element on 𝑍. Volume elements are related by the following relations: 𝑑𝑣𝑜 = 𝐽𝑜 𝑑𝑣𝑦̄ = 𝐽𝑜 𝐽𝑦̄ 𝑑𝑉 ,

𝑑𝑣𝑥̄ = 𝐽𝑥̄ 𝑑𝑉 ,

𝑑𝑣𝑥 = 𝐽𝑥 𝑑𝑉 .

𝐅𝑥̄ = 𝐑𝑥̄ 𝐔𝑥̄ ,

(18)

with 𝐔𝑘 ∈ Sym and positive definite, and 𝐑𝑘 ∈ Rot. The right stretch tensor 𝐔𝑥 and 𝐔𝑥̄ filter out the rigid part of the tangent map, while retaining all the metric information of a deformation. From (18), it follows

(15)

𝐂𝑥 = 𝐔2𝑥 ,

Our approach is based on: (i) morphing from a target metric; (ii) transporting a metric between body manifolds. The algorithm proceeds as follows:

−1 −1 𝑇 𝐂𝑥𝑥 ̄ = 𝐑𝑥̄ 𝐔𝑥̄ 𝐂𝑥 𝐔𝑥̄ 𝐑𝑥̄ ;

(19)

The representation (19)2 deserves one more comment; let us define a right stretch tensor 𝐔𝑥𝑥 ̄ as follows 2 −1 −1 −1 𝑇 −1 (𝐔𝑥𝑥 ̄ ) ∶= 𝐔𝑥̄ 𝐂𝑥 𝐔𝑥̄ = (𝐔𝑥 𝐔𝑥̄ ) (𝐔𝑥 𝐔𝑥̄ ) .

̄ 𝑋, that is, • Compute the metric tensor associated to the pair 𝑋, the metric 𝐂𝑥𝑥 of the deformation 𝑓 with respect to 𝑓𝑥̄ . This is ̄ 𝑥 the metric associated to the map 𝑓𝑥𝑥 . ̄ ̄ • Transport the metric info contained in 𝐂𝑥𝑥 ̄ onto 𝑌 so to define a ̄ distortion tensor 𝐅𝑜 for 𝑌 . • Define the elastic metric 𝐂𝑒𝑦̄ associated to the deformation 𝑓𝑦𝑦 ̄ , basing on the distortion 𝐅𝑜 ; • Morph 𝑌̄ by minimizing an elastic energy density 𝜑 = 𝜑(𝐂𝑒 ) with respect to the unknown function 𝑓 .

(20)

It turns out that 𝐔𝑥𝑥 ̄ is an endomorphism of 𝑇 𝑍, and measures the ̄ We use now a noteworthy relative stretch of 𝑋 with respect to 𝑋. property involving√the square root a rotated tensor: given 𝐀 ∈ Lin and √ 𝐑 ∈ Rot, it holds: 𝐑 𝐀𝐑𝑇 = 𝐑 𝐀𝐑𝑇 . From (19), (20), it follows √ 𝑇 𝐂𝑥𝑥 (21) ̄ = 𝐑𝑥̄ 𝐔𝑥𝑥 ̄ 𝐑𝑥̄ ; √ thus, the right stretch 𝐂𝑥𝑥 ̄ associated to the metric tensor 𝐂𝑥𝑥 ̄ may be interpreted as the relative stretch 𝐔𝑥𝑥 ̄ , rotated by 𝐑𝑥̄ .

All the relevant tangent maps are shown in Fig. 3.

3.2. The distortion 𝐅𝑜 for 𝑌̄ , based on the metric 𝐂𝑥𝑥 ̄

̄ 3.1. The metric 𝐂𝑥𝑥 ̄ of the pair 𝑋, 𝑋

Our goal is to transport metric information from 𝑇 𝑋̄ to 𝑇 𝑌̄ , and use this information to define a distortion 𝐅𝑜 for 𝑌̄ . As before, to filter out the rigid part of the tangent maps, we consider the polar decomposition of 𝐅𝑦̄ as well: 𝐅𝑦̄ = 𝐑𝑦̄ 𝐔𝑦̄ . By using √ ̄ ̄ only the rotational parts 𝐑𝑥̄ , 𝐑𝑦̄ , we can transport 𝐂𝑥𝑥 ̄ ∈ 𝑇 𝑋 to 𝑇 𝑌 in two steps: √ √ 𝑇 𝑇 𝑋̄ ∋ 𝐂𝑥𝑥 𝐂𝑥𝑥 ̄ ↦ 𝐑𝑥̄ ̄ 𝐑𝑥̄ = 𝐔𝑥𝑥 ̄ ∈ 𝑇𝑍 ;

We start by assuming that the relevant information about the deformation 𝑓𝑥𝑥 ̄ are given by the so-called right √ stretch, defined as the square root of the metric tensor, that is, by √ 𝐂𝑥𝑥 ̄ . Thus, we first define ̄ ̄ 𝐂𝑥𝑥 𝐂𝑥𝑥 ̄ , and then we transport its square root ̄ from 𝑇 𝑋 to 𝑇 𝑌 . ̄ Given the tangent map 𝐅𝑥𝑥 ∶ 𝑇 𝑋 → 𝑇 𝑋, defined by 𝐅 = ̄ 𝑥𝑥 ̄ 𝐅𝑥 𝐅𝑥−1 , the metric tensor 𝐂 associated to the deformation 𝑓 has the 𝑥𝑥 ̄ 𝑥𝑥 ̄ ̄ following representation 𝑇 −𝑇 −1 𝐂𝑥𝑥 ̄ = 𝐅𝑥𝑥 ̄ = 𝐅𝑥̄ 𝐂𝑥 𝐅𝑥̄ , ̄ 𝐅𝑥𝑥

with

𝐂𝑥 = 𝐅𝑇𝑥 𝐅𝑥 .

We note that 𝐂𝑥 is a field defined on 𝑍, while 𝐂𝑥𝑥 ̄ 𝐂𝑥 ∶ 𝑇 𝑍 → 𝑇 𝑍 ;

̄ ̄ 𝐂𝑥𝑥 ̄ ∶ 𝑇𝑋 → 𝑇𝑋 .

𝑇 ̄ 𝑇 𝑍 ∋ 𝐔𝑥𝑥 (22) ̄ ↦ 𝐑𝑦̄ 𝐔𝑥𝑥 ̄ 𝐑𝑦̄ ∈ 𝑇 𝑌 ; √ We note that 𝐂𝑥𝑥 ̄ has been transported by the composite rotation 𝐑𝑥̄ 𝑦̄ = 𝐑𝑦̄ 𝐑𝑥𝑇̄ , representing the relative rotation of 𝑇 𝑌̄ with respect to ̄ From (22), we define the distortion field 𝐅𝑜 as follows: 𝑇 𝑋.

(16) ̄ is defined on 𝑋;

𝑇 𝐅𝑜 = 𝐑𝑦̄ 𝐔𝑥𝑥 ̄ 𝐑𝑦̄ .

(17) 3

(23)

L. Teresi, F. Milicchio, S. Gabriele et al.

International Journal of Non-Linear Mechanics 119 (2020) 103326

Thus, given the initial configuration 𝑌̄ , and the pair B𝑥 , B𝑦 , the sought configuration 𝑌 = 𝑓𝑦 (𝑍) can be found by minimizing with respect to 𝑓𝑦 the elastic energy (36).

3.3. The elastic metric 𝐂𝑒𝑦̄ based on the distortion 𝐅𝑜 We follow the framework of Section 2; given the elastic deformation 𝐅𝑒𝑦̄ = 𝐅 𝐅−1 𝑜 , and using (23), we may represent 𝐅𝑒𝑦̄ as −1 𝑇 𝐅𝑒𝑦̄ = 𝐅 𝐅−1 𝑜 = 𝐅 𝐑𝑦̄ 𝐔𝑥𝑥 ̄ 𝐑𝑦̄ .

𝑌 = 𝑓𝑦 (𝑍) = min 𝜑𝑒 (𝐔𝑦 , 𝐔𝑦̄ , 𝐔𝑥 , 𝐔𝑥̄ ) . 𝑓𝑦

(24)

The deformation of 𝑌 with respect to 𝑌̄ will be as close as possible ̄ the closeness of the two to the deformation of 𝑋 with respect to 𝑋; deformations is measured by the energy (36) in terms of the elastic metric 𝐂𝑒 . We note that, as the new configuration 𝑌 is found by minimizing an elastic energy, it will be stressed (due to self-stresses) if 𝐂𝑜𝑦 is not compatible, and the solution of (37) will depend on how the elastic energy gauges the incompatibility. This implies that also the choice of the elastic moduli may be important; see discussion in Appendix E In general, being a target metric not compatible [8,14,18], it can be useful define an error function 𝜌, measuring a posteriori the transport error: ‖𝐔𝑦𝑦 ̄ ‖ 𝜌= − 1, (38) ‖𝐔𝑥𝑥 ̄ ‖

Thus, the elastic metric 𝐂𝑒𝑦̄ is given by −1 −1 𝑇 −1 𝑇 𝐂𝑒𝑦̄ = 𝐅𝑇𝑒𝑦̄ 𝐅𝑒𝑦̄ = 𝐅−𝑇 𝑜 𝐂 𝐅𝑜 = ( 𝐑𝑦̄ 𝐔𝑥𝑥 ̄ 𝐑𝑦̄ ) 𝐂 ( 𝐑𝑦̄ 𝐔𝑥𝑥 ̄ 𝐑𝑦̄ ) ,

(25)

𝐅𝑇

with 𝐂 = 𝐅; in the last equality the parentheses are used only to emphasize the structure of 𝐂𝑒𝑦̄ . It is then convenient to pull-back the elastic metric on the common reference configuration 𝑍 ( ) −1 𝐂𝑒 = 𝐑𝑦𝑇̄ 𝐂𝑒𝑦̄ 𝐑𝑦̄ = 𝐔𝑥𝑥 𝐑𝑦𝑇̄ 𝐂 𝐑𝑦̄ 𝐔−1 (26) ̄ 𝑥𝑥 ̄ . The term among parentheses in (26) may be worked out to obtain a representation having the same structure of 𝐂𝑥𝑥 ̄ as in (19); given that 𝐅 = 𝐅𝑦 𝐅𝑦−1 ̄ , and using the polar decomposition of 𝐅𝑦̄ , the metric 𝐂 can be rewritten as −1 −1 −1 𝑇 𝐂 = 𝐅𝑇 𝐅 = 𝐅−𝑇 𝑦̄ 𝐂𝑦 𝐅𝑦̄ = 𝐑𝑦̄ 𝐔𝑦̄ 𝐂𝑦 𝐔𝑦̄ 𝐑𝑦̄ ,

(27)

where ‖𝐔‖ = (𝐔 ⋅ 𝐔)1∕2 is the norm of the tensor 𝐔. In practical terms, because the assigned distortion is a change of metric, and we transport this change onto a shape with a different intrinsic metric, then the error (38) gives a quantification of the difference in the change of metric. We note that 𝜌 must not be interpreted as a computational error but rather as an intrinsic geometrical error.

𝐅𝑇𝑦

with 𝐂𝑦 = 𝐅𝑦 . This prompts the definition of a right stretch tensor 𝐔𝑦𝑦 ̄ , similar to (20): 2 −1 −1 (𝐔𝑦𝑦 ̄ ) = 𝐔𝑦̄ 𝐂𝑦 𝐔𝑦̄ ,

(28)

for which the same remarks done for (20) hold. Eventually, the elastic metric 𝐂𝑒 can be given a representation which highlights its geometric meaning: −1 2 −1 −1 𝑇 −1 𝑇 𝐂𝑒 = 𝐔𝑥𝑥 ̄ 𝐔𝑥𝑥 ̄ 𝐔𝑥𝑥 ̄ 𝐔𝑦𝑦 ̄ 𝐔𝑥𝑥 ̄ = (𝐔𝑦𝑦 ̄ ) (𝐔𝑦𝑦 ̄ ) = 𝐅𝑒 𝐅𝑒 .

4. The shell-like bodies

(29) The algorithm for the metric’s transport developed in the previous section may be used to the full hierarchy of geometric manifolds, that is, for 1D curves, 2D surfaces, and 3D solids, provided the appropriate metric information and elastic energy is used. In this section, starting from the results for 3D manifolds, we derive the transport of metric between 2D manifolds; the key point is the construction of the first and the second fundamental forms of a surface by using the information of the 3D elastic metric 𝐂𝑒 . In computer graphics bodies are often represented as surfaces in a 3D space. In continuum mechanics, a surface whose mechanical response is sensitive to both its first and second fundamental form, is named a shell. In particular, for the so called Koiter–Shell model [19], the mechanical response of a shell is determined by the sum a membrane energy (𝜑𝑚 ) and a bending energy (𝜑𝑏 ). The first one gauges the local changes in the tangent plane metric (first fundamental form), while the second one gauges the local changes in the curvatures (second fundamental form). For very small thickness and very flexible materials, the bending energy can be neglected, and the shell is named a membrane.

The tensor 𝐅𝑒 , defined by −1 𝐅𝑒 = 𝐔𝑦𝑦 ̄ 𝐔𝑥𝑥 ̄ ,

(30)

measures, on the common reference domain 𝑍, the difference between the actual stretch 𝐔𝑦𝑦 ̄ and the target stretch 𝐔𝑥𝑥 ̄ . From (29) it follows that condition (6) rewrites as 𝐂𝑒 = 𝐈

⇔

2 𝐔2𝑦𝑦 ̄ = 𝐔𝑥𝑥 ̄ .

(31)

As done in Section 2, we introduce the notion of target metric 𝐂𝑜 for 𝑌̄ , an endomorphism of 𝑇 𝑌̄ ; from (23), it follows 𝐂𝑜 = 𝐅𝑇𝑜 𝐅𝑜 = 𝐑𝑦̄ 𝐔𝑥𝑥 ̄ 𝐑𝑦̄ ;

(32)

then, we transport 𝐂𝑜 on 𝑇 𝑍, obtaining the target metric 𝐂𝑜𝑦 , an endomorphism of 𝑇 𝑍: −1 −1 𝐂𝑜𝑦 = 𝐅𝑦𝑇̄ 𝐂𝑜 𝐅−𝑇 𝑦̄ = 𝐔𝑦̄ 𝐔𝑥̄ 𝐂𝑥 𝐔𝑥̄ 𝐔𝑦̄ .

(33)

Given (33), we may write 𝐂𝑒 = 𝐈

⇔

−1 𝐂𝑦 = 𝐂𝑜𝑦 = 𝐔𝑦̄ 𝐔−1 𝑥̄ 𝐂𝑥 𝐔𝑥̄ 𝐔𝑦̄ .

(34)

When 𝐂𝑦 = 𝐂𝑜𝑦 , the metric transport is called exact, and there exist the deformation 𝑓𝑦 such that (∇𝑓𝑦 )𝑇 ∇𝑓𝑦 = 𝐂𝑦 = 𝐂𝑜𝑦 . As noted in Section 2, the necessary condition for the existence of a deformation 𝑓𝑦 which realizes a given target metric 𝐂𝑜𝑦 is that the Riemannian curvature associated to 𝐂𝑜𝑦 be null. Upon 𝐂𝑒 , we define the Green strain tensor 𝐄𝑒 :

4.1. The shell-like deformation Let {𝑜; 𝐞1 , 𝐞2 , 𝐞3 } be an orthonormal frame of E ; to denote the Cartesian components, we shall use the subscripts 𝛼, 𝛽 = 1, 2 and 𝑖, 𝑗 = 1, 2, 3. We assume as reference configuration 𝛺 of a shell-like region the Cartesian product 𝑍 × ℎ, with 𝑍 ∈ span(𝐞1 , 𝐞2 ) and ℎ ∈ span(𝐞3 ), whose characteristic diameter is much larger than the thickness ℎ: √ area(𝑍)∕ℎ ≫ 1. Let (𝑧1 , 𝑧2 , 𝜁) be the coordinates of a point ∈ 𝛺; following [19], we represent the shell-like deformation 𝑓𝑘 ∶ 𝛺 → E as the sum of two fields defined on 𝑍, describing the mid-surface 𝑓̂𝑘 and its normal vector 𝐧𝑘 ; we have:

) 1( 𝐂𝑒 − 𝐈 2 (35) ( ) 1 −1 2 − 𝐔2 = 𝐔𝑥𝑥 𝐔𝑦𝑦 𝐔−1 ; ̄ ̄ 𝑥𝑥 ̄ 𝑥𝑥 ̄ 2 moreover, we define the elastic energy density as: ) 1( 𝜑𝑒 = C𝐄𝑒 ⋅ 𝐄𝑒 𝐽𝑜 𝐽𝑦 . (36) 2 where the elastic tensor C ∶ Sym → Sym is a linear map describing the stiffness of the body. It is worth noting that when 𝐂𝑒 = 𝐈, we have 𝐄𝑒 = 0 and the elastic energy achieves a minimum at zero.

𝐄𝑒

(37)

=

𝑓𝑘 (𝑧1 , 𝑧2 , 𝜁) = 𝑓̂𝑘 (𝑧1 , 𝑧2 ) + 𝜁 𝐧𝑘 (𝑧1 , 𝑧2 ) ,

(39)

where, as before, 𝑘 = 𝑥, ̄ 𝑥, 𝑦, ̄ 𝑦, refers to the four configurations we consider. 4

L. Teresi, F. Milicchio, S. Gabriele et al.

International Journal of Non-Linear Mechanics 119 (2020) 103326

The tangent vectors 𝐚𝑘1 , 𝐚𝑘2 to the surface 𝑓̂𝑘 are defined in terms of its partial derivatives: 𝐚𝑘1 = 𝑓̂𝑘,1 = 𝜕 𝑓̂𝑘 ∕𝜕𝑧1 ;

𝐚𝑘2 = 𝑓̂𝑘,2 = 𝜕 𝑓̂𝑘 ∕𝜕𝑧2 .

2 , and the terms of 𝑜(𝜁). At first, we need a linear expansion of 𝐔2𝑦𝑦 ̄ , 𝐔𝑥𝑥 ̄ −1 𝐔𝑥𝑥 ̄ ; we have

(40)

2 = 𝐔−1 𝐂 𝐔−1 ≐ 𝐀 + 2 𝜁 𝐁 𝐔𝑦𝑦 𝑦 𝑦̄ 𝑦𝑦 ̄ 𝑦𝑦 ̄ ̄ 𝑦̄

The vector field 𝐧𝑘 , the normal to the surface, is defined by 𝐚 × 𝐚𝑘2 . 𝐧𝑘 = 𝑘1 ‖𝐚𝑘1 × 𝐚𝑘2 ‖

𝑓̂𝑘1,2 𝑓̂𝑘2,2 𝑓̂𝑘3,2

𝑛𝑘1 ⎤ ⎥ 𝑛𝑘2 ⎥ + 𝜁 ⎥ 𝑛𝑘2 ⎦

⎡𝑛𝑘1,1 ⎢ ⎢𝑛𝑘2,1 ⎢ ⎣𝑛𝑘3,1

𝑛𝑘1,2 𝑛𝑘2,2 𝑛𝑘3,2

0⎤ ⎥ 0⎥ ⎥ 0⎦

=

𝐅̂ 𝑇𝑘 𝐅̂ 𝑘 + 𝐅̂ 𝑇𝑘 (𝐧𝑘 ⊗ 𝐞3 ) + (𝐞3 ⊗ 𝐧𝑘 ) 𝐅̂ 𝑘

+

𝜁 ( 𝐅̂ 𝑇𝑘 ∇𝐧𝑘 + ∇𝐧𝑇𝑘 𝐅̂ 𝑘 )

+

𝜁 2 ∇𝐧𝑇𝑘 ∇𝐧𝑘 .

(42)

−1 𝐔−1 ̄ , 𝑥𝑥 ̄ = 𝐔𝑥𝑥𝑜 ̄ − 𝜁 𝐃𝑥𝑥

(43)

is the 1st fund. form; is null, being 𝐚𝑘𝛼 ⊥ 𝐧;

𝜁 2 ∇𝐧𝑇𝑘 ∇𝐧𝑘 = 𝑜(𝜁)

−1 − 𝜁 𝐃 ) (𝐀 + 2 𝜁 𝐁 ) (𝐔−1 − 𝜁 𝐃 ) = (𝐔𝑥𝑥𝑜 𝑥𝑥 ̄ 𝑦𝑦 ̄ 𝑦𝑦 ̄ 𝑥𝑥 ̄ ̄ 𝑥𝑥𝑜 ̄ −1 𝐀 𝐔−1 + 2 𝜁 ( 𝐔−1 𝐁 𝐔−1 − sym(𝐃 𝐀 𝐔−1 ) ) = 𝐔𝑥𝑥𝑜 𝑦𝑦 ̄ 𝑦𝑦 ̄ 𝑥𝑥 ̄ 𝑦𝑦 ̄ ̄ 𝑥𝑥𝑜 ̄ 𝑥𝑥𝑜 ̄ 𝑥𝑥𝑜 ̄ 𝑥𝑥𝑜 ̄ [ ] −1 −1 −1 𝐔 = 𝐔𝑥𝑥𝑜 𝐀 + 2 𝜁 ( 𝐁 − sym(𝐀 𝐔 )) 𝐔 𝑦𝑦 ̄ 𝑦𝑦 ̄ 𝑦𝑦 ̄ 𝑥𝑥1 ̄ 𝑥𝑥𝑜 ̄ 𝑥𝑥𝑜 ̄ ̄

At this point, we can identify the elastic fundamental forms 𝐀𝑒 and 𝐁𝑒 : −1 𝐀 𝐔−1 , = 𝐔𝑥𝑥𝑜 𝑦𝑦 ̄ ̄ 𝑥𝑥𝑜 ̄ ) −1 ( −1 −1 𝐔 = 𝐔𝑥𝑥𝑜 𝐁 − sym(𝐀𝑦𝑦 𝑥𝑥1 ̄ ) 𝐔𝑥𝑥𝑜 𝑦𝑦 ̄ ̄ 𝐔𝑥𝑥𝑜 ̄ ̄ ̄ .

𝐀𝑒

(44)

(45)

where the symbol ≐ denotes a truncation of 𝑜(𝜁). Being 𝐂𝑘 = assume a similar representation for 𝐔𝑘 : 𝐔𝑘 ≐ 𝐔𝑘𝑜 + 𝜁 𝐔𝑘1 .

=

we can write √ 𝐔𝑘𝑜 = 𝐀𝑘 , 𝐁𝑘 = sym(𝐔𝑘𝑜 𝐔𝑘1 ) .

1 2

(56)

(𝐀𝑒 − 𝐈) + 𝜁 𝐁𝑒 = 𝐄𝑒𝑜 + 𝜁 𝐄𝑒1 .

1 2 1 2

−1 − 𝜁 𝐃 ) (𝐀 − 𝐀 + 2 𝜁 (𝐁 − 𝐁 )) (𝐔−1 − 𝜁 𝐃 ) (𝐔𝑥𝑥𝑜 𝑥𝑥 ̄ 𝑦𝑦 ̄ 𝑥𝑥 ̄ 𝑦𝑦 ̄ 𝑥𝑥 ̄ 𝑥𝑥 ̄ ̄ 𝑥𝑥𝑜 ̄

(57)

−1 (𝐀 − 𝐀 ) 𝐔−1 𝐔𝑥𝑥𝑜 𝑦𝑦 ̄ 𝑥𝑥 ̄ ̄ 𝑥𝑥𝑜 ̄ ) ] −1 [ ( −1 𝐔 𝐔𝑥𝑥𝑜 +𝜁 𝐔−1 (𝐁 − 𝐁 ) − sym (𝐀𝑦𝑦 𝑦𝑦 ̄ 𝑥𝑥 ̄ 𝑥𝑥1 ̄ ̄ − 𝐀𝑥𝑥 ̄ ) 𝐔𝑥𝑥𝑜 𝑥𝑥𝑜 ̄ ̄ ̄

=

Comparing (57) with (56), we have (46)

1 −1 , (𝐀𝑒 − 𝐈) = 12 𝐔−1 𝐄𝑒1 = 𝐁𝑒 , ̄ − 𝐀𝑥𝑥 ̄ ) 𝐔𝑥𝑥𝑜 𝑥𝑥𝑜 ̄ (𝐀𝑦𝑦 ̄ 2 [ ( )] −1 −1 𝐔𝑥𝑥𝑜 (𝐁𝑦𝑦 ̄ − 𝐁𝑥𝑥 ̄ ) − sym (𝐀𝑦𝑦 ̄ − 𝐀𝑥𝑥 ̄ ) 𝐔𝑥𝑥𝑜 ̄ ̄ ̄ 𝐔𝑥𝑥1

𝐄𝑒𝑜 = 𝐁𝑒 =

𝐔−1 𝑥𝑥𝑜 ̄

(58)

The apparent diversity of the representations (55) and (58) of 𝐁𝑒 can 2 be worked out by considering that 𝐀𝑥𝑥 and using the analogous ̄ = 𝐔𝑥𝑥𝑜 ̄ of (50)2 for 𝐁𝑥𝑥 , from which it follows that ̄ ( ( ) ) −1 sym 𝐀𝑥𝑥 = sym 𝐔𝑥𝑥𝑜 = 𝐁𝑥𝑥 (59) ̄ 𝐔𝑥𝑥𝑜 ̄ 𝐔𝑥𝑥1 ̄ ̄ ; ̄ ̄ 𝐔𝑥𝑥1

(47) 𝐔2𝑘 ,

we

(48)

thus, the two representations coincide. Remembering that 𝐔𝑥𝑥𝑜 = ̄ √ 𝐀𝑥𝑥 , we can further manipulate 𝐄 and 𝐄 ; we have ̄ 𝑒𝑜 𝑒1

Then, from the equality 2 2 (𝐔𝑘0 ̄ + 𝜁 𝐔𝑘1 ̄ ) = 𝐔 ̄ + 𝜁 sym(𝐔𝑘𝑜 𝐔𝑘1 ) = 𝐀𝑘̄ + 2 𝜁 𝐁𝑘̄ , 𝑘0

(𝐂𝑒 − 𝐈) =

Then, we use the representation (35), together with (52) and (53) which hold for the 2D case, to highlights the geometric meaning of the two summands in 𝐄𝑒 : ( ) 2 𝐔2𝑦𝑦 𝐔−1 𝐄𝑒 = 12 𝐔−1 𝑥𝑥 ̄ ̄ − 𝐔𝑥𝑥 ̄ 𝑥𝑥 ̄

Given (45), we shall represent 𝐂𝑘 as a linear function in 𝜁: 𝐂𝑘 ≐ 𝐀𝑘 + 2 𝜁 𝐁𝑘 ,

1 2

𝐄𝑒 =

negligible for small thickness.

(𝐧 ⋅ 𝐚𝛽 ),𝛼 = 0 ⇒ 𝐧,𝛼 ⋅ 𝐚𝛽 = −𝐧 ⋅ 𝐚𝛽,𝛼 ; 𝐧,𝛼 ⋅ 𝐚𝛽 = −𝐧 ⋅ 𝐚𝛽,𝛼 = −𝐧 ⋅ 𝑓̂, 𝛽 𝛼 = −𝐧 ⋅ 𝑓̂, 𝛼 𝛽 = 𝐧,𝛽 ⋅ 𝐚𝛼 .

(55)

Using the same procedure, we compute 𝐄𝑒 ; at first, we exploit equation (54) and write

is the 2nd fund. form;

⇒

(54)

= 𝐀𝑒 + 2 𝜁 𝐁𝑒 .

To prove the previous relations we shall drop the subscript 𝑘 to simplify notation; the proof of (45)1 is straightforward, as the components of 𝐀 are given by 𝐴𝛼𝛽 = 𝐚𝛼 ⋅ 𝐚𝛽 . To prove (45)3 , we start from the definition, ̂ and we note that (𝐅̂ 𝑇 ∇𝐧)𝑇 = ∇𝐧𝑇 𝐅̂ = 𝐁. Let us remember 𝐁 = ∇𝐧𝑇 𝐅, that the symmetry of 𝐁 follows from the fact that the second derivatives of 𝑓̂ commute: 𝐧 ⋅ 𝐚𝛼 = 0

(53)

−1 𝐔2 𝐔−1 = 𝐔𝑥𝑥 ̄ 𝑦𝑦 ̄ 𝑥𝑥 ̄

𝐂𝑒

The different summands in (44) deserve further remarks:

𝐅̂ 𝑇𝑘 (𝐧𝑘 ⊗ 𝐞3 ) + (𝐞3 ⊗ 𝐧𝑘 ) 𝐅̂ 𝑘 = 0 ( ) 1 ̂𝑇 𝐁𝑘 ∶= 𝐅𝑘 ∇𝐧𝑘 + ∇𝐧𝑇𝑘 𝐅̂ 𝑘 2

−1 −1 with 𝐃𝑥𝑥 ̄ = 𝐔𝑥𝑥𝑜 ̄ 𝐔𝑥𝑥𝑜 ̄ 𝐔𝑥𝑥1 ̄ .

Then, we can compute the linear expansion of 𝐂𝑒 ; we have:

𝐁𝑒

𝐀𝑘 ∶= 𝐅̂ 𝑇𝑘 𝐅̂ 𝑘

≐ 𝐔𝑥𝑥𝑜 ̄ + 𝜁 𝐔𝑥𝑥1 ̄

The identification of the six tensors 𝐀𝑦𝑦 from ̄ , 𝐀𝑥𝑥 ̄ , 𝐁𝑦𝑦 ̄ , 𝐁𝑥𝑥 ̄ , 𝐔𝑥𝑥𝑜 ̄ , 𝐔𝑥𝑥1 ̄ the four tensor 𝐂𝑘 as given in (47) is reported in the Appendix. It is worth noting that the maps 𝑓𝑘 → 𝐂𝑘 are straightforward, while the −1 maps 𝐂𝑥̄ → 𝐔−1 𝑥̄ and 𝐂𝑦̄ → 𝐔𝑦̄ involve the square root and the inversion of a tensor. From (52)3 , it follows

with 𝑓̂𝑘𝑖 = 𝑓̂𝑘 ⋅ 𝐞𝑖 , and 𝑛𝑘𝑖 = 𝐧𝑘 ⋅ 𝐞𝑖 the Cartesian components of 𝑓̂𝑘 , 𝐧𝑘 , respectively. It is worth noting that both 𝐅̂ 𝑘 and ∇𝐧𝑘 may be represented by 3 × 2 matrices; moreover, given (41), the gradient 𝐅𝑘 contains both first and second derivatives of the surface 𝑓̂𝑘 . We also note that the first two columns of the first summands in (43) are the tangent vectors 𝐚𝑘1 , 𝐚𝑘2 . The metric 𝐂𝑘 associated to 𝐅𝑘 is represented by 𝐂𝑘 = 𝐅𝑇𝑘 𝐅𝑘

1∕2 𝐂𝑥 𝐔−1 𝑥̄ )

𝐔𝑥𝑥 ̄ =

with 𝐅̂ 𝑘 = ∇𝑓̂𝑘 . The matrix-like representation |𝐅𝑘 | of 𝐅𝑘 is given by ⎡𝑓̂𝑘1,1 ⎢ |𝐅𝑘 | = ⎢𝑓̂𝑘2,1 ⎢̂ ⎣𝑓𝑘3,1

(𝐔−1 𝑥̄

(41)

Given (39), its gradient 𝐅𝑘 = ∇𝑓𝑘 is represented by 𝐅𝑘 = 𝐅̂ 𝑘 + 𝐧𝑘 ⊗ 𝐞3 + 𝜁 ∇𝐧𝑘 .

(52)

−1 −1 𝐔2𝑥𝑥 ̄ + 2 𝜁 𝐁𝑥𝑥 ̄ ̄ = 𝐔𝑥̄ 𝐂𝑥 𝐔𝑥̄ ≐ 𝐀𝑥𝑥

(49)

𝐄𝑒𝑜 𝐄𝑒1

(50)

The identity (50)2 can be solved for 𝐔𝑘1 , and yields this noteworthy expression ( ) 1 −1 (51) 𝐔𝑘1 = 𝐁𝑘 + det(𝐔𝑘𝑜 ) 𝐔−1 𝑘𝑜 𝐁𝑘 𝐔𝑘𝑜 . tr(𝐔𝑘𝑜 )

1 −1∕2 −1∕2 𝐀̄ (𝐀𝑦𝑦 , ̄ − 𝐀𝑥𝑥 ̄ ) 𝐀𝑥𝑥 ̄ 2 𝑥𝑥 −1∕2 −1∕2 = 𝐀𝑥𝑥 (𝐁𝑦𝑦 ̄ − 𝐁𝑥𝑥 ̄ ) 𝐀𝑥𝑥 ̄ ̄ ( ) −1∕2 1∕2 −1∕2 −1∕2 − 2 𝐀𝑥𝑥 sym 𝐀 𝐄 𝐔𝑥𝑥1 𝐀𝑥𝑥 𝑒𝑜 𝐀𝑥𝑥 ̄ ̄ 𝑥𝑥 ̄ ̄ ̄

=

−1∕2

−1∕2

= 𝐀𝑥𝑥 (𝐁𝑦𝑦 ̄ − 𝐁𝑥𝑥 ̄ ) 𝐀𝑥𝑥 ̄ ̄ ( ) 2 −1∕2 −1∕2 − sym 𝐄𝑒𝑜 𝐀𝑥𝑥 𝐁𝑥𝑥 ̄ 𝐀𝑥𝑥 ̄ ̄ tr(𝐔𝑥𝑥𝑜 ̄ ) ( ) 2 det(𝐔𝑥𝑥𝑜 ̄ ) −1 . − sym 𝐄𝑒𝑜 𝐀−1 ̄ 𝐀𝑥𝑥 𝑥𝑥 ̄ 𝐁𝑥𝑥 ̄ tr(𝐔𝑥𝑥𝑜 ̄ )

We now replicate the construction of the elastic metric 𝐂𝑒 and the Green strain 𝐄𝑒 , following the definitions (29), (35), and neglecting all 5

(60)

L. Teresi, F. Milicchio, S. Gabriele et al.

International Journal of Non-Linear Mechanics 119 (2020) 103326

4.2. Transport error We propose some norms to measure the error done in transporting the metric. A metric transport is exact when 𝐂𝑒 = 𝐈; from (54), we have 𝐂𝑒 = 𝐈

⇔

𝐔𝑥𝑥 ̄ = 𝐔𝑦𝑦 ̄ , and 𝐀𝑒 = 𝐈 , 𝐁𝑒 = 𝟎 .

(61)

Then, from (55)1 and (58)2 , it follows 𝐀𝑒 = 𝐈

⇔

𝐀𝑥𝑥 ̄ = 𝐀𝑦𝑦 ̄ ,

𝐁𝑒 = 𝟎

⇔

𝐁𝑥𝑥 ̄ = 𝐁𝑦𝑦 ̄ , and 𝐀𝑥𝑥 ̄ = 𝐀𝑦𝑦 ̄ ,

(62)

−1 𝐔 or, equivalently 𝐁𝑦𝑦 ̄ = sym(𝐀𝑦𝑦 ̄ 𝐔𝑥𝑥𝑜 𝑥𝑥1 ̄ ). ̄

Given (62), we define the following functions to measure the error done in transporting the metric: | ‖𝐀𝑦𝑦 | ̄ ‖ | | 𝜌𝐴 = | − 1| | ‖𝐀𝑥𝑥 | ‖ ̄ | | (( )2 )1∕2 −1 𝐔 ‖𝐁𝑦𝑦 ‖ − sym(𝐀𝑦𝑦 𝜌𝐵 = ℎ ̄ ̄ 𝐔𝑥𝑥𝑜 𝑥𝑥1 ̄ )∥ ̄

stretching error , (63) bending error ,

with ℎ a characteristic thickness of the shell. Definition (63)2 avoids a division by zero as the second fundamental form may be zero. 5. Worked examples: Shell-like bodies

Fig. 4. Square-like body deformation transported onto annular geometry. The 𝑍 domain is scaled down; the 𝑌 is rotated to fit into the figure. All the domains have been rendered in wireframe to highlight the change of metric between the different configurations. In particular, it can be appreciated how the orthogonal gridlines on 𝑋̄ and 𝑌̄ get distorted on 𝑋 and 𝑌 .

In this section we present some worked examples of the proposed algorithm, where it is shown the transfer of the deformation between different shell-like bodies. At first, we consider two examples of shell embedded in a flat, 2D Euclidean space, undergoing only large planar stretching; then, we consider a shell in a 3D space, having both stretching and bending deformations. We use the shell model presented in the Appendix D, and we solve all the examples by using the finite element method, implementing the weak form of the shell model into the software COMSOL Multiphysics. We discuss the effects of the elastic moduli in Appendix E. 5.1. Shell in a 2D ambient space We consider shells embedded in a 2D space, and we assume as reference domain 𝑍 ⊂ R2 a rectangular region (𝑧1, 𝑧2) ∈ [−𝐿∕2, 𝐿∕2] × [−𝐻∕2, 𝐻∕2]. 5.1.1. A square body In this example both 𝑍 and 𝑋̄ are square (𝐿 = 𝐻), while 𝑌̄ is a quarter of an annulus; the configuration 𝑋 is a non uniform scaling, whose intensity depends on the distance from the origin, see Fig. 4. ̄ 𝑋 and 𝑌̄ are defined by The three configurations 𝑋, { { 𝑥1 = 𝑧1 , 𝑥̄ 1 = (1 + 𝑔(𝑟)) 𝑧1 , 𝑋̄ ∶ 𝑋∶ 𝑥2 = 𝑧2 , 𝑥̄ 2 = (1 + 𝑔(𝑟)) 𝑧2 , { 𝑦 ̄ = 2 sin[𝑧 𝜋∕(2 𝐿)] + 𝑧2 sin[𝑧1 𝜋∕(2 𝐿)], 1 1 𝑌̄ ∶ 𝑦̄2 = 2 cos[𝑧1 𝜋∕(2 𝐿)] + 𝑧2 cos[𝑧1 𝜋∕(2 𝐿)], (64) with ( ) 2 ⎧ 𝑙 ⎪ 2 2 𝑟 = (𝑧21 + 𝑧22 )1∕2 , 𝑔(𝑟) = ⎨𝛼 𝑒 𝑟 − 𝑙 −𝑙 < 𝑟 < 𝑙, ⎪0 elsewhere. ⎩

Fig. 5. Square-like example. The stretching error 𝜌𝐴 for 𝛼 = 1 is shown both on 𝑌 (top) and on 𝑋 (bottom). Note that 𝜌𝐴 is a reference field defined on 𝑍; here we plot its push forward on both 𝑌 and 𝑋. In the background, wireframe shows the distorted gridlines.

by the stretching error 𝜌𝐴 , see (63)1 ; for this example 𝜌𝐵 = 0. Fig. 5 shows 𝜌𝐴 for the maximum value of the parameter 𝛼 = 1; it may be noted that the error is confined where the deformation gradient is very high. Fig. 6 shows the average value of 𝜌𝐴 versus 𝛼. The value of 𝜌𝐴 is monotone increasing with 𝛼: our opinion is that the parameter 𝛼 increases the amount of incompatibility of the transported metric. The proposed algorithm yields the best embedding, measured with the energy norm (37); thus, the more a metric is incompatible, the more the error is large.

The non dimensional parameter 𝛼 is used to amplify the deformation ̄ the following results are for 𝑙 = 𝐿∕2 and of 𝑋 with respect to 𝑋; 𝛼 ∈ (0, 1). Fig. 4 shows how the proposed method is able to gauge the metric change from 𝑋̄ to 𝑋 and to transport this metric on the annulus 𝑌̄ , which is quite different from the square, to obtain the sought configuration 𝑌 . This result is interesting from a qualitative point of view, since the deformation 𝑋̄ → 𝑋 of the square turns out to be similar to the one of the annulus. The similarity of the two deformations is measured

5.1.2. A rectangular body In this example both 𝑍 and 𝑋̄ are rectangle (𝐿 = 10 𝐻), while 𝑌̄ is a bended rectangle; the configuration 𝑋 is obtained from 𝑋̄ by modū lating the height along the axis, see Fig. 7. The three configurations 𝑋, 6

L. Teresi, F. Milicchio, S. Gabriele et al.

International Journal of Non-Linear Mechanics 119 (2020) 103326

Fig. 6. Square-like example. Average value of 𝜌𝐴 versus 𝛼. Being this embedding 2D, the bending error 𝜌𝐵 is null.

Fig. 8. Rectangle-like example. The stretching error 𝜌𝐴 for 𝛼 = 1 is shown both on 𝑌 (left) and on 𝑋 (right). Note that 𝜌𝐴 is a reference field defined on 𝑍; here we plot its push forward on both 𝑌 and 𝑋. In the background, wireframe shows the distorted gridlines.

Fig. 9. Rectangle-like example. Average value of 𝜌𝐴 versus 𝛼. Being this embedding 2D, the bending error 𝜌𝐵 is null.

̄ 𝑋 and 𝑌̄ are defined by 𝑋, ⎧𝑥 = 𝑧 , 1 ⎪ 1 ̄ 𝑋 ∶ ⎨𝑥2 = 𝑧2 , ⎪𝑥3 = 0, ⎩

Fig. 7. Rectangle-like body deformation transported onto the annular geometry. The domains have been rotated to fit into the figure. All the domains have been rendered in wireframe to highlight the change of metric between the different configurations. In particular, it can be appreciated how the orthogonal gridlines on 𝑋̄ and 𝑌̄ get distorted on 𝑋 and 𝑌 .

⎧𝑦̄1 = 𝑧1 , ⎪ ̄ 𝑌 ∶ ⎨𝑦̄2 = 𝑧2 , ⎪𝑦̄ = − 4 𝑅 (𝑧2 + 𝑧2 ) + 𝑅, 2 ⎩ 3 𝐿2 1

(66)

with 𝑅 = 0.5.

The function 𝑏(., .) describes the bumps in 𝑋, and is defined by means of the piecewise exponential 𝑔(.) given by

𝑋 and 𝑌̄ are defined by { { 𝑥1 = 𝑧1 , 𝑥̄ 1 = 𝑧1 + 𝛼 (𝑥1 − 𝑎(𝑧1 ) 𝑧1 𝑧2 ), 𝑋̄ ∶ 𝑋∶ 𝑥2 = 𝑧2 , 𝑥̄ 2 = 𝑧2 + 𝛼 ( 𝑎(𝑧1 ) 𝑧2 (1 − 𝑧2 )), and 𝑎(𝑥1 ) = 𝐻 + 𝛼 (10 𝑔(𝑥1 ) − 𝑔(𝑥1 + 4)) , 𝛼 ∈ (0, 1) ⎧𝑦̄ = 𝑧 + 2 𝐻 𝑧 8 𝑅 𝑧1 , 1 2 ⎪ 1 𝐿2 𝑌̄ ∶ ⎨ 2 4 𝑅 𝑧 1 ⎪𝑦̄ = − + 𝑅 + 2 𝐻 𝑧2 , with 𝑅 = 2, ⎩ 2 𝐿2

⎧𝑥̄ = 𝑧 1 ⎪ 1 𝑋 ∶ ⎨𝑥̄ 2 = 𝑧2 , ⎪𝑥̄ 3 = 𝑏(𝑧1 , 𝑧2 ), ⎩

) ( ⎧ 1 ⎪10 exp , 𝑔(𝑠) = ⎨ (𝑠∕𝑠𝑜 )2 − 1 ⎪0, ⎩

(65)

−𝑠𝑜 < 𝑠 < 𝑠𝑜 ,

(67)

elsewhere.

We set 𝑠𝑜 = 0.3. Then, the function 𝑏(𝑧1 , 𝑧2 ) is the sum of four bumps, with the parameter 𝛼 used to amplify the bump amplitude; here 𝛼 ∈ (0, 0.1). The bumps are defined by:

The function 𝑔(.) is a Gaussian pulse located at 2 with standard deviation = 0.5. Fig. 7 shows how the proposed method is able to gauge the large metric change from 𝑋̄ to 𝑋 and to transport this information the annulus 𝑌̄ , which deforms to configuration 𝑌 . In Fig. 8 we show the stretching error 𝜌𝐴 for 𝛼 = 1.

𝑏 = 𝛼 (𝑏1 + 𝑏2 + 𝑏3 + 𝑏4 ), ( ) 𝑏1 (𝑧1 , 𝑧2 ) = 𝑔 (𝑧21 + 𝑧22 )1∕2 , [( )1∕2 ] 𝑏2 (𝑧1 , 𝑧2 ) = 𝑔 (𝑧1 − 0.5)2 + (𝑧2 − 0.5)2 , [( )1∕2 ] 1 𝑏3 (𝑧1 , 𝑧2 ) = 𝑔 (𝑧1 + 0.5)2 + (𝑧2 + 0.5)2 , 2 ( ) 𝑏4 (𝑧1 , 𝑧2 ) = 𝑔 (𝑧21𝑟 + 𝑧22𝑟 )1∕2 ,

As done in the previous example, we estimate the similarity of the two deformations plotting the average value of the stretching error 𝜌𝐴 versus the parameter 𝛼, see Fig. 9.

(68)

where 𝑧1𝑟 , 𝑧2𝑟 , defining the rotated bump 𝑏4 , are given by 5.2. Shell in a 3D ambient space

𝑧1𝑟 = cos(𝜋∕4) (𝑧1 − 0.5) − sin(𝜋∕4) (𝑧2 + 0.5), 𝑧2𝑟 = sin(𝜋∕4) (𝑧1 − 0.5) − cos(𝜋∕4) (𝑧2 + 0.5).

In this example both 𝑍 and 𝑋̄ are square (𝐿 = 𝐻), while 𝑌̄ is a circular parabolic surface; the configuration 𝑋 is obtained from 𝑋̄ by modulating the height with bumps, see Fig. 7. The three configurations

(69)

Fig. 10 shows the results of the proposed algorithm; the bumps on the flat surface 𝑋 are reproduced on the curved surface 𝑌 by the metric transport. Figs. 11 and 12 show the stretching error 𝜌𝐴 , and the bending 7

L. Teresi, F. Milicchio, S. Gabriele et al.

International Journal of Non-Linear Mechanics 119 (2020) 103326

Fig. 12. Dome example: Bending error 𝜌𝐵 for 𝛼 = 0.05, on both 𝑌 (left) and 𝑋 (right).

Fig. 13. Dome example: Average values of 𝜌𝐴 and 𝜌𝐵 versus 𝛼. Fig. 10. Dome example: The bumps at 𝑋 have been transported onto 𝑌 . The domain 𝑍 has been scaled down to fit into the figure.

metric realized by the embedding. The numerical experiments clearly show good results, and the method is able to transport compatible metric changes, even in presence of very different geometries. The error function shows that the lack of compatibility is localized in small portion of the manifold, where very large deformations need to be realized to meet the demands of the target metric. The presence of this errors seems to not affect the global performance of the method.

Acknowledgment The authors acknowledge the support of the Italian Group of Mathematical Physics (GNFM-INdAM).

Fig. 11. Dome example: Stretching error 𝜌𝐴 for 𝛼 = 0.05, on both 𝑌 (left) and 𝑋 (right).

Appendix A. 3D flowchart error 𝜌𝐵 , respectively, for 𝛼 = 0.05. Fig. 13 shows the average values of 𝜌𝐴 (100 x) and 𝜌𝐵 versus 𝛼.

𝑓𝑥̄ , 𝑓𝑦̄ , 𝑓𝑥 given, 𝑓𝑦 unknown ↓

6. Conclusions

gradient ↓

This work presents a new method for transporting the deformations gauged by comparing two different configurations of a same body, onto a different body. The aim of this new transport is to obtain an automatic morphing procedure, that is able to replicate, as close as possible, the same metric change of the source deformation in the transported deformation. The morphing is obtained by assigning a distortion field, which in turns, induces a target metric: the actual configuration will be the one that minimizes the elastic energy, measuring the distance between the target metric and the actual one. This procedure could be useful in computer graphics and shape analysis, as example for physics-based modeling and character animation for the movie industry. We note that the proposed technique may be used to the full hierarchy of geometric manifolds, that is, for 1D curves, 2D surfaces, and 3D solids, provided the appropriate metric information, and elastic energy is used. We provided some simple, yet non trivial, examples to demonstrate the effectiveness of the procedure; we also defined an error function to measure the metric gap between the assigned target metric and the

𝐅𝑘 = ∇𝑓𝑘 ↓ tensor product ↓ 𝐂𝑘 = 𝐅𝑇𝑘 𝐅𝑘 ↓ square root of 𝐂𝑥̄ , 𝐂𝑦̄ 𝐔𝑥̄ =

↓ √ √ 𝐂𝑥̄ , 𝐔𝑦̄ = 𝐂𝑦̄ ↓

inverse of 𝐔𝑥̄ , 𝐔𝑦̄ ↓ −1 −1 2 (𝐔𝑦𝑦 ̄ ) = 𝐔𝑦̄ 𝐂𝑦 𝐔𝑦̄ −1 −1 2 (𝐔𝑥𝑥 ̄ ) = 𝐔𝑥̄ 𝐂𝑥 𝐔𝑥̄

↓ 8

L. Teresi, F. Milicchio, S. Gabriele et al.

International Journal of Non-Linear Mechanics 119 (2020) 103326

Appendix C. Linearization square root and inverse of 𝐔𝑥𝑥 ̄

Identification of the four tensors 𝐀𝑦𝑦 ̄ , 𝐀𝑥𝑥 ̄ , 𝐁𝑦𝑦 ̄ , 𝐁𝑥𝑥 ̄ , from the four tensor 𝐂𝑘 as given in (47)

↓ 𝐔−1 𝑥𝑥 ̄

=

(𝐔𝑥−1 ̄

−1∕2 𝐂𝑥 𝐔𝑥−1 ̄ )

𝐂𝑘 = 𝐀𝑘 + 2 𝜁 𝐁𝑘 , 𝜁-expansion of 𝐂𝑘 , 𝐀𝑘 = 𝐅̂ 𝑇 𝐅̂ 𝑘 , 𝐁𝑘 = sym(𝐅̂ 𝑇 ∇𝐧𝑘 ) .

↓ tensor product

𝑘

𝑘

↓

↓ 2 𝐔−1 𝐂𝑒 = 𝐔−1 ̄ 𝐔𝑦𝑦 ̄ 𝑥𝑥 ̄ ) (𝑥𝑥 1 2 − 𝐔2 −1 𝐄𝑒 = 𝐔−1 𝐔𝑦𝑦 𝐔𝑥𝑥 𝑥𝑥 ̄ ̄ 𝑥𝑥 ̄ ̄ 2 ↓

square root of 𝐂𝑥̄ , 𝐂𝑦̄ √ 𝐔𝑘̄ = 𝐂𝑘̄ = 𝐔𝑘0 ̄ + 𝜁 𝐔𝑘1 ̄ , only for 𝑘 = 𝑥 , 𝑦 √ ( ) 𝐀𝑘̄ = 1𝑡 𝐀𝑘̄ + 𝑠 𝐈 , 𝐔𝑘0 ̄ = √ √ with 𝑠 = det(𝐀𝑘̄ ) , 𝑡 = tr(𝐀𝑘̄ ) + 2 𝑠 , ( ) 1 −1 −1 𝐔𝑘1 𝐁𝑘̄ + det(𝐔𝑘0 , ̄ = ̄ ) 𝐔 ̄ 𝐁𝑘̄ 𝐔 ̄ 𝑘0 𝑘0 tr(𝐔𝑘𝑜 ̄ ) ↓

(A.1)

assemble energy density ↓ 𝜑 = 𝜑(𝐂𝑒 ) ↓

inverse of 𝐔𝑥̄ , 𝐔𝑦̄ −1 𝐔−1 ̄ = 𝐔 ̄ − 𝜁 𝐃𝑘̄ ,

minimization

𝑘

↓

𝑘𝑜

−1 𝐃𝑘̄ = 𝐔−1 ̄ 𝐔̄ ̄ 𝐔𝑘1

(C.1)

𝑘𝑜

𝑘𝑜

↓

𝑓𝑦

2 𝜁 expansion of 𝐔2𝑥𝑥 ̄ , 𝐔𝑦𝑦 ̄ −1 𝐔2̄ = 𝐔−1 ̄ 𝐂𝑘 𝐔 ̄

Appendix B. 2D flowchart

𝑘𝑘

𝑘

𝑘

−1 = (𝐔−1 ̄ − 𝜁 𝐃𝑘̄ ) (𝐀𝑘 + 2 𝜁 𝐁𝑘 ) (𝐔 ̄ − 𝜁 𝐃𝑘̄ ) 𝑘𝑜

𝑘𝑜

= 𝐀𝑘𝑘 ̄ ̄ + 2 𝜁 𝐁𝑘𝑘

𝑓𝑘 = 𝑓̂𝑘 + 𝜁 𝐧𝑘 𝜁-expansion

↓

↓

−1 𝐀 𝐔−1 𝐀𝑥𝑥 ̄ = 𝐔𝑥𝑜 𝑥 𝑥𝑜 ̄ ̄ ( ) −1 𝐁 − sym(𝐀 𝐔−1 𝐔 ) 𝐔−1 𝐁𝑥𝑥 ̄ = 𝐔𝑥𝑜 𝑥 𝑥 𝑥𝑜 𝑥1 ̄ 𝑥𝑜 ̄ ̄ ̄

gradient ↓

−1 𝐀 𝐔−1 𝐀𝑦𝑦 ̄ = 𝐔𝑦𝑜 ̄ (𝑦 𝑦𝑜 ̄ ) −1 −1 𝐁𝑦𝑦 𝐁𝑦 − sym(𝐀𝑦 𝐔−1 ̄ = 𝐔𝑦𝑜 ̄ ) 𝐔𝑦𝑜 ̄ 𝑦𝑜 ̄ 𝐔𝑦1 ̄

𝐅𝑘 = 𝐅̂ 𝑘 + 𝐧𝑘 ⊗ 𝐞3 + 𝜁 ∇𝐧𝑘 ↓ tensor product

Consequences of small membranal deformation:

↓

𝐀𝑘̄ = 𝐈

𝐂𝑘 = 𝐀𝑘 + 2 𝜁 𝐁𝑘 , with 𝐀𝑘 = 𝐅̂ 𝑇𝑘 𝐅̂ 𝑘 , 𝐁𝑘 = sym(𝐅̂ 𝑇𝑘 ∇𝐧𝑘 ) ↓

𝑘

↓

↓

√ 𝐀𝑘 + 2 𝜁 𝐁𝑘 = 𝐔𝑘0 + 𝜁 𝐔𝑘1 , only for 𝑘 = 𝑥̄ , 𝑦̄

−1 𝐔2̄ = 𝐔−1 ̄ 𝐂𝑘 𝐔 ̄ = (𝐈 − 𝜁 𝐁𝑘̄ ) (𝐀𝑘 + 2 𝜁 𝐁𝑘 ) (𝐈 − 𝜁 𝐁𝑘̄ ) 𝑘𝑘

𝑘

↓

↓

𝐀𝑥𝑥 ̄ = 𝐀𝑥 ,

−1 −1 2 (𝐔𝑦𝑦 ̄ ) = 𝐔𝑦̄ 𝐂𝑦 𝐔𝑦̄ = 𝐀𝑦𝑦 ̄ + 2 𝜁 𝐁𝑦𝑦 ̄

𝐀𝑦𝑦 ̄ = 𝐀𝑦 ,

(𝐔𝑥𝑥 ̄

=

𝐔𝑥−1 ̄

𝐂𝑥 𝐔𝑥−1 ̄

(C.2)

𝑘

= 𝐀𝑘𝑘 ̄ ̄ + 2 𝜁 𝐁𝑘𝑘

↓ inverse of 𝐔𝑥̄ , 𝐔𝑦̄

)2

𝐔𝑘1 ̄ = 𝐁𝑘̄ ,

𝐔−1 ̄ = 𝐈 − 𝜁 𝐁𝑘̄ ,

square root of 𝐂𝑥̄ , 𝐂𝑦̄ 𝐔𝑘 =

𝐔𝑘0 ̄ = 𝐈,

⇒

↓

= 𝐀𝑥𝑥 ̄ + 2 𝜁 𝐁𝑥𝑥 ̄

(B.1)

( ) 𝐁𝑥𝑥 ̄ = 𝐁𝑥 − sym(𝐀𝑥 𝐁𝑥̄ ) , ( ) 𝐁𝑦𝑦 ̄ = 𝐁𝑦 − sym(𝐀𝑦 𝐁𝑦̄ ) ,

Identification of the three tensors 𝐀𝑒 , 𝐁𝑒 , 𝐄𝑒𝑜 from the two tensors 𝐔𝑥𝑥 ̄ , 𝐔𝑦𝑦 ̄ . We follow the same procedure as before:

↓ square root and inverse of 𝐔𝑥𝑥 ̄

2 𝐔2̄ = 𝐀𝑘𝑘 ̄ + 2 𝜁 𝐁𝑘𝑘 ̄ , 𝜁-expansion of 𝐔 ̄ 𝑘𝑘

↓

𝑘𝑘

↓

−1∕2 = 𝐔−1 − 𝜁 𝐃 𝐔−1 ̄ + 2 𝜁 𝐁𝑥𝑥 ̄ ) 𝑥𝑥 ̄ 𝑥𝑥 ̄ = (𝐀𝑥𝑥 𝑥𝑥𝑜 ̄

2 square root of 𝐔𝑥𝑥 ̄ √ √ 2 𝐔𝑥𝑥 = 𝐔 = 𝐀𝑥𝑥 ̄ ̄ + 2 𝜁 𝐁𝑥𝑥 ̄ = 𝐔𝑥𝑥𝑜 ̄ + 𝜁 𝐔𝑥𝑥1 ̄ 𝑥𝑥 ̄ √ 𝐔𝑥𝑥𝑜 = 𝐀𝑥𝑥 ̄ ̄ , ( ) 1 −1 𝐁 𝐔−1 𝐔𝑥𝑥1 = 𝐁𝑥𝑥 ̄ ̄ + det(𝐔𝑥𝑥𝑜 ̄ ) 𝐔𝑥𝑥0 𝑥𝑥 ̄ 𝑥𝑥𝑜 ̄ ̄ tr(𝐔𝑥𝑥𝑜 ) ̄ ↓

↓ tensor product ↓ 2 𝐔−1 = 𝐀 + 2 𝜁 𝐁 𝐂𝑒 = 𝐔−1 ̄ ) ) 𝑒 𝑒 ̄ (𝐔𝑦𝑦 𝑥𝑥 ̄ ( 𝑥𝑥 1 −1 2 2 −1 = 1 (𝐀 − 𝐈) + 𝜁 𝐁 𝐄𝑒 = 𝐔𝑥𝑥 𝐔 − 𝐔 𝐔 𝑒 𝑒 𝑦𝑦 ̄ 𝑥𝑥 ̄ 𝑥𝑥 ̄ 2 2 ̄ ↓

inverse of 𝐔𝑥𝑥 ̄ , −1 − 𝜁 𝐃 , 𝐔−1 𝑥𝑥 ̄ 𝑥𝑥 ̄ = 𝐔𝑥𝑥𝑜 ̄

assemble energy density

−1 𝐔 −1 𝐃𝑥𝑥 ̄ = 𝐔𝑥𝑥𝑜 𝑥𝑥1 ̄ 𝐔𝑥𝑥𝑜 ̄ ̄

↓

↓

𝜁 expansion of 𝐂𝑒 (same as expansion of 𝐔2̄ )

𝜑 = 𝜑(𝐂𝑒 ) = 𝜑(𝐀𝑒 , 𝐁𝑒 )

𝑘𝑘

2 𝐔−1 𝐂𝑒 = 𝐔−1 𝑥𝑥 ̄ 𝐔𝑦𝑦 ̄ 𝑥𝑥 ̄

↓

−1 − 𝜁 𝐃 ) (𝐀 + 2 𝜁 𝐁 ) (𝐔−1 − 𝜁 𝐃 ) = (𝐔𝑥𝑥𝑜 𝑥𝑥 ̄ 𝑦𝑦 ̄ 𝑦𝑦 ̄ 𝑥𝑥 ̄ ̄ 𝑥𝑥𝑜 ̄

minimization

= 𝐀𝑒 + 2 𝜁 𝐁𝑒

↓ 𝑓̂𝑦

↓ 9

L. Teresi, F. Milicchio, S. Gabriele et al.

International Journal of Non-Linear Mechanics 119 (2020) 103326

−1 𝐀𝑒 = 𝐔−1 ̄ 𝐔𝑥𝑥𝑜 𝑥𝑥𝑜 ̄ 𝐀𝑦𝑦 ̄ ( ) −1 −1 𝐔 𝐁𝑒 = 𝐔−1 𝐁𝑦𝑦 ̄ − sym(𝐀𝑦𝑦 ̄ 𝐔𝑥𝑥𝑜 𝑥𝑥1 ̄ ) 𝐔𝑥𝑥𝑜 𝑥𝑥𝑜 ̄ ̄ ̄

↓ 𝜁 expansion of 𝐄𝑒 = 1∕2 (𝐂𝑒 − 𝐈) 𝐄𝑒 = 21 (𝐂𝑒 − 𝐈) = 12 (𝐀𝑒 − 𝐈) + 𝜁 𝐁𝑒 ( ) 1 −1 = 𝐔−1 𝐔2𝑦𝑦 − 𝐔2𝑥𝑥 𝐔𝑥𝑥 𝑥𝑥 ̄ ̄ ̄ ̄ 2 = 𝐄𝑒𝑜 + 𝜁 𝐄𝑒1 𝐄𝑒1 = 𝐁𝑒

(C.3) Fig. E.14. Rectangle-like example. Deformed configurations resulting from metric’s transport for different values of 𝑘∕𝜇. Left (blue): highly compressible material with 𝑘∕𝜇 = 10−2 ; the change of shape is transported very well, but the area change is overshoot. Center (black): 𝑘∕𝜇 = 1, our benchmark; change of shape and area are well balanced. Right (red): highly incompressible material with 𝑘∕𝜇 = 102 ; the incompressibility is so high that our algorithm fails. It may be noticed in the central part a mesh inversion. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

↓ −1 , (𝐀𝑒 − 𝐈) = 12 𝐔−1 ̄ − 𝐀𝑥𝑥 ̄ ) 𝐔𝑥𝑥𝑜 𝑥𝑥𝑜 ̄ (𝐀𝑦𝑦 ̄ ) ] −1 [ ( −1 −1 𝐔 𝐔𝑥𝑥𝑜 = 𝐔𝑥𝑥𝑜 (𝐁𝑦𝑦 ̄ − 𝐁𝑥𝑥 ̄ ) − sym (𝐀𝑦𝑦 ̄ − 𝐀𝑥𝑥 ̄ ) 𝐔𝑥𝑥𝑜 𝑥𝑥1 ̄ ̄ ̄ ̄ ( ) −1 −1 𝐔 −1 = 𝐔𝑥𝑥𝑜 𝐁 − sym(𝐀 𝐔 ) 𝐔 𝑦𝑦 ̄ 𝑦𝑦 ̄ 𝑥𝑥1 ̄ ̄ 𝑥𝑥𝑜 ̄ 𝑥𝑥𝑜 ̄ 1 2

𝐄𝑒𝑜 = 𝐄𝑒1

Appendix D. Shell model The variation of 𝜓𝑟 identifies the dynamical actions We summarize the Naghdi–Shell model [20]. The state variables of the model are 𝐮 = (𝑢1 , 𝑢2 , 𝑢3 )

= (𝛽1 , 𝛽2 )

= (𝛾1 , 𝛾2 )  = (𝜏1 , 𝜏2 )

𝜓 ̃𝑟

displacements;

𝐚𝑦1 ̄ ‖𝐚𝑦1 ̄ ‖

(

+

(D.1)

shear stretch;

𝐞2 =

𝐚𝑦2 ̄ ‖𝐚𝑦2 ̄ ‖

,

shear stress

(D.2)

𝐧𝑦̄ = 𝐞1 × 𝐞2 ;

𝜇1 = ℎ

𝑓̂𝑦 = 𝑓̂𝑦̄ + 𝑢1 𝐞1 + 𝑢2 𝐞2 + 𝑢3 𝐧𝑦̄ .

𝑏1 =

‖𝐧𝑦̄ + 𝛽1 𝐞1 + 𝛽2 𝐞2 ‖

,



𝑑 = 𝐅̂ 𝑇𝑦 𝐝 .

= 𝜇1 tr(𝐄𝑒𝑜 )2 + 𝜇2 tr(𝐄2𝑒𝑜 )

𝜓𝑏 (𝐄𝑒1 )

= 𝑏1 tr(𝐄𝑒1 )2 + 𝑏2 tr(𝐄2𝑒1 )

𝜓𝛾 ( )

=

1 2

𝜇2

(D.5)

(D.6)

⋅ .

In the framework of mixed methods, we require the shear variable to be equal to the shear deformation 𝑑 , that is, we introduce the shear constraint = 𝑑 ; it follows that the shear stress must be considered as the reaction enforcing such a kinematical constraint. Then, we relax the energy 𝜓 by considering the work done by the shear stress on the shear constraint





 ⋅ ( − 𝑑 )) .

𝜓𝑟 = (𝜓𝑚 (𝐄𝑒𝑜 ) + 𝜓𝑏 (𝐄𝑒1 ) + 𝜓𝛾 ( ) −

ℎ3 𝜆 𝜇 , 12 𝜆 + 2 𝜇

𝑏2 =

ℎ3 𝜇. 12

(D.10)

(D.11)

As aforementioned, the final results of metric’s transport can be affected by the choice of the elastic moduli; in particular, for isotropic materials, it is important the ratio 𝑘∕𝜇, between the bulk and the shear modulus. In all our examples, we used 𝑘∕𝑚𝑢 ∼ 1. For the example in Section 5.1.2 we investigated the influence of the bulk modulus by solving the problem with 𝑘∕𝜇 = 10−2 , 10−1 , 1, 10, 102 , that is, with 𝑘∕𝜇 ranging 5 orders of magnitude. Fig. E.14 shows the configurations 𝑌 resulting from the metric transport discussed in 5.1.2 for 𝑘∕𝜇 = 1 (black, center), representing our benchmark, for a highly compressible material having 𝑘∕𝜇 = 10−2 (blue, left), and a highly incompressible material with 𝑘∕𝜇 = 102 (red, right). With reference to the benchmark at the center, it may be noticed an area overshoot (left) or a mesh inversion (right). This problem could be avoided by using solution strategies for high incompressible materials, such as the use of a weak constraint to enforce incompressibility; a discussion of this issue is beyond the scope of the paper. Fig. E.15 shows the average values of the stretching error 𝜌𝐴 versus 𝛼, for different values of the ratio 𝑘∕𝜇. We note that the curve 𝜌𝐴 = 𝜌𝐴 (𝛼) has an asymptotic curve as 𝑘∕𝜇 growths above 1. The down step at 𝛼 ∼ 0.9 for 𝑘∕𝜇 = 10, 100 indicates the occurrence of a mesh inversion; in such a case the solution in not valid any more.

where 𝐽𝑜 = det(𝐔𝑥 )∕ det(𝐔𝑥𝑜 ), and 𝐽𝑦̄ = det(𝐔𝑦𝑜 ̄ ). The three energy densities are defined as: 𝜓𝑚 (𝐄𝑒𝑜 )

(D.9)

Appendix E. The effect of the elastic moduli

The second fundamental form 𝐁𝑦 is then computed by using the gradient of 𝐝: 𝐁𝑦 = ∇𝐝𝑇 𝐅̂ 𝑦 . We assume the elastic energy density 𝜓, a density per unit reference area, to be the sum of three contributions, a membrane energy 𝜓𝑚 , a bending energy 𝜓𝑏 and a shear energy 𝜓𝛾 :

𝜇2 = ℎ 𝜇 .

(3 𝑘 − 2 𝜇) 𝜇 𝜆𝜇 = . 𝜆 + 2𝜇 3𝑘 + 4𝜇

(D.4)

𝜓 = (𝜓𝑚 (𝐄𝑒𝑜 ) + 𝜓𝑏 (𝐄𝑒1 ) + 𝜓𝛾 ( )) 𝐽𝑜 𝐽𝑦̄ ,

𝜆𝜇 , 𝜆 + 2𝜇

We note that introducing the bulk modulus 𝑘 = 𝜆 + 2∕3 𝜇, we can write the stiffness in 𝜇1 and 𝑏1 as a function of 𝜇, 𝑘

We proceed as follows; we define the normal 𝐝 to the actual surface in terms of the state variable , and the shear deformation 𝑑 in terms of 𝐝 as 𝐧𝑦̄ + 𝛽1 𝐞1 + 𝛽2 𝐞2

constraint

the bending stiffnesses 𝑏𝑖 are given by

(D.3)



(D.8)

The dimension of the shell-like region 𝛺 = 𝑍 × ℎ, with 𝑍 = 𝐿 × 𝐻, are 𝐻 = 1 m, 𝐿 = 𝐻, or 𝐿 = 10 𝐻, and ℎ = 𝐻∕100. We use the following values for the Lamé parameters: 𝜆 = 1.5 N/m2 , and 𝜇 = 1 N/m2 ; the membrane stiffnesses 𝜇𝑖 are given by

then, the target shell 𝑓𝑦 = 𝑓̂𝑦 + 𝜁 𝐧𝑦 is described by displacing 𝑓̂𝑦̄ along 𝐞1 , 𝐞2 , 𝐧𝑦̄ with

𝐝=

𝜕𝜓𝛾 ( )

⏟⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏟

shear stress.

,

bending stress

− ) ⋅ ̃ +  ⋅ ̃ − ̃ ⋅ ( − ) . 𝑑 𝑑 𝜕 ⏟⏞⏞⏟⏞⏞⏟

membrane stress

rotation parameters;

Given the source shell 𝑓𝑦̄ = 𝑓̂𝑦̄ + 𝜁 𝐧𝑦̄ , using (40) we define the unit tangent vectors 𝐞𝛼 and the normal vector 𝐧𝑦̄ as 𝐞1 =

𝜕𝜓𝑚 (𝐄𝑒𝑜 ) ̃ 𝜕𝜓𝑏 (𝐊) ̃ ⋅𝐄𝑒𝑜 + ⋅𝐊 𝜕𝐄𝑒𝑜 𝜕𝐊 ⏟⏞⏟⏞⏟ ⏟⏞⏞⏟⏞⏞⏟

=



(D.7) 10

L. Teresi, F. Milicchio, S. Gabriele et al.

International Journal of Non-Linear Mechanics 119 (2020) 103326 [5] R. Sumner, J. Popović, Deformation transfer for triangle meshes, ACM Trans. Graph. 23 (3) (2004) 399–405. [6] V. Varano, S. Gabriele, L. Teresi, I. Dryden, P. Puddu, C. Torromeo, P. Piras, The tps direct transport: A new method for transporting deformations in the size-and-shape space, Int. J. Comput. Vis. 124 (3) (2017) 384–408. [7] Y. Klein, E. Efrati, E. Sharon, Shaping of elastic sheets by prescription of non-euclidean metrics, Science 315 (5815) (2007) 1116–1120. [8] P. Nardinocchi, L. Teresi, V. Varano, The elastic metric: A review of elasticity with large distortions, Int. J. Non-Linear Mech. 56 (2013) 34–42. [9] P.G. Ciarlet, An Introduction to Differential Geometry with Applications to Elasticity, Springer, 2005. [10] Q. Xie, S. Kurtek, H. Le, A. Srivastava, Parallel transport of deformations in shape space of elastic surfaces, in: Proceedings of the IEEE International Conference on Computer Vision, 2013, pp. 865–872. [11] C. Twining, S. Marsland, C. Taylor, Metrics, connections, and correspondence: the setting for groupwise shape analysis, in: International Workshop on Energy Minimization Methods, 2011, pp. 399–412. [12] A. Trouvé, L. Younes, Metamorphoses through lie group action, Found. Comput. Math. 5 (2) (2005) 173–198. [13] V. Varano, S. Gabriele, L. Teresi, I. Dryden, P. Puddu, C. Torromeo, P. Piras, Comparing shape trajectories of biological soft tissues in the size-andshape space, in: Biomat 2014 International Symposium on Mathematical and Computational Biology, 2015, 2014, pp. 351–365. [14] F. Milicchio, V. Varano, S. Gabriele, L. Teresi, P. Puddu, P. Piras, Parallel transport of local strains, Comput. Methods Biomech. Biomed. Eng.: Imaging & Visualization (2018) 1–9. [15] A. Dicarlo, S. Quiligotti, Growth and balance, Mech. Res. Commun. 29 (2002) 449–456. [16] E. Kröner, Allgemeine kontinuumstheorie der versetzungen und eigenspannungen, Arch. Ration. Mech. Anal. 4 (1959). [17] E. Lee, Elastic–plastic deformation at finite strains, J. Appl. Mech. 36 (1969). [18] C. Davini, Some remarks on the continuum theory of defects in solids, Int. J. Solids Struct. 38 (2001) 1169–1182. [19] W. Koiter, On the nonlinear theory of thin elastic shells, in: Proc. K. Ned. Akad. Wet. B, Vol. 69, 1966, pp. 1–54. [20] I.N. Sneddon, R. Hill, P.M. Naghdi, H. Ziegler, Progress in Solid Mechanics, Vol. 4, North-Holland Publishing Company, John Wiley & Sons, Amsterdam, New York, 1963.

Fig. E.15. Rectangle-like example. Average value of 𝜌𝐴 versus 𝛼, for different values of the ratio 𝑘∕𝜇. The curve 𝜌𝐴 = 𝜌𝐴 (𝛼) has an asymptotic curve as 𝑘∕𝜇 growths above 1. The down step at 𝛼 ∼ 0.9 for 𝑘∕𝜇 = 10, 100 indicates the occurrence of a mesh inversion.

References [1] A. Chern, F. Knöppel, U. Pinkall, P. Schröder, Shape from metric, ACM Trans. Graph. 37 (4) (2018) 63:1–63:17. [2] D. Boscaini, D. Eynard, D. Kourounis, M.M. Bronstein, Shape-from-operator: recovering shapes from intrinsic operators, Comput. Graph. Forum 34 (2) (2015) 265–274. [3] M. Lembo, On the determination of deformation from strain, Meccanica 52 (9) (2017) 2111–2125. [4] W. Pietraszkiewicz, M. Szwabowicz, C. Vallée, Determination of the midsurface of a deformed shell from prescribed surface strains and bendings via the polar decomposition, Int. J. Non-Linear Mech. 43 (7) (2008) 579–587.

11