Skip to content

Fix inconsistent strain convention: implement proper Voigt notation with engineering shear strains and consistent variable naming #80

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
Copilot wants to merge 7 commits into
base: master
Choose a base branch
from
Draft

Fix inconsistent strain convention: implement proper Voigt notation with engineering shear strains and consistent variable naming #80

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
23e4b1f
Initial plan
Copilot Jan 4, 2026
d596d20
Fix strain convention to use proper Voigt notation with engineering s…
Copilot Jan 4, 2026
0875a54
Update documentation to clarify Voigt notation convention for strains
Copilot Jan 4, 2026
6afb604
Fix monoclinic and triclinic coupling terms
Copilot Jan 4, 2026
d2ccd99
Fix RST syntax and add clarification about trigonal symmetry
Copilot Jan 4, 2026
339a1be
Rename variables to use Voigt notation consistently throughout symmet…
Copilot Jan 5, 2026
47b4aed
Fix tetragonal symmetry function: correct e4/e5 ordering
Copilot Jan 5, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 12 additions & 2 deletions docs/source/intro.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -159,8 +159,18 @@ symmetry relations and strains :math:`s_{j}^{a}` create a new equation matrix
:math:`S`. :math:`S_{ju}(s^{a})C_{u}=\sigma_{j}^{a}`. The :math:`S(s)` matrix is
a linear function of the strain vector s with all symmetry relations taken into
account. The index a runs over all data sets we have in the calculation while
index u runs over all independent components of the :math:`C_{ij}` matrix. For
the cubic crystal the above equation takes explicit form:
index u runs over all independent components of the :math:`C_{ij}` matrix.

.. note::
In the mathematical formulation below, the strain vector :math:`s` uses
tensor components (i.e., :math:`s_4 = \epsilon_{23}`) and the factor of 2 appears
explicitly in the matrices. However, the Elastic module's ``get_strain()``
function returns strains in standard Voigt notation with engineering shear strains
(i.e., :math:`\epsilon_4 = 2\epsilon_{23}`), and the symmetry functions account
for this convention. Both formulations are mathematically equivalent and yield
identical elastic constants.

For the cubic crystal the above equation takes explicit form:

.. math::
\left[\begin{array}{ccc}
Expand Down
159 changes: 87 additions & 72 deletions elastic/elastic.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,9 +28,12 @@
Elastic is a module for calculation of :math:`C_{ij}` components of elastic
tensor from the strain-stress relation.

The strain components here are ordered in standard way which is different
to ordering in previous versions of the code (up to 4.0).
The ordering is: :math:`u_{xx}, u_{yy}, u_{zz}, u_{yz}, u_{xz}, u_{xy}`.
The strain components here are in standard Voigt notation with engineering
shear strains (also known as the standard Voigt notation for strain).
The ordering is: :math:`\\epsilon_{1}, \\epsilon_{2}, \\epsilon_{3}, \\epsilon_{4}, \\epsilon_{5}, \\epsilon_{6}`,
corresponding to :math:`\\epsilon_{xx}, \\epsilon_{yy}, \\epsilon_{zz}, 2\\epsilon_{yz}, 2\\epsilon_{xz}, 2\\epsilon_{xy}`,
where the shear components (:math:`\\epsilon_{4}, \\epsilon_{5}, \\epsilon_{6}`) are engineering strains
(twice the corresponding tensor shear components).

The general ordering of :math:`C_{ij}` components is (except for triclinic
symmetry and taking into account customary names of constants - e.g.
Expand All @@ -48,11 +51,6 @@
down as a quadratic form of the deformations with respect to elastic constant
and deformation.

*Note:*
The elements for deformations :math:`u_{xy}, u_{xz}, u_{yz}`
have to be divided by 2 to properly match the usual definition
of elastic constants.

See: [LL]_ L.D. Landau, E.M. Lifszyc, "Theory of elasticity"

There is some usefull summary also at:
Expand Down Expand Up @@ -99,19 +97,20 @@ def regular(u):
.. math::
C_{11}, C_{12}, C_{44}

:param u: vector of deformations:
[ :math:`u_{xx}, u_{yy}, u_{zz}, u_{yz}, u_{xz}, u_{xy}` ]
:param u: vector of deformations in Voigt notation:
[ :math:`\\epsilon_{1}, \\epsilon_{2}, \\epsilon_{3}, \\epsilon_{4}, \\epsilon_{5}, \\epsilon_{6}` ]
where shear strains are engineering strains (2× tensor components)

:returns: Symmetry defined stress-strain equation matrix
'''
uxx, uyy, uzz, uyz, uxz, uxy = u[0], u[1], u[2], u[3], u[4], u[5]
e1, e2, e3, e4, e5, e6 = u[0], u[1], u[2], u[3], u[4], u[5]
return array(
[[uxx, uyy + uzz, 0],
[uyy, uxx + uzz, 0],
[uzz, uxx + uyy, 0],
[0, 0, 2*uyz],
[0, 0, 2*uxz],
[0, 0, 2*uxy]])
[[e1, e2 + e3, 0],
[e2, e1 + e3, 0],
[e3, e1 + e2, 0],
[0, 0, e4],
[0, 0, e5],
[0, 0, e6]])


def tetragonal(u):
Expand All @@ -122,20 +121,21 @@ def tetragonal(u):
.. math::
C_{11}, C_{33}, C_{12}, C_{13}, C_{44}, C_{66}

:param u: vector of deformations:
[ :math:`u_{xx}, u_{yy}, u_{zz}, u_{yz}, u_{xz}, u_{xy}` ]
:param u: vector of deformations in Voigt notation:
[ :math:`\\epsilon_{1}, \\epsilon_{2}, \\epsilon_{3}, \\epsilon_{4}, \\epsilon_{5}, \\epsilon_{6}` ]
where shear strains are engineering strains (2× tensor components)

:returns: Symmetry defined stress-strain equation matrix
'''

uxx, uyy, uzz, uyz, uxz, uxy = u[0], u[1], u[2], u[3], u[4], u[5]
e1, e2, e3, e4, e5, e6 = u[0], u[1], u[2], u[3], u[4], u[5]
return array(
[[uxx, 0, uyy, uzz, 0, 0],
[uyy, 0, uxx, uzz, 0, 0],
[0, uzz, 0, uxx+uyy, 0, 0],
[0, 0, 0, 0, 2*uxz, 0],
[0, 0, 0, 0, 2*uyz, 0],
[0, 0, 0, 0, 0, 2*uxy]])
[[e1, 0, e2, e3, 0, 0],
[e2, 0, e1, e3, 0, 0],
[0, e3, 0, e1+e2, 0, 0],
[0, 0, 0, 0, e4, 0],
[0, 0, 0, 0, e5, 0],
[0, 0, 0, 0, 0, e6]])


def orthorombic(u):
Expand All @@ -147,20 +147,21 @@ def orthorombic(u):
C_{11}, C_{22}, C_{33}, C_{12}, C_{13}, C_{23},
C_{44}, C_{55}, C_{66}

:param u: vector of deformations:
[ :math:`u_{xx}, u_{yy}, u_{zz}, u_{yz}, u_{xz}, u_{xy}` ]
:param u: vector of deformations in Voigt notation:
[ :math:`\\epsilon_{1}, \\epsilon_{2}, \\epsilon_{3}, \\epsilon_{4}, \\epsilon_{5}, \\epsilon_{6}` ]
where shear strains are engineering strains (2× tensor components)

:returns: Symmetry defined stress-strain equation matrix
'''

uxx, uyy, uzz, uyz, uxz, uxy = u[0], u[1], u[2], u[3], u[4], u[5]
e1, e2, e3, e4, e5, e6 = u[0], u[1], u[2], u[3], u[4], u[5]
return array(
[[uxx, 0, 0, uyy, uzz, 0, 0, 0, 0],
[0, uyy, 0, uxx, 0, uzz, 0, 0, 0],
[0, 0, uzz, 0, uxx, uyy, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 2*uyz, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 2*uxz, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 2*uxy]])
[[e1, 0, 0, e2, e3, 0, 0, 0, 0],
[0, e2, 0, e1, 0, e3, 0, 0, 0],
[0, 0, e3, 0, e1, e2, 0, 0, 0],
[0, 0, 0, 0, 0, 0, e4, 0, 0],
[0, 0, 0, 0, 0, 0, 0, e5, 0],
[0, 0, 0, 0, 0, 0, 0, 0, e6]])


def trigonal(u):
Expand All @@ -174,22 +175,26 @@ def trigonal(u):
.. math::
C_{11}, C_{33}, C_{12}, C_{13}, C_{44}, C_{14}

:param u: vector of deformations:
[ :math:`u_{xx}, u_{yy}, u_{zz}, u_{yz}, u_{xz}, u_{xy}` ]
:param u: vector of deformations in Voigt notation:
[ :math:`\\epsilon_{1}, \\epsilon_{2}, \\epsilon_{3}, \\epsilon_{4}, \\epsilon_{5}, \\epsilon_{6}` ]
where shear strains are engineering strains (2× tensor components)

:returns: Symmetry defined stress-strain equation matrix
'''

# TODO: Not tested yet.
# TODO: There is still some doubt about the :math:`C_{14}` constant.
uxx, uyy, uzz, uyz, uxz, uxy = u[0], u[1], u[2], u[3], u[4], u[5]
# Note: Some coupling terms retain factors of 2 due to trigonal symmetry
# relationships, not due to Voigt notation. These are part of the
# crystallographic symmetry and are distinct from the engineering strain factors.
e1, e2, e3, e4, e5, e6 = u[0], u[1], u[2], u[3], u[4], u[5]
return array(
[[ uxx, 0, uyy, uzz, 0, 2*uxz ],
[ uyy, 0, uxx, uzz, 0, -2*uxz ],
[ 0, uzz, 0, uxx+uyy, 0, 0 ],
[ 0, 0, 0, 0, 2*uyz, -4*uxy ],
[ 0, 0, 0, 0, 2*uxz, 2*(uxx-uyy)],
[ 2*uxy, 0, -2*uxy, 0, 0, -4*uyz ]])
[[ e1, 0, e2, e3, 0, e5 ],
[ e2, 0, e1, e3, 0, -e5 ],
[ 0, e3, 0, e1+e2, 0, 0 ],
[ 0, 0, 0, 0, e4, -2*e6 ],
[ 0, 0, 0, 0, e5, (e1-e2) ],
[ e6, 0, -e6, 0, 0, -2*e4 ]])


def hexagonal(u):
Expand All @@ -203,21 +208,22 @@ def hexagonal(u):
.. math::
C_{11}, C_{33}, C_{12}, C_{13}, C_{44}

:param u: vector of deformations:
[ :math:`u_{xx}, u_{yy}, u_{zz}, u_{yz}, u_{xz}, u_{xy}` ]
:param u: vector of deformations in Voigt notation:
[ :math:`\\epsilon_{1}, \\epsilon_{2}, \\epsilon_{3}, \\epsilon_{4}, \\epsilon_{5}, \\epsilon_{6}` ]
where shear strains are engineering strains (2× tensor components)

:returns: Symmetry defined stress-strain equation matrix
'''

# TODO: Still needs good verification
uxx, uyy, uzz, uyz, uxz, uxy = u[0], u[1], u[2], u[3], u[4], u[5]
e1, e2, e3, e4, e5, e6 = u[0], u[1], u[2], u[3], u[4], u[5]
return array(
[[ uxx, 0, uyy, uzz, 0 ],
[ uyy, 0, uxx, uzz, 0 ],
[ 0, uzz, 0, uxx+uyy, 0 ],
[ 0, 0, 0, 0, 2*uyz ],
[ 0, 0, 0, 0, 2*uxz ],
[ uxy, 0, -uxy, 0, 0 ]])
[[ e1, 0, e2, e3, 0 ],
[ e2, 0, e1, e3, 0 ],
[ 0, e3, 0, e1+e2, 0 ],
[ 0, 0, 0, 0, e4 ],
[ 0, 0, 0, 0, e5 ],
[ e6, 0, -e6, 0, 0 ]])


def monoclinic(u):
Expand All @@ -229,20 +235,21 @@ def monoclinic(u):
C_{11}, C_{22}, C_{33}, C_{12}, C_{13}, C_{23},
C_{44}, C_{55}, C_{66}, C_{16}, C_{26}, C_{36}, C_{45}

:param u: vector of deformations:
[ :math:`u_{xx}, u_{yy}, u_{zz}, u_{yz}, u_{xz}, u_{xy}` ]
:param u: vector of deformations in Voigt notation:
[ :math:`\\epsilon_{1}, \\epsilon_{2}, \\epsilon_{3}, \\epsilon_{4}, \\epsilon_{5}, \\epsilon_{6}` ]
where shear strains are engineering strains (2× tensor components)

:returns: Symmetry defined stress-strain equation matrix
'''

uxx, uyy, uzz, uyz, uxz, uxy = u[0], u[1], u[2], u[3], u[4], u[5]
e1, e2, e3, e4, e5, e6 = u[0], u[1], u[2], u[3], u[4], u[5]
return array(
[[uxx, 0, 0,uyy,uzz, 0, 0, 0, 0,uxy, 0, 0, 0],
[ 0,uyy, 0,uxx, 0,uzz, 0, 0, 0, 0,uxy, 0, 0],
[ 0, 0,uzz, 0,uxx,uyy, 0, 0, 0, 0, 0,uxy, 0],
[ 0, 0, 0, 0, 0, 0,2*uyz, 0, 0, 0, 0, 0,uxz],
[ 0, 0, 0, 0, 0, 0, 0,2*uxz, 0, 0, 0, 0,uyz],
[ 0, 0, 0, 0, 0, 0, 0, 0,2*uxy,uxx,uyy,uzz, 0]])
[[e1, 0, 0,e2,e3, 0, 0, 0, 0,e6, 0, 0, 0 ],
[ 0,e2, 0,e1, 0,e3, 0, 0, 0, 0,e6, 0, 0 ],
[ 0, 0,e3, 0,e1,e2, 0, 0, 0, 0, 0,e6, 0 ],
[ 0, 0, 0, 0, 0, 0, e4, 0, 0, 0, 0, 0,e5 ],
[ 0, 0, 0, 0, 0, 0, 0, e5, 0, 0, 0, 0,e4 ],
[ 0, 0, 0, 0, 0, 0, 0, 0, e6,e1,e2,e3, 0 ]])


