added rotation conversions

modified versions from 3Drotations code (available on GitHub) by Marc De Graef
This commit is contained in:
Martin Diehl 2018-12-08 08:02:55 +01:00
parent 590f83a944
commit 40d38ebf55
6 changed files with 1757 additions and 6 deletions

View File

@ -49,18 +49,30 @@ add_library(FEsolving OBJECT "FEsolving.f90")
add_dependencies(FEsolving DEBUG)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:FEsolving>)
add_library(DAMASK_MATH OBJECT "math.f90")
add_dependencies(DAMASK_MATH FEsolving)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:DAMASK_MATH>)
add_library(MATH OBJECT "math.f90")
add_dependencies(MATH DEBUG)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:MATH>)
add_library(QUATERNIONS OBJECT "quaternions.f90")
add_dependencies(QUATERNIONS MATH)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:QUATERNIONS>)
add_library(LAMBERT OBJECT "Lambert.f90")
add_dependencies(LAMBERT MATH)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:LAMBERT>)
add_library(ROTATIONS OBJECT "rotations.f90")
add_dependencies(ROTATIONS LAMBERT QUATERNIONS)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:ROTATIONS>)
# SPECTRAL solver and FEM solver use different mesh files
if (PROJECT_NAME STREQUAL "DAMASK_spectral")
add_library(MESH OBJECT "mesh.f90")
add_dependencies(MESH DAMASK_MATH)
add_dependencies(MESH ROTATIONS FEsolving)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:MESH>)
elseif (PROJECT_NAME STREQUAL "DAMASK_FEM")
add_library(FEZoo OBJECT "FEM_zoo.f90")
add_dependencies(FEZoo DAMASK_MATH)
add_dependencies(FEZoo ROTATIONS FEsolving)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:FEZoo>)
add_library(MESH OBJECT "meshFEM.f90")
add_dependencies(MESH FEZoo)

213
src/Lambert.f90 Normal file
View File

