domain: phys spi

text:
C Tannous

Tensor calculus and anisotropic elasticity for the impatient

Keywords: numbers: 46.25.-y, 46.35.+z, 62.20.D-Keywords: Elasticity in continuum mechanics of solids, mechanical properties of solids Contents I. Introduction Covariant derivative C. Voigt

Tensor calculus is introduced to Physics and Mechanical engineering students in 2D and 3D and applied to anisotropic elasticity such as in condensed matter physics approaching the subject from the practical tool aspect point of view. It provides powerful mathematical techniques to tackle many aspects of Vector Calculus, Continuum mechanics, Solid State Physics, Electromagnetism...

I. INTRODUCTION

Tensor calculus is usually introduced at the graduate/undergraduate level in Special Relativity or General Relativity/Gravitation/Cosmology courses. It has been used by Einstein to tackle the daunting mathematical operations needed in performing many required differential geometry tasks in curved space.

Prior to relativity courses, physics students are typically exposed to Vector calculus [START_REF] Spiegel | Vector analysis and an Introduction to Tensor Analysis, Schaum's outline series in mathematics[END_REF] in Classical Mechanics, Electromagnetism or in Mathematics for Physicists since they employ, in these courses, vector differential operators like the gradient, divergence, curl, Laplacian...

Tensor calculus provides the students a powerful geometrical extension of Vector calculus to arbitrary dimension and to mathematical objects they were not familiar with to tackle problems in fluids, elasticity, continuum mechanics... Another field of application of Tensor calculus [START_REF] Kay | Tensor calculus, Schaum's outline series in Mathematics[END_REF] is Solid State Physics where students become acquainted with geometry of anisotropic materials and recently with mathematics of Topological materials.

It is natural to suppose that a quantitative description of physical phenomena cannot depend on the coordinate system in which they are expressed. The argument may be returned: since physical phenomena are independent of the coordinate system, what are the possible implications on the nature of the quantities describing these phenomena?

Einstein covariance principle (ECP) states that Laws of Physics are same in every Lorentz reference (observer) frame, implying a geometric view of physical laws and calling for the mathematical formulation of any change in a physical quantity under a change of reference frame [START_REF]When Einstein worked at the Federal Office for Intellectual Property[END_REF].

The study of implications and the resulting classification of physical quantities constitutes the theory of tensor calculus which was used by Einstein. In this work we introduce some calculation rules used in tensor manipulation and show how the tensors are transformed by change of axis. We apply these results to the elasticity case.

This work is not only useful to Physics and Mechanical engineering students but also to those endeavoring in Fluid and Continuum Mechanics. It is organized as follows: In section 2, we introduce tensor fields and coordinate systems in compliance with ECP and in Section 3 Tensor methods for treating Linear Algebra and Vector Calculus [START_REF] Spiegel | Vector analysis and an Introduction to Tensor Analysis, Schaum's outline series in mathematics[END_REF][START_REF] Kay | Tensor calculus, Schaum's outline series in Mathematics[END_REF] is described. We treat three distinct cases: Cartesian, orthogonal curvilinear and non-orthogonal curvilinear coordinates. Section 4 is dedicated to Tensor extension and contraction of physical quantities and Section 5 is dedicated to Elasticity problems as an application. Section 6 contains a conclusion and perspectives.

Appendix A describes Tensor etymology and its evolution, Appendix B is about covariant derivative needed in curvilinear coordinates and Appendix C discusses Voigt notation heavily used in Elasticity.

II. PHYSICAL FIELDS AND TENSOR FIELDS IN GENERALIZED COORDINATES

Physical fields are function of position in space thus it is imperative to introduce coordinate systems that will allow us to express accurately these fields. A field is a quantity that takes values "instantaneously" at some arbitrary position in space. It is obvious this can't be true according to Feynman and Special Relativity since we know that any change a physical quantity undergoes a signal is needed to transmit the information. As an example, if we have a charge q at a point considered in R 3 as the origin r = 0 and if we want know the value of the electric field at any point r = 0, it is necessary that a photon (called "longitudinal") is transmitted from the origin to the point r. Since the speed of the photon and of light are finite, a certain time for this information to travel.

Let us contemplate some fields at some point r with different nature such as scalar, vector, tensorial.... Temperature field in a room is a scalar field T (r) while the speed within a fluid is a vector field V (r). When we consider the electromagnetic field, we have two vector quantities E(r) and H(r) knowing that the electric field E(r) is intrinsically related to the magnetic field H(r) via Maxwell's equations. One relevant question that we answer later in subsection IV A is: since E(r) H(r) fields are intimately related, is it possible to encapsulate them compactly in Maxwell's equations?

A. 2D covariant and contravariant coordinates ECP being our guide, we start with a general description of coordinates before describing reference (observer) frame transformation. This word is distinct from the covariant adjective that means "form invariant" when it concerns physical laws (pertinent to ECP) or "transforms as standard basis vectors" when it concerns basis vectors. We use this notion in order to describe the different ways of expressing a 2D position vector in a given set of basis vectors e 1 , e 2 (cf fig. 1). Covariance principle has been debated extensively and several workers questioned its validity and pondered about its consequences (see for instance references [START_REF] Misner | Gravitation[END_REF][START_REF] Norton | [END_REF]).

We assume that e 1 , e 2 are normalized (||e 1 || = ||e 2 || = 1) and making an angle θ (cf fig. 1).

Let us express the components of some arbitrary 2D vector r : r = x 1 e 1 + x 2 e 2 where x i are the ordinary components (contravariant). Geometrically, the covariant components x i are obtained by perpendicular projection on to the standard basis axes (cf fig. 1). Thus x i = r • e i , i = 1, 2. Substituting r in the latter yields:

x i = (x 1 e 1 + x 2 e 2 ) • e i , i.e: x 1 = (x 1 e 1 + x 2 e 2 ) • e 1 = x 1 + x 2 e 2 • e 1 = x 1 + x 2 cos θ, x 2 = (x 1 e 1 + x 2 e 2 ) • e 2 = x 2 + x 1 e 2 • e 1 = x 1 cos θ + x 2
Assuming a symmetry principle between covariance and contravariance, let us rewrite r = i x i e i (cf fig. 1).

Taking the scalar product with e j yields: e j •r = i x i e j •e i . Assuming the covariant coordinate x j is the projection of r over e j we get: x j = e j • r = i x i e j • e i then the scalar product e j • e i = δ i j the Kronecker symbol written in a mixed fashion but meaning same as generally known i.e. δ i j = 0 if i = j and δ i j = 1 if i = j. Thus, the contravariant basis is orthogonal to the covariant and is termed dual basis in Mathematics or Reciprocal basis in Crystallography [START_REF] Kittel | Introduction to Solid State Physics[END_REF].

We can express this mathematically by introducing a e 3 unit vector perpendicular to the e 1 , e 2 plane. An explicit expression of the contravariant basis vectors can be derived by exploiting the reciprocal of the standard basis formula [START_REF] Kittel | Introduction to Solid State Physics[END_REF]:

e 1 =
e 2 × e 3 (e 1 , e 2 , e 3 )

, e 2 = e 3 × e 1 (e 1 , e 2 , e 3 )

(1)

The triple product (e 1 , e 2 , e 3 ) = sin θ confirming the geometrical construct in fig. 1 that e 1 is the projection of e 1 and e 2 is the projection of e 2 with the values of the norms ||e 1 || = ||e 2 || = 1/ sin θ.

Let us rewrite the position vector as r = x i e i with Einstein summation rule stating that whenever we have an index appearing twice as covariant and contravariant, a summation over it is implied.

