Geologic utilities about deformation

PLEASE DO NOT REDISTRIBUTE WITHOUT PERMISSION
2014 February 14 by Joshua R. Davis, jdavis@carleton.edu

Deformations

Here is the matrixLog function. When applied to a position gradient tensor, it produces a canonical choice of velocity gradient tensor. For computing the logarithm of a rotation or ellipsoid, see the Rotations or Ellipsoids sections below.

Clear[jordanBlocks, directSum, blockLog, matrixLog];

(* Returns a list of the Jordan blocks of a Jordan matrix. *)
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}]]];
(* Returns the direct sum of two matrices. *)
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]]];
(* Returns the matrix logarithm of a Jordan block. *)
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];

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

Here are velocity gradient tensors for some common steady homogeneous deformations, along with some functions for changing coordinates with deformation tensors. I’m not sure that I like the latter. The striking function takes an angle (in radians,in geology coordinates) and a velocity gradient tensor for a deformation that strikes 270 degrees,and produces a velocity gradient tensor with the indicated strike. The dipping function takes does the same thing for a 90-degree-striking vertical deformation,using the right-hand rule. Ordinarily,you make a deformation dip,and then make it strike.

Clear[striking, unstriking, dipping, undipping, simpleShear, monoclinicTranspression, inclinedTranspression, triclinicTranspression];

striking[str_, mat_] :=
	Module[{s},
		s = -str - Pi / 2;
		s = {{Cos[s], -Sin[s], 0}, {Sin[s], Cos[s], 0}, {0, 0, 1}};
		s.mat.Transpose[s]];

unstriking[str_, mat_] :=
	striking[-str - Pi, mat];

dipping[dip_, mat_] :=
	Module[{d},
		d = Pi / 2 - dip;
		d = {{1, 0, 0}, {0, Cos[d], -Sin[d]}, {0, Sin[d], Cos[d]}};
		d.mat.Transpose[d]];

undipping[dip_, mat_] :=
	dipping[-dip + Pi, mat];

(* Returns the VGT for 3D simple shear. *)
simpleShear[gamma_] :=
	{{0, gamma, 0}, {0, 0, 0}, {0, 0, 0}};

(* Returns the VGT for monoclinic transpression. Assumes k > 0. *)
monoclinicTranspression[gamma_, k_] :=
	{{0, gamma, 0}, {0, Log[k], 0}, {0, 0, -Log[k]}};

(* Returns the VGT for inclined transpression (Jones et al., 2004). *)
(* Assumes 0 Degree < dip < 180 Degree, 0 Degree < obl < 180 Degree, sho < 1. *)
(* Strike is 270 degrees, according to right-handed strike-dip convention. *)
inclinedTranspression[dip_, obl_, sho_] :=
	Module[{delta, v1, v2, v3, l},
		delta = Pi / 2 - dip;
		v1 = sho Cot[obl] / Cos[delta];
		v2 = -sho;
		v3 = sho Tan[delta];
		l = {{0, Log[1 + v2] v1 / v2, 0}, {0, Log[1 + v2], 0}, {0, Log[1 + v2] v3 / v2, -Log[1 + v2]}};
		dipping[dip, l]];

(* Returns the VGT for triclinic transpression with block velocity v. *)
(* See Davis and Titus (2011; Section 3.2). Assumes v[[2]] > -1. *)
triclinicTranspression[v_] :=
	If[v[[2]] == 0,
		{{0, v[[1]], 0}, {0, 0, 0}, {0, v[[3]], 0}},
		{{0, Log[1 + v[[2]]] v[[1]] / v[[2]], 0}, {0, Log[1 + v[[2]]], 0}, {0, Log[1 + v[[2]]] v[[3]] / v[[2]], -Log[1 + v[[2]]]}}];

Here are some functions for analyzing deformations.

Clear[shearPlanes, shortening, monoclinicAngleOblConv, inclinedTranspressionOffset, dexterity, kinematicVorticity];

(* Returns the kinematic vorticity w_k of Truesdell (1953) for the given VGT. *)
kinematicVorticity[l_] := 
Norm[Flatten[l - Transpose[l]]] / Norm[Flatten[l + Transpose[l]]];

(* Returns one to three possible shear planes for the given VGT. *)
shearPlanes[l_] :=
	Map[unit, Select[Eigenvectors[Transpose[l]], Function[v, Re[v] == v && Norm[v] > 0]]];

(* Returns the shortening perpendicular to the plane with the given pole. *)
shortening[l_, pole_] :=
	1 - Dot[MatrixExp[l].pole, pole] / Dot[pole, pole];

(* Returns the angle of oblique convergence, in radians. *)
monoclinicAngleOblConv[gamma_, k_] :=
	If[gamma == 0, Pi / 2, -ArcTan[Log[k] / gamma]];