@ -0,0 +1,213 @@
! ###################################################################
! Copyright (c) 2013-2015, Marc De Graef/Carnegie Mellon University
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification, are
! permitted provided that the following conditions are met:
!
! - Redistributions of source code must retain the above copyright notice, this list
! of conditions and the following disclaimer.
! - Redistributions in binary form must reproduce the above copyright notice, this
! list of conditions and the following disclaimer in the documentation and/or
! other materials provided with the distribution.
! - Neither the names of Marc De Graef, Carnegie Mellon University nor the names
! of its contributors may be used to endorse or promote products derived from
! this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
! USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! ###################################################################
!--------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief everything that has to do with the modified Lambert projections
!
!> @details This module contains a number of projection functions for the modified
!> Lambert projection between square lattice and 2D hemisphere, hexagonal lattice
!> and 2D hemisphere, as well as the more complex mapping between a 3D cubic grid
!> and the unit quaternion hemisphere with positive scalar component. In addition, there
!> are some other projections, such as the stereographic one. Each function is named
!> by the projection, the dimensionality of the starting grid, and the forward or inverse
!> character. For each function, there is also a single precision and a double precision
!> version, but we use the interface formalism to have only a single call. The Forward
!> mapping is taken to be the one from the simple grid to the curved grid. Since the module
!> deals with various grids, we also add a few functions/subroutines that apply symmetry
!> operations on those grids.
!> References:
!> D. Rosca, A. Morawiec, and M. De Graef. A new method of constructing a grid
!> in the space of 3D rotations and its applications to texture analysis.
!> Modeling and Simulations in Materials Science and Engineering 22, 075013 (2014).
!--------------------------------------------------------------------------
module Lambert
use math
use prec
implicit none
real(pReal), private :: &
sPi = sqrt(PI), &
pref = sqrt(6.0_pReal/PI), &
! the following constants are used for the cube to quaternion hemisphere mapping
ap = PI**(2.0_pReal/3.0_pReal), &
sc = 0.897772786961286_pReal, & ! a/ap
beta = 0.962874509979126_pReal, & ! pi^(5/6)/6^(1/6)/2
R1 = 1.330670039491469_pReal, & ! (3pi/4)^(1/3)
r2 = sqrt(2.0_pReal), &
pi12 = PI/12.0_pReal, &
prek = 1.643456402972504_pReal, & ! R1 2^(1/4)/beta
r24 = sqrt(24.0_pReal)
private
public :: &
LambertCubeToBall, &
LambertBallToCube
private :: &
GetPyramidOrder
contains
!--------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief map from 3D cubic grid to 3D ball
!--------------------------------------------------------------------------
function LambertCubeToBall(cube) result(ball)
use, intrinsic :: IEEE_ARITHMETIC
implicit none
real(pReal), intent(in), dimension(3) :: cube
real(pReal), dimension(3) :: ball, LamXYZ, XYZ
real(pReal) :: T(2), c, s, q
real(pReal), parameter :: eps = 1.0e-8_pReal
integer(pInt), dimension(3) :: p
integer(pInt), dimension(2) :: order
if (maxval(abs(cube)) > ap/2.0+eps) then
ball = IEEE_value(cube,IEEE_positive_inf)
return
end if
! transform to the sphere grid via the curved square, and intercept the zero point
center: if (all(dEq0(cube))) then
ball = 0.0_pReal
else center
! get pyramide and scale by grid parameter ratio
p = GetPyramidOrder(cube)
XYZ = cube(p) * sc
! intercept all the points along the z-axis
special: if (all(dEq0(XYZ(1:2)))) then
LamXYZ = [ 0.0_pReal, 0.0_pReal, pref * XYZ(3) ]
else special
order = merge( [2,1], [1,2], abs(XYZ(2)) <= abs(XYZ(1))) ! order of absolute values of XYZ
q = pi12 * XYZ(order(1))/XYZ(order(2)) ! smaller by larger
c = cos(q)
s = sin(q)
q = prek * XYZ(order(2))/ sqrt(r2-c)
T = [ (r2*c - 1.0), r2 * s] * q
! transform to sphere grid (inverse Lambert)
! [note that there is no need to worry about dividing by zero, since XYZ(3) can not become zero]
c = sum(T**2)
s = Pi * c/(24.0*XYZ(3)**2)
c = sPi * c / r24 / XYZ(3)
q = sqrt( 1.0 - s )
LamXYZ = [ T(order(2)) * q, T(order(1)) * q, pref * XYZ(3) - c ]
endif special
! reverse the coordinates back to the regular order according to the original pyramid number
ball = LamXYZ(p)
endif center
end function LambertCubeToBall
!--------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief map from 3D ball to 3D cubic grid
!--------------------------------------------------------------------------
pure function LambertBallToCube(xyz) result(cube)
use, intrinsic :: IEEE_ARITHMETIC
implicit none
real(pReal), intent(in), dimension(3) :: xyz
real(pReal), dimension(3) :: cube, xyz1, xyz3
real(pReal), dimension(2) :: Tinv, xyz2
real(pReal) :: rs, qxy, q2, sq2, q, tt
integer(pInt) , dimension(3) :: p
rs = norm2(xyz)
if (rs > R1) then
cube = IEEE_value(cube,IEEE_positive_inf)
return
endif
center: if (all(dEq0(xyz))) then
cube = 0.0_pReal
else center
p = GetPyramidOrder(xyz)
xyz3 = xyz(p)
! inverse M_3
xyz2 = xyz3(1:2) * sqrt( 2.0*rs/(rs+abs(xyz3(3))) )
! inverse M_2
qxy = sum(xyz2**2)
special: if (dEq0(qxy)) then
Tinv = 0.0
else special
q2 = qxy + maxval(abs(xyz2))**2
sq2 = sqrt(q2)
q = (beta/r2/R1) * sqrt(q2*qxy/(q2-maxval(abs(xyz2))*sq2))
tt = (minval(abs(xyz2))**2+maxval(abs(xyz2))*sq2)/r2/qxy
Tinv = q * sign(1.0,xyz2) * merge([ 1.0_pReal, acos(math_clip(tt,-1.0_pReal,1.0_pReal))/pi12], &
[ acos(math_clip(tt,-1.0_pReal,1.0_pReal))/pi12, 1.0_pReal], &
abs(xyz2(2)) <= abs(xyz2(1)))
endif special
! inverse M_1
xyz1 = [ Tinv(1), Tinv(2), sign(1.0,xyz3(3)) * rs / pref ] /sc
! reverst the coordinates back to the regular order according to the original pyramid number
cube = xyz1(p)
endif center
end function LambertBallToCube
!--------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief determine to which pyramid a point in a cubic grid belongs
!--------------------------------------------------------------------------
pure function GetPyramidOrder(xyz)
implicit none
real(pReal),intent(in),dimension(3) :: xyz
integer(pInt), dimension(3) :: GetPyramidOrder
if (((abs(xyz(1)) <= xyz(3)).and.(abs(xyz(2)) <= xyz(3))) .or. &
((abs(xyz(1)) <= -xyz(3)).and.(abs(xyz(2)) <= -xyz(3)))) then
GetPyramidOrder = [1,2,3]
else if (((abs(xyz(3)) <= xyz(1)).and.(abs(xyz(2)) <= xyz(1))) .or. &
((abs(xyz(3)) <= -xyz(1)).and.(abs(xyz(2)) <= -xyz(1)))) then
GetPyramidOrder = [2,3,1]
else if (((abs(xyz(1)) <= xyz(2)).and.(abs(xyz(3)) <= xyz(2))) .or. &
((abs(xyz(1)) <= -xyz(2)).and.(abs(xyz(3)) <= -xyz(2)))) then
GetPyramidOrder = [3,1,2]
else
GetPyramidOrder = -1 ! should be impossible, but might simplify debugging
end if
end function GetPyramidOrder
end module Lambert

View File

@ -11,6 +11,9 @@
#include "debug.f90"
#include "config.f90"
#include "math.f90"
#include "quaternions.f90"
#include "Lambert.f90"
#include "rotations.f90"
#include "FEsolving.f90"
#include "mesh.f90"
#include "material.f90"

View File

@ -30,7 +30,9 @@ module prec
integer(pInt), allocatable, dimension(:) :: realloc_lhs_test
type, public :: group_float !< variable length datatype used for storage of state
real(pReal), parameter, public :: epsijk = -1.0_pReal !< parameter for orientation conversion. ToDo: Better place?
type, public :: group_float !< variable length datatype used for storage of state
real(pReal), dimension(:), pointer :: p
end type group_float

433
src/quaternions.f90 Normal file
View File

@ -0,0 +1,433 @@
! ###################################################################
! Copyright (c) 2013-2015, Marc De Graef/Carnegie Mellon University
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification, are
! permitted provided that the following conditions are met:
!
! - Redistributions of source code must retain the above copyright notice, this list
! of conditions and the following disclaimer.
! - Redistributions in binary form must reproduce the above copyright notice, this
! list of conditions and the following disclaimer in the documentation and/or
! other materials provided with the distribution.
! - Neither the names of Marc De Graef, Carnegie Mellon University nor the names
! of its contributors may be used to endorse or promote products derived from
! this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
! USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! ###################################################################
module quaternions
use prec
implicit none
public
type, public :: quaternion
real(pReal) :: w = 0.0_pReal
real(pReal) :: x = 0.0_pReal
real(pReal) :: y = 0.0_pReal
real(pReal) :: z = 0.0_pReal
contains
procedure, private :: add__
procedure, private :: pos__
generic, public :: operator(+) => add__,pos__
procedure, private :: sub__
procedure, private :: neg__
generic, public :: operator(-) => sub__,neg__
procedure, private :: mul_quat__
procedure, private :: mul_scal__
generic, public :: operator(*) => mul_quat__, mul_scal__
procedure, private :: div_quat__
procedure, private :: div_scal__
generic, public :: operator(/) => div_quat__, div_scal__
procedure, private :: eq__
generic, public :: operator(==) => eq__
procedure, private :: neq__
generic, public :: operator(/=) => neq__
procedure, private :: pow_quat__
procedure, private :: pow_scal__
generic, public :: operator(**) => pow_quat__, pow_scal__
procedure, private :: abs__
procedure, private :: dot_product__
procedure, private :: conjg__
procedure, private :: exp__
procedure, private :: log__
procedure, public :: homomorphed => quat_homomorphed
!procedure,private :: quat_write
!generic :: write(formatted) => quat_write
end type
interface assignment (=)
module procedure assign_quat__
module procedure assign_vec__
end interface assignment (=)
interface quaternion
module procedure init__
end interface quaternion
interface abs
procedure abs__
end interface abs
interface dot_product
procedure dot_product__
end interface dot_product
interface conjg
module procedure conjg__
end interface conjg
interface exp
module procedure exp__
end interface exp
interface log
module procedure log__
end interface log
contains
!--------------------------------------------------------------------------
!> constructor for a quaternion from a 4-vector
!--------------------------------------------------------------------------
type(quaternion) pure function init__(array)
implicit none
real(pReal), intent(in), dimension(4) :: array
init__%w=array(1)
init__%x=array(2)
init__%y=array(3)
init__%z=array(4)
end function init__
!--------------------------------------------------------------------------
!> assing a quaternion
!--------------------------------------------------------------------------
elemental subroutine assign_quat__(self,other)
implicit none
type(quaternion), intent(out) :: self
type(quaternion), intent(in) :: other
self%w = other%w
self%x = other%x
self%y = other%y
self%z = other%z
end subroutine assign_quat__
!--------------------------------------------------------------------------
!> assing a 4-vector
!--------------------------------------------------------------------------
pure subroutine assign_vec__(self,other)
implicit none
type(quaternion), intent(out) :: self
real(pReal), intent(in), dimension(4) :: other
self%w = other(1)
self%x = other(2)
self%y = other(3)
self%z = other(4)
end subroutine assign_vec__
!--------------------------------------------------------------------------
!> addition of two quaternions
!--------------------------------------------------------------------------
type(quaternion) elemental function add__(self,other)
implicit none
class(quaternion), intent(in) :: self,other
add__%w = self%w + other%w
add__%x = self%x + other%x
add__%y = self%y + other%y
add__%z = self%z + other%z
end function add__
!--------------------------------------------------------------------------
!> unary positive operator
!--------------------------------------------------------------------------
type(quaternion) elemental function pos__(self)
implicit none
class(quaternion), intent(in) :: self
pos__%w = self%w
pos__%x = self%x
pos__%y = self%y
pos__%z = self%z
end function pos__
!--------------------------------------------------------------------------
!> subtraction of two quaternions
!--------------------------------------------------------------------------
type(quaternion) elemental function sub__(self,other)
implicit none
class(quaternion), intent(in) :: self,other
sub__%w = self%w - other%w
sub__%x = self%x - other%x
sub__%y = self%y - other%y
sub__%z = self%z - other%z
end function sub__
!--------------------------------------------------------------------------
!> unary positive operator
!--------------------------------------------------------------------------
type(quaternion) elemental function neg__(self)
implicit none
class(quaternion), intent(in) :: self
neg__%w = -self%w
neg__%x = -self%x
neg__%y = -self%y
neg__%z = -self%z
end function neg__
!--------------------------------------------------------------------------
!> multiplication of two quaternions
!--------------------------------------------------------------------------
type(quaternion) elemental function mul_quat__(self,other)
implicit none
class(quaternion), intent(in) :: self, other
mul_quat__%w = self%w*other%w - self%x*other%x - self%y*other%y - self%z*other%z
mul_quat__%x = self%w*other%x + self%x*other%w + epsijk * (self%y*other%z - self%z*other%y)
mul_quat__%y = self%w*other%y + self%y*other%w + epsijk * (self%z*other%x - self%x*other%z)
mul_quat__%z = self%w*other%z + self%z*other%w + epsijk * (self%x*other%y - self%y*other%x)
end function mul_quat__
!--------------------------------------------------------------------------
!> multiplication of quaternions with scalar
!--------------------------------------------------------------------------
type(quaternion) elemental function mul_scal__(self,scal)
implicit none
class(quaternion), intent(in) :: self
real(pReal), intent(in) :: scal
mul_scal__%w = self%w*scal
mul_scal__%x = self%x*scal
mul_scal__%y = self%y*scal
mul_scal__%z = self%z*scal
end function mul_scal__
!--------------------------------------------------------------------------
!> division of two quaternions
!--------------------------------------------------------------------------
type(quaternion) elemental function div_quat__(self,other)
implicit none
class(quaternion), intent(in) :: self, other
div_quat__ = self * (conjg(other)/(abs(other)**2.0_pReal))
end function div_quat__
!--------------------------------------------------------------------------
!> divisiont of quaternions by scalar
!--------------------------------------------------------------------------
type(quaternion) elemental function div_scal__(self,scal)
implicit none
class(quaternion), intent(in) :: self
real(pReal), intent(in) :: scal
div_scal__ = [self%w,self%x,self%y,self%z]/scal
end function div_scal__
!--------------------------------------------------------------------------
!> equality of two quaternions
!--------------------------------------------------------------------------
logical elemental function eq__(self,other)
implicit none
class(quaternion), intent(in) :: self,other
eq__ = all(dEq([ self%w, self%x, self%y, self%z], &
[other%w,other%x,other%y,other%z]))
end function eq__
!--------------------------------------------------------------------------
!> inequality of two quaternions
!--------------------------------------------------------------------------
logical elemental function neq__(self,other)
implicit none
class(quaternion), intent(in) :: self,other
neq__ = .not. self%eq__(other)
end function neq__
!--------------------------------------------------------------------------
!> quaternion to the power of a scalar
!--------------------------------------------------------------------------
type(quaternion) elemental function pow_scal__(self,expon)
implicit none
class(quaternion), intent(in) :: self
real(pReal), intent(in) :: expon
pow_scal__ = exp(log(self)*expon)
end function pow_scal__
!--------------------------------------------------------------------------
!> quaternion to the power of a quaternion
!--------------------------------------------------------------------------
type(quaternion) elemental function pow_quat__(self,expon)
implicit none
class(quaternion), intent(in) :: self
type(quaternion), intent(in) :: expon
pow_quat__ = exp(log(self)*expon)
end function pow_quat__
!--------------------------------------------------------------------------
!> exponential of a quaternion
!> ToDo: Lacks any check for invalid operations
!--------------------------------------------------------------------------
type(quaternion) elemental function exp__(self)
implicit none
class(quaternion), intent(in) :: self
real(pReal) :: absImag
absImag = norm2([self%x, self%y, self%z])
exp__ = exp(self%w) * [ cos(absImag), &
self%x/absImag * sin(absImag), &
self%y/absImag * sin(absImag), &
self%z/absImag * sin(absImag)]
end function exp__
!--------------------------------------------------------------------------
!> logarithm of a quaternion
!> ToDo: Lacks any check for invalid operations
!--------------------------------------------------------------------------
type(quaternion) elemental function log__(self)
implicit none
class(quaternion), intent(in) :: self
real(pReal) :: absImag
absImag = norm2([self%x, self%y, self%z])
log__ = [log(abs(self)), &
self%x/absImag * acos(self%w/abs(self)), &
self%y/absImag * acos(self%w/abs(self)), &
self%z/absImag * acos(self%w/abs(self))]
end function log__
!--------------------------------------------------------------------------
!> norm of a quaternion
!--------------------------------------------------------------------------
real(pReal) elemental function abs__(a)
implicit none
class(quaternion), intent(in) :: a
abs__ = norm2([a%w,a%x,a%y,a%z])
end function abs__
!--------------------------------------------------------------------------
!> dot product of two quaternions
!--------------------------------------------------------------------------
real(pReal) elemental function dot_product__(a,b)
implicit none
class(quaternion), intent(in) :: a,b
dot_product__ = a%w*b%w + a%x*b%x + a%y*b%y + a%z*b%z
end function dot_product__
!--------------------------------------------------------------------------
!> conjugate complex of a quaternion
!--------------------------------------------------------------------------
type(quaternion) elemental function conjg__(a)
implicit none
class(quaternion), intent(in) :: a
conjg__ = quaternion([a%w, -a%x, -a%y, -a%z])
end function conjg__
!--------------------------------------------------------------------------
!> homomorphed quaternion of a quaternion
!--------------------------------------------------------------------------
type(quaternion) elemental function quat_homomorphed(a)
implicit none
class(quaternion), intent(in) :: a
quat_homomorphed = quaternion(-[a%w,a%x,a%y,a%z])
end function quat_homomorphed
end module quaternions

1088
src/rotations.f90 Normal file

File diff suppressed because it is too large Load Diff