Skip to content

Fix handedness for bingham average - #7297

Open
Theo-Geo wants to merge 2 commits into
geodynamics:mainfrom
Theo-Geo:fix_handedness_for_bingham_average
Open

Fix handedness for bingham average#7297
Theo-Geo wants to merge 2 commits into
geodynamics:mainfrom
Theo-Geo:fix_handedness_for_bingham_average

Conversation

@Theo-Geo

Copy link
Copy Markdown
Contributor

In this pull request i address problems with outputing the bingham average as euler angles.

In the cpo bingham average particle property if Use rotation matrix is set to false, the eigenvectors with the highest eigenvalue are combined as column vectors to form a 3x3 matrix. The problem here is that the bingham average is an antipodal average i.e. the eigenvectors can be either + or - and are still eigenvectors of the second order structure tensor.
When the principal eigenvectors are combined they can form a left-handed coordinate system. This corresponds to a rotation matrix that is not in the group of rotations as it does not preserve orientation. Similarly a reflection also does not preserve orientation and is therefore not in the group of rotations.

Euler angles are a parametrization of the group of rotations, but are not a parametrization of all orthogonal matrices, e.g. they can not parametrize a reflection. This has bad consequences if one wants to recover the principle eigenvector of the first axis also called the mean 100 axis or mean a-axis.
To see the issue i attached a .zip file with a "minimum" working example (it only contains a .prm).
euler_angle_handedness.zip
The parameter file runs a 2d convection box, the crystal preferred orientation plugin and the bingham average.
The steady state convection pattern and bingham average output before applying the changes (No Check for Handness) and after applying the changes (Check for Handedness) can be seen below.

euler_angle_handedness

In this figure the black lines indicate the bingham average a-axis. Usually one would expect that it aligns with the velocity direction in a steady state regime as is created in this convection box. To recreate the a-axis orientation from the euler angles in paraview one can use the calculator and put in:
( cos("cpo mineral 0 phi2")*cos("cpo mineral 0 phi1") - cos("cpo mineral 0 theta")*sin("cpo mineral 0 phi1")*sin("cpo mineral 0 phi2") )*iHat - ( cos("cpo mineral 0 phi2")*sin("cpo mineral 0 phi1") + cos("cpo mineral 0 theta")*cos("cpo mineral 0 phi1")*sin("cpo mineral 0 phi2") )*jHat - sin("cpo mineral 0 phi2")*sin("cpo mineral 0 theta")*kHat

The solution in the right panel was created by checking for handedness before orthogonalizing the rotation matrix.
The handedness $h$ can be computed by taking the cross product between the first two vectors and multiplying it by the last vector,

$$h =( \mathbf{v}^1 \times \mathbf{v}^2) \cdot \mathbf{v}^3 = \sum_{i,j,k=0}^3 \epsilon_{ijk}v^1_i v^2_j v^3_k \text{ .}$$

If h is smaller than 0, $\mathbf{v}^3$ approximately points in the opposite direction of $( \mathbf{v}^1 \times \mathbf{v}^2) $. To ensure the correct handedness one can then multiply $\mathbf{v}^3$ by -1 before constructing the rotation matrix.

TODO I still need to update all tests where euler angles are output, but feel welcome to read and comment on the changes before i do this so i do not have to do it twice.

Before your first pull request:

For all pull requests:

For new features/models or changes of existing features:

  • I have tested my new feature locally to ensure it is correct.
  • I have created a testcase for the new feature/benchmark in the tests/ directory.
  • I have added a changelog entry in the doc/modules/changes directory that will inform other users of my change.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant