APSG tutorial - Part 2

from apsg import *

Matrix like classes and tensors

APSG provides matrix-like classes to work with tensor quantities used commonly in structural geology analysis. It includes DeformationGradient3 and VelocityGradient3 for deformation and velocity gradient, Stress3 for stress tensor, Ellipsoid for quadratic forms and OrientationTensor3 for orientation tensor. All these classes support common matrix mathematical operations and provide basic methods and properties.

All matrix-like objects could be created either by passing nested list or tuple or providing individual components to class method from_comp

F = defgrad([[2, 0, 1], [0, 1, 0], [0, 0, 0.5]])
F

DeformationGradient3
[[2.  0.  1. ]
 [0.  1.  0. ]
 [0.  0.  0.5]]

F = defgrad.from_comp(xx=2, zz=0.5, xz=1)
F

DeformationGradient3
[[2.  0.  1. ]
 [0.  1.  0. ]
 [0.  0.  0.5]]

For multiplications of matrix or vectors we have to use matmul @ operator

v = vec('z') # unit-length vector in direction af axis z
u = F @ v
u

Vector3(1, 0, 0.5)

The I property returns the inverse matrix

F.I @ u

Vector3(0, 0, 1)

To transpose matrix, we can use T property and for multiplication we have to use @ operator

F.T @ F

DeformationGradient3
[[4.   0.   2.  ]
 [0.   1.   0.  ]
 [2.   0.   1.25]]

v @ F.T @ F @ v

1.25

Eigenvalues and eigenvectors are obtained with the methods eigenvalues and eigenvectors. Individual eigenvalues and eigenvectors could be accessed by properties E1, E2, E3 and V1, V2, V3

Deformation gradient and rotations

The deformation gradient DeformationGradient3 could describe distortion, dilation, and rigid-body rotation. All APSG features provide a transform method that applies the provided deformation gradient to them.

The rigid-body rotation Rotation3 could be either extracted from the deformation gradient using the R method:

R = F.R
R

Rotation3
[[ 0.928  0.     0.371]
 [ 0.     1.     0.   ]
 [-0.371  0.     0.928]]

or could be created by one of the Rotation3 class methods like from_axisangle, defining axis of rotation and angle:

R = rotation.from_axisangle(lin(120, 50), 60)
R

Rotation3
[[ 0.552 -0.753  0.359]
 [ 0.574  0.655  0.492]
 [-0.605 -0.065  0.793]]

from_two_vectors, where axis of rotation is perpendicular to both vectors and angle is angle of vectors

R = rotation.from_two_vectors(lin(120, 50), lin(270, 80))
R

Rotation3
[[ 0.938  0.074  0.339]
 [ 0.186  0.718 -0.671]
 [-0.294  0.692  0.66 ]]

lin(120, 50).transform(R)

L:270/80

or by from_vectors_axis, where axis does not need to be perpendicular to vectors. Note that rotation axis needs to be adjusted to provide correct rotation of vector.

R = rotation.from_vectors_axis(lin(45,30), lin(135, 30), lin(90, 70))
R

Rotation3
[[-0.393 -0.864  0.315]
 [ 0.864 -0.23   0.448]
 [-0.315  0.448  0.837]]

lin(45,30).transform(R)

L:135/30

a, ang = R.axisangle()
print(lin(a), ang)

L:90/70 113.11571469196133

from_two_pairs method, to describe rotation between two coordinate systems. Note that pair defines the X axis as lineation vector and Z axis as foliation vector.

p1 = pair(150, 60, 90, 40)
p2 = pair(45, 30, 10, 25)
R = rotation.from_two_pairs(p1, p2)
R

Rotation3
[[-0.071  0.97   0.234]
 [-0.874 -0.174  0.453]
 [ 0.48  -0.173  0.86 ]]

p1.transform(R)

P:45/30-10/25

Quaternion and Euler angles representation could be obtained by quat and euler methods:

R.quat()

array([-0.24628826, -0.09662418, -0.72537752,  0.63547881])

R.euler('zxz')

array([109.79936603,  30.68210123, 152.64165336])

Similarly, rotation could be defined by class methods from_quat and from_euler:

R = rotation.from_euler('zxz', [90, -20, -50])
lin(0, 50).transform(R)

L:40/30

R = rotation.from_quat([0, 0, 1, 0.15])
lin(0, 50).transform(R)

L:163/50

Ellipsoid

In deformation analysis, the quadratic forms are represented by Ellipsoid class. It can be used to represent either ellipsoid objects or finite strain ellipsoid.

It provides additional methods and properties including S1, S2 and S3 for principal stretches, Woodcock’s shape and strength, k, K, d and D for Flinn’s and Ramsay symmetries and intensities, lode for Lode’s parameter etc. For more check documentation. Eigenvectors could be also represented by linear or planar features using properties eigenlins and eigenfols.

We can create Ellipsoid object similarly to Matrix3 (note that only components of upper triangular part are available in from_comp method due to matrix symmetry), or you can use additional class methods from_defgrad and from_stretch.

B = ellipsoid.from_defgrad(F)  # Finger deformation tensor
B

Ellipsoid
[[5.   0.   0.5 ]
 [0.   1.   0.  ]
 [0.5  0.   0.25]]
(S1:2.25, S2:1, S3:0.445)

In the above example, the left Cauchy-Green deformation tensor (Finger tensor) B represents finite strain ellipsoid resulting from deformation described by deformation gradient F. We can explore several parameters:

print(f'Principal stretches: Sx={B.S1}, Sy={B.S2}, Sz={B.S3}')
print(f'Principal strain ratios: Rxy={B.Rxy}, Ryz={B.Ryz}')
print(f"Flinn's finite strain parameters: d={B.d}, k={B.k}")
print(f"Ramsay's finite strain parameters: d={B.D}, k={B.K}")
print(f"Woodcock's parameters: strength={B.strength}, shape={B.shape}")
print(f"Watterson's strain intesity: s{B.r}")
print(f"Nadai's natural octahedral unit shear: {B.goct}")
print(f"Nadai's natural octahedral unit strain: {B.eoct}")
print(f"Lode's parameter: {B.lode}")

Principal stretches: Sx=2.247679020649623, Sy=1.0, Sz=0.44490338291762865
Principal strain ratios: Rxy=2.247679020649623, Ryz=2.2476790206496235
Flinn's finite strain parameters: d=1.7644845924910784, k=0.9999999999999997
Ramsay's finite strain parameters: d=1.311869986019497, k=0.9999999999999998
Woodcock's parameters: strength=1.6197962748565, shape=0.9999999999999998
Watterson's strain intesity: s3.495358041299246
Nadai's natural octahedral unit shear: 1.3225581202197993
Nadai's natural octahedral unit strain: 1.1453689300917398
Lode's parameter: 1.3708180983729115e-16

C = ellipsoid.from_defgrad(F, 'right')  # Green's deformation tensor
C

Ellipsoid
[[4.   0.   2.  ]
 [0.   1.   0.  ]
 [2.   0.   1.25]]
(S1:2.25, S2:1, S3:0.445)

v @ C @ v

1.25

Orientation tensor

OrientationTensor3 class represents orientation tensor of set of vectors, linear or planar features. In addition to Ellipsoid methods and properties, it provides properties to describe orientation distribution, e.g. Vollmer’s P, G, R and B indexes, Intensity for Lisle intensity index and MAD for approximate angular deviation.

l = linset.random_fisher(position=lin(120,40))
ot = l.ortensor()
# or
ot = ortensor.from_features(l)
ot

OrientationTensor3
[[ 0.209 -0.258 -0.231]
 [-0.258  0.434  0.342]
 [-0.231  0.342  0.358]]
(S1:0.954, S2:0.227, S3:0.194)

ot.eigenvalues()

array([0.910695  , 0.05168365, 0.03762135])

ot.eigenvectors()

(Vector3(0.443, -0.668, -0.598),
 Vector3(0.011, -0.663, 0.749),
 Vector3(-0.897, -0.338, -0.286))

ot.kind

'L'

The instances of Stress3, Ellipsoid and OrientationTensor3 also provides eigenlins and eigenfols properties to represent principal axes and planes

ot.eigenlins()

(L:124/37, L:271/48, L:21/17)

ot.eigenfols()

(S:304/53, S:91/42, S:201/73)

ot.strength, ot.shape

(1.5933181276385737, 9.034445996834622)

ot.k, ot.d

(18.581928840995992, 3.2023124578445437)

ot.K, ot.D

(9.034445996834622, 2.083098321407397)

ot.P, ot.G, ot.R

(0.8590113445091124, 0.028124596911186656, 0.11286405857970182)

ot.MAD

17.387909861533412