(* Returns the dextral offset. *)
inclinedTranspressionOffset[dip_, obl_, sho_, finalWidth_] :=
	Module[{l, f12, initialWidth},
		l = undipping[dip, inclinedTranspression[dip, obl, sho]];
		f12 = MatrixExp[l][[1, 2]];
		initialWidth = 1 / (1 - sho);
		f12 initialWidth];

(* Returns the dextral offset, along the plane with the given pole? *)
dexterity[l_, pole_] :=
	Module[{up},
		up = mostUpwardOrthogonal[pole];
	Dot[MatrixExp[l].pole, Cross[pole, up]] / (Norm[Cross[pole, up]] Norm[pole])];

Euler and Runge-Kutta methods

In this notebook, an "error record" is a list of at least one element, such that the first element is a string. Several functions in this notebook Throw error records to mark an error's happening. The passError function is a helper function for processing error records. The function argument is any Mathematica function. The arguments argument is a list of arguments to be passed to the function. The header argument is an error record. The passError function applies the function to its arguments. If no errors or messages occur, then it returns the result. If an error or message occurs, then it returns an error record, built from the header and the arguments or error record.
	These are general Euler and Runge-Kutta methods on Euclidean spaces and Lie groups. The Lie group methods work on a matrix Lie group with its left- and right-invariant affine connections. In all cases, the input g is the starting point in the group. The input x is a function that takes in a g and outputs a tangent vector (velocity vector). The input n is the number of steps to use. Time runs from 0 to 1. If x is constant, then the Euclidean methods output g + x and the Lie group methods output g (exp x) and (exp x) g. The Lie group Runge-Kutta methods come from Munthe-Kaas (1998) and Munthe-Kaas (1999).

Clear[passError, euclideanEuler, euclideanRungeKutta, leftLieEuler, leftLieRungeKutta, rightLieEuler, dExpInv, rightLieRungeKutta];

(* Helper function for Euler and Runge-Kutta methods. *)
passError[function_, arguments_, header_] :=
	Module[{value},
		value = Catch[Quiet[Check[Apply[function, arguments], {"message:", function, arguments}]]];
		If[ListQ[value] && StringQ[value[[1]]],
			Throw[Append[header, value]],
			value]];

(* Euler and RK4 methods in any Euclidean space. *)
(* g0 = g(0) is the initial condition, at time t = 0. *)
(* x is a function that returns the velocity for any given g. *)
(* n is the number of steps to use. *)
(* Returns g(1), the value of g at time t = 1. *)
euclideanEuler[g0_, x_, n_] :=
	Module[{h, last, i, xval},
		h = 1.0 / n;
		last = g0;
		For[i = 0, i < n, i += 1,
			(* Compute xval = x[last] in a way that lets errors be reported. *)
			xval = passError[x, {last}, {"euclideanEuler: g0, n, i:", g0, n, i}];
			last += h xval];
		last];
euclideanRungeKutta[g0_, x_, n_] :=
	Module[{h, last, i, k1, k2, k3, k4},
		h = 1.0 / n;
		last = g0;
		For[i = 0, i < n, i += 1,
			(* Compute x[last], etc. in a way that lets errors be reported. *)
			k1 = passError[x, {last}, {"euclideanRungeKutta: k1: g, n, i:", g0, n, i}];
			k2 = passError[x, {last + 0.5 h k1}, {"euclideanRungeKutta: k2: g, n, i:", g0, n, i}];
			k3 = passError[x, {last + 0.5 h k2}, {"euclideanRungeKutta: k3: g, n, i:", g0, n, i}];
			k4 = passError[x, {last + h k3}, {"euclideanRungeKutta: k4: g, n, i:", g0, n, i}];
			last += (h / 6.0) (k1 + 2.0 k2 + 2.0 k3 + k4)];
		last];

(* Euler and RK4 methods for the left-invariant affine connection on *)
(* any matrix Lie group. g0 = g(0) is the initial condition, at time t = 0. *)
(* x is a function that returns the velocity for any given g. *)
(* n is the number of steps to use. *)
(* Returns g(1), the value of g at time t = 1. *)
leftLieEuler[g0_, x_, n_] :=
	Module[{h, last, i, xval},
		h = 1.0 / n;
		last = g0;
		For[i = 0, i < n, i += 1,
			(* Compute xval = x[last] in a way that lets errors be reported. *)
			xval = passError[x, {last}, {"leftLieEuler: g, n, i:", g0, n, i}];
			last = last.MatrixExp[h xval]];
		last];
