223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge branch 'integrate-lambert' into 'development' 

Integrate lambert

See merge request damask/DAMASK!166
This commit is contained in:
Karo 2020-05-12 23:36:08 +02:00
parent ae408baf27 5f3f87cd68
commit 90f93d2399
6 changed files with 489 additions and 533 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

164
python/damask/_Lambert.py
View File

@ -1,164 +0,0 @@
####################################################################################################
# Code below available according to the following conditions on
# https://github.com/MarDiehl/3Drotations
####################################################################################################
# Copyright (c) 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH
# Copyright (c) 2013-2014, 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.
####################################################################################################
import numpy as np
sc = np.pi**(1./6.)/6.**(1./6.)
beta = np.pi**(5./6.)/6.**(1./6.)/2.
R1 = (3.*np.pi/4.)**(1./3.)
def cube_to_ball(cube):
"""
Map a point in a uniform refinable cubical grid to a point on a uniform refinable grid on a ball.
Parameters
----------
cube : numpy.ndarray
coordinates of a point in a uniform refinable cubical grid.
References
----------
D. Roşca et al., Modelling and Simulation in Materials Science and Engineering 22:075013, 2014
https://doi.org/10.1088/0965-0393/22/7/075013
"""
cube_ = np.clip(cube,None,np.pi**(2./3.) * 0.5) if np.isclose(np.abs(np.max(cube)),np.pi**(2./3.) * 0.5,atol=1e-6) else cube
# transform to the sphere grid via the curved square, and intercept the zero point
if np.allclose(cube_,0.0,rtol=0.0,atol=1.0e-16):
ball = np.zeros(3)
else:
# get pyramide and scale by grid parameter ratio
p = _get_order(cube_)
XYZ = cube_[p[0]] * sc
# intercept all the points along the z-axis
if np.allclose(XYZ[0:2],0.0,rtol=0.0,atol=1.0e-16):
ball = np.array([0.0, 0.0, np.sqrt(6.0/np.pi) * XYZ[2]])
else:
order = [1,0] if np.abs(XYZ[1]) <= np.abs(XYZ[0]) else [0,1]
q = np.pi/12.0 * XYZ[order[0]]/XYZ[order[1]]
c = np.cos(q)
s = np.sin(q)
q = R1*2.0**0.25/beta * XYZ[order[1]] / np.sqrt(np.sqrt(2.0)-c)
T = np.array([ (np.sqrt(2.0)*c - 1.0), np.sqrt(2.0) * s]) * q
# transform to sphere grid (inverse Lambert)
# note that there is no need to worry about dividing by zero, since XYZ[2] can not become zero
c = np.sum(T**2)
s = c * np.pi/24.0 /XYZ[2]**2
c = c * np.sqrt(np.pi/24.0)/XYZ[2]
q = np.sqrt( 1.0 - s )
ball = np.array([ T[order[1]] * q, T[order[0]] * q, np.sqrt(6.0/np.pi) * XYZ[2] - c ])
# reverse the coordinates back to the regular order according to the original pyramid number
ball = ball[p[1]]
return ball
def ball_to_cube(ball):
"""
Map a point on a uniform refinable grid on a ball to a point in a uniform refinable cubical grid.
Parameters
----------
ball : numpy.ndarray
coordinates of a point on a uniform refinable grid on a ball.
References
----------
D. Roşca et al., Modelling and Simulation in Materials Science and Engineering 22:075013, 2014
https://doi.org/10.1088/0965-0393/22/7/075013
"""
ball_ = ball/np.linalg.norm(ball)*R1 if np.isclose(np.linalg.norm(ball),R1,atol=1e-6) else ball
rs = np.linalg.norm(ball_)
if np.allclose(ball_,0.0,rtol=0.0,atol=1.0e-16):
cube = np.zeros(3)
else:
p = _get_order(ball_)
xyz3 = ball_[p[0]]
# inverse M_3
xyz2 = xyz3[0:2] * np.sqrt( 2.0*rs/(rs+np.abs(xyz3[2])) )
# inverse M_2
qxy = np.sum(xyz2**2)
if np.isclose(qxy,0.0,rtol=0.0,atol=1.0e-16):
Tinv = np.zeros(2)
else:
q2 = qxy + np.max(np.abs(xyz2))**2
sq2 = np.sqrt(q2)
q = (beta/np.sqrt(2.0)/R1) * np.sqrt(q2*qxy/(q2-np.max(np.abs(xyz2))*sq2))
tt = np.clip((np.min(np.abs(xyz2))**2+np.max(np.abs(xyz2))*sq2)/np.sqrt(2.0)/qxy,-1.0,1.0)
Tinv = np.array([1.0,np.arccos(tt)/np.pi*12.0]) if np.abs(xyz2[1]) <= np.abs(xyz2[0]) else \
np.array([np.arccos(tt)/np.pi*12.0,1.0])
Tinv = q * np.where(xyz2<0.0,-Tinv,Tinv)
# inverse M_1
cube = np.array([ Tinv[0], Tinv[1], (-1.0 if xyz3[2] < 0.0 else 1.0) * rs / np.sqrt(6.0/np.pi) ]) /sc
# reverse the coordinates back to the regular order according to the original pyramid number
cube = cube[p[1]]
return cube
def _get_order(xyz):
"""
Get order of the coordinates.
Depending on the pyramid in which the point is located, the order need to be adjusted.
Parameters
----------
xyz : numpy.ndarray
coordinates of a point on a uniform refinable grid on a ball or
in a uniform refinable cubical grid.
References
----------
D. Roşca et al., Modelling and Simulation in Materials Science and Engineering 22:075013, 2014
https://doi.org/10.1088/0965-0393/22/7/075013
"""
if (abs(xyz[0])<= xyz[2]) and (abs(xyz[1])<= xyz[2]) or \
(abs(xyz[0])<=-xyz[2]) and (abs(xyz[1])<=-xyz[2]):
return [[0,1,2],[0,1,2]]
elif (abs(xyz[2])<= xyz[0]) and (abs(xyz[1])<= xyz[0]) or \
(abs(xyz[2])<=-xyz[0]) and (abs(xyz[1])<=-xyz[0]):
return [[1,2,0],[2,0,1]]
elif (abs(xyz[0])<= xyz[1]) and (abs(xyz[2])<= xyz[1]) or \
(abs(xyz[0])<=-xyz[1]) and (abs(xyz[2])<=-xyz[1]):
return [[2,0,1],[1,2,0]]

189
python/damask/_rotation.py
View File

@ -1,10 +1,14 @@
import numpy as np import numpy as np
from ._Lambert import ball_to_cube, cube_to_ball
from . import mechanics from . import mechanics
_P = -1 _P = -1
# parameters for conversion from/to cubochoric
_sc = np.pi**(1./6.)/6.**(1./6.)
_beta = np.pi**(5./6.)/6.**(1./6.)/2.
_R1 = (3.*np.pi/4.)**(1./3.)
def iszero(a): def iszero(a):
return np.isclose(a,0.0,atol=1.0e-12,rtol=0.0) return np.isclose(a,0.0,atol=1.0e-12,rtol=0.0)
@ -237,7 +241,7 @@ class Rotation:
"""Homochoric vector: (h_1, h_2, h_3).""" """Homochoric vector: (h_1, h_2, h_3)."""
return Rotation.qu2ho(self.quaternion) return Rotation.qu2ho(self.quaternion)
def asCubochoric(self): def as_cubochoric(self):
"""Cubochoric vector: (c_1, c_2, c_3).""" """Cubochoric vector: (c_1, c_2, c_3)."""
return Rotation.qu2cu(self.quaternion) return Rotation.qu2cu(self.quaternion)
@ -261,6 +265,7 @@ class Rotation:
asMatrix = as_matrix asMatrix = as_matrix
asRodrigues = as_Rodrigues asRodrigues = as_Rodrigues
asHomochoric = as_homochoric asHomochoric = as_homochoric
asCubochoric = as_cubochoric
################################################################################################ ################################################################################################
# Static constructors. The input data needs to follow the conventions, options allow to # Static constructors. The input data needs to follow the conventions, options allow to
@ -382,7 +387,7 @@ class Rotation:
return Rotation(Rotation.ho2qu(ho)) return Rotation(Rotation.ho2qu(ho))
@staticmethod @staticmethod
def fromCubochoric(cubochoric, def from_cubochoric(cubochoric,
P = -1): P = -1):
cu = np.array(cubochoric,dtype=float) cu = np.array(cubochoric,dtype=float)
@ -457,6 +462,7 @@ class Rotation:
fromMatrix = from_matrix fromMatrix = from_matrix
fromRodrigues = from_Rodrigues fromRodrigues = from_Rodrigues
fromHomochoric = from_homochoric fromHomochoric = from_homochoric
fromCubochoric = from_cubochoric
fromRandom = from_random fromRandom = from_random
#################################################################################################### ####################################################################################################
@ -1047,12 +1053,71 @@ class Rotation:
@staticmethod @staticmethod
def ho2cu(ho): def ho2cu(ho):
"""Homochoric vector to cubochoric vector.""" """
if len(ho.shape) == 1: Homochoric vector to cubochoric vector.
return ball_to_cube(ho)
else:
raise NotImplementedError('Support for multiple rotations missing')
References
----------
D. Roşca et al., Modelling and Simulation in Materials Science and Engineering 22:075013, 2014
https://doi.org/10.1088/0965-0393/22/7/075013
"""
if len(ho.shape) == 1:
rs = np.linalg.norm(ho)
if np.allclose(ho,0.0,rtol=0.0,atol=1.0e-16):
cu = np.zeros(3)
else:
xyz3 = ho[Rotation._get_pyramid_order(ho,'forward')]
# inverse M_3
xyz2 = xyz3[0:2] * np.sqrt( 2.0*rs/(rs+np.abs(xyz3[2])) )
# inverse M_2
qxy = np.sum(xyz2**2)
if np.isclose(qxy,0.0,rtol=0.0,atol=1.0e-16):
Tinv = np.zeros(2)
else:
q2 = qxy + np.max(np.abs(xyz2))**2
sq2 = np.sqrt(q2)
q = (_beta/np.sqrt(2.0)/_R1) * np.sqrt(q2*qxy/(q2-np.max(np.abs(xyz2))*sq2))
tt = np.clip((np.min(np.abs(xyz2))**2+np.max(np.abs(xyz2))*sq2)/np.sqrt(2.0)/qxy,-1.0,1.0)
Tinv = np.array([1.0,np.arccos(tt)/np.pi*12.0]) if np.abs(xyz2[1]) <= np.abs(xyz2[0]) else \
np.array([np.arccos(tt)/np.pi*12.0,1.0])
Tinv = q * np.where(xyz2<0.0,-Tinv,Tinv)
# inverse M_1
cu = np.array([ Tinv[0], Tinv[1], (-1.0 if xyz3[2] < 0.0 else 1.0) * rs / np.sqrt(6.0/np.pi) ]) /_sc
cu = cu[Rotation._get_pyramid_order(ho,'backward')]
else:
rs = np.linalg.norm(ho,axis=-1,keepdims=True)
xyz3 = np.take_along_axis(ho,Rotation._get_pyramid_order(ho,'forward'),-1)
with np.errstate(invalid='ignore',divide='ignore'):
# inverse M_3
xyz2 = xyz3[...,0:2] * np.sqrt( 2.0*rs/(rs+np.abs(xyz3[...,2:3])) )
qxy = np.sum(xyz2**2,axis=-1,keepdims=True)
q2 = qxy + np.max(np.abs(xyz2),axis=-1,keepdims=True)**2
sq2 = np.sqrt(q2)
q = (_beta/np.sqrt(2.0)/_R1) * np.sqrt(q2*qxy/(q2-np.max(np.abs(xyz2),axis=-1,keepdims=True)*sq2))
tt = np.clip((np.min(np.abs(xyz2),axis=-1,keepdims=True)**2\
+np.max(np.abs(xyz2),axis=-1,keepdims=True)*sq2)/np.sqrt(2.0)/qxy,-1.0,1.0)
T_inv = np.where(np.abs(xyz2[...,1:2]) <= np.abs(xyz2[...,0:1]),
np.block([np.ones_like(tt),np.arccos(tt)/np.pi*12.0]),
np.block([np.arccos(tt)/np.pi*12.0,np.ones_like(tt)]))*q
T_inv[xyz2<0.0] *= -1.0
T_inv[np.broadcast_to(np.isclose(qxy,0.0,rtol=0.0,atol=1.0e-12),T_inv.shape)] = 0.0
cu = np.block([T_inv, np.where(xyz3[...,2:3]<0.0,-np.ones_like(xyz3[...,2:3]),np.ones_like(xyz3[...,2:3])) \
* rs/np.sqrt(6.0/np.pi),
])/ _sc
cu[np.isclose(np.sum(np.abs(ho),axis=-1),0.0,rtol=0.0,atol=1.0e-16)] = 0.0
cu = np.take_along_axis(cu,Rotation._get_pyramid_order(ho,'backward'),-1)
return cu
#---------- Cubochoric ---------- #---------- Cubochoric ----------
@staticmethod @staticmethod
@ -1082,8 +1147,110 @@ class Rotation:
@staticmethod @staticmethod
def cu2ho(cu): def cu2ho(cu):
"""Cubochoric vector to homochoric vector.""" """
Cubochoric vector to homochoric vector.
References
----------
D. Roşca et al., Modelling and Simulation in Materials Science and Engineering 22:075013, 2014
https://doi.org/10.1088/0965-0393/22/7/075013
"""
if len(cu.shape) == 1: if len(cu.shape) == 1:
return cube_to_ball(cu) # transform to the sphere grid via the curved square, and intercept the zero point
if np.allclose(cu,0.0,rtol=0.0,atol=1.0e-16):
ho = np.zeros(3)
else: else:
raise NotImplementedError('Support for multiple rotations missing') # get pyramide and scale by grid parameter ratio
XYZ = cu[Rotation._get_pyramid_order(cu,'forward')] * _sc
# intercept all the points along the z-axis
if np.allclose(XYZ[0:2],0.0,rtol=0.0,atol=1.0e-16):
ho = np.array([0.0, 0.0, np.sqrt(6.0/np.pi) * XYZ[2]])
else:
order = [1,0] if np.abs(XYZ[1]) <= np.abs(XYZ[0]) else [0,1]
q = np.pi/12.0 * XYZ[order[0]]/XYZ[order[1]]
c = np.cos(q)
s = np.sin(q)
q = _R1*2.0**0.25/_beta * XYZ[order[1]] / np.sqrt(np.sqrt(2.0)-c)
T = np.array([ (np.sqrt(2.0)*c - 1.0), np.sqrt(2.0) * s]) * q
# transform to sphere grid (inverse Lambert)
# note that there is no need to worry about dividing by zero, since XYZ[2] can not become zero
c = np.sum(T**2)
s = c * np.pi/24.0 /XYZ[2]**2
c = c * np.sqrt(np.pi/24.0)/XYZ[2]
q = np.sqrt( 1.0 - s )
ho = np.array([ T[order[1]] * q, T[order[0]] * q, np.sqrt(6.0/np.pi) * XYZ[2] - c ])
ho = ho[Rotation._get_pyramid_order(cu,'backward')]
else:
with np.errstate(invalid='ignore',divide='ignore'):
# get pyramide and scale by grid parameter ratio
XYZ = np.take_along_axis(cu,Rotation._get_pyramid_order(cu,'forward'),-1) * _sc
order = np.abs(XYZ[...,1:2]) <= np.abs(XYZ[...,0:1])
q = np.pi/12.0 * np.where(order,XYZ[...,1:2],XYZ[...,0:1]) \
/ np.where(order,XYZ[...,0:1],XYZ[...,1:2])
c = np.cos(q)
s = np.sin(q)
q = _R1*2.0**0.25/_beta/ np.sqrt(np.sqrt(2.0)-c) \
* np.where(order,XYZ[...,0:1],XYZ[...,1:2])
T = np.block([ (np.sqrt(2.0)*c - 1.0), np.sqrt(2.0) * s]) * q
# transform to sphere grid (inverse Lambert)
c = np.sum(T**2,axis=-1,keepdims=True)
s = c * np.pi/24.0 /XYZ[...,2:3]**2
c = c * np.sqrt(np.pi/24.0)/XYZ[...,2:3]
q = np.sqrt( 1.0 - s)
ho = np.where(np.isclose(np.sum(np.abs(XYZ[...,0:2]),axis=-1,keepdims=True),0.0,rtol=0.0,atol=1.0e-16),
np.block([np.zeros_like(XYZ[...,0:2]),np.sqrt(6.0/np.pi) *XYZ[...,2:3]]),
np.block([np.where(order,T[...,0:1],T[...,1:2])*q,
np.where(order,T[...,1:2],T[...,0:1])*q,
np.sqrt(6.0/np.pi) * XYZ[...,2:3] - c])
)
ho[np.isclose(np.sum(np.abs(cu),axis=-1),0.0,rtol=0.0,atol=1.0e-16)] = 0.0
ho = np.take_along_axis(ho,Rotation._get_pyramid_order(cu,'backward'),-1)
return ho
@staticmethod
def _get_pyramid_order(xyz,direction=None):
"""
Get order of the coordinates.
Depending on the pyramid in which the point is located, the order need to be adjusted.
Parameters
----------
xyz : numpy.ndarray
coordinates of a point on a uniform refinable grid on a ball or
in a uniform refinable cubical grid.
References
----------
D. Roşca et al., Modelling and Simulation in Materials Science and Engineering 22:075013, 2014
https://doi.org/10.1088/0965-0393/22/7/075013
"""
order = {'forward':np.array([[0,1,2],[1,2,0],[2,0,1]]),
'backward':np.array([[0,1,2],[2,0,1],[1,2,0]])}
if len(xyz.shape) == 1:
if np.maximum(abs(xyz[0]),abs(xyz[1])) <= xyz[2] or \
np.maximum(abs(xyz[0]),abs(xyz[1])) <=-xyz[2]:
p = 0
elif np.maximum(abs(xyz[1]),abs(xyz[2])) <= xyz[0] or \
np.maximum(abs(xyz[1]),abs(xyz[2])) <=-xyz[0]:
p = 1
elif np.maximum(abs(xyz[2]),abs(xyz[0])) <= xyz[1] or \
np.maximum(abs(xyz[2]),abs(xyz[0])) <=-xyz[1]:
p = 2
else:
p = np.where(np.maximum(np.abs(xyz[...,0]),np.abs(xyz[...,1])) <= np.abs(xyz[...,2]),0,
np.where(np.maximum(np.abs(xyz[...,1]),np.abs(xyz[...,2])) <= np.abs(xyz[...,0]),1,2))
return order[direction][p]

98
python/tests/test_Rotation.py
View File

@ -78,9 +78,9 @@ def default():
specials_scatter /= np.linalg.norm(specials_scatter,axis=1).reshape(-1,1) specials_scatter /= np.linalg.norm(specials_scatter,axis=1).reshape(-1,1)
specials_scatter[specials_scatter[:,0]<0]*=-1 specials_scatter[specials_scatter[:,0]<0]*=-1
return [Rotation.fromQuaternion(s) for s in specials] + \ return [Rotation.from_quaternion(s) for s in specials] + \
[Rotation.fromQuaternion(s) for s in specials_scatter] + \ [Rotation.from_quaternion(s) for s in specials_scatter] + \
[Rotation.fromRandom() for _ in range(n-len(specials)-len(specials_scatter))] [Rotation.from_random() for _ in range(n-len(specials)-len(specials_scatter))]
@pytest.fixture @pytest.fixture
def reference_dir(reference_dir_base): def reference_dir(reference_dir_base):
@ -92,41 +92,41 @@ class TestRotation:
def test_Eulers(self,default): def test_Eulers(self,default):
for rot in default: for rot in default:
m = rot.asQuaternion() m = rot.as_quaternion()
o = Rotation.fromEulers(rot.asEulers()).asQuaternion() o = Rotation.from_Eulers(rot.as_Eulers()).as_quaternion()
ok = np.allclose(m,o,atol=atol) ok = np.allclose(m,o,atol=atol)
if np.isclose(rot.asQuaternion()[0],0.0,atol=atol): if np.isclose(rot.as_quaternion()[0],0.0,atol=atol):
ok = ok or np.allclose(m*-1.,o,atol=atol) ok = ok or np.allclose(m*-1.,o,atol=atol)
print(m,o,rot.asQuaternion()) print(m,o,rot.as_quaternion())
assert ok and np.isclose(np.linalg.norm(o),1.0) assert ok and np.isclose(np.linalg.norm(o),1.0)
def test_AxisAngle(self,default): def test_AxisAngle(self,default):
for rot in default: for rot in default:
m = rot.asEulers() m = rot.as_Eulers()
o = Rotation.fromAxisAngle(rot.asAxisAngle()).asEulers() o = Rotation.from_axis_angle(rot.as_axis_angle()).as_Eulers()
u = np.array([np.pi*2,np.pi,np.pi*2]) u = np.array([np.pi*2,np.pi,np.pi*2])
ok = np.allclose(m,o,atol=atol) ok = np.allclose(m,o,atol=atol)
ok = ok or np.allclose(np.where(np.isclose(m,u),m-u,m),np.where(np.isclose(o,u),o-u,o),atol=atol) ok = ok or np.allclose(np.where(np.isclose(m,u),m-u,m),np.where(np.isclose(o,u),o-u,o),atol=atol)
if np.isclose(m[1],0.0,atol=atol) or np.isclose(m[1],np.pi,atol=atol): if np.isclose(m[1],0.0,atol=atol) or np.isclose(m[1],np.pi,atol=atol):
sum_phi = np.unwrap([m[0]+m[2],o[0]+o[2]]) sum_phi = np.unwrap([m[0]+m[2],o[0]+o[2]])
ok = ok or np.isclose(sum_phi[0],sum_phi[1],atol=atol) ok = ok or np.isclose(sum_phi[0],sum_phi[1],atol=atol)
print(m,o,rot.asQuaternion()) print(m,o,rot.as_quaternion())
assert ok and (np.zeros(3)-1.e-9 <= o).all() and (o <= np.array([np.pi*2.,np.pi,np.pi*2.])+1.e-9).all() assert ok and (np.zeros(3)-1.e-9 <= o).all() and (o <= np.array([np.pi*2.,np.pi,np.pi*2.])+1.e-9).all()
def test_Matrix(self,default): def test_Matrix(self,default):
for rot in default: for rot in default:
m = rot.asAxisAngle() m = rot.as_axis_angle()
o = Rotation.fromAxisAngle(rot.asAxisAngle()).asAxisAngle() o = Rotation.from_axis_angle(rot.as_axis_angle()).as_axis_angle()
ok = np.allclose(m,o,atol=atol) ok = np.allclose(m,o,atol=atol)
if np.isclose(m[3],np.pi,atol=atol): if np.isclose(m[3],np.pi,atol=atol):
ok = ok or np.allclose(m*np.array([-1.,-1.,-1.,1.]),o,atol=atol) ok = ok or np.allclose(m*np.array([-1.,-1.,-1.,1.]),o,atol=atol)
print(m,o,rot.asQuaternion()) print(m,o,rot.as_quaternion())
assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi++1.e-9 assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi++1.e-9
def test_Rodrigues(self,default): def test_Rodrigues(self,default):
for rot in default: for rot in default:
m = rot.asMatrix() m = rot.as_matrix()
o = Rotation.fromRodrigues(rot.asRodrigues()).asMatrix() o = Rotation.from_Rodrigues(rot.as_Rodrigues()).as_matrix()
ok = np.allclose(m,o,atol=atol) ok = np.allclose(m,o,atol=atol)
print(m,o) print(m,o)
assert ok and np.isclose(np.linalg.det(o),1.0) assert ok and np.isclose(np.linalg.det(o),1.0)
@ -134,27 +134,27 @@ class TestRotation:
def test_Homochoric(self,default): def test_Homochoric(self,default):
cutoff = np.tan(np.pi*.5*(1.-1e-4)) cutoff = np.tan(np.pi*.5*(1.-1e-4))
for rot in default: for rot in default:
m = rot.asRodrigues() m = rot.as_Rodrigues()
o = Rotation.fromHomochoric(rot.asHomochoric()).asRodrigues() o = Rotation.from_homochoric(rot.as_homochoric()).as_Rodrigues()
ok = np.allclose(np.clip(m,None,cutoff),np.clip(o,None,cutoff),atol=atol) ok = np.allclose(np.clip(m,None,cutoff),np.clip(o,None,cutoff),atol=atol)
ok = ok or np.isclose(m[3],0.0,atol=atol) ok = ok or np.isclose(m[3],0.0,atol=atol)
print(m,o,rot.asQuaternion()) print(m,o,rot.as_quaternion())
assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) assert ok and np.isclose(np.linalg.norm(o[:3]),1.0)
def test_Cubochoric(self,default): def test_Cubochoric(self,default):
for rot in default: for rot in default:
m = rot.asHomochoric() m = rot.as_homochoric()
o = Rotation.fromCubochoric(rot.asCubochoric()).asHomochoric() o = Rotation.from_cubochoric(rot.as_cubochoric()).as_homochoric()
ok = np.allclose(m,o,atol=atol) ok = np.allclose(m,o,atol=atol)
print(m,o,rot.asQuaternion()) print(m,o,rot.as_quaternion())
assert ok and np.linalg.norm(o) < (3.*np.pi/4.)**(1./3.) + 1.e-9 assert ok and np.linalg.norm(o) < (3.*np.pi/4.)**(1./3.) + 1.e-9
def test_Quaternion(self,default): def test_Quaternion(self,default):
for rot in default: for rot in default:
m = rot.asCubochoric() m = rot.as_cubochoric()
o = Rotation.fromQuaternion(rot.asQuaternion()).asCubochoric() o = Rotation.from_quaternion(rot.as_quaternion()).as_cubochoric()
ok = np.allclose(m,o,atol=atol) ok = np.allclose(m,o,atol=atol)
print(m,o,rot.asQuaternion()) print(m,o,rot.as_quaternion())
assert ok and o.max() < np.pi**(2./3.)*0.5+1.e-9 assert ok and o.max() < np.pi**(2./3.)*0.5+1.e-9
@pytest.mark.parametrize('function',[Rotation.from_quaternion, @pytest.mark.parametrize('function',[Rotation.from_quaternion,
@ -185,9 +185,11 @@ class TestRotation:
Rotation.qu2eu, Rotation.qu2eu,
Rotation.qu2ax, Rotation.qu2ax,
Rotation.qu2ro, Rotation.qu2ro,
Rotation.qu2ho]) Rotation.qu2ho,
Rotation.qu2cu
])
def test_quaternion_vectorization(self,default,conversion): def test_quaternion_vectorization(self,default,conversion):
qu = np.array([rot.asQuaternion() for rot in default]) qu = np.array([rot.as_quaternion() for rot in default])
conversion(qu.reshape(qu.shape[0]//2,-1,4)) conversion(qu.reshape(qu.shape[0]//2,-1,4))
co = conversion(qu) co = conversion(qu)
for q,c in zip(qu,co): for q,c in zip(qu,co):
@ -199,9 +201,10 @@ class TestRotation:
Rotation.om2ax, Rotation.om2ax,
Rotation.om2ro, Rotation.om2ro,
Rotation.om2ho, Rotation.om2ho,
Rotation.om2cu
]) ])
def test_matrix_vectorization(self,default,conversion): def test_matrix_vectorization(self,default,conversion):
om = np.array([rot.asMatrix() for rot in default]) om = np.array([rot.as_matrix() for rot in default])
conversion(om.reshape(om.shape[0]//2,-1,3,3)) conversion(om.reshape(om.shape[0]//2,-1,3,3))
co = conversion(om) co = conversion(om)
for o,c in zip(om,co): for o,c in zip(om,co):
@ -213,9 +216,10 @@ class TestRotation:
Rotation.eu2ax, Rotation.eu2ax,
Rotation.eu2ro, Rotation.eu2ro,
Rotation.eu2ho, Rotation.eu2ho,
Rotation.eu2cu
]) ])
def test_Euler_vectorization(self,default,conversion): def test_Euler_vectorization(self,default,conversion):
eu = np.array([rot.asEulers() for rot in default]) eu = np.array([rot.as_Eulers() for rot in default])
conversion(eu.reshape(eu.shape[0]//2,-1,3)) conversion(eu.reshape(eu.shape[0]//2,-1,3))
co = conversion(eu) co = conversion(eu)
for e,c in zip(eu,co): for e,c in zip(eu,co):
@ -227,9 +231,10 @@ class TestRotation:
Rotation.ax2eu, Rotation.ax2eu,
Rotation.ax2ro, Rotation.ax2ro,
Rotation.ax2ho, Rotation.ax2ho,
Rotation.ax2cu
]) ])
def test_axisAngle_vectorization(self,default,conversion): def test_axisAngle_vectorization(self,default,conversion):
ax = np.array([rot.asAxisAngle() for rot in default]) ax = np.array([rot.as_axis_angle() for rot in default])
conversion(ax.reshape(ax.shape[0]//2,-1,4)) conversion(ax.reshape(ax.shape[0]//2,-1,4))
co = conversion(ax) co = conversion(ax)
for a,c in zip(ax,co): for a,c in zip(ax,co):
@ -242,9 +247,10 @@ class TestRotation:
Rotation.ro2eu, Rotation.ro2eu,
Rotation.ro2ax, Rotation.ro2ax,
Rotation.ro2ho, Rotation.ro2ho,
Rotation.ro2cu
]) ])
def test_Rodrigues_vectorization(self,default,conversion): def test_Rodrigues_vectorization(self,default,conversion):
ro = np.array([rot.asRodrigues() for rot in default]) ro = np.array([rot.as_Rodrigues() for rot in default])
conversion(ro.reshape(ro.shape[0]//2,-1,4)) conversion(ro.reshape(ro.shape[0]//2,-1,4))
co = conversion(ro) co = conversion(ro)
for r,c in zip(ro,co): for r,c in zip(ro,co):
@ -256,11 +262,41 @@ class TestRotation:
Rotation.ho2eu, Rotation.ho2eu,
Rotation.ho2ax, Rotation.ho2ax,
Rotation.ho2ro, Rotation.ho2ro,
Rotation.ho2cu
]) ])
def test_homochoric_vectorization(self,default,conversion): def test_homochoric_vectorization(self,default,conversion):
ho = np.array([rot.asHomochoric() for rot in default]) ho = np.array([rot.as_homochoric() for rot in default])
conversion(ho.reshape(ho.shape[0]//2,-1,3)) conversion(ho.reshape(ho.shape[0]//2,-1,3))
co = conversion(ho) co = conversion(ho)
for h,c in zip(ho,co): for h,c in zip(ho,co):
print(h,c) print(h,c)
assert np.allclose(conversion(h),c) assert np.allclose(conversion(h),c)
@pytest.mark.parametrize('conversion',[Rotation.cu2qu,
Rotation.cu2om,
Rotation.cu2eu,
Rotation.cu2ax,
Rotation.cu2ro,
Rotation.cu2ho
])
def test_cubochoric_vectorization(self,default,conversion):
cu = np.array([rot.as_cubochoric() for rot in default])
conversion(cu.reshape(cu.shape[0]//2,-1,3))
co = conversion(cu)
for u,c in zip(cu,co):
print(u,c)
assert np.allclose(conversion(u),c)
@pytest.mark.parametrize('direction',['forward',
'backward'])
def test_pyramid_vectorization(self,direction):
p = np.random.rand(n,3)
o = Rotation._get_pyramid_order(p,direction)
for i,o_i in enumerate(o):
assert np.all(o_i==Rotation._get_pyramid_order(p[i],direction))
def test_pyramid_invariant(self):
a = np.random.rand(n,3)
f = Rotation._get_pyramid_order(a,'forward')
b = Rotation._get_pyramid_order(a,'backward')
assert np.all(np.take_along_axis(np.take_along_axis(a,f,-1),b,-1) == a)

201
src/Lambert.f90
View File

@ -1,201 +0,0 @@
! ###################################################################
! Copyright (c) 2013-2015, Marc De Graef/Carnegie Mellon University
! Modified 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH
! 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
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Mapping homochoric <-> cubochoric
!
!> @details
!> 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 prec
use math
implicit none
private
real(pReal), parameter :: &
SPI = sqrt(PI), &
PREF = sqrt(6.0_pReal/PI), &
A = PI**(5.0_pReal/6.0_pReal)/6.0_pReal**(1.0_pReal/6.0_pReal), &
AP = PI**(2.0_pReal/3.0_pReal), &
SC = A/AP, &
BETA = A/2.0_pReal, &
R1 = (3.0_pReal*PI/4.0_pReal)**(1.0_pReal/3.0_pReal), &
R2 = sqrt(2.0_pReal), &
PI12 = PI/12.0_pReal, &
PREK = R1 * 2.0_pReal**(1.0_pReal/4.0_pReal)/BETA
public :: &
Lambert_CubeToBall, &
Lambert_BallToCube
contains
!--------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief map from 3D cubic grid to 3D ball
!--------------------------------------------------------------------------
pure function Lambert_CubeToBall(cube) result(ball)
real(pReal), intent(in), dimension(3) :: cube
real(pReal), dimension(3) :: ball, LamXYZ, XYZ
real(pReal), dimension(2) :: T
real(pReal) :: c, s, q
real(pReal), parameter :: eps = 1.0e-8_pReal
integer, dimension(3,2) :: p
integer, 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(:,1)) * 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 / sqrt(24.0_pReal) / 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 order according to the original pyramid number
ball = LamXYZ(p(:,2))
endif center
end function Lambert_CubeToBall
!--------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief map from 3D ball to 3D cubic grid
!--------------------------------------------------------------------------
pure function Lambert_BallToCube(xyz) result(cube)
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, dimension(3,2) :: 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(:,1))
! 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_pReal
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_pReal,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_pReal,xyz3(3)) * rs / pref ] /sc
! reverse the coordinates back to order according to the original pyramid number
cube = xyz1(p(:,2))
endif center
end function Lambert_BallToCube
!--------------------------------------------------------------------------
!> @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
!--------------------------------------------------------------------------
pure function GetPyramidOrder(xyz)
real(pReal),intent(in),dimension(3) :: xyz
integer, dimension(3,2) :: 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 = reshape([[1,2,3],[1,2,3]],[3,2])
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 = reshape([[2,3,1],[3,1,2]],[3,2])
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 = reshape([[3,1,2],[2,3,1]],[3,2])
else
GetPyramidOrder = -1 ! should be impossible, but might simplify debugging
end if
end function GetPyramidOrder
end module Lambert

1
src/commercialFEM_fileList.f90
View File

@ -13,7 +13,6 @@
#include "LAPACK_interface.f90" #include "LAPACK_interface.f90"
#include "math.f90" #include "math.f90"
#include "quaternions.f90" #include "quaternions.f90"
#include "Lambert.f90"
#include "rotations.f90" #include "rotations.f90"
#include "FEsolving.f90" #include "FEsolving.f90"
#include "element.f90" #include "element.f90"

141
src/rotations.f90
View File

@ -50,7 +50,6 @@ module rotations
use prec use prec
use IO use IO
use math use math
use Lambert
use quaternions use quaternions
implicit none implicit none
@ -81,6 +80,18 @@ module rotations
procedure, public :: standardize procedure, public :: standardize
end type rotation end type rotation
real(pReal), parameter :: &
SPI = sqrt(PI), &
PREF = sqrt(6.0_pReal/PI), &
A = PI**(5.0_pReal/6.0_pReal)/6.0_pReal**(1.0_pReal/6.0_pReal), &
AP = PI**(2.0_pReal/3.0_pReal), &
SC = A/AP, &
BETA = A/2.0_pReal, &
R1 = (3.0_pReal*PI/4.0_pReal)**(1.0_pReal/3.0_pReal), &
R2 = sqrt(2.0_pReal), &
PI12 = PI/12.0_pReal, &
PREK = R1 * 2.0_pReal**(1.0_pReal/4.0_pReal)/BETA
public :: & public :: &
rotations_init, & rotations_init, &
eu2om eu2om
@ -516,7 +527,7 @@ pure function qu2ho(qu) result(ho)
omega = 2.0 * acos(math_clip(qu(1),-1.0_pReal,1.0_pReal)) omega = 2.0 * acos(math_clip(qu(1),-1.0_pReal,1.0_pReal))
if (dEq0(omega)) then if (dEq0(omega,tol=1.e-5_pReal)) then
ho = [ 0.0_pReal, 0.0_pReal, 0.0_pReal ] ho = [ 0.0_pReal, 0.0_pReal, 0.0_pReal ]
else else
ho = qu(2:4) ho = qu(2:4)
@ -1120,16 +1131,56 @@ pure function ho2ro(ho) result(ro)
end function ho2ro end function ho2ro
!--------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University !> @author Marc De Graef, Carnegie Mellon University
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief convert homochoric to cubochoric !> @brief convert homochoric to cubochoric
!--------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------
pure function ho2cu(ho) result(cu) pure function ho2cu(ho) result(cu)
real(pReal), intent(in), dimension(3) :: ho real(pReal), intent(in), dimension(3) :: ho
real(pReal), dimension(3) :: cu real(pReal), dimension(3) :: cu, xyz1, xyz3
real(pReal), dimension(2) :: Tinv, xyz2
real(pReal) :: rs, qxy, q2, sq2, q, tt
integer, dimension(3,2) :: p
cu = Lambert_BallToCube(ho) rs = norm2(ho)
if (rs > R1+1.e-6_pReal) then
cu = IEEE_value(cu,IEEE_positive_inf)
return
endif
center: if (all(dEq0(ho))) then
cu = 0.0_pReal
else center
p = GetPyramidOrder(ho)
xyz3 = ho(p(:,1))
! 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_pReal
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_pReal,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_pReal,xyz3(3)) * rs / pref ] /sc
! reverse the coordinates back to order according to the original pyramid number
cu = xyz1(p(:,2))
endif center
end function ho2cu end function ho2cu
@ -1204,20 +1255,88 @@ pure function cu2ro(cu) result(ro)
end function cu2ro end function cu2ro
!--------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University !> @author Marc De Graef, Carnegie Mellon University
!> @brief convert cubochoric to homochoric !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!--------------------------------------------------------------------------------------------------- !> @brief map from 3D cubic grid to 3D ball
!--------------------------------------------------------------------------
pure function cu2ho(cu) result(ho) pure function cu2ho(cu) result(ho)
real(pReal), intent(in), dimension(3) :: cu real(pReal), intent(in), dimension(3) :: cu
real(pReal), dimension(3) :: ho real(pReal), dimension(3) :: ho, LamXYZ, XYZ
real(pReal), dimension(2) :: T
real(pReal) :: c, s, q
real(pReal), parameter :: eps = 1.0e-8_pReal
integer, dimension(3,2) :: p
integer, dimension(2) :: order
ho = Lambert_CubeToBall(cu) if (maxval(abs(cu)) > AP/2.0+eps) then
ho = IEEE_value(cu,IEEE_positive_inf)
return
end if
! transform to the sphere grid via the curved square, and intercept the zero point
center: if (all(dEq0(cu))) then
ho = 0.0_pReal
else center
! get pyramide and scale by grid parameter ratio
p = GetPyramidOrder(cu)
XYZ = cu(p(:,1)) * 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 / sqrt(24.0_pReal) / 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 order according to the original pyramid number
ho = LamXYZ(p(:,2))
endif center
end function cu2ho 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
!--------------------------------------------------------------------------
pure function GetPyramidOrder(xyz)
real(pReal),intent(in),dimension(3) :: xyz
integer, dimension(3,2) :: 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 = reshape([[1,2,3],[1,2,3]],[3,2])
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 = reshape([[2,3,1],[3,1,2]],[3,2])
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 = reshape([[3,1,2],[2,3,1]],[3,2])
else
GetPyramidOrder = -1 ! should be impossible, but might simplify debugging
end if
end function GetPyramidOrder
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief check correctness of some rotations functions !> @brief check correctness of some rotations functions
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------