Using laws relating covariant and contravariant components, we define a "metric matrix" g ij in the following way

g ij = 1 cos θ cos θ 1 , i = 1, 2.
g ij allows us to go between covariant and contravariant coordinates since:

x i = g ij x j .
g ij makes the indices "move" down implying by symmetry the existence of another contravariant metric matrix g ij that makes the indices "move" up. Note that g ij = e i • e j (respectively g ij = e i • e j ) relates two vectors e i and e j (resp. e i and e j ) and if the standard basis e i is orthonormal (θ = π/2), g ij = δ ij the 2D unit matrix and x i = x i making covariant and contravariant coordinates same. Evaluating g ij from the scalar products of e i , e j using the previous definitions 1 and fig

1 yields: [g ij ] = 1 sin 2 θ 1 -cos θ -cos θ 1 , i = 1, 2.
[g ij ] and [g ij ] are inverse of one another as consequence of their respective definitions. We provide an alternate proof based on their determinants: det[g ij ] = sin 2 θ = (e 1 × e 2 ) 2 and det[g ij ] = (e 1 × e 2 ) 2 moreover (e 1 × e 2 ) e 3 and (e 1 × e 2 ) e 3 . Thus (e 1 × e 2 ) 2 (e

1 × e 2 ) 2 = [(e 1 × e 2 ) • (e 1 × e 2 )] 2 yielding det[g ij ]det[g ij ] = [(e 1 × e 2 ) • (e 1 × e 2 )] 2 .
Using the dot product of two cross products [START_REF] Arfken | Mathematical methods for physicists[END_REF]:

(a × b) • (c × d) = a • c a • d b • c b • d (2)
we get:

(e 1 × e 2 ) • (e 1 × e 2 ) = e 1 • e 1 e 1 • e 2 e 2 • e 1 e 2 • e 2 = 1 0 0 1 (3)
proving once again that [g ij ] and [g ij ] are inverse of one another. We prove further below that [g ij ] and [g ij ] are rank-2 tensors allowing to write compactly index lowering x i = g ij x j and lifting x i = g ij x j operations. 2: 3D covariant and contravariant curvilinear coordinates: contravariant vectors are perpendicular to u i (x, y, z) = constant surfaces whereas covariant basis vectors ei are tangent to u i curves that are intersections of u i surfaces. Position vector r = (x, y, z) can be expressed with u i coordinates as r(u 1 , u 2 , u 3 ).

3D Cartesian orthonormal coordinates (x, y, z) are not appropriate to use in physical problems with special symmetry. It is beneficial to move to curvilinear coordinates {x, y, z} → {u 1 , u 2 , u 3 } defined by extending the previous 2D case to contravariant basis vectors e i perpendicular to u i (x, y, z) = constant surface as displayed in fig. 2 whereas covariant basis vectors e i are tangent to u i curves lying at the intersection between respective u i surfaces.

Defining e i = ∂r ∂u i and e i = ∇u i = ( ∂u i ∂x , ∂u i ∂y , ∂u i ∂z ), where the operator ∇ = ( ∂ ∂x , ∂ ∂y , ∂ ∂z ), it is possible to show that their properties are similar to the 2D case and that many results obtained previously can be extended to the 3D case.

Extending the previous 2D case to 3D, we write explicit expressions of the contravariant basis vectors versus the covariant one {e 1 , e 2 , e 3 } as [START_REF] Kittel | Introduction to Solid State Physics[END_REF]:

e 1 = e 2 × e 3 (e 1 , e 2 , e 3 )
, e 2 = e 3 × e 1 (e 1 , e 2 , e 3 )

,

e 3 = e 1 × e 2 (e 1 , e 2 , e 3 ) (4)
The orthogonality relations e i • e j = 0, i = j can be proven directly from the line-element dr obtained from Taylor expanding the position vector r(u 1 , u 2 , u 3 ) in the u i basis as:

dr = ∂r ∂u i du i (5)
Taking the dot product with e j = ∇u j yields ∇u j • dr = ∇u j • ∂r ∂u i du i and using the gradient definition ∇u j • dr ≡ du j , we conclude e i • e j = 0, i = j.

Similarly to the 2D case, a covariant "metric matrix" g ij can be defined from: g ij = e i • e j , i, j = 1, 2, 3 as well as a contravariant one g ij = e i • e j , i, j = 1, 2, 3 with the properties of raising/lowering indices x i = g ij x j , x i = g ij x j and mutual orthogonality which translates into [g ij ][g ij ] = 1 1 and indicially as: g ik g kj = δ j i . Note that ds 2 = dx 2 + dy 2 + dz 2 in a Cartesian orthonormal frame {x, y, z} with g ij = g ij = δ ij . In a frame characterized by cylindrical coordinates {ρ, φ, z} the line-element squared is given by: ds 2 = dρ 2 + ρ 2 dφ 2 + dz 2 . Using the correspondence {1, 2, 3} → {ρ, φ, z}, the metric matrix components are:

g 11 = 1, g 22 = ρ 2 , g 33 = 1, g ij = 0 if i = j, i, j = 1, 2, 3.
In spherical coordinates {r, θ, φ} the line-element squared is given by: ds 2 = dr 2 + r 2 dθ 2 + r 2 sin 2 θdφ 2 . Using the correspondence {r, θ, φ} → {1, 2, 3}, the metric matrix components are:

g 11 = 1, g 22 = r 2 , g 33 = r 2 sin 2 θ, g ij = 0 if i = j, i, j = 1, 2, 3.
In general (non-orthogonal) curvilinear coordinates, the line-element squared is given by ds 2 = dr•dr = g ij du i du j = g ij du i du j yielding non-diagonal "metric matrices".

C. Einstein covariance principle and tensor transformation

Another aspect of ECP is that a physical quantity when it undergoes a change of reference frame transforms in a fashion revealing its true nature. Thus, it suffices to study the mathematics of reference frame change in order to understand the underlying nature of any physical quantity.

Using the classical (contravariant) coordinates of a vector x i in a reference frame R, we call its transform xi in the new reference R. Note that x i can be considered as a function of xj the new coordinates and mathematically express it as

x i = x i (x j ).
Performing a Taylor expansion in n dimensions allow us to find a small displacement as:

dx i = n j=1 ∂x i ∂ xj dx j (6)
It is expressed with Einstein summation rule as:

dx i = ∂x i ∂ xj dx j (7)
This shows that we have a "Jacobian" factor of the style ∂x i ∂ xj which expresses the laws of transformation. If we start from a scalar field (such as temperature or electric scalar potential), we obviously have:

φ(x i ) = φ(x j ) (8)
meaning that for a scalar field, there is no "Jacobian" factor while that for a vector field like x i , we have only one factor. Generalizing, we might say a "Jacobian" factor indicates the rank of the tensor.

Let us take the derivatives of the two members of the equation 8 using the compound derivative since we can consider x i = x i (x j ) just like xj = xj (x i ):

∂φ ∂x i = ∂ φ ∂ xj ∂ xj ∂x i (9)
We infer an "inverse Jacobian" factor ∂ xj ∂x i appears in the gradient transformation. We conclude that there are two possibilities during a coordinate transformation: a "direct Jacobian" ∂x i ∂ xj is present in the contravariant case whereas an "inverse Jacobian" ∂ xj ∂x i appears in the covariant case. Thus

∂x i ∂ xj ∂ xj ∂x k = δ i k .
For a scalar field (rank-0) there is no Jacobian factor, whereas for a vector field (rank-1) there is a single factor and so on.