leftLieRungeKutta[g0_, x_, n_] :=
	Module[{h, last, i, k1, k2, k3, k4, v, w},
		h = 1.0 / n;
		last = g0;
		For[i = 0, i < n, i += 1,
			(* Compute xval = x[last], etc. in a way that lets errors be reported. *)
			k1 = passError[x, {last}, {"leftLieRungeKutta: k1: g, n, i:", g0, n, i}];
			k2 = passError[x, {last.MatrixExp[(h / 2) k1]}, {"leftLieRungeKutta: k2: g, n, i:", g0, n, i}];
			k3 = passError[x, {last.MatrixExp[(h / 2) k2 + (h^2 / 24) (k1.k2 - k2.k1)]}, {"leftLieRungeKutta: k3: g, n, i:", g0, n, i}];
			k4 = passError[x, {last.MatrixExp[h k3 + (h^2 / 6) (k1.k3 - k3.k1)]}, {"leftLieRungeKutta: k4: g, n, i:", g0, n, i}];
			v = (h / 6) (k1 + 2 k2 + 2 k3 + k4);
			w = (h / 24) (3 k1 + 2 k2 + 2 k3 - k4);
			last = last.MatrixExp[v + (w.v - v.w)]];
		last];

(* Euler and RK4 methods for the right-invariant affine connection on *)
(* any matrix Lie group. g0 = g(0) is the initial condition, at time t = 0. *)
(* x is a function that returns the velocity for any given g. *)
(* n is the number of steps to use. *)
(* Returns g(1), the value of g at time t = 1. *)
rightLieEuler[g0_, x_, n_] :=
	Module[{h, last, i, xval},
		h = 1.0 / n;
		last = g0;
		For[i = 0, i < n, i += 1,
			(* Compute xval = x[last] in a way that lets errors be reported. *)
			xval = passError[x, {last}, {"rightLieEuler: g, n, i:", g0, n, i}];
			last = MatrixExp[h xval].last];
		last];
dExpInv[u_, v_] :=
	Module[{commutator},
		commutator = u.v - v.u;
		v - commutator / 2.0 + (u.commutator - commutator.u) / 24.0];
rightLieRungeKutta[g0_, x_, n_] :=
	Module[{h, last, i, k1, k2, k3, k4},
		h = 1.0 / n;
		last = g0;
		For[i = 0, i < n, i += 1,
			(* Compute xval = x[last], etc. in a way that lets errors be reported. *)
			k1 = passError[x, {last}, {"rightLieRungeKutta: k1: g, n, i:", g0, n, i}];
			k2 = passError[dExpInv, {(h / 2) k1, x[MatrixExp[(h / 2) k1].last]}, {"rightLieRungeKutta: k2: g, n, i:", g0, n, i}];
			k3 = passError[dExpInv, {(h / 2) k2, x[MatrixExp[(h / 2) k2].last]}, {"rightLieRungeKutta: k3: g, n, i:", g0, n, i}];
			k4 = passError[dExpInv, {h k3, x[MatrixExp[h k3].last]}, {"rightLieRungeKutta: k4: g, n, i:", g0, n, i}];
			last = MatrixExp[(h / 6) (k1 + 2 k2 + 2 k3 + k4)].last];
		last];

Visualizing homogeneous deformation

Here is a general user interface construct for deforming stuff in 2D or 3D. The given matrix must be 2x2 or 3x3, and it must contain t as a free parameter, for this deformer to be useful.

Clear[deformedCurvesPoints, deformer];

(* Helper function for deformer. *)
deformedCurvesPoints[bound_, curves_, points_, matrix_] :=
	Module[{objects},
		objects = {
			RGBColor[0.0, 0.0, 0.0],
			PointSize[Large],
			Point[Map[Function[point, Chop[point.Transpose[matrix]]], points]],
			AbsoluteThickness[2],
			Line[Map[Function[curve, Chop[curve.Transpose[matrix]]], curves]]};
		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]]]];
(* Returns a manipulator for a 2D or 3D deformation. *)
(* bound determines the length of the coordinate axes and often the view. *)
(* curves and points are some curves and points to be deformed. *)
(* matrix is a t-dependent matrix used to deform the curves and points. *)
deformer[bound_, curves_, points_, matrix_] :=
	Manipulate[deformedCurvesPoints[bound, curves, points, matrix], {t, 0, 1}];

The hemisphericalDeformer function is similar to the 3D version of the deformer function, but hemisphericalDeformer projects everything to the unit lower hemisphere and plots it using a hemispherical projection. The user must make sure that the origin does not appear in the points to deform or the curves to deform.

Clear[hemisphericalDeformer];

(* Like deformer, but with a lower hemispherical view instead of a 3D view. *)
Options[hemisphericalDeformer] = {pointSize -> 0.02, projection -> equalAreaProjection};
hemisphericalDeformer[curvesToDeform_, pointsToDeform_, matrix_, OptionsPattern[]] :=
	Manipulate[Module[{points, curves},
			points = Map[Function[point, point.Transpose[matrix]], pointsToDeform];
			curves = Map[Function[curve, curve.Transpose[matrix]], curvesToDeform];
			hemisphericalPlot[{}, points, curves, 
				pointSize -> OptionValue[pointSize], projection -> OptionValue[projection]]],
		{t, 0, 1}];