Merge branch 'no-crystallite_subF' into no-subXXX
This commit is contained in:
Martin Diehl
commit
942bd8932e
|
|
@ -21,6 +21,9 @@ class NiceDumper(yaml.SafeDumper):
|
|||
return self.represent_data(dict(data)) if isinstance(data, dict) and type(data) != dict else \
|
||||
super().represent_data(data)
|
||||
|
||||
def ignore_aliases(self, data):
|
||||
"""No references."""
|
||||
return True
|
||||
|
||||
class Config(dict):
|
||||
"""YAML-based configuration."""
|
||||
|
|
|
|||
|
|
@ -9,6 +9,16 @@ from . import Orientation
|
|||
class ConfigMaterial(Config):
|
||||
"""Material configuration."""
|
||||
|
||||
_defaults = {'material': [],
|
||||
'homogenization': {},
|
||||
'phase': {}}
|
||||
|
||||
def __init__(self,d={}):
|
||||
"""Initialize object with default dictionary keys."""
|
||||
super().__init__(d)
|
||||
for k,v in self._defaults.items():
|
||||
if k not in self: self[k] = v
|
||||
|
||||
def save(self,fname='material.yaml',**kwargs):
|
||||
"""
|
||||
Save to yaml file.
|
||||
|
|
@ -75,6 +85,8 @@ class ConfigMaterial(Config):
|
|||
fraction: 1.0
|
||||
phase: Steel
|
||||
homogenization: SX
|
||||
homogenization: {}
|
||||
phase: {}
|
||||
|
||||
"""
|
||||
constituents_ = {k:table.get(v) for k,v in constituents.items()}
|
||||
|
|
@ -261,6 +273,8 @@ class ConfigMaterial(Config):
|
|||
fraction: 1.0
|
||||
phase: Aluminum
|
||||
homogenization: SX
|
||||
homogenization: {}
|
||||
phase: {}
|
||||
|
||||
"""
|
||||
length = -1
|
||||
|
|
@ -274,6 +288,7 @@ class ConfigMaterial(Config):
|
|||
|
||||
c = [{} for _ in range(length)] if constituents is None else \
|
||||
[{'constituents':u} for u in ConfigMaterial._constituents(**constituents)]
|
||||
|
||||
if len(c) == 1: c = [copy.deepcopy(c[0]) for _ in range(length)]
|
||||
|
||||
if length != 1 and length != len(c):
|
||||
|
|
|
|||
|
|
@ -226,9 +226,9 @@ class Orientation(Rotation):
|
|||
|
||||
"""
|
||||
return super().__eq__(other) \
|
||||
and self.family == other.family \
|
||||
and self.lattice == other.lattice \
|
||||
and self.parameters == other.parameters
|
||||
and hasattr(other, 'family') and self.family == other.family \
|
||||
and hasattr(other, 'lattice') and self.lattice == other.lattice \
|
||||
and hasattr(other, 'parameters') and self.parameters == other.parameters
|
||||
|
||||
|
||||
def __matmul__(self,other):
|
||||
|
|
|
|||
|
|
@ -144,6 +144,11 @@ class Rotation:
|
|||
return self.copy(rotation=Rotation(np.block([np.cos(pwr*phi),np.sin(pwr*phi)*p]))._standardize())
|
||||
|
||||
|
||||
def __mul__(self,other):
|
||||
"""Standard multiplication is not implemented."""
|
||||
raise NotImplementedError('Use "R@b", i.e. matmul, to apply rotation "R" to object "b"')
|
||||
|
||||
|
||||
def __matmul__(self,other):
|
||||
"""
|
||||
Rotation of vector, second or fourth order tensor, or rotation object.
|
||||
|
|
@ -199,8 +204,16 @@ class Rotation:
|
|||
|
||||
|
||||
def append(self,other):
|
||||
"""Extend rotation array along first dimension with other array."""
|
||||
return self.copy(rotation=np.vstack((self.quaternion,other.quaternion)))
|
||||
"""
|
||||
Extend rotation array along first dimension with other array(s).
|
||||
|
||||
Parameters
|
||||
----------
|
||||
other : Rotation or list of Rotations.
|
||||
|
||||
"""
|
||||
return self.copy(rotation=np.vstack(tuple(map(lambda x:x.quaternion,
|
||||
[self]+other if isinstance(other,list) else [self,other]))))
|
||||
|
||||
|
||||
def flatten(self,order = 'C'):
|
||||
|
|
@ -258,7 +271,7 @@ class Rotation:
|
|||
"""Intermediate representation supporting quaternion averaging."""
|
||||
return np.einsum('...i,...j',quat,quat)
|
||||
|
||||
if not weights:
|
||||
if weights is None:
|
||||
weights = np.ones(self.shape,dtype=float)
|
||||
|
||||
eig, vec = np.linalg.eig(np.sum(_M(self.quaternion) * weights[...,np.newaxis,np.newaxis],axis=-3) \
|
||||
|
|
|
|||
|
|
@ -1,10 +1,10 @@
|
|||
homogenization:
|
||||
SX:
|
||||
N_constituents: 2
|
||||
mech: {type: none}
|
||||
N_constituents: 1
|
||||
mechanics: {type: none}
|
||||
Taylor:
|
||||
N_constituents: 2
|
||||
mech: {type: isostrain}
|
||||
mechanics: {type: isostrain}
|
||||
|
||||
material:
|
||||
- constituents:
|
||||
|
|
@ -34,11 +34,11 @@ material:
|
|||
phase:
|
||||
Aluminum:
|
||||
lattice: cF
|
||||
mech:
|
||||
mechanics:
|
||||
output: [F, P, F_e, F_p, L_p]
|
||||
elasticity: {C_11: 106.75e9, C_12: 60.41e9, C_44: 28.34e9, type: hooke}
|
||||
Steel:
|
||||
lattice: cI
|
||||
mech:
|
||||
mechanics:
|
||||
output: [F, P, F_e, F_p, L_p]
|
||||
elasticity: {C_11: 233.3e9, C_12: 135.5e9, C_44: 118.0e9, type: hooke}
|
||||
|
|
|
|||
|
|
@ -800,6 +800,14 @@ class TestRotation:
|
|||
print(f'append 2x {shape} --> {s.shape}')
|
||||
assert s[0,...] == r[0,...] and s[-1,...] == p[-1,...]
|
||||
|
||||
@pytest.mark.parametrize('shape',[None,1,(1,),(4,2),(3,3,2)])
|
||||
def test_append_list(self,shape):
|
||||
r = Rotation.from_random(shape=shape)
|
||||
p = Rotation.from_random(shape=shape)
|
||||
s = r.append([r,p])
|
||||
print(f'append 3x {shape} --> {s.shape}')
|
||||
assert s[0,...] == r[0,...] and s[-1,...] == p[-1,...]
|
||||
|
||||
@pytest.mark.parametrize('quat,standardized',[
|
||||
([-1,0,0,0],[1,0,0,0]),
|
||||
([-0.5,-0.5,-0.5,-0.5],[0.5,0.5,0.5,0.5]),
|
||||
|
|
|
|||
|
|
@ -11,7 +11,6 @@
|
|||
#include "config.f90"
|
||||
#include "LAPACK_interface.f90"
|
||||
#include "math.f90"
|
||||
#include "quaternions.f90"
|
||||
#include "rotations.f90"
|
||||
#include "element.f90"
|
||||
#include "HDF5_utilities.f90"
|
||||
|
|
|
|||
|
|
@ -840,10 +840,11 @@ subroutine crystallite_init
|
|||
co, & !< counter in integration point component loop
|
||||
ip, & !< counter in integration point loop
|
||||
el, & !< counter in element loop
|
||||
so, &
|
||||
cMax, & !< maximum number of integration point components
|
||||
iMax, & !< maximum number of integration points
|
||||
eMax, & !< maximum number of elements
|
||||
so
|
||||
eMax !< maximum number of elements
|
||||
|
||||
|
||||
class(tNode), pointer :: &
|
||||
num_crystallite, &
|
||||
|
|
@ -875,7 +876,6 @@ subroutine crystallite_init
|
|||
|
||||
allocate(crystallite_orientation(cMax,iMax,eMax))
|
||||
|
||||
|
||||
num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)
|
||||
|
||||
num%subStepMinCryst = num_crystallite%get_asFloat ('subStepMin', defaultVal=1.0e-3_pReal)
|
||||
|
|
|
|||
|
|
@ -1015,6 +1015,7 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul
|
|||
|
||||
enddo iteration
|
||||
|
||||
|
||||
contains
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
|
|||
|
|
@ -1,534 +0,0 @@
|
|||
!---------------------------------------------------------------------------------------------------
|
||||
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
|
||||
!> @author Philip Eisenlohr, Michigan State University
|
||||
!> @brief general quaternion math, not limited to unit quaternions
|
||||
!> @details w is the real part, (x, y, z) are the imaginary parts.
|
||||
!> @details https://en.wikipedia.org/wiki/Quaternion
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
module quaternions
|
||||
use prec
|
||||
|
||||
implicit none
|
||||
private
|
||||
|
||||
real(pReal), parameter, public :: P = -1.0_pReal !< parameter for orientation conversion.
|
||||
|
||||
type, public :: quaternion
|
||||
real(pReal), private :: w = 0.0_pReal
|
||||
real(pReal), private :: x = 0.0_pReal
|
||||
real(pReal), private :: y = 0.0_pReal
|
||||
real(pReal), private :: 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, public :: abs => abs__
|
||||
procedure, public :: conjg => conjg__
|
||||
procedure, public :: real => real__
|
||||
procedure, public :: aimag => aimag__
|
||||
|
||||
procedure, public :: homomorphed
|
||||
procedure, public :: asArray
|
||||
procedure, public :: inverse
|
||||
|
||||
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
|
||||
|
||||
interface real
|
||||
module procedure real__
|
||||
end interface real
|
||||
|
||||
interface aimag
|
||||
module procedure aimag__
|
||||
end interface aimag
|
||||
|
||||
public :: &
|
||||
quaternions_init, &
|
||||
assignment(=), &
|
||||
conjg, aimag, &
|
||||
log, exp, &
|
||||
abs, dot_product, &
|
||||
inverse, &
|
||||
real
|
||||
|
||||
contains
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Do self test.
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine quaternions_init
|
||||
|
||||
print'(/,a)', ' <<<+- quaternions init -+>>>'; flush(6)
|
||||
|
||||
call selfTest
|
||||
|
||||
end subroutine quaternions_init
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief construct a quaternion from a 4-vector
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
type(quaternion) pure function init__(array)
|
||||
|
||||
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__
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief assign a quaternion
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
elemental pure subroutine assign_quat__(self,other)
|
||||
|
||||
type(quaternion), intent(out) :: self
|
||||
type(quaternion), intent(in) :: other
|
||||
|
||||
self = [other%w,other%x,other%y,other%z]
|
||||
|
||||
end subroutine assign_quat__
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief assign a 4-vector
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
pure subroutine assign_vec__(self,other)
|
||||
|
||||
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__
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief add a quaternion
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
type(quaternion) elemental pure function add__(self,other)
|
||||
|
||||
class(quaternion), intent(in) :: self,other
|
||||
|
||||
add__ = [ self%w, self%x, self%y ,self%z] &
|
||||
+ [other%w, other%x, other%y,other%z]
|
||||
|
||||
end function add__
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief return (unary positive operator)
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
type(quaternion) elemental pure function pos__(self)
|
||||
|
||||
class(quaternion), intent(in) :: self
|
||||
|
||||
pos__ = self * (+1.0_pReal)
|
||||
|
||||
end function pos__
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief subtract a quaternion
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
type(quaternion) elemental pure function sub__(self,other)
|
||||
|
||||
class(quaternion), intent(in) :: self,other
|
||||
|
||||
sub__ = [ self%w, self%x, self%y ,self%z] &
|
||||
- [other%w, other%x, other%y,other%z]
|
||||
|
||||
end function sub__
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief negate (unary negative operator)
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
type(quaternion) elemental pure function neg__(self)
|
||||
|
||||
class(quaternion), intent(in) :: self
|
||||
|
||||
neg__ = self * (-1.0_pReal)
|
||||
|
||||
end function neg__
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief multiply with a quaternion
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
type(quaternion) elemental pure function mul_quat__(self,other)
|
||||
|
||||
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 + P * (self%y*other%z - self%z*other%y)
|
||||
mul_quat__%y = self%w*other%y + self%y*other%w + P * (self%z*other%x - self%x*other%z)
|
||||
mul_quat__%z = self%w*other%z + self%z*other%w + P * (self%x*other%y - self%y*other%x)
|
||||
|
||||
end function mul_quat__
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief multiply with a scalar
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
type(quaternion) elemental pure function mul_scal__(self,scal)
|
||||
|
||||
class(quaternion), intent(in) :: self
|
||||
real(pReal), intent(in) :: scal
|
||||
|
||||
mul_scal__ = [self%w,self%x,self%y,self%z]*scal
|
||||
|
||||
end function mul_scal__
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief divide by a quaternion
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
type(quaternion) elemental pure function div_quat__(self,other)
|
||||
|
||||
class(quaternion), intent(in) :: self, other
|
||||
|
||||
div_quat__ = self * (conjg(other)/(abs(other)**2.0_pReal))
|
||||
|
||||
end function div_quat__
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief divide by a scalar
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
type(quaternion) elemental pure function div_scal__(self,scal)
|
||||
|
||||
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__
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief test equality
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
logical elemental pure function eq__(self,other)
|
||||
|
||||
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__
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief test inequality
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
logical elemental pure function neq__(self,other)
|
||||
|
||||
class(quaternion), intent(in) :: self,other
|
||||
|
||||
neq__ = .not. self%eq__(other)
|
||||
|
||||
end function neq__
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief raise to the power of a quaternion
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
type(quaternion) elemental pure function pow_quat__(self,expon)
|
||||
|
||||
class(quaternion), intent(in) :: self
|
||||
type(quaternion), intent(in) :: expon
|
||||
|
||||
pow_quat__ = exp(log(self)*expon)
|
||||
|
||||
end function pow_quat__
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief raise to the power of a scalar
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
type(quaternion) elemental pure function pow_scal__(self,expon)
|
||||
|
||||
class(quaternion), intent(in) :: self
|
||||
real(pReal), intent(in) :: expon
|
||||
|
||||
pow_scal__ = exp(log(self)*expon)
|
||||
|
||||
end function pow_scal__
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief take exponential
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
type(quaternion) elemental pure function exp__(a)
|
||||
|
||||
class(quaternion), intent(in) :: a
|
||||
real(pReal) :: absImag
|
||||
|
||||
absImag = norm2(aimag(a))
|
||||
|
||||
exp__ = merge(exp(a%w) * [ cos(absImag), &
|
||||
a%x/absImag * sin(absImag), &
|
||||
a%y/absImag * sin(absImag), &
|
||||
a%z/absImag * sin(absImag)], &
|
||||
IEEE_value(1.0_pReal,IEEE_SIGNALING_NAN), &
|
||||
dNeq0(absImag))
|
||||
|
||||
end function exp__
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief take logarithm
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
type(quaternion) elemental pure function log__(a)
|
||||
|
||||
class(quaternion), intent(in) :: a
|
||||
real(pReal) :: absImag
|
||||
|
||||
absImag = norm2(aimag(a))
|
||||
|
||||
log__ = merge([log(abs(a)), &
|
||||
a%x/absImag * acos(a%w/abs(a)), &
|
||||
a%y/absImag * acos(a%w/abs(a)), &
|
||||
a%z/absImag * acos(a%w/abs(a))], &
|
||||
IEEE_value(1.0_pReal,IEEE_SIGNALING_NAN), &
|
||||
dNeq0(absImag))
|
||||
|
||||
end function log__
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief return norm
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
real(pReal) elemental pure function abs__(self)
|
||||
|
||||
class(quaternion), intent(in) :: self
|
||||
|
||||
abs__ = norm2([self%w,self%x,self%y,self%z])
|
||||
|
||||
end function abs__
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief calculate dot product
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
real(pReal) elemental pure function dot_product__(a,b)
|
||||
|
||||
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__
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief take conjugate complex
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
type(quaternion) elemental pure function conjg__(self)
|
||||
|
||||
class(quaternion), intent(in) :: self
|
||||
|
||||
conjg__ = [self%w,-self%x,-self%y,-self%z]
|
||||
|
||||
end function conjg__
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief homomorph
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
type(quaternion) elemental pure function homomorphed(self)
|
||||
|
||||
class(quaternion), intent(in) :: self
|
||||
|
||||
homomorphed = - self
|
||||
|
||||
end function homomorphed
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief return as plain array
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
pure function asArray(self)
|
||||
|
||||
real(pReal), dimension(4) :: asArray
|
||||
class(quaternion), intent(in) :: self
|
||||
|
||||
asArray = [self%w,self%x,self%y,self%z]
|
||||
|
||||
end function asArray
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief real part (scalar)
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
pure function real__(self)
|
||||
|
||||
real(pReal) :: real__
|
||||
class(quaternion), intent(in) :: self
|
||||
|
||||
real__ = self%w
|
||||
|
||||
end function real__
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief imaginary part (3-vector)
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
pure function aimag__(self)
|
||||
|
||||
real(pReal), dimension(3) :: aimag__
|
||||
class(quaternion), intent(in) :: self
|
||||
|
||||
aimag__ = [self%x,self%y,self%z]
|
||||
|
||||
end function aimag__
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief inverse
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
type(quaternion) elemental pure function inverse(self)
|
||||
|
||||
class(quaternion), intent(in) :: self
|
||||
|
||||
inverse = conjg(self)/abs(self)**2.0_pReal
|
||||
|
||||
end function inverse
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief check correctness of some quaternions functions
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine selfTest
|
||||
|
||||
real(pReal), dimension(4) :: qu
|
||||
type(quaternion) :: q, q_2
|
||||
|
||||
if(dNeq(abs(P),1.0_pReal)) error stop 'P not in {-1,+1}'
|
||||
|
||||
call random_number(qu)
|
||||
qu = (qu-0.5_pReal) * 2.0_pReal
|
||||
q = quaternion(qu)
|
||||
|
||||
q_2= qu
|
||||
if(any(dNeq(q%asArray(),q_2%asArray()))) error stop 'assign_vec__'
|
||||
|
||||
q_2 = q + q
|
||||
if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) error stop 'add__'
|
||||
|
||||
q_2 = q - q
|
||||
if(any(dNeq0(q_2%asArray()))) error stop 'sub__'
|
||||
|
||||
q_2 = q * 5.0_pReal
|
||||
if(any(dNeq(q_2%asArray(),5.0_pReal*qu))) error stop 'mul__'
|
||||
|
||||
q_2 = q / 0.5_pReal
|
||||
if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) error stop 'div__'
|
||||
|
||||
q_2 = q * 0.3_pReal
|
||||
if(dNeq0(abs(q)) .and. q_2 == q) error stop 'eq__'
|
||||
|
||||
q_2 = q
|
||||
if(q_2 /= q) error stop 'neq__'
|
||||
|
||||
if(dNeq(abs(q),norm2(qu))) error stop 'abs__'
|
||||
if(dNeq(abs(q)**2.0_pReal, real(q*q%conjg()),1.0e-14_pReal)) &
|
||||
error stop 'abs__/*conjg'
|
||||
|
||||
if(any(dNeq(q%asArray(),qu))) error stop 'eq__'
|
||||
if(dNeq(q%real(), qu(1))) error stop 'real()'
|
||||
if(any(dNeq(q%aimag(), qu(2:4)))) error stop 'aimag()'
|
||||
|
||||
q_2 = q%homomorphed()
|
||||
if(q /= q_2* (-1.0_pReal)) error stop 'homomorphed'
|
||||
if(dNeq(q_2%real(), qu(1)* (-1.0_pReal))) error stop 'homomorphed/real'
|
||||
if(any(dNeq(q_2%aimag(),qu(2:4)*(-1.0_pReal)))) error stop 'homomorphed/aimag'
|
||||
|
||||
q_2 = conjg(q)
|
||||
if(dNeq(abs(q),abs(q_2))) error stop 'conjg/abs'
|
||||
if(q /= conjg(q_2)) error stop 'conjg/involution'
|
||||
if(dNeq(q_2%real(), q%real())) error stop 'conjg/real'
|
||||
if(any(dNeq(q_2%aimag(),q%aimag()*(-1.0_pReal)))) error stop 'conjg/aimag'
|
||||
|
||||
if(abs(q) > 0.0_pReal) then
|
||||
q_2 = q * q%inverse()
|
||||
if( dNeq(real(q_2), 1.0_pReal,1.0e-15_pReal)) error stop 'inverse/real'
|
||||
if(any(dNeq0(aimag(q_2), 1.0e-15_pReal))) error stop 'inverse/aimag'
|
||||
|
||||
q_2 = q/abs(q)
|
||||
q_2 = conjg(q_2) - inverse(q_2)
|
||||
if(any(dNeq0(q_2%asArray(),1.0e-15_pReal))) error stop 'inverse/conjg'
|
||||
endif
|
||||
if(dNeq(dot_product(qu,qu),dot_product(q,q))) error stop 'dot_product'
|
||||
|
||||
#if !(defined(__GFORTRAN__) && __GNUC__ < 9)
|
||||
if (norm2(aimag(q)) > 0.0_pReal) then
|
||||
if (dNeq0(abs(q-exp(log(q))),1.0e-13_pReal)) error stop 'exp/log'
|
||||
if (dNeq0(abs(q-log(exp(q))),1.0e-13_pReal)) error stop 'log/exp'
|
||||
endif
|
||||
#endif
|
||||
|
||||
end subroutine selfTest
|
||||
|
||||
|
||||
end module quaternions
|
||||
|
|
@ -47,16 +47,16 @@
|
|||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
module rotations
|
||||
use prec
|
||||
use IO
|
||||
use math
|
||||
use quaternions
|
||||
|
||||
implicit none
|
||||
private
|
||||
|
||||
real(pReal), parameter :: P = -1.0_pReal !< parameter for orientation conversion.
|
||||
|
||||
type, public :: rotation
|
||||
type(quaternion) :: q
|
||||
real(pReal), dimension(4) :: q
|
||||
contains
|
||||
procedure, public :: asQuaternion
|
||||
procedure, public :: asEulers
|
||||
|
|
@ -103,7 +103,6 @@ contains
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine rotations_init
|
||||
|
||||
call quaternions_init
|
||||
print'(/,a)', ' <<<+- rotations init -+>>>'; flush(IO_STDOUT)
|
||||
|
||||
print*, 'Rowenhorst et al., Modelling and Simulation in Materials Science and Engineering 23:083501, 2015'
|
||||
|
|
@ -122,7 +121,7 @@ pure function asQuaternion(self)
|
|||
class(rotation), intent(in) :: self
|
||||
real(pReal), dimension(4) :: asQuaternion
|
||||
|
||||
asQuaternion = self%q%asArray()
|
||||
asQuaternion = self%q
|
||||
|
||||
end function asQuaternion
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
|
@ -131,7 +130,7 @@ pure function asEulers(self)
|
|||
class(rotation), intent(in) :: self
|
||||
real(pReal), dimension(3) :: asEulers
|
||||
|
||||
asEulers = qu2eu(self%q%asArray())
|
||||
asEulers = qu2eu(self%q)
|
||||
|
||||
end function asEulers
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
|
@ -140,7 +139,7 @@ pure function asAxisAngle(self)
|
|||
class(rotation), intent(in) :: self
|
||||
real(pReal), dimension(4) :: asAxisAngle
|
||||
|
||||
asAxisAngle = qu2ax(self%q%asArray())
|
||||
asAxisAngle = qu2ax(self%q)
|
||||
|
||||
end function asAxisAngle
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
|
@ -149,7 +148,7 @@ pure function asMatrix(self)
|
|||
class(rotation), intent(in) :: self
|
||||
real(pReal), dimension(3,3) :: asMatrix
|
||||
|
||||
asMatrix = qu2om(self%q%asArray())
|
||||
asMatrix = qu2om(self%q)
|
||||
|
||||
end function asMatrix
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
|
@ -158,7 +157,7 @@ pure function asRodrigues(self)
|
|||
class(rotation), intent(in) :: self
|
||||
real(pReal), dimension(4) :: asRodrigues
|
||||
|
||||
asRodrigues = qu2ro(self%q%asArray())
|
||||
asRodrigues = qu2ro(self%q)
|
||||
|
||||
end function asRodrigues
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
|
@ -167,7 +166,7 @@ pure function asHomochoric(self)
|
|||
class(rotation), intent(in) :: self
|
||||
real(pReal), dimension(3) :: asHomochoric
|
||||
|
||||
asHomochoric = qu2ho(self%q%asArray())
|
||||
asHomochoric = qu2ho(self%q)
|
||||
|
||||
end function asHomochoric
|
||||
|
||||
|
|
@ -259,7 +258,7 @@ pure elemental function rotRot__(self,R) result(rRot)
|
|||
type(rotation) :: rRot
|
||||
class(rotation), intent(in) :: self,R
|
||||
|
||||
rRot = rotation(self%q*R%q)
|
||||
rRot = rotation(multiply_quaternion(self%q,R%q))
|
||||
call rRot%standardize()
|
||||
|
||||
end function rotRot__
|
||||
|
|
@ -272,14 +271,14 @@ pure elemental subroutine standardize(self)
|
|||
|
||||
class(rotation), intent(inout) :: self
|
||||
|
||||
if (real(self%q) < 0.0_pReal) self%q = self%q%homomorphed()
|
||||
if (self%q(1) < 0.0_pReal) self%q = - self%q
|
||||
|
||||
end subroutine standardize
|
||||
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @author Marc De Graef, Carnegie Mellon University
|
||||
!> @brief rotate a vector passively (default) or actively
|
||||
!> @brief Rotate a vector passively (default) or actively.
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
pure function rotVector(self,v,active) result(vRot)
|
||||
|
||||
|
|
@ -288,9 +287,8 @@ pure function rotVector(self,v,active) result(vRot)
|
|||
real(pReal), intent(in), dimension(3) :: v
|
||||
logical, intent(in), optional :: active
|
||||
|
||||
real(pReal), dimension(3) :: v_normed
|
||||
type(quaternion) :: q
|
||||
logical :: passive
|
||||
real(pReal), dimension(4) :: v_normed, q
|
||||
logical :: passive
|
||||
|
||||
if (present(active)) then
|
||||
passive = .not. active
|
||||
|
|
@ -301,13 +299,13 @@ pure function rotVector(self,v,active) result(vRot)
|
|||
if (dEq0(norm2(v))) then
|
||||
vRot = v
|
||||
else
|
||||
v_normed = v/norm2(v)
|
||||
v_normed = [0.0_pReal,v]/norm2(v)
|
||||
if (passive) then
|
||||
q = self%q * (quaternion([0.0_pReal, v_normed(1), v_normed(2), v_normed(3)]) * conjg(self%q) )
|
||||
q = multiply_quaternion(self%q, multiply_quaternion(v_normed, conjugate_quaternion(self%q)))
|
||||
else
|
||||
q = conjg(self%q) * (quaternion([0.0_pReal, v_normed(1), v_normed(2), v_normed(3)]) * self%q )
|
||||
q = multiply_quaternion(conjugate_quaternion(self%q), multiply_quaternion(v_normed, self%q))
|
||||
endif
|
||||
vRot = q%aimag()*norm2(v)
|
||||
vRot = q(2:4)*norm2(v)
|
||||
endif
|
||||
|
||||
end function rotVector
|
||||
|
|
@ -315,8 +313,8 @@ end function rotVector
|
|||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!> @author Marc De Graef, Carnegie Mellon University
|
||||
!> @brief rotate a rank-2 tensor passively (default) or actively
|
||||
!> @details: rotation is based on rotation matrix
|
||||
!> @brief Rotate a rank-2 tensor passively (default) or actively.
|
||||
!> @details: Rotation is based on rotation matrix
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
pure function rotTensor2(self,T,active) result(tRot)
|
||||
|
||||
|
|
@ -403,7 +401,7 @@ pure elemental function misorientation(self,other)
|
|||
type(rotation) :: misorientation
|
||||
class(rotation), intent(in) :: self, other
|
||||
|
||||
misorientation%q = other%q * conjg(self%q)
|
||||
misorientation%q = multiply_quaternion(other%q, conjugate_quaternion(self%q))
|
||||
|
||||
end function misorientation
|
||||
|
||||
|
|
@ -1338,7 +1336,7 @@ end function cu2ho
|
|||
!--------------------------------------------------------------------------
|
||||
!> @author Marc De Graef, Carnegie Mellon University
|
||||
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
|
||||
!> @brief determine to which pyramid a point in a cubic grid belongs
|
||||
!> @brief Determine to which pyramid a point in a cubic grid belongs.
|
||||
!--------------------------------------------------------------------------
|
||||
pure function GetPyramidOrder(xyz)
|
||||
|
||||
|
|
@ -1362,7 +1360,39 @@ end function GetPyramidOrder
|
|||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief check correctness of some rotations functions
|
||||
!> @brief Multiply two quaternions.
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
pure function multiply_quaternion(qu1,qu2)
|
||||
|
||||
real(pReal), dimension(4), intent(in) :: qu1, qu2
|
||||
real(pReal), dimension(4) :: multiply_quaternion
|
||||
|
||||
|
||||
multiply_quaternion(1) = qu1(1)*qu2(1) - qu1(2)*qu2(2) - qu1(3)*qu2(3) - qu1(4)*qu2(4)
|
||||
multiply_quaternion(2) = qu1(1)*qu2(2) + qu1(2)*qu2(1) + P * (qu1(3)*qu2(4) - qu1(4)*qu2(3))
|
||||
multiply_quaternion(3) = qu1(1)*qu2(3) + qu1(3)*qu2(1) + P * (qu1(4)*qu2(2) - qu1(2)*qu2(4))
|
||||
multiply_quaternion(4) = qu1(1)*qu2(4) + qu1(4)*qu2(1) + P * (qu1(2)*qu2(3) - qu1(3)*qu2(2))
|
||||
|
||||
end function multiply_quaternion
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Calculate conjugate complex of a quaternion.
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
pure function conjugate_quaternion(qu)
|
||||
|
||||
real(pReal), dimension(4), intent(in) :: qu
|
||||
real(pReal), dimension(4) :: conjugate_quaternion
|
||||
|
||||
|
||||
conjugate_quaternion = [qu(1), -qu(2), -qu(3), -qu(4)]
|
||||
|
||||
|
||||
end function conjugate_quaternion
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Check correctness of some rotations functions.
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine selfTest
|
||||
|
||||
|
|
@ -1374,7 +1404,8 @@ subroutine selfTest
|
|||
real :: A,B
|
||||
integer :: i
|
||||
|
||||
do i=1,10
|
||||
|
||||
do i = 1, 10
|
||||
|
||||
#if defined(__GFORTRAN__) && __GNUC__<9
|
||||
if(i<7) cycle
|
||||
|
|
|
|||
Loading…
Reference in New Issue