It is important to note that for rank-2 tensors, we have a matrix representation by letting its coefficients depend on position r. However one should recall that a matrix and a tensor are different mathematical objects.

We summarize below the various transformations of a rank-1 and rank-2 tensors:

• Transformation law for a contravariant field:

A i = ∂x i ∂ xj Āj ( 10
)
• Transformation law for a covariant field:

A i = ∂ xj ∂x i Āj (11)
• Transformation law for a twice contravariant tensor:

A ij = ∂x i ∂ xk ∂x j ∂ xl Ākl (12)
• Transformation law for a twice covariant tensor:

A ij = ∂ xk ∂x i ∂ xl ∂x j Ākl (13)
• Transformation law for a mixed singly covariant, singly contravariant tensor:

A k i = ∂ xm ∂x i ∂x k ∂ xl Āl m (14)

D. Tensor construction

After reckoning the 2D metric matrix is built from two basis vectors, a mechanical question arises concerning the pressure involving also two vectors (force and area). Thus, what is the underlying nature of the pressure field?

In fact, pressure is a tensor field of rank-2 since it connects two vector fields as described in the elasticity section V. Essentially, there are several methods to build a tensor:

• From association of several vector (rank-1) components:

The "metric tensor" we called previously "metric matrix" is obtained from a scalar product operation between two vectors (rank-1 tensors): g ij = e i • e j or g ij = e i • e j , which are both rank-2 tensors. Since ds 2 the squared line-element is a scalar, it can be expressed invariantly in two different coordinate systems {x i } → {x i } where the set {x i } is orthonormal and {x i } is not. ds 2 = δ ij dx i dx j = g ij dx i dx j . Performing a change of reference frame gives:

dx i = ∂x i ∂ xj dx j . Thus ds 2 = δ ij dx i dx j = δ ij ∂x i ∂ xk dx k ∂x j ∂ xl dx l = g kl dx k dx l .
Covariant metric tensor components are given by:

g kl = δ ij ∂x i ∂ xk ∂x j
∂ xl where coordinates xi are e.g. curvilinear. The contravariant metric tensor g ij = e i • e j is the inverse of the covariant g ij as proven previously in subsection II B. Another proof is that e i being the dual (contravariant, reciprocal) basis, the two "matrices" are inverse of one another [START_REF] Kittel | Introduction to Solid State Physics[END_REF].

• From direct multiplication of vector components (outer product):

As an example, take two vector fields A i and B j and make the component product by component to obtain a tensor of rank-2: T ij = A i B j which can be formally written like T = A ⊗ B. This is of course generalizable to any dimension. For example with three vector fields A, B, C on can do: 3). This implies that a tensor of rank-r is constructed from vector fields in D dimensions, we have D r components.

T = A ⊗ B ⊗ C so T ijk = A i B j C k (cf fig.
• From spatial derivation [9] of a scalar field:

The gradient of a scalar field φ(r) is in fact a means to transform a rank-0 field into a vector field (rank-1):

[gradφ] i = ∂φ ∂x i . Note however that the gradient behaves as a covariant vector field (It is shown further below that there are contravariant and covariant varieties of the gradient). As a consequence, derivation (respectively integration) of another tensor field allow to change the rank by 1 (resp. -1).

• From contraction of another tensor:

Contraction is an operation consisting of equating two indices (say i = j) in a tensorial expression and summing over (say i). One example is the scalar product of two covariant vectors: a • b = a i b i that transforms a rank-2 tensor a i b j after a contraction into a scalar (rank-0). Another example drawn from linear algebra is to apply a matrix M with coefficients M ij to a vector a = (a 1 , a 2 , a 3 ) resulting into another vector b such that b = M a, meaning indicially b i = 3 j=1 M ij a j . Using Einstein summation rule, we write b i = M ij a j which in fact might be viewed as a contraction (j = k) from a rank-3 tensor M ij a k to a (rank-1) vector b i . Rephrasing in Tensor language writes: Suppose T i1,i2,i3,,im transforms as a tensor of rank-m then T j,j,i3,,im transforms as a tensor of rank-(m -2). This operation corresponds to a general contraction of a tensor.

• By combining all previous operations In Cartesian (Affine) coordinates, the basis vectors are same independent of position. Moreover we do not make any distinction, in this paragraph between contravariant and covariant components since basis vectors are also orthonormal e i • e j = δ ij . This implies the Kroenecker symbol [START_REF]There are several possibilities for the rank-2 Kroenecker tensor: a-Doubly covariant: δij = ei • ej, b-Doubly contravariant: δ ij = e i • e j and c-Mixed covariant-contravariant: δ j i = ei • e j assuming the right properties of orthonormality of the bases ei and e j[END_REF] is in fact a rank-2 tensor built from the scalar product of two vectors.

T ijk = A i B j C k , i, j, k = 1, 2,
Let us introduce rank-3 Levi-Civita antisymmetric tensor: ijk = e i • (e j × e k ) = (e i , e j , e k ), i, j, k = 1..3 with the following properties: ijk = 0 for i = j, j = k, i = k whereas ijk = ±1 depending on whether i, j, k is an even or odd permutation of 1,2,3.

Levi-Civita ijk is in fact a pseudo-tensor [START_REF]as: ijk = g il gjmg kn lmn[END_REF], since built from a triple-product of vectors (rank-1 tensors). The product of two Levi-Civita tensors is found from a Linear Algebra theorem stating that the determinant of the product of two triple products is the determinant of scalar products: ijk lmn = (e i , e j , e k )(e l , e m , e n ).

Since a triple product is a determinant we use det(A)det(B) = det(AB), thus:

ijk klm = δ il δ im δ jl δ jm = δ il δ jm -δ jl δ im
Using the former identity yields ijk ijl = 2δ kl . Contracting over all indices yields: ijk ijk = 2δ kk = 6

Let us apply to the scalar product, vector product and double vector product.

• Scalar product: Given two vectors A = A i e i and B = B i e i , we write:

A • B = A i e i • B j e j = A i B j e i • e j = A i B j δ ij = A i B i • Vector product: A × B = A i e i × B j e j = A i B j e i × e j .
Taking the scalar product with a basis vector:

(A × B) • e k = A i B j (e i × e j ) • e k = (e i , e j , e k )A i B j = ijk A i B j • Double vector product: A × (B × C). Taking the i-th component we get: ijk A j (B × C) k = ijk A j klm B l C m This yields: (δ il δ jm -δ jl δ im )A j B l C m = A j B i C j -A j B j C i = (A • C)B i -(A • B)C i Thus: A × (B × C) = (A • C)B -(A • B)C B.

Differential vector calculus with tensors in Cartesian coordinates

Differential vector calculus is based on the use of the nabla vector operator ∇ yielding the gradient, divergence, curl and Laplacian. In (Affine) Cartesian coordinates, we do not make any distinction between contravariant and covariant components (A i = A i ) since we have e i • e j = δ ij . Moreover, the e i are independent of coordinates giving:

∇ = e i ∂ ∂xi = e i ∂ i = ∂ i e i .
• Scalar field gradient:

Given a scalar field φ(r), the i-th component of the gradient is [gradφ(r)] i = ∇ i φ(r) = ∂φ(r) ∂xi = ∂ i φ. • Vector field gradient:

For a vector field A(r), the i-th component of the gradient of j-th

component of A is [∇(A) j ] i = ∇ i A j = ∂ i A j
• Vector field divergence: Given a vector field A(r), its divergence is:

∇ • A(r) = e i ∂ i • A j e j = ∂ i A j e i • e j = ∂ i A j δ ij = ∂ i A i .
• Vector field curl: Given a vector field A(r), its curl is:

∇ × A(r) = e i ∂ i × A j e j = ∂ i A j (e i × e j ). The i-th component is [curlA] i = ijk ∂ j A k .
• Scalar field Laplacian:

The Laplacian of a scalar field φ(r

) is ∆φ(r) = ∇ • ∇φ(r) = ∂ i ∂ i φ. • Curl divergence: div curlA(r) = ∂ i ijk ∂ j A k = ijk ∂ i ∂ j A k = 0.
This stems from the fact a contraction of two tensors with one symmetric (∂ i ∂ j ) [12] and the other antisymmetric ijk gives zero.

• Gradient curl: [curl gradφ(r)] i = ijk ∂ j ∂ k φ = 0.
This originates as before from the contraction of a symmetric (∂ j ∂ k ) and antisymmetric ijk tensors.

• curl of curl:

[curl curlA(r)] i = ijk ∂ j [curlA(r)] k = ijk ∂ j klm ∂ l A m = ijk klm ∂ j ∂ l A m = (δ il δ jm -δ jl δ im )∂ j ∂ l A m = ∂ i ∂ m A m -∂ j ∂ j A i = ∂ i (divA) -[∆A] i (15)
We end up with: curl curlA = grad(divA) -∆A

• Vector field Laplacian: In case of a vector field A(r) we apply the Laplacian operator on every component: [∆A] i = ∂ j ∂ j A i . Note that this is not true for curvilinear coordinates case as seen further below.

C. Differential operators in orthogonal curvilinear coordinates

The gradient of a scalar field φ(r(u i )) is evaluated by comparing its differential dφ = ∂φ ∂u i du i to the corresponding gradient expression dφ = gradφ • dr. The line-element expression: dr = h i e i du i uses h i coefficients that are scale factors [START_REF]Scale factors (also called Lamé coefficients) are introduced to make the basis vectors ei normalized. In orthogonal curvilinear coordinates we obtain from expansion eq. 5 a scale factor hi = | ∂r ∂u i | yielding ∂r ∂u i = hiei such that: dr = hieidu i (summation over i being assumed despite its presence three times) with ei normalized. The scale factors can be defined with the metric tensor hi = √ ei • ei = √ gii ([END_REF] rendering e i normalized. An important difference between this case and the Cartesian is that e i depend on coordinates u j . Moreover any vector is expressed with the scale factors as:

A = h i A i e i .
The ∇ operator in orthogonal curvilinear coordinates is defined by: ∇ = e i ∂ hi∂u i . Note that we keep using the contravariant notation (such as A i , u i ...) despite the fact they are same in orthogonal curvilinear coordinates. The expressions for the gradient, div, curl and Laplacian are obtained from the ∇ operator as done in Cartesian coordinates. Thus:

gradφ = e i ∂ h i ∂u i φ div A = e i ∂ h i ∂u i • A = e i ∂ h i ∂u i • (h j A j e j ) = 1 h 1 h 2 h 3 ∂ ∂u i (h j h k A i ) , i, j, k = 1, 2, 3 curl A = e i ∂ h i ∂u i × A = e i ∂ h i ∂u i × (h j A j e j ) = 1 h 1 h 2 h 3 h 1 e 1 h 2 e 2 h 3 e 3 ∂ ∂u 1 ∂ ∂u 2 ∂ ∂u 3 h 1 A 1 h 2 A 2 h 3 A 3 ∆φ = ∇ • ∇φ = e i ∂ h i ∂u i • e j ∂ h j ∂u j φ = 1 h 1 h 2 h 3 ∂ ∂u i h j h k h i ∂φ ∂u i , i, j, k = 1, 2, 3
In order to evaluate the Laplacian of a vector field [START_REF] Moon | [END_REF]15] (called the Beltrami operator) in coordinates other than Cartesian, one resorts to the double curl formula derived above: ∆A = ∇(∇ • A) -∇ × (∇ × A).

D. Differential operators in non-orthogonal curvilinear coordinates

In non-orthogonal curvilinear coordinates, one distinguishes covariant from contravariant quantities as well as dependence of e i , e i on u j , u j . The line-element squared versus metric tensor is expressed as:

ds 2 = dr • dr = ∂r ∂u i du i • ∂r ∂u j du j ≡ g ij du i du j (16)
The ∇ operator in curvilinear coordinates is written as: ∇ = e i ∂ ∂u i without using any scale factors [START_REF]Scale factors (also called Lamé coefficients) are introduced to make the basis vectors ei normalized. In orthogonal curvilinear coordinates we obtain from expansion eq. 5 a scale factor hi = | ∂r ∂u i | yielding ∂r ∂u i = hiei such that: dr = hieidu i (summation over i being assumed despite its presence three times) with ei normalized. The scale factors can be defined with the metric tensor hi = √ ei • ei = √ gii ([END_REF] and consequently possessing un-normalized e i , e i .

The gradient [START_REF]∂φ ∂u i e i with all e i normalized[END_REF] of a scalar field φ(r) is ∇φ = e i ∂φ ∂u i corresponding to the same expression as in the orthogonal case (with the scale factors h i = 1). This is not true for the divergence, curl and Laplacian where the covariant derivative intervenes because of the appearance of the Jacobian terms (cf Appendix B).

The divergence of a vector field A is written as:

∇ • A = e j ∂ ∂u j • (A i e i ) = e j • [( ∂A i ∂u j e i ) + A i Γ k ij e k ]. The term [( ∂A i ∂u j e i ) + A i Γ k ij e k ] is transformed into [( ∂A i ∂u j ) + A k Γ i kj ]e i (cf Appendix B) containing A covariant derivative D j A i = [( ∂A i ∂u j ) + A k Γ i kj ]
. Thus the divergence of a vector field A is e j • D j A i e i is similar to the Cartesian case with the covariant derivative replacing the ordinary one. This can be further expressed compactly using Γ i ij = 1 2 g lk ∂g lk ∂u j (cf Appendix B) as:

1 √ g ∂ ∂u k ( √ gA k ) where g = det[g ij ].
Similarly, the curl of a vector field A is defined by: ∇ × A = e i ∂ ∂u i × A j e j = ( ∂A k ∂u i + A i Γ k ij )(e k × e j ). Again it is similar to the Cartesian case, the covariant derivative replacing the ordinary one.

The scalar Laplacian is obtained from a double application of the covariant derivative following the Cartesian case: ∆φ = ∂ i ∂ i φ. Starting from ∆φ = D i D i φ, we introduce the metric tensor [START_REF] Arfken | Mathematical methods for physicists[END_REF] such that:

∆φ = D i (g ij D j φ) = [ ∂(g ij Dj φ) ∂u i + (g ij D j φ)Γ k ij ] = 1 √ g ∂ ∂u j ( √ gg jk ∂φ ∂u k ).
In order to evaluate the vector Laplacian (Beltrami operator) we follow Hirota et al. [15] work who have explicitly shown that the operators ∆ = ∇ • ∇ and ∇(∇• ) -∇ × (∇× ) are same, since the equivalence taken for granted from the Cartesian case is not acceptable. Thus we have:

∆A = ∇ • ∇A = e i ∂ ∂u i • e j ∂ ∂u j A = g ij ∂ 2 ∂u i ∂u j -Γ k ij ∂ ∂u k A.

IV. TENSOR EXTENSION AND REDUCTION OF PHYSICAL QUANTITIES A. Tensor extension of physical quantities

Many possibilities exist to extend the attributes [START_REF] Gallego | [END_REF][START_REF] Newnham | Properties of Materials -Anisotropy, Symmetry, Structure[END_REF][START_REF] Nye | Physical properties of crystals and their representation by tensors and matrices[END_REF] of any physical quantity initially attached to isotropic media (non-crystalline) to some corresponding anisotropic medium (crystalline).

For example, in electrostatics of isotropic media, the relation:

D = E with D E when extended to crystals becomes D i = ij E j with D ∦ E.
Similarly, in magnetostatics:

B = µH translates to B i = µ ij H j .
Electric conduction in materials is described by Ohm's law J = σE that translates into J i = σ ij E j in the anisotropic case. The resistivity tensor ρ ij inverse of the conductivity tensor σ ij allows us to rewrite Ohm's law in the following form:

E i = ρ ij J j .
In magnetic (anisotropic) crystals, the susceptibility defined by M i = χ ij H j where M is the magnetization induced by a magnetic field H is extended from the simple case where M = χH with M H and χ a scalar.

In metallic, semiconducting and insulating crystals, anisotropic effective mass that originated from ordinary mass of Newton's law F = mγ would be translated into F i = m ij γ j where γ is acceleration.

In piezoelectric crystals [START_REF] Landau | Electrodynamics of Continuous Media[END_REF][START_REF] Cady | Piezoelectricity: An introduction to the theory and applications of electromechanical phenomena in crystals[END_REF][START_REF] Kao | Dielectric phenomena in solids[END_REF], a stress σ creates a polarization P such that P = dσ with d the piezoelectric modulus. In the anisotropic case, a stress σ jk creates a polarization P i such that P i = d ijk σ jk with d ijk a rank-3 tensor of piezoelectric moduli since it relates a rank-2 (stress) and a rank-1 (polarization) tensors.

When a magnetic field H is applied to a metal or a semiconductor traversed by an electric current density J , the Hall effect arises as a voltage, transverse to J , described by an electric field E i = ρ ijk J j H k where ρ ijk is a rank-3 tensor relating three vectors E, J and H.

It is tempting to generalize Ohm and Hall results by introducing [START_REF] Gallego | [END_REF] an ordered field expansion:

E i = ρ (0) ij J j + ρ (1) ijk J j H k + ρ (2) ijkl J j H k H l ... where ρ (0) ij is zero order (no field present) resistivity (Ohm), ρ (1)
ijk is first-order (field present) Hall effect, ρ

ijkl is second-order (field present) Hall effect... Note that only ρ (0) ij has dimension of resistivity whereas higher order tensors ρ (n) ijk... , n ≥ 1 do not, because of the presence of magnetic field powers. Maxwell equations can be rewritten in a tensorial covariant fashion by encapsulating the scalar Φ and vector A potentials into a single 4-potential A α = (Φ, A) where α = 0, 1, 2, 3. Time corresponds to α = 0 whereas 3D spatial degrees of freedom are represented by α = 1, 2, 3. This leads to write the charge continuity equation ∂ρ ∂t + ∇ • J = 0 as ∂ α J α = 0 with J α = (cρ, J ) where c is light velocity.

The antisymmetric rank-2 electromagnetic tensor given by F α,β = ∂ α A β -∂ β A α allows to encapsulate the inhomogeneous Maxwell equations:

∇ • E = 4πρ, ∇ × B -1 c ∂E ∂t = 4π
c J into a covariant form given [START_REF] Jackson | Classical Electromagnetism[END_REF] by: ∂ α F α,β = 4π c J β .

B. Reduction to a scalar by projection along a single or two orthogonal directions

A physical property written by a rank-2 tensor like ↔ , for example, associated with the dielectric constant appearing in D i = ij E j , can be estimated according to a direction given by the unit vector n like n in the following way:

n = ij n i n j .
We contract the tensor ij with the components n i , n j direction n to make it a representative scalar n . To prove this, let D n be the projection of D following a field E parallel to n. This gives

D n = D • E/|E| with D n = n |E|. We then obtain n |E| = D • E/|E|.
We get:

n |E| = ↔ E • E/|E| and n = ↔ E • E/(|E||E|) = ↔ (E/|E|) • (E/|E|
). Now n = E/|E| and therefore:

n = ij n i n j .
A physical property expressed as a rank-2 tensor like the dielectric constant ij or Shear elastic modulus G(n, q) can be probed along two different directions (generally orthogonal) carried by two unit vectors n and q. Projecting along those two directions, we obtain: n,q = ij n i q j .

C. Symmetry reduction

Starting from the fact that a component tensor T ij... ∼ x i x j x k ... according to Fumi rule [START_REF] Fumi | [END_REF], we can easily simplify tensors by exploiting crystal symmetry. Let us consider a rank-2 tensor in 3D and carry on simplification of its components as we we go along triclinic to monoclinic to orthorhombic, tetragonal and finally cubic symmetry.

• In the triclinic case (no symmetry):

We have 3x3=9 components:

  T xx T xy T xz T yx T yy T yz T zx T zy T zz   .
For any physical tensor (conductivity, permittivity, etc...) the symmetry T ij = T ji is true in general since corresponding energy is akin to a quadratic form (for example, electric energy is σ ij E i E j with σ ij the conductivity tensor and E i the electric field). This symmetry leads to only six non-zero components.

• Monoclinic symmetry:

We have a mirror plane perpendicular to the elementary mesh: z → -z, which causes all components containing z once to be removed by symmetry. There are four components left:

  T xx T xy 0 T yx T yy 0 0 0 T zz  
• Orthorhombic symmetry:

We have 3 possible mirror planes: x → -x, y → -y, z → -z. The non-zero components of the tensor are those containing even combinations of x, y, z. We then have three components:

  T xx 0 0 0 T yy 0 0 0 T zz  
• Tetragonal symmetry:

We have the equivalence x y :

  T xx 0 0 0 T xx 0 0 0 T zz   • Cubic symmetry:
We have the equivalence x y z :

  T xx 0 0 0 T xx 0 0 0 T xx   = T xx 1
with 1 the unit matrix.

NB: All matrix representations of the tensors mentioned above are made in the orthonormal basis {xyz}.

V. APPLICATION TO ANISOTROPIC ELASTICITY

Pressure possesses a tensor character since it involves a force dF and a surface dS = ndS that are both represented mathematically by vectors. Pression is therefore a rank-2 tensor since it associates two vectors (rank-1 tensors). It is represented by the stress tensor σ ij that originates from force i component i "divided" by surface n component j.

Mathematically it is written to bypass division operation as dF

i = σ ij dS j .
Stress can be applied to an object in various fashions:

1. Uniaxial stress: σ ij = σn i n j with σ applied along n orthogonal to a surface element.

2. Hydrostatic stress: σ ij = -σδ ij with σ applied equally along three directions.

3. Simple shear stress: σ ij = σ(n i q j + n j q i ) where σ is applied along direction n belonging to a surface element orthogonal to q. A displacement field δ(r) represents local geometrical alterations of a continuum (isotropic or anisotropic), when subjected to internal or external mechanical efforts.

Adopting the same notation as in ECP we consider x i in a reference frame and call its correspondent xi in the deformed reference (cf fig. 5). Deformation alters basis vectors in a way such that they are ēi in the deformed object whereas they are given by e i in the undeformed one (cf fig. 5). Extending the deformation picture illustrated in fig. 5 from 1D to 3D we define the i-th component of the displacement field as [δ(r)] i = x i (ē i -e i ) (no summation involved). Thus it suffices to define a deformation matrix E ij , i, j = 1, 2, 3 transforming the basis vectors e j into ēi such that ēi = E ij e j , allowing us to define the displacement vector field as: δ(r) = x i (ē i -e i ) (with summation).

Stress efforts are represented by rank-2 tensor σ ij (r) and displacement γ ij a rank-2 tensor is defined by γ ij = ∂δ i ∂x j = E ji . It is decomposable, like any tensor, into a symmetric ij (strain/deformation tensor) and an antisymmetric part ω ij (small rotations tensor):

ij = 1 2 ( ∂δ i ∂x j + ∂δ j ∂x i ), ω ij = 1 2 ( ∂δ i ∂x j - ∂δ j ∂x i ), γ ij = ij + ω ij (17)
In linear elasticity the extension of Hooke's law F = -kx yields a linear relation between σ and , σ ij = C ij,kl kl obtained from the correspondence: F ↔ σ, -k ↔ C ij,kl , x ↔ kl . For a general solid, the 1D spring elastic constant k transforms into a rank-4 tensor C ij,kl containing elastic constants linking σ and that are both rank-2 tensors.

Elastic energy writes: U E = 1 2 C ij,kl ij kl extending the 1D spring energy definition U E = 1 2 kx 2 . Thus:

σ ij = ∂U E ∂ ij
and since is symmetric ( ij = ji ) we infer σ is as well (σ ij = σ ji ). Hooke's law is extended to:

C ij,kl = ∂ 2 U E ∂ ij ∂ kl .
Instead of using Hooke's law, we might use the relationship ij = S ij,kl σ kl with S ij,kl the compliance tensor inverse of C ij,kl . S ij,kl possesses the same symmetry properties as C ij,kl .

B. Elastic moduli

Elastic moduli [START_REF] Royer | Elastic Waves in Solids I: Free and Guided Propagation[END_REF] (Young, shear G, bulk K) and Poisson ratio are given generally along a single direction n or a pair of orthogonal directions n, q with the following operations:

• Strain represented by ij tensor along single direction n :

(∆l/l) n = n = ij n i n j (18)
with l the initial and ∆l its stress induced extra lengths.

• Strain represented by ij tensor along two orthogonal directions n, q :

(∆l/l) nq = nq = ij n i q j [START_REF] Nye | Physical properties of crystals and their representation by tensors and matrices[END_REF] with l the initial and ∆l nq its stress induced extra lengths.

1. Young modulus (see for instance example [START_REF]S3311σ11. Poisson coefficient ν is the ratio of transverse to longitudinal strains: ν = -22 11 = -S 2211 σ 11 S 1111 σ 11 = -ES[END_REF]) along a given direction E(n) Applying a stress σ n i.e. a pressure p along direction n induces a strain n along same direction yielding Young modulus E(n) :

1 E(n) = n σ n (20)
Relating n to σ n is done using σ ij = C ij,kl kl , ij = S ijkl σ kl :

σ n = σ ij n i n j , n = ij n i n j (21)
Introducing σ ij = pn i n j with p the pressure into σ n = σ ij n i n j , we get Young modulus as:

1 E(n) = n σ n = pS ijkl n i n j n k n l p = S ijkl n i n j n k n l (22)
Using Voigt notation (see Appendix) we get the expression of E(n) in the triclinic case for any orientation [START_REF] Ting | [END_REF] as:

1 E(n) =
The example of Silver is displayed in Fig. 6.

Poisson ratio

Applying a pressure p along n and measuring deformation along q perpendicularly to n yields Poisson coefficient (see for instance example [START_REF]S3311σ11. Poisson coefficient ν is the ratio of transverse to longitudinal strains: ν = -22 11 = -S 2211 σ 11 S 1111 σ 11 = -ES[END_REF]) from deformation ratio along directions q and n :

ν(n, q) = - (∆l/l) q (∆l/l) n (25)
Uniaxial pressure p is related to stress tensor by: σ kl = pn k n l whereas strains are given by: (∆l/l) n = ij n i n j , (∆l/l) q = ij q i q j (26)

Using compliance tensor relations:

ij = S ijkl σ kl , σ kl = pn k n l , (27)
we obtain: ij = S ijkl pn k n l .

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5

[100] -1.5 -0.5 0.5 1.5

[010] -3 -2 -1 0 1 2 [001]
FIG. 6: Inverse Young modulus 1 E(n) of Silver versus direction n. Silver compliances [28] are s11 = 2.29, s12 = -0.983, s44 = 2.17 in 1/[100GPa]. Calculated [START_REF] Press | Numerical Recipes in C: The Art of Scientific Computing[END_REF] values along all directions are in 1/[100GPa] units.

Using the above relations, the Poisson ratio is given by:

ν(n, q) = - ij q i q j ij n i n j = - S ijkl q i q j n k n l S ijkl n i n j n k n l (28)
In the triclinic case we get for the numerator (21 terms):

S ijkl q i q j n k n l = s 11 n 2 1 q 2 1 + s 22 n 2 2 q 2 2 + s 33 n 2 3 q 2 3
+ s 44 n 2 n 3 q 2 q 3 + s 55 n 1 n 3 q 1 q 3 + s 66 n 1 n 2 q 1 q 2 + s 12 n 2 1 q 2 2 + s 13 n 2 1 q 2 3 + s 14 n 2 1 q 2 q 3 + s 15 n 2 1 q 1 q 3 + s 16 n 2 1 q 1 q 2 + s 23 n 2 2 q 2 3 + s 24 n 2 2 q 2 q 3 + s 25 n 2 2 q 1 q 3 + s 26 n 2 2 q 1 q 2 + s 34 n 2 3 q 2 q 3 + s 35 n 2 3 q 1 q 3 + s 36 n 2 3 q 1 q 2 + s 45 n 2 n 3 q 1 q 3 + s 46 n 2 n 3 q 1 q 2 + s 56 n 1 n 3 q 1 q 2 [START_REF] Press | Numerical Recipes in C: The Art of Scientific Computing[END_REF] whereas the denominator is E(n) (see eq. 23).

In the cubic case [START_REF] Brantley | [END_REF], the Poisson ratio is (using the conversion rules in Appendix C):

ν(n, q) = - s 12 + (s 11 -s 12 -1 2 s 44 )(n 2 1 q 2 1 + n 2 2 q 2 2 + n 2 3 q 2 3 ) s 11 -2((s 11 -s 12 -1 2 s 44 )(n 2 1 n 2 2 + n 2 2 n 2 3 + n 2 1 n 2 3 ) (30)
The ν(n, q) example of Silver with n in a plane orthogonal to q = [001] is displayed in Fig. 7.

Shear modulus

It is defined by G = Shear stress Angular shear strain (see for instance example [START_REF]Shear modulus is obtained by applying a single stress σ12 inducing an angular deformation 2 12[END_REF]) with two efforts are applied along n and q with n⊥q as required by shear stress. Like Poisson ratio it is a symmetric function of n and q given by G(n, q) defined as inverse ratio of angular shear deformation over shear stress:

1 G(n, q) = 2 ij σ ij nq (31)
This is a generalization of formula 18 to the nq symmetric case to express stress and strain in a symmetric fashion. nq symmetrized angular shear strain is given by: We get the symmetrized stress as:

[2 ij ] nq = 2 ij n i q j = 2S ijkl σ kl n i q j (32)
[σ ij ] nq = σ ij n i q j = σ(n i q j + n j q i )n i q j = σ (33)
Shear modulus is given by:

1 G(n, q) = 2 ij n i q j σ ij n i q j = 2S ijkl σ kl n i q j σ = 2S ijkl (n k q l + n l q k ) n i q j (34)
Using Tensor symmetry S ijkl = S ijlk , we finally obtain:

1 4G(n, q) = S ijkl n i q j n k q l (35)
Using Voigt notation (cf. Appendix B) we get the expression of 1 4G(n,q) in the triclinic case for any orientation [START_REF] Ting | [END_REF] as:

1 G(n, q) = 4[2s 12 -(s 11 + s 22 -s 66 )]n 1 q 1 n 2 q 2 + 4[2s 23 -(s 22 + s 33 -s 44 )]n 2 q 2 n 3 q 3 + 4[2s 31 -(s 33 + s 11 -s 55 )]n 3 q 3 n 1 q 1 + 4(n 1 q 2 + n 2 q 1 )[(s 16 -s 36 )n 1 q 1 + (s 26 -s 36 )n 2 q 2 ] + 4(n 2 q 3 + n 3 q 2 )[(s 24 -s 14 )n 2 q 2 + (s 34 -s 14 )n 3 q 3 ] + 4(n 3 q 1 + n 1 q 3 )[(s 35 -s 25 )n 3 q 3 + (s 15 -s 25 )n 1 q 1 ] + s 44 (n 2 q 3 -n 3 q 2 ) 2 + s 55 (n 3 q 1 -n 1 q 3 ) 2 + s 66 (n 1 q 2 -n 2 q 1 ) 2 + 2s 45 (n 2 q 3 + n 3 q 2 )(n 3 q 1 + n 1 q 3 ) + 2s 56 (n 3 q 1 + n 1 q 3 )(n 1 q 2 + n 2 q 1 ) + 2s 64 (n 1 q 2 + n 2 q 1 )(n 2 q 3 + n 3 q 2 ) (36)
Note that we used above the orthogonality of n, q translating indicially into n i q i = 0. We have when considered as a single elastic coefficient ... In the cubic case, we get (using conversion rules of Appendix C):

1 G(n, q) = s 44 + 4(s 11 -s 12 - 1 2 s 44 )(n 2 1 q 2 1 + n 2 2 q 2 2 + n 2 3 q 2 3 ) (37)
From data compute element point of view, a CPU processes (1 × 1) data units (such as product of two scalars), a GPU a (1 × N ) data unit (such as the scalar product of two N dimensional vectors) whereas TPU processes (N × N ) data units (such as the product of two N × N matrices). A CPU may perform tens of instructions per cycle (IPS), a GPU ∼ 10 4 IPS whereas a TPU may crunch several ∼ 10 5 IPS.

Appendix B: Covariant derivative

In general, the derivative of a tensor (with rank = 0) is not a tensor, since the derivative of the Jacobian terms intervene violating the tensor character as prescribed by ECP. As an example let us take the derivative [START_REF] Weinberg | Gravitation and Cosmology[END_REF][START_REF] Moore | A General Relativity Workbook[END_REF] of a a contravariant vector field (eq. 10) and examine its transformation according to ECP:

∂A i ∂u k = ∂ 2 u i ∂u k ∂ ūj Āj + ∂u i ∂ ūj ∂ Āj ∂u k (B1)
The first term in the RHS breaks ECP since we were expecting a single term with Jacobian factors. Thus we have to introduce a special derivative (called covariant derivative) that absorbs the derivative terms of the Jacobian factors. This mathematical problem disappears in the Cartesian case since the basis vectors do not depend on local coordinates. Let us define the derivative of the covariant basis vector ∂ei ∂u j = Γ k ij e k such that the vector differential dA = d(A i e i ) = ( ∂A i ∂u j du j )e i + A i ∂ei ∂u j du j = [( ∂A i ∂u j e i ) + A i Γ k ij e k ]du j . Note that the coefficients Γ k ij that are called Christoffel symbols are not tensors.

The covariant derivative [( ∂A i ∂u j e i ) + A i Γ k ij e k ] can be rewritten by exchanging the dummy indices i, k in the second term to yield: [( ∂A i ∂u j )+A k Γ i kj ]e i . This allows us to rewrite dA = d(A i e i ) = (D j A i )e i du j where the covariant derivative D j A i = [( ∂A i ∂u j ) + A k Γ i kj ] is a rank-2 tensor containing the Jacobian factor derivative. The evaluation of the Christoffel symbol consists of taking the scalar product with e l of ∂ei ∂u j = Γ k ij e k . Thus e l • ∂ei ∂u j = Γ k ij e l • e k = Γ k ij δ l k = Γ l ij using the orthogonality of e l and e k . The relation Γ l ij = e l • ∂ei ∂u j implies symmetry of the Christoffel symbol Γ l ij = Γ l ji since Γ l ij = e l • ∂ 2 r ∂u i ∂u j after using the definition e i = ∂r ∂u i and allows to relate Γ l ij to the metric tensor by taking its spatial derivative: ∂u i + ∂g ki ∂u j -∂gij ∂u k . Contracting Γ l ij yields: Γ i ij = 1 2 g lk ∂g lk ∂u j since the last two terms cancel by symmetry and index transformation. Using the derivative of determinant g = det[g ij ] formula [START_REF] Golberg | The Amer. Math. Monthly[END_REF] given by ∂g ∂u j = gg ik ∂g ik ∂u j , the contracted Christoffel symbol is rewritten as:

Γ i ij = 1 √ g ∂ √ g ∂u j .
Contrary to simple intuition (like in the rank-2 case) leading to a single elastic constant in a cubic solid, three elastic constants C 11 , C

1 FIG. 1 :

 11 FIG.1:(Color on-line) 2D Covariant and contravariant coordinates of r = OM in the "standard basis" made of covariant vectors (e1, e2). Ordinary (contravariant) coordinates: OA = x 1 and OC = x 2 are made of projections parallel to basis axes whereas covariant coordinates OB = x1 and OD = x2 are made from perpendicular projections. e3 is a unit vector perpendicular to the (e1, e2) plane. Contravariant basis vectors e 1 , e 2 (in red) are such that e1⊥e 2 and e2⊥e 1 . Note that e1 is the projection of e 1 while e2 is the projection of e 2 .

  FIG.2: 3D covariant and contravariant curvilinear coordinates: contravariant vectors are perpendicular to u i (x, y, z) = constant surfaces whereas covariant basis vectors ei are tangent to u i curves that are intersections of u i surfaces. Position vector r = (x, y, z) can be expressed with u i coordinates as r(u 1 , u 2 , u 3 ).

FIG. 3 :

 3 FIG. 3: Pictorial representation of tensor componentsT ijk = A i B j C k , i, j, k = 1, 2, 3at nodes of a 3D lattice. Adapted from ref.[START_REF] Dullemond | Introduction to Tensor Calculus[END_REF]

  ijk lmn = e i • e l e i • e m e i • e n e j • e l e j • e m e j • e n e k • e l e k • e m e k • e n = δ il δ im δ in δ jl δ jm δ jn δ kl δ km δ kn The contraction of two Levi-Civita product:

FIG. 4 :

 4 FIG.4: (Color online) Applied pressure with non collinear force F and normal n to (blue) surface element. When α = 0, we get the ordinary scalar pressure P = dF dS . Moving on to arbitrary α = 0 yields the stress tensor: σij ∼ dF i dSn j rigorously written as dFi = σijdSnj

FIG. 5 :

 5 FIG.5:(Color online) Applying stress to an object results in displacement. This 1D diagram shows an undeformed spring (at left) with a given red marker at xex displaced to a new position xēx after applying stretch. The value x does not change since we are defining positions with respect to corresponding basis vectors ex (initial) and ēx (stretched) yielding the 1D displacement field as δ(x) = x(ēx -ex). One may assume that ēx = Eex where E is the stretching coefficient.

FIG. 7 :

 7 FIG.7: Calculated[START_REF] Press | Numerical Recipes in C: The Art of Scientific Computing[END_REF] Poisson ratio ν(n, q) of Silver versus n angle in [100], [010] plane orthogonal to q = [001]. Silver compliances[28] are s11 = 2.29, s12 = -0.983, s44 = 2.17 in 1/[100GPa].

  e j • ∂ei ∂u k . Combining the spatial derivatives we finally get: Γ l ij = g lk 2 ∂g jk

  12 and C 44 (using Voigt notation explained in Appendix B) are needed in this case yielding: C 11 = C xx,xx = C yy,yy = C zz,zz , C 12 = C xx,yy = C xx,zz = C yy,zz , C 44 = C xy,xy = C xz,xz = C yz,yz , C xy,xz = C xy,yz = C xz,yz = 0. This means one has to distinguish between geometrical symmetry and mechanical symmetry. Consequently we have a hierarchy of symmetries: 1. Cubic symmetry (3 constants: C 11 , C 12 and C 44 ). 2. Isotropic elasticity with rotational symmetry about some reference direction (2 constants: Lamé {λ, µ} coefficients or Young and Poisson {E, ν}). 3. Full isotropy in all directions like in a Newtonian fluid (1 constant). Taking x, y, z along cubic axes, we write C xy = C 44 and C zz = 1 2 (C 11 -C 12 ). Consequently, the elastic constant tensor writes: in the cubic case are: σ xx = C 11 xx + C 12 yy + C 12 zz , σ yy = C 12 xx + C 11 yy + C 12 zz σ zz = C 12 xx + C 12 yy + C 11 zz σ xy = 2C 44 xy , σ xz = 2C 44 xz , σ yz = 2C 44 yz Particular case: In the isotropic elastic case we have only two elastic constants with C xy = C zz , leading to a couple of Lamé coefficients λ = C 12 and µ = C xy = C zz yielding: σ ij = λ kk δ ij + 2µ ij . This means a stress σ 11 induces a longitudinal strain 11 and transverse identical values 22 = 33 resulting in: σ 12 = 2µ 12 and G = µ.

  Note that we have 15 terms instead of 21 due to the mixed coefficient terms such as (s 44 +2s 23 ) and (s 55 +2s 31 )... In the cubic case[START_REF] Ting | [END_REF], using the conversion rules from triclinic to cubic (see Appendix C): s 22 = s 33 = s 11 , s 55 = s 66 = s 44 , s 13 = s 23 = s 12 , s 14 = s 15 = s 16 = s 24 = s 25 = s 26 = s 34 = s 35 = s 36 = s 45 = s 46 = s 56 = 0, we get:

		s 11 n 4 1 + s 22 n 4 2 + s 33 n 4 3 + (s 44 + 2s 23 )n 2 2 n 2 3
		+ (s 55 + 2s 31 )n 2 3 n 2 1 + (s 66 + 2s 12 )n 2 1 n 2 2
		+ 2n 2 n 3 [(s 14 + s 56 )n 2 1 + s 24 n 2 2 + s 34 n 2 3 ]
		+ 2n 3 n 1 [s 15 n 2 1 + (s 25 + s 46 )n 2 2 + s 35 n 2 3 ]
		+ 2n 1 n 2 [s 16 n 2 1 + s 26 n 2 2 + (s 36 + s 45 )n 2 3 ]
				(23)
	1 E(n)	= s 11 (n 4 1 + n 4 2 + n 4 3 ) + 2(s 12 +	1 2	s 44 )(n 2 1 n 2 2 + n 2 2 n 2 3 + n 2 3 n 2 1 )

  15 terms again like in the Young modulus case because of the mixed terms such as [2s 12 -(s 11 +s 22 -s 66 )] or [2s 23 -(s 22 +s 33 -s 44 )]

FIG. 8: Calculated [START_REF] Press | Numerical Recipes in C: The Art of Scientific Computing[END_REF] inverse shear modulus 1 G(n,q) of Silver versus n angle in [100], [010] plane orthogonal to q = [001]. Silver compliances [28] are s11 = 2.29, s12 = -0.983, s44 = 2.17 in 1/[100GPa]. Values along all directions are in 1/[100GPa] units.

The 1 G(n,q) example of Silver versus n in a plane orthogonal to q = [100] is displayed in Fig. 8.

Bulk modulus

The definition K T,S = -1 V ∂V ∂P T,S for isothermal or adiabatic bulk modulus implies a change of volume for an applied hydrostatic pressure.

and using ii = ∆V V , we get K ≈ii ∆P . Using ij = S ijkl σ kl with σ kl = -∆P δ kl , we finally get: K = S iijj .

Thus the most general expression for the bulk modulus exploiting symmetry and Voigt notation (see Appendix B) contains 9 terms since i, j = 1, 2, 3 :

VI. CONCLUSION AND PERSPECTIVES

Tensor calculus is powerful and very useful for Physics students since it provides them not only with elegant procedures to simplify complicated algebraic, vectorial expressions but also to describe anisotropic materials by extending simple scalar physical properties to more general appropriate mathematical expressions.

Appendix A: Tensor etymology and its evolution

The word tensor originally pertains to muscular taxonomy (extensor is a type of muscle implicated in stretch effort as opposed to flexor implicated in angular effort) then moved to differential geometry and more recently to computer science and neural networks. Recently Google built TPU (Tensor Processing Units) for dealing with deep learning problems of Artificial Intelligence deals with cognitive problems such as chess and go games, face recognition, speaker recognition... TPU-based processors are fast because of the massive parallelism of their architecture specially tailored for tackling stacked layers in deep neural network problems. They are faster that traditional CPU (Von Neumann architecture or with some modifications for speedup) and GPU (Graphic Processing Units) dedicated to graphic (Video) operations possessing a Harvard architecture like DSP (Digital signal processor) chips targeted for real-time applications.

Google TPU are made of parallel (systolic) arrays that contain each 65,536 (256 x 256) circuits called "Tensor Processing Elements (TPE)" that perform matrix multiplication with 8-bit based multiply-and-accumulate (MAC) operations in a single clock cycle. The TPU runs at 700 MHz, thus it can compute 65,536x7x10 8 = 46x10 12 MAC per second [START_REF] Jouppi | In-Data-Center Performance Analysis of a Tensor Processing Unit[END_REF].

Moreover it is invariant under {ij} ↔ {kl} interchange. Thus we infer C ij,kl = C ji,kl = C ij,lk = C kl,ij . Thus we may replace a couple of indices ij with a single index I and replace rank-4 tensor C ij,kl = with its matrix representation C IJ with I, J = 1...6.

The index replacement is done according to the recipe: 11 → 1, 22 → 2, 33 → 3, 23 → 4, 13 → 5, 12 → 6. This can be written in a more compact way as: ii → i, ij → 9 -(i + j) when i = j

The 6 × 6 matrix (36 components) represents completely C ij,kl with its 3 4 = 81 components taking account of symmetry. This matrix is in fact symmetric C ij,kl = C kl,ij originating from the property of elastic energy 1 2 C ij,kl ij kl providing another justification for practicality of Voigt notation. We end up with 21 components (since [36-6]/2 +6, i.e. 15+6=21) for the triclinic crystal.

• Application to the cubic crystal case

The number of components in the cubic case can be done along the same lines we treated rank-2 tensors previously.