pygplates.Strain

class pygplates.Strain

Bases: instance

The strain at a particular location (parcel of crust) is tracked over time and represents the accumulated deformation undergone by that parcel of crust.

The strain is represented internally by the deformation gradient \(\boldsymbol F\), which is a tensor that transforms a direction \(\hat{\boldsymbol N}\) in the initial configuration (eg, at a time before deformation began) to a direction \(\hat{\boldsymbol n}\) in the current configuration (eg, at a time during or after deformation), and stretches it by a factor of \(\Lambda\):

\[\hat{\boldsymbol n} \Lambda = \boldsymbol F \cdot \hat{\boldsymbol N}\]

Note

\(\hat{\boldsymbol N}\) is a 2D unit vector direction in the local South-East coordinate system at a location (\(\Theta, \Phi\)) that a point occupies at the initial time. \(\hat{\boldsymbol n}\) is a 2D unit vector direction in the local South-East coordinate system at a new location (\(\theta, \phi\)) that the same point now occupies at the current time (eg, it has moved as its tectonic plate rotates/deforms). For example, if an initial South direction transforms to an Eastwards direction then \(\hat{\boldsymbol N} = (1, 0)\) and \(\hat{\boldsymbol n} = (0, 1)\).

The Lagrangian deformation tensor \(\boldsymbol C\) is defined in terms of the deformation gradient tensor \(\boldsymbol F\):

\[\boldsymbol C = \boldsymbol{F}^T \cdot \boldsymbol F\]

…and determines the stretch factor \(\Lambda_{(\hat{\boldsymbol N})}\) given a direction \(\hat{\boldsymbol N}\) in the initial configuration:

\[\Lambda_{(\hat{\boldsymbol N})} = \sqrt{\hat{\boldsymbol N} \cdot \boldsymbol C \cdot \hat{\boldsymbol N}}\]

The Eulerian deformation tensor \(\boldsymbol c\) is also defined in terms of the deformation gradient tensor \(\boldsymbol F\):

\[\boldsymbol c = (\boldsymbol{F}^{-1})^T \cdot \boldsymbol{F}^{-1}\]

…and determines the stretch factor \(\Lambda_{(\hat{\boldsymbol n})}\) given a direction \(\hat{\boldsymbol n}\) in the current configuration:

\[\frac{1}{\Lambda_{(\hat{\boldsymbol n})}} = \sqrt{\hat{\boldsymbol n} \cdot \boldsymbol c \cdot \hat{\boldsymbol n}}\]

Note

These stretch factors (\(\Lambda_{(\hat{\boldsymbol N})}\) and \(\Lambda_{(\hat{\boldsymbol n})}\)) are stretches along the surface direction, not the depth direction (as is the case with crustal thinning). Although surface stretch is used to determine crustal thinning (obtained by passing pygplates.ScalarType.gpml_crustal_stretching_factor to ReconstructedGeometryTimeSpan.get_scalar_values() returned by TopologicalModel.reconstruct_geometry()).

The normal strain is the change in length per unit initial length in a given direction, and is the stretch factor minus one. For a direction \(\hat{\boldsymbol N}\) in the initial configuration, the normal strain is:

\[\begin{split}e_{(\hat{\boldsymbol N})} &= \Lambda_{(\hat{\boldsymbol N})} - 1\\ &= \sqrt{\hat{\boldsymbol N} \cdot \boldsymbol C \cdot \hat{\boldsymbol N}} - 1\end{split}\]

And for a direction \(\hat{\boldsymbol n}\) in the current configuration, the normal strain is:

\[\begin{split}e_{(\hat{\boldsymbol n})} &= \Lambda_{(\hat{\boldsymbol n})} - 1\\ &= \frac{1}{\sqrt{\hat{\boldsymbol n} \cdot \boldsymbol c \cdot \hat{\boldsymbol n}}} - 1\end{split}\]

The Lagrangian finite strain tensor \(\boldsymbol E\) is defined in terms of the Lagrangian deformation tensor \(\boldsymbol C\):