def triclinic(u):
Expand All @@ -259,22 +266,23 @@ def triclinic(u):
C_{16}, C_{26}, C_{36}, C_{46}, C_{56},
C_{14}, C_{15}, C_{25}, C_{45}

:param u: vector of deformations:
[ :math:`u_{xx}, u_{yy}, u_{zz}, u_{yz}, u_{xz}, u_{xy}` ]
:param u: vector of deformations in Voigt notation:
[ :math:`\\epsilon_{1}, \\epsilon_{2}, \\epsilon_{3}, \\epsilon_{4}, \\epsilon_{5}, \\epsilon_{6}` ]
where shear strains are engineering strains (2× tensor components)

:returns: Symmetry defined stress-strain equation matrix
'''

# Based on the monoclinic matrix and not tested on real case.
# If you have test cases for this symmetry send them to the author.
uxx, uyy, uzz, uyz, uxz, uxy = u[0], u[1], u[2], u[3], u[4], u[5]
e1, e2, e3, e4, e5, e6 = u[0], u[1], u[2], u[3], u[4], u[5]
return array(
[[uxx, 0, 0,uyy,uzz, 0, 0, 0, 0,uxy, 0, 0, 0, 0,uyz,uxz, 0, 0],
[ 0,uyy, 0,uxx, 0,uzz, 0, 0, 0, 0,uxy, 0, 0, 0, 0, 0,uxz, 0],
[ 0, 0,uzz, 0,uxx,uyy, 0, 0, 0, 0, 0,uxy, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0,2*uyz, 0, 0, 0, 0, 0,uxy, 0,uxx, 0, 0,uxz],
[ 0, 0, 0, 0, 0, 0, 0,2*uxz, 0, 0, 0, 0, 0,uxy, 0,uxx,uyy,uyz],
[ 0, 0, 0, 0, 0, 0, 0, 0,2*uxy,uxx,uyy,uzz,uyz,uxz, 0, 0, 0, 0]])
[[e1, 0, 0,e2,e3, 0, 0, 0, 0,e6, 0, 0, 0, 0,e4,e5, 0, 0],
[ 0,e2, 0,e1, 0,e3, 0, 0, 0, 0,e6, 0, 0, 0, 0, 0,e5, 0],
[ 0, 0,e3, 0,e1,e2, 0, 0, 0, 0, 0,e6, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, e4, 0, 0, 0, 0, 0,e6, 0,e1, 0, 0,e5],
[ 0, 0, 0, 0, 0, 0, 0, e5, 0, 0, 0, 0, 0,e6, 0,e1,e2,e4],
[ 0, 0, 0, 0, 0, 0, 0, 0, e6,e1,e2,e3,e4,e5, 0, 0, 0, 0]])


def get_cij_order(cryst):
Expand Down Expand Up @@ -705,11 +713,17 @@ def get_strain(cryst, refcell=None):
Computes the strain tensor in the Voight notation as a conventional
6-vector. The calculation is done with respect to the crystal
geometry passed in refcell parameter.

The shear strain components (ε₄, ε₅, ε₆) are returned as engineering
strains, which are twice the corresponding tensor shear components:
ε₄ = 2ε₂₃, ε₅ = 2ε₁₃, ε₆ = 2ε₁₂. This is the standard Voigt notation
convention for strain used in elasticity theory.

:param cryst: deformed structure
:param refcell: reference, undeformed structure

:returns: 6-vector of strain tensor in the Voight notation
[ε₁, ε₂, ε₃, ε₄, ε₅, ε₆] = [ε₁₁, ε₂₂, ε₃₃, 2ε₂₃, 2ε₁₃, 2ε₁₂]
'''
if refcell is None:
refcell = cryst
Expand All @@ -718,7 +732,8 @@ def get_strain(cryst, refcell=None):
m = inv(m)
u = dot(m, du)
u = (u+u.T)/2
return array([u[0, 0], u[1, 1], u[2, 2], u[2, 1], u[2, 0], u[1, 0]])
# Return Voigt notation: [ε₁₁, ε₂₂, ε₃₃, 2ε₂₃, 2ε₁₃, 2ε₁₂]
return array([u[0, 0], u[1, 1], u[2, 2], 2*u[2, 1], 2*u[2, 0], 2*u[1, 0]])


if __name__ == '__main__':
Expand Down