Homogeneous steady deformation: A review of computational techniques, Appendix F

Joshua R. Davis and Sarah J. Titus

Matrix logarithm

Here is the matrixLog function which takes a finite deformation tensor F to a canonical choice of velocity gradient tensor L. Several helper functions are required.

jordanBlocks[jordan_] :=
	Module[{blocks, n, oldStop, newStop},
		blocks = {};
		n = Length[jordan];
		oldStop = 0;
		newStop = 1;
		While[newStop < n,
			If[jordan[[newStop, newStop + 1]] == 0,
				blocks = Append[blocks, Take[jordan, {oldStop + 1, newStop}, {oldStop + 1, newStop}]];
				oldStop=newStop];
	newStop += 1];
		Append[blocks, Take[jordan, {oldStop + 1, n}, {oldStop + 1, n}]]];

directSum[mat1_, mat2_] :=
	Module[{sumMat, width, height, i},
		If[Length[mat1] == 0, mat2,
			If[Length[mat2] == 0, mat1,
				sumMat = {};
				width = Length[mat1[[1]]] + Length[mat2[[1]]];
				height = Length[mat1];
				For[i = 1, i <= height, i += 1,
					sumMat = Append[sumMat, PadRight[mat1[[i]], width]]];
				height = Length[mat2];
				For[i = 1, i <= height, i += 1,
					sumMat = Append[sumMat, PadLeft[mat2[[i]], width]]];
				sumMat]]];

blockLog[block_] :=
Module[{result, n, i, j},
		result = block;
		n = Length[result];
		For[i = 1, i <= n, i += 1,
			result[[i, i]] = Log[block[[i, i]]];
			For[j = i + 1, j <= n, j += 1,
				result[[i, j]] = (-1)^(j - i + 1)/((j - i) block[[i, i]]^(j - i))]];
		result];

matrixLog[mat_] :=
Module[{p, j, blockList, blockLogList},
		{p, j} = JordanDecomposition[mat];
		blockList = jordanBlocks[j];
		blockLogList = Map[blockLog, blockList];
		p.(Fold[directSum, {}, blockLogList]).Inverse[p]];

To simultaneously superimpose two or more finite deformations, one takes their logs, adds them, and exponentiates. Here is the corrected version of the inclined transpression of Jones and Holdsworth (1998).

Clear[gxy, gzy, az, t];
a = {{1, gxy, 0}, {0, 1, 0}, {0, 0, 1}};
b = {{1, 0, 0}, {0, 1, 0},{0, gzy, 1}};
c = {{1, 0, 0}, {0, 1 / az, 0}, {0, 0, az}};
l = matrixLog[a] + matrixLog[b] + matrixLog[c];
l = Simplify[l, az > 0];
MatrixExp[t l] // MatrixForm

Visualizing progressive deformation

Here is a general user interface construct for visualizing homogeneous steady deformation in 2D or 3D. The argument bound is a positive number indicating how far to extend the coordinate axes. The argument curves to deform is a list of 2D or 3D curves, where each curve is a list of points (to be joined by straight lines). The pointsToDeform argument is a list of 2D or 3D points. The matrix argument must be 2x2 or 3x3, and it must contain t as a free parameter, for this deformer to be useful.

deformer[bound_, curvesToDeform_, pointsToDeform_, matrix_] :=
	Manipulate[Module[{objects},
			objects = {
				RGBColor[0.0, 0.0, 0.0],
				PointSize[Large],
				Point[Map[Function[point, point.Transpose[matrix]], pointsToDeform]],
				AbsoluteThickness[2],
				Line[Map[Function[curve, curve.Transpose[matrix]], curvesToDeform]]};
			If[Length[matrix] == 2,
				Graphics[Join[{
					RGBColor[1.0, 0.0, 0.0],
					Line[{{bound, 0}, {0, 0}}],
					Line[{{-bound, 0}, {-bound/2, 0}}],
					RGBColor[0.0, 1.0, 0.0],
					Line[{{0, bound}, {0, 0}}],
					Line[{{0, -bound}, {0, -bound/2}}]}, objects]],
				Graphics3D[Join[{
					RGBColor[1.0, 0.0, 0.0],
					Line[{{bound, 0, 0}, {0, 0, 0}}],
					Line[{{-bound, 0, 0}, {-bound/2, 0, 0}}],
					RGBColor[0.0, 1.0, 0.0],
					Line[{{0, bound, 0}, {0, 0, 0}}],
					Line[{{0, -bound, 0}, {0, -bound/2, 0}}],
					RGBColor[0.0, 0.0, 1.0],
					Line[{{0, 0, bound}, {0, 0, 0}}],
					Line[{{0, 0, -bound}, {0, 0, -bound/2}}]}, objects]]]],
		{t, 0, 1}];

box2D = {{{1, 1}, {-1, 1}, {-1, -1}, {1, -1}, {1, 1}}, 
{{1, 0}, {0, 0}}, {{0, 1}, {0, 0}}};
circle2D = Table[{Cos[t], Sin[t]}, {t, 0, 2 Pi, 2 Pi / 100}];
rock2D = Append[box2D, circle2D];
box3D = {{{-1,-1,-1},{-1,-1,1}},{{-1,1,-1},{-1,1,1}},{{1,-1,-1},{1,-1,1}},{{1,1,-1},{1,1,1}},{{-1,-1,-1},{-1,1,-1}},{{-1,-1,1},{-1,1,1}},{{1,-1,-1},{1,1,-1}},{{1,-1,1},{1,1,1}},{{-1,-1,-1},{1,-1,-1}},{{-1,-1,1},{1,-1,1}},{{-1,1,-1},{1,1,-1}},{{-1,1,1},{1,1,1}},
{{0,0,0},{0,0,1}},{{0,1,0},{0,1,1}},{{1,0,0},{1,0,1}},{{0,0,0},{0,1,0}},{{0,0,1},{0,1,1}},{{1,0,0},{1,1,0}},{{0,0,0},{1,0,0}},{{0,0,1},{1,0,1}},{{0,1,0},{1,1,0}}};

Here are some examples of 2D finite deformations: simple shear, pure shear, and pure shear about a rotated shear plane.

f = {{1, 1}, {0, 1}};
Clear[t];
deformer[3, rock2D, {}, MatrixExp[t matrixLog[f]]]

f = {{2, 0}, {0, 0.5}};
Clear[t];
deformer[3, rock2D, {}, MatrixExp[t matrixLog[f]]]

f = {{2, 1}, {1, 1}};
Clear[t];
deformer[3, rock2D, {}, MatrixExp[t matrixLog[f]]]

Here are some examples of 3D finite deformations: simple shear and inclined transpression.

f = {{1, 1, 0}, {0, 1, 0}, {0, 0, 1}};
Clear[t];
deformer[3, box3D, {}, MatrixExp[t matrixLog[f]]]

Clear[gxy, gzy, az, t];
gxy = 0.7;
gzy = 0.1;
az = 1.1;
a = {{1, gxy, 0}, {0, 1, 0}, {0, 0, 1}};
b = {{1, 0, 0}, {0, 1, 0},{0, gzy, 1}};
c = {{1, 0, 0}, {0, 1 / az, 0}, {0, 0, az}};
l = matrixLog[a] + matrixLog[b] + matrixLog[c];
l = Simplify[l, az > 0];
deformer[3, box3D, {}, MatrixExp[t l]]

Here is an example of coaxial deformation and coaxial deformation with extraneous rotation. The parameters k and lambda can be adjusted to vary the amount of rotation and stretching. The kinematic vorticity, which is
    (2 Pi / Sqrt[3]) Abs[k / lambda],
can be made arbitrarily large or small. See Figure 4, Section 4.7, and Eq. (C.1).

Clear[k, t, lambda];
lambda = 0.4;
l = {{-lambda, -2 Pi k, 0}, {2 Pi k, -lambda, 0}, {0, 0, 2 lambda}};
deformer[3, box3D, {}, MatrixExp[t l /. {k -> 0}]]
deformer[3, box3D, {}, MatrixExp[t l /. {k -> 1}]]

This final example involves a great deal of shearing and flattening, but it preserves volume and is not physically unreasonable. We have singled out a particular initial position x for special consideration. The 2D plot below the deformer shows time on the horizontal axis and the length squared of the vector (exp t L) x on the vertical axis. The 2D plot shows that the material line that initially passes through the origin and x exhibits an unusual behavior: it elongates, then shortens, then elongates again.

l = {{1.2, 2.3, 0.0}, {0.0, 0.6, 0.0}, {0.0, 0.0, -1.8}};
x = {-0.88,0.48,0.00};
Clear[t];
deformer[10, box3D, {x}, MatrixExp[t l]]
Plot[x.Transpose[MatrixExp[t l]].MatrixExp[t l].x, {t, 0, 1}, AxesOrigin->{0, 0.7}]

Non-uniqueness of logarithms

Let's dig deeper into the logarithms of a finite deformation tensor F with a two-dimensional eigenspace (Appendix C.2). The following code creates a velocity gradient tensor L that depends on two parameters lambda and k, from its complex diagonalization L = P D P^-1. Notice that exp L does not depend on k (as long as k is an integer). Therefore exp L has infinitely many logarithms.

Clear[k, lambda];
d = {{lambda + I k 2 Pi, 0}, {0, lambda - I k 2 Pi}};
d // MatrixForm
p = {{I, -I}, {1, 1}};
l = Simplify[p.d.Inverse[p], Element[k, Integers]];
l // MatrixForm
Simplify[MatrixExp[l], Element[k, Integers]] // MatrixForm

So far, this example is essentially the same as the coaxial deformation with extraneous rotation above. But Gantmacher (1959, p. 239) describes how to use D to create even more real logarithms of exp L. Namely, let S be any complex matrix such that S D S^-1 is real and not equal to D. Then
    exp (S D S^-1) = S (exp D) S^-1 = exp D = exp L,
so S D S^-1 is also a logarithm of exp L. This is powerful in that the matrices S are quite abundant. Here we demonstrate a circle's worth of them, depending on a parameter a. When a = -Pi / 2 we recover S D S^-1 = L.

Clear[a];
s = {{1, Cos[a] + I Sin[a]}, {Cos[a] + I Sin[a], 1}};
Simplify[s.d.Inverse[s], Element[k, Integers]] // MatrixForm
Simplify[s.d.Inverse[s] /. {a -> -Pi / 2}, Element[k, Integers]] // MatrixForm

Here are two examples, that differ only in the parameter a, but that exhibit surprisingly different behavior.

Clear[t];
vgt = Chop[s.d.Inverse[s] /. {lambda -> 0.25, k -> 1, a -> -Pi / 2}];
vgt // MatrixForm
deformer[3, rock2D, {}, Chop[Expand[MatrixExp[t vgt]]]]

Clear[t];
vgt = Chop[s.d.Inverse[s] /. {lambda -> 0.25, k -> 1, a -> -1.1}];
vgt // MatrixForm
deformer[3, rock2D, {}, Chop[Expand[MatrixExp[t vgt]]]]

Averaging finite deformation tensors

The set of finite deformation tensors forms a curved subspace of the space of all 3x3 matrices. If we wish to average two finite deformation tensors, we should not simply average them as matrices; this does not respect their geometric properties such as determinant. On the other hand, the velocity gradient tensors form a flat vector space, where averaging does make sense. So to average two finite deformations we first take their logarithms, then average those, and then exponentiate to obtain the average finite deformation. This can be viewed as taking a geometric mean of finite deformations, rather than an arithmetic mean. This is a first step toward a Lie-group-based statistics of deformation tensors. Notice that this averaging process is nearly identical to simultaneous deformation; the only difference is a factor of 1/2.

Here is an example. We create two finite deformation tensors F1 and F2. (They happen to be rotations, but that is not important.)

zRotation = {{Cos[angle], -Sin[angle], 0}, {Sin[angle], Cos[angle], 0}, {0, 0, 1}};
f1 = N[zRotation /. {angle -> 20 Degree}];
f2 = N[zRotation /. {angle -> 170 Degree}];
MatrixForm[f1]
MatrixForm[f2]

Here's what happens if we simply compute the arithmetic mean of F1 and F2. We are averaging two determinant-1 tensors, and we get a tensor of determinant 0.06. That is undesirable.

arithmetic = (f1 + f2) / 2;
arithmetic // MatrixForm
Det[arithmetic]

In contrast, here's the geometric mean of F1 and F2. It's a rotation by 95 degrees, which is what we expect from the average of a rotation of 20 degrees and a rotation of 170 degrees. Its determinant is 1.

geometric = MatrixExp[(matrixLog[f1] + matrixLog[f2]) / 2];
MatrixForm[Chop[geometric]]
MatrixForm[N[zRotation /. {angle -> 95 Degree}]]