\[\boldsymbol E = \frac{1}{2} (\boldsymbol C - \boldsymbol I)\]

…and it’s called a strain tensor because for very small strains the normal strain \(e_{(\hat{\boldsymbol N})}\) can be approximated as:

\[e_{(\hat{\boldsymbol N})} \simeq \hat{\boldsymbol N} \cdot \boldsymbol E \cdot \hat{\boldsymbol N}\]

The Eulerian finite strain tensor \(\boldsymbol e\) is defined in terms of the Eulerian deformation tensor \(\boldsymbol c\):

\[\boldsymbol e = \frac{1}{2} (\boldsymbol I - \boldsymbol c)\]

…and for very small strains the normal strain \(e_{(\hat{\boldsymbol n})}\) can be approximated as:

\[e_{(\hat{\boldsymbol n})} \simeq \frac{\hat{\boldsymbol n} \cdot \boldsymbol e \cdot \hat{\boldsymbol n}}{\hat{\boldsymbol n} \cdot \boldsymbol c \cdot \hat{\boldsymbol n}}\]

The angle \(\alpha\) between two directions in the current configuration that were originally in the directions \(\hat{\boldsymbol N}_1\) and \(\hat{\boldsymbol N}_2\) in the initial configuration is given by:

\[\cos{\alpha} = \frac{\hat{\boldsymbol N}_1 \cdot \boldsymbol C \cdot \hat{\boldsymbol N}_2}{\Lambda_{(\hat{\boldsymbol N}_1)} \Lambda_{(\hat{\boldsymbol N}_2)}}\]

…and the angle \(A\) between two directions in the initial configuration that are currently in the directions \(\hat{\boldsymbol n}_1\) and \(\hat{\boldsymbol n}_2\) in the current configuration is given by:

\[\cos{A} = \Lambda_{(\hat{\boldsymbol n}_1)} \Lambda_{(\hat{\boldsymbol n}_2)} (\hat{\boldsymbol n}_1 \cdot \boldsymbol c \cdot \hat{\boldsymbol n}_2)\]

References:

Strains are equality (==, !=) comparable (but not hashable - cannot be used as a key in a dict).

Convenience class static data is available for the identity strain:

  • pygplates.Strain.identity

A Strain can also be pickled.

Added in version 0.46.

__init__(...)

A Strain object can be constructed in more than one way…

__init__()

Construct an identity strain (no deformation).

identity_strain = pygplates.Strain()

Note

Alternatively you can use identity_strain = pygplates.Strain.identity.

__init__(deformation_gradient_theta_theta,, deformation_gradient_theta_phi, deformation_gradient_phi_theta, deformation_gradient_phi_phi)

Create from the deformation gradient \(\boldsymbol F\) in spherical polar coordinates (ignoring radial dimension).

\[\begin{split}\boldsymbol F = \begin{bmatrix} F_{\theta\theta} & F_{\theta\phi} \\ F_{\phi\theta} & F_{\phi\phi} \end{bmatrix}\end{split}\]
param deformation_gradient_theta_theta:

\(F_{\theta\theta}\)

type deformation_gradient_theta_theta:

float

param deformation_gradient_theta_phi:

\(F_{\theta\phi}\)

type deformation_gradient_theta_phi:

float

param deformation_gradient_phi_theta:

\(F_{\phi\theta}\)

type deformation_gradient_phi_theta:

float

param deformation_gradient_phi_phi:

\(F_{\phi\phi}\)

type deformation_gradient_phi_phi:

float

See also

get_deformation_gradient()

Note

Typically you wouldn’t need this since you can start with no deformation (pygplates.Strain.identity) at the initial time and use accumulate() to incrementally update strain over time using your own calculations of strain rate.

Methods

__init__(...)

A Strain object can be constructed in more than one way...

accumulate(previous_strain, ...)

[staticmethod] Accumulate previous strain using both previous and current strain rates (in units of \(second^{-1}\)) over a time increment (in units of \(second\)).

get_deformation_gradient()

Return the deformation gradient tensor \(\boldsymbol F\) in spherical polar coordinates (ignoring radial dimension).

get_dilatation()

Return the change in crustal area with respect to the original area.

get_principal_strain(...)

Return the maximum and minimum strains (along principal axes), and the angle of the major principal axis.

Attributes

identity

static accumulate(previous_strain, previous_strain_rate, current_strain_rate, time_increment)

[staticmethod] Accumulate previous strain using both previous and current strain rates (in units of \(second^{-1}\)) over a time increment (in units of \(second\)).

Parameters:

  • previous_strain (Strain) – the previous strain

  • previous_strain_rate (StrainRate) – the previous strain rate

  • current_strain_rate (StrainRate) – the current strain rate

  • time_increment (float) – the time increment to accumulate strain over (in units of \(second^{-1}\))

Returns:

the current strain (accumulated from previous strain)

Return type:

Strain

To accumulate strain from an initial undeformed state at 100Ma to its final deformed strain at present day:

time_increment_1myr_in_seconds = 1e6 * 365 * 24 * 60 * 60
previous_strain = pygplates.Strain.identity
previous_strain_rate = pygplates.StrainRate.zero

for time in range(100, -1, -1):
    current_strain_rate = pygplates.StrainRate(...)
    current_strain = pygplates.Strain.accumulate(previous_strain,
        previous_strain_rate, current_strain_rate, time_increment_1myr_in_seconds)

    previous_strain = current_strain
    previous_strain_rate = current_strain_rate

Strain is accumulated by integrating the ordinary differential equation defining the rate of change of the deformation gradient \(\boldsymbol F\) in terms of the spatial gradients of velocity \(\boldsymbol L\):

\[\dot{\boldsymbol F} = \boldsymbol L \cdot \boldsymbol F\]

…which is approximated using the central differencing scheme over a time increment \(\Delta t\):

\[\frac{\boldsymbol{F}_{t + \Delta t} - \boldsymbol{F}}{\Delta t} = \frac{\boldsymbol{L}_{t + \Delta t} \boldsymbol{F}_{t + \Delta t} + \boldsymbol{L}_{t} \boldsymbol{F}_{t}}{2}\]

…which rearranges to become (in matrix form, with identity matrix \(\boldsymbol I\)):

\[\boldsymbol{F}_{t + \Delta t} = (\boldsymbol I - \boldsymbol{L}_{t + \Delta t} \frac{\Delta t}{2})^{-1} (\boldsymbol I + \boldsymbol{L}_{t} \frac{\Delta t}{2}) \boldsymbol{F}_{t}\]

…where \(\Delta t\) is time_increment, \(\boldsymbol{L}_{t}\) is previous_strain_rate, \(\boldsymbol{L}_{t + \Delta t}\) is current_strain_rate, \(\boldsymbol{F}_{t}\) is previous_strain and \(\boldsymbol{F}_{t + \Delta t}\) is returned by this function.

get_deformation_gradient()

Return the deformation gradient tensor \(\boldsymbol F\) in spherical polar coordinates (ignoring radial dimension).

\[\begin{split}\boldsymbol F = \begin{bmatrix} F_{\theta\theta} & F_{\theta\phi} \\ F_{\phi\theta} & F_{\phi\phi} \end{bmatrix}\end{split}\]
Returns:

the tuple of \((F_{\theta\theta}, F_{\theta\phi}, F_{\phi\theta}, F_{\phi\phi})\)

Return type:

tuple (float, float, float, float)

Note

\(\theta\) is co-latitude and hence increases from North to South (and \(\phi\) increases from West to East, as expected).

get_dilatation()

Return the change in crustal area with respect to the original area.

Return type:

float

The dilatation is the increase (if positive) or decrease (if negative) of crustal area with respect to the original area in the initial configuration (eg, at a time before deformation began).

The dilatation is defined as the determinant of the deformation gradient tensor \(\boldsymbol F\) minus one. So if we define \(A\) as the original area of a parcel of crust in the initial configuration, then the dilatation is given by:

\[\frac{\Delta A}{A} = \det \boldsymbol F - 1\]

See also

Chapter 4.11 in Continuum Mechanics for Engineers for a derivation of the change in volume. We ignore the radial dimension, hence volume becomes area.

Note

The dilatation is invariant with respect to the local coordinate axes (South and East) since \(\det \boldsymbol F\) is the third invariant of \(\boldsymbol F\) (see Chapter 3.6 in Continuum Mechanics for Engineers).

get_principal_strain([principal_angle_type=PrincipalAngleType.major_south])

Return the maximum and minimum strains (along principal axes), and the angle of the major principal axis.

Parameters:

principal_angle_type (PrincipalAngleType.major_south, PrincipalAngleType.major_east or PrincipalAngleType.major_azimuth) – how the angle of the major principal axis is defined relative to the local coordinate system (defaults to PrincipalAngleType.major_south)

Returns:

the tuple of maximum strain, minimum strain and major axis angle \((e_{(1)}, e_{(2)}, \alpha)\)

Return type:

tuple (float, float, float)

principal_angle_type supports the following enumeration types:

Value

Description

PrincipalAngleType.major_south

The major principal axis points South when the angle is zero. The angle ranges from \(-\pi\) to \(\pi\) radians anti-clockwise (observed from above the globe).

PrincipalAngleType.major_east

The major principal axis points East when the angle is zero. The angle ranges from \(-\pi\) to \(\pi\) radians anti-clockwise (observed from above the globe). This is equiavlent to MajorAngle in the GPlates deformation export.

PrincipalAngleType.major_azimuth

The major principal axis points North when the angle is zero. The angle ranges from \(0\) to \(2\pi\) radians clockwise (observed from above the globe). This is equiavlent to MajorAzimuth in the GPlates deformation export.

Note

Regardless of the value of principal_angle_type, the direction of the minimum principal axis is always an anti-clockwise rotation of \(\frac{\pi}{2}\) radians (90 degrees) of the major principal axis (observed from above the globe).

The principal strains are the maximum and minimum strains that occur along the principal axes (where shear strain is zero). The principal axes are the coordinate axes rotated anti-clockwise (when observed from above the globe) by an angle \(\alpha\) which is defined in terms of the Eulerian deformation tensor \(\boldsymbol c\) (see Strain):

\[\tan{2 \alpha} = \frac{2 c_{\theta\phi}}{c_{\theta\theta} - c_{\phi\phi}}\]

Then the two perpendicular principal strain directions are (in the local South-East coordinate system):

\[\begin{split}\hat{\boldsymbol n}_{(1)} &= (\cos \alpha, \sin \alpha)\\ \hat{\boldsymbol n}_{(2)} &= (-\sin \alpha, \cos \alpha)\end{split}\]

…and the two principal strains are determined by (see Strain):

\[\begin{split}e_{(1)} &= \Lambda_{(1)} - 1 &= \frac{1}{\sqrt{\hat{\boldsymbol n}_{(1)} \cdot \boldsymbol c \cdot \hat{\boldsymbol n}_{(1)}}} - 1\\ e_{(2)} &= \Lambda_{(2)} - 1 &= \frac{1}{\sqrt{\hat{\boldsymbol n}_{(2)} \cdot \boldsymbol c \cdot \hat{\boldsymbol n}_{(2)}}} - 1\end{split}\]

Note

This function returns principal strains. To get the principal stretches simply add one since \(\Lambda_{(i)} = 1 + e_{(i)}\) (see Strain).

To get the principal stretches (ie, strains plus one) with the angle (in degrees) of the major principal axis clockwise relative to North (ie, azimuth):

import math
...
max_strain, min_strain, major_azimuth_radians = strain.get_principal_strain(
    principal_angle_type=pygplates.PrincipalAngleType.major_azimuth)

max_stretch = 1 + max_strain
min_stretch = 1 + min_strain
major_azimuth_degrees = math.degrees(major_azimuth_radians)