diff --git a/python/damask/_Lambert.py b/python/damask/_Lambert.py deleted file mode 100644 index a823764e9..000000000 --- a/python/damask/_Lambert.py +++ /dev/null @@ -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]] diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 716977a79..a13a9af0b 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -1,10 +1,14 @@ import numpy as np -from ._Lambert import ball_to_cube, cube_to_ball from . import mechanics _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): 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).""" return Rotation.qu2ho(self.quaternion) - def asCubochoric(self): + def as_cubochoric(self): """Cubochoric vector: (c_1, c_2, c_3).""" return Rotation.qu2cu(self.quaternion) @@ -261,6 +265,7 @@ class Rotation: asMatrix = as_matrix asRodrigues = as_Rodrigues asHomochoric = as_homochoric + asCubochoric = as_cubochoric ################################################################################################ # Static constructors. The input data needs to follow the conventions, options allow to @@ -382,7 +387,7 @@ class Rotation: return Rotation(Rotation.ho2qu(ho)) @staticmethod - def fromCubochoric(cubochoric, + def from_cubochoric(cubochoric, P = -1): cu = np.array(cubochoric,dtype=float) @@ -457,6 +462,7 @@ class Rotation: fromMatrix = from_matrix fromRodrigues = from_Rodrigues fromHomochoric = from_homochoric + fromCubochoric = from_cubochoric fromRandom = from_random #################################################################################################### @@ -1047,12 +1053,71 @@ class Rotation: @staticmethod def ho2cu(ho): - """Homochoric vector to cubochoric vector.""" - if len(ho.shape) == 1: - return ball_to_cube(ho) - else: - raise NotImplementedError('Support for multiple rotations missing') + """ + Homochoric vector to cubochoric 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(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 ---------- @staticmethod @@ -1082,8 +1147,110 @@ class Rotation: @staticmethod 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: - 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: + # 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: - raise NotImplementedError('Support for multiple rotations missing') + 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] diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 302e895ee..b7442035f 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -78,9 +78,9 @@ def default(): specials_scatter /= np.linalg.norm(specials_scatter,axis=1).reshape(-1,1) specials_scatter[specials_scatter[:,0]<0]*=-1 - return [Rotation.fromQuaternion(s) for s in specials] + \ - [Rotation.fromQuaternion(s) for s in specials_scatter] + \ - [Rotation.fromRandom() for _ in range(n-len(specials)-len(specials_scatter))] + return [Rotation.from_quaternion(s) for s in specials] + \ + [Rotation.from_quaternion(s) for s in specials_scatter] + \ + [Rotation.from_random() for _ in range(n-len(specials)-len(specials_scatter))] @pytest.fixture def reference_dir(reference_dir_base): @@ -92,41 +92,41 @@ class TestRotation: def test_Eulers(self,default): for rot in default: - m = rot.asQuaternion() - o = Rotation.fromEulers(rot.asEulers()).asQuaternion() + m = rot.as_quaternion() + o = Rotation.from_Eulers(rot.as_Eulers()).as_quaternion() 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) - print(m,o,rot.asQuaternion()) + print(m,o,rot.as_quaternion()) assert ok and np.isclose(np.linalg.norm(o),1.0) def test_AxisAngle(self,default): for rot in default: - m = rot.asEulers() - o = Rotation.fromAxisAngle(rot.asAxisAngle()).asEulers() + m = rot.as_Eulers() + o = Rotation.from_axis_angle(rot.as_axis_angle()).as_Eulers() u = np.array([np.pi*2,np.pi,np.pi*2]) 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) 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]]) 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() def test_Matrix(self,default): for rot in default: - m = rot.asAxisAngle() - o = Rotation.fromAxisAngle(rot.asAxisAngle()).asAxisAngle() + m = rot.as_axis_angle() + o = Rotation.from_axis_angle(rot.as_axis_angle()).as_axis_angle() ok = np.allclose(m,o,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) - 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 def test_Rodrigues(self,default): for rot in default: - m = rot.asMatrix() - o = Rotation.fromRodrigues(rot.asRodrigues()).asMatrix() + m = rot.as_matrix() + o = Rotation.from_Rodrigues(rot.as_Rodrigues()).as_matrix() ok = np.allclose(m,o,atol=atol) print(m,o) assert ok and np.isclose(np.linalg.det(o),1.0) @@ -134,27 +134,27 @@ class TestRotation: def test_Homochoric(self,default): cutoff = np.tan(np.pi*.5*(1.-1e-4)) for rot in default: - m = rot.asRodrigues() - o = Rotation.fromHomochoric(rot.asHomochoric()).asRodrigues() + m = rot.as_Rodrigues() + 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 = 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) def test_Cubochoric(self,default): for rot in default: - m = rot.asHomochoric() - o = Rotation.fromCubochoric(rot.asCubochoric()).asHomochoric() + m = rot.as_homochoric() + o = Rotation.from_cubochoric(rot.as_cubochoric()).as_homochoric() 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 def test_Quaternion(self,default): for rot in default: - m = rot.asCubochoric() - o = Rotation.fromQuaternion(rot.asQuaternion()).asCubochoric() + m = rot.as_cubochoric() + o = Rotation.from_quaternion(rot.as_quaternion()).as_cubochoric() 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 @pytest.mark.parametrize('function',[Rotation.from_quaternion, @@ -185,9 +185,11 @@ class TestRotation: Rotation.qu2eu, Rotation.qu2ax, Rotation.qu2ro, - Rotation.qu2ho]) + Rotation.qu2ho, + Rotation.qu2cu + ]) 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)) co = conversion(qu) for q,c in zip(qu,co): @@ -199,9 +201,10 @@ class TestRotation: Rotation.om2ax, Rotation.om2ro, Rotation.om2ho, + Rotation.om2cu ]) 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)) co = conversion(om) for o,c in zip(om,co): @@ -213,9 +216,10 @@ class TestRotation: Rotation.eu2ax, Rotation.eu2ro, Rotation.eu2ho, + Rotation.eu2cu ]) 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)) co = conversion(eu) for e,c in zip(eu,co): @@ -227,9 +231,10 @@ class TestRotation: Rotation.ax2eu, Rotation.ax2ro, Rotation.ax2ho, + Rotation.ax2cu ]) 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)) co = conversion(ax) for a,c in zip(ax,co): @@ -242,9 +247,10 @@ class TestRotation: Rotation.ro2eu, Rotation.ro2ax, Rotation.ro2ho, + Rotation.ro2cu ]) 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)) co = conversion(ro) for r,c in zip(ro,co): @@ -256,11 +262,41 @@ class TestRotation: Rotation.ho2eu, Rotation.ho2ax, Rotation.ho2ro, + Rotation.ho2cu ]) 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)) co = conversion(ho) for h,c in zip(ho,co): print(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) diff --git a/src/Lambert.f90 b/src/Lambert.f90 deleted file mode 100644 index 932fe221b..000000000 --- a/src/Lambert.f90 +++ /dev/null @@ -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 diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index f13778f01..f8a5ad96a 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -13,7 +13,6 @@ #include "LAPACK_interface.f90" #include "math.f90" #include "quaternions.f90" -#include "Lambert.f90" #include "rotations.f90" #include "FEsolving.f90" #include "element.f90" diff --git a/src/rotations.f90 b/src/rotations.f90 index c20b9e905..fea642170 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -3,27 +3,27 @@ ! Modified 2017-2020, 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 +! 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 +! - 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 +! - 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 +! - 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 +! 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. ! ################################################################### @@ -31,7 +31,7 @@ !> @author Marc De Graef, Carnegie Mellon University !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief rotation storage and conversion -!> @details: rotation is internally stored as quaternion. It can be inialized from different +!> @details: rotation is internally stored as quaternion. It can be inialized from different !> representations and also returns itself in different representations. ! ! All methods and naming conventions based on Rowenhorst_etal2015 @@ -50,9 +50,8 @@ module rotations use prec use IO use math - use Lambert use quaternions - + implicit none private @@ -80,7 +79,19 @@ module rotations procedure, public :: misorientation procedure, public :: standardize 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 :: & rotations_init, & eu2om @@ -106,16 +117,16 @@ pure function asQuaternion(self) class(rotation), intent(in) :: self real(pReal), dimension(4) :: asQuaternion - + asQuaternion = self%q%asArray() end function asQuaternion !--------------------------------------------------------------------------------------------------- pure function asEulers(self) - + class(rotation), intent(in) :: self real(pReal), dimension(3) :: asEulers - + asEulers = qu2eu(self%q%asArray()) end function asEulers @@ -124,16 +135,16 @@ pure function asAxisAngle(self) class(rotation), intent(in) :: self real(pReal), dimension(4) :: asAxisAngle - + asAxisAngle = qu2ax(self%q%asArray()) end function asAxisAngle !--------------------------------------------------------------------------------------------------- pure function asMatrix(self) - + class(rotation), intent(in) :: self real(pReal), dimension(3,3) :: asMatrix - + asMatrix = qu2om(self%q%asArray()) end function asMatrix @@ -142,20 +153,20 @@ pure function asRodrigues(self) class(rotation), intent(in) :: self real(pReal), dimension(4) :: asRodrigues - + asRodrigues = qu2ro(self%q%asArray()) - + end function asRodrigues !--------------------------------------------------------------------------------------------------- pure function asHomochoric(self) class(rotation), intent(in) :: self real(pReal), dimension(3) :: asHomochoric - + asHomochoric = qu2ho(self%q%asArray()) end function asHomochoric - + !--------------------------------------------------------------------------------------------------- ! Initialize rotation from different representations !--------------------------------------------------------------------------------------------------- @@ -207,7 +218,7 @@ subroutine fromAxisAngle(self,ax,degrees,P) else angle = merge(ax(4)*INRAD,ax(4),degrees) endif - + if (.not. present(P)) then axis = ax(1:3) else @@ -217,7 +228,7 @@ subroutine fromAxisAngle(self,ax,degrees,P) if(dNeq(norm2(axis),1.0_pReal) .or. angle < 0.0_pReal .or. angle > PI) & call IO_error(402,ext_msg='fromAxisAngle') - + self%q = ax2qu([axis,angle]) end subroutine fromAxisAngle @@ -240,10 +251,10 @@ end subroutine fromMatrix !> @brief: Rotate a rotation !--------------------------------------------------------------------------------------------------- pure elemental function rotRot__(self,R) result(rRot) - + type(rotation) :: rRot class(rotation), intent(in) :: self,R - + rRot = rotation(self%q*R%q) call rRot%standardize() @@ -251,12 +262,12 @@ end function rotRot__ !--------------------------------------------------------------------------------------------------- -!> @brief quaternion representation with positive q +!> @brief quaternion representation with positive q !--------------------------------------------------------------------------------------------------- pure elemental subroutine standardize(self) class(rotation), intent(inout) :: self - + if (real(self%q) < 0.0_pReal) self%q = self%q%homomorphed() end subroutine standardize @@ -267,22 +278,22 @@ end subroutine standardize !> @brief rotate a vector passively (default) or actively !--------------------------------------------------------------------------------------------------- pure function rotVector(self,v,active) result(vRot) - + real(pReal), dimension(3) :: vRot class(rotation), intent(in) :: self real(pReal), intent(in), dimension(3) :: v logical, intent(in), optional :: active - + real(pReal), dimension(3) :: v_normed type(quaternion) :: q logical :: passive - + if (present(active)) then passive = .not. active else passive = .true. endif - + if (dEq0(norm2(v))) then vRot = v else @@ -304,12 +315,12 @@ end function rotVector !> @details: rotation is based on rotation matrix !--------------------------------------------------------------------------------------------------- pure function rotTensor2(self,T,active) result(tRot) - + real(pReal), dimension(3,3) :: tRot class(rotation), intent(in) :: self real(pReal), intent(in), dimension(3,3) :: T logical, intent(in), optional :: active - + logical :: passive if (present(active)) then @@ -317,7 +328,7 @@ pure function rotTensor2(self,T,active) result(tRot) else passive = .true. endif - + if (passive) then tRot = matmul(matmul(self%asMatrix(),T),transpose(self%asMatrix())) else @@ -339,7 +350,7 @@ pure function rotTensor4(self,T,active) result(tRot) class(rotation), intent(in) :: self real(pReal), intent(in), dimension(3,3,3,3) :: T logical, intent(in), optional :: active - + real(pReal), dimension(3,3) :: R integer :: i,j,k,l,m,n,o,p @@ -370,7 +381,7 @@ pure function rotTensor4sym(self,T,active) result(tRot) class(rotation), intent(in) :: self real(pReal), intent(in), dimension(6,6) :: T logical, intent(in), optional :: active - + if (present(active)) then tRot = math_sym3333to66(rotTensor4(self,math_66toSym3333(T),active)) else @@ -384,10 +395,10 @@ end function rotTensor4sym !> @brief misorientation !--------------------------------------------------------------------------------------------------- pure elemental function misorientation(self,other) - + type(rotation) :: misorientation class(rotation), intent(in) :: self, other - + misorientation%q = other%q * conjg(self%q) end function misorientation @@ -401,7 +412,7 @@ pure function qu2om(qu) result(om) real(pReal), intent(in), dimension(4) :: qu real(pReal), dimension(3,3) :: om - + real(pReal) :: qq qq = qu(1)**2-sum(qu(2:4)**2) @@ -431,13 +442,13 @@ pure function qu2eu(qu) result(eu) real(pReal), intent(in), dimension(4) :: qu real(pReal), dimension(3) :: eu - + real(pReal) :: q12, q03, chi - + q03 = qu(1)**2+qu(4)**2 q12 = qu(2)**2+qu(3)**2 chi = sqrt(q03*q12) - + degenerated: if (dEq0(q12)) then eu = [atan2(-P*2.0_pReal*qu(1)*qu(4),qu(1)**2-qu(4)**2), 0.0_pReal, 0.0_pReal] elseif (dEq0(q03)) then @@ -460,7 +471,7 @@ pure function qu2ax(qu) result(ax) real(pReal), intent(in), dimension(4) :: qu real(pReal), dimension(4) :: ax - + real(pReal) :: omega, s if (dEq0(sum(qu(2:4)**2))) then @@ -481,13 +492,13 @@ end function qu2ax !> @brief convert unit quaternion to Rodrigues vector !--------------------------------------------------------------------------------------------------- pure function qu2ro(qu) result(ro) - + real(pReal), intent(in), dimension(4) :: qu real(pReal), dimension(4) :: ro - + real(pReal) :: s real(pReal), parameter :: thr = 1.0e-8_pReal - + if (abs(qu(1)) < thr) then ro = [qu(2), qu(3), qu(4), IEEE_value(1.0_pReal,IEEE_positive_inf)] else @@ -497,7 +508,7 @@ pure function qu2ro(qu) result(ro) else ro = [qu(2)/s,qu(3)/s,qu(4)/s, tan(acos(math_clip(qu(1),-1.0_pReal,1.0_pReal)))] endif - + end if end function qu2ro @@ -511,12 +522,12 @@ pure function qu2ho(qu) result(ho) real(pReal), intent(in), dimension(4) :: qu real(pReal), dimension(3) :: ho - + real(pReal) :: omega, f 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 ] else ho = qu(2:4) @@ -532,7 +543,7 @@ end function qu2ho !> @brief convert unit quaternion to cubochoric !--------------------------------------------------------------------------------------------------- pure function qu2cu(qu) result(cu) - + real(pReal), intent(in), dimension(4) :: qu real(pReal), dimension(3) :: cu @@ -565,18 +576,18 @@ pure function om2eu(om) result(eu) real(pReal), intent(in), dimension(3,3) :: om real(pReal), dimension(3) :: eu real(pReal) :: zeta - + if (abs(om(3,3)) < 1.0_pReal) then zeta = 1.0_pReal/sqrt(1.0_pReal-om(3,3)**2.0_pReal) eu = [atan2(om(3,1)*zeta,-om(3,2)*zeta), & acos(om(3,3)), & atan2(om(1,3)*zeta, om(2,3)*zeta)] - else + else eu = [atan2(om(1,2),om(1,1)), 0.5_pReal*PI*(1.0_pReal-om(3,3)),0.0_pReal ] end if where(eu<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI]) - + end function om2eu @@ -588,19 +599,19 @@ function om2ax(om) result(ax) real(pReal), intent(in), dimension(3,3) :: om real(pReal), dimension(4) :: ax - + real(pReal) :: t real(pReal), dimension(3) :: Wr, Wi real(pReal), dimension((64+2)*3) :: work real(pReal), dimension(3,3) :: VR, devNull, om_ integer :: ierr, i - + om_ = om - + ! first get the rotation angle t = 0.5_pReal * (math_trace33(om) - 1.0_pReal) ax(4) = acos(math_clip(t,-1.0_pReal,1.0_pReal)) - + if (dEq0(ax(4))) then ax(1:3) = [ 0.0_pReal, 0.0_pReal, 1.0_pReal ] else @@ -674,7 +685,7 @@ pure function eu2qu(eu) result(qu) real(pReal) :: cPhi, sPhi ee = 0.5_pReal*eu - + cPhi = cos(ee(2)) sPhi = sin(ee(2)) @@ -692,15 +703,15 @@ end function eu2qu !> @brief Euler angles to orientation matrix !--------------------------------------------------------------------------------------------------- pure function eu2om(eu) result(om) - + real(pReal), intent(in), dimension(3) :: eu real(pReal), dimension(3,3) :: om - - real(pReal), dimension(3) :: c, s - + + real(pReal), dimension(3) :: c, s + c = cos(eu) s = sin(eu) - + om(1,1) = c(1)*c(3)-s(1)*s(3)*c(2) om(1,2) = s(1)*c(3)+c(1)*s(3)*c(2) om(1,3) = s(3)*s(2) @@ -710,7 +721,7 @@ pure function eu2om(eu) result(om) om(3,1) = s(1)*s(2) om(3,2) = -c(1)*s(2) om(3,3) = c(2) - + where(dEq0(om)) om = 0.0_pReal end function eu2om @@ -721,19 +732,19 @@ end function eu2om !> @brief convert euler to axis angle !--------------------------------------------------------------------------------------------------- pure function eu2ax(eu) result(ax) - + real(pReal), intent(in), dimension(3) :: eu real(pReal), dimension(4) :: ax - + real(pReal) :: t, delta, tau, alpha, sigma - + t = tan(eu(2)*0.5_pReal) sigma = 0.5_pReal*(eu(1)+eu(3)) delta = 0.5_pReal*(eu(1)-eu(3)) tau = sqrt(t**2+sin(sigma)**2) - + alpha = merge(PI, 2.0_pReal*atan(tau/cos(sigma)), dEq(sigma,PI*0.5_pReal,tol=1.0e-15_pReal)) - + if (dEq0(alpha)) then ! return a default identity axis-angle pair ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ] else @@ -741,7 +752,7 @@ pure function eu2ax(eu) result(ax) ax(4) = alpha if (alpha < 0.0_pReal) ax = -ax ! ensure alpha is positive end if - + end function eu2ax @@ -753,7 +764,7 @@ pure function eu2ro(eu) result(ro) real(pReal), intent(in), dimension(3) :: eu real(pReal), dimension(4) :: ro - + ro = eu2ax(eu) if (ro(4) >= PI) then ro(4) = IEEE_value(ro(4),IEEE_positive_inf) @@ -762,7 +773,7 @@ pure function eu2ro(eu) result(ro) else ro(4) = tan(ro(4)*0.5_pReal) end if - + end function eu2ro @@ -799,7 +810,7 @@ end function eu2cu !> @brief convert axis angle pair to quaternion !--------------------------------------------------------------------------------------------------- pure function ax2qu(ax) result(qu) - + real(pReal), intent(in), dimension(4) :: ax real(pReal), dimension(4) :: qu @@ -825,7 +836,7 @@ pure function ax2om(ax) result(om) real(pReal), intent(in), dimension(4) :: ax real(pReal), dimension(3,3) :: om - + real(pReal) :: q, c, s, omc c = cos(ax(4)) @@ -839,11 +850,11 @@ pure function ax2om(ax) result(om) q = omc*ax(1)*ax(2) om(1,2) = q + s*ax(3) om(2,1) = q - s*ax(3) - + q = omc*ax(2)*ax(3) om(2,3) = q + s*ax(1) om(3,2) = q - s*ax(1) - + q = omc*ax(3)*ax(1) om(3,1) = q + s*ax(2) om(1,3) = q - s*ax(2) @@ -875,12 +886,12 @@ pure function ax2ro(ax) result(ro) real(pReal), intent(in), dimension(4) :: ax real(pReal), dimension(4) :: ro - + real(pReal), parameter :: thr = 1.0e-7_pReal - + if (dEq0(ax(4))) then ro = [ 0.0_pReal, 0.0_pReal, P, 0.0_pReal ] - else + else ro(1:3) = ax(1:3) ! we need to deal with the 180 degree case ro(4) = merge(IEEE_value(ro(4),IEEE_positive_inf),tan(ax(4)*0.5_pReal),abs(ax(4)-PI) < thr) @@ -897,9 +908,9 @@ pure function ax2ho(ax) result(ho) real(pReal), intent(in), dimension(4) :: ax real(pReal), dimension(3) :: ho - + real(pReal) :: f - + f = 0.75_pReal * ( ax(4) - sin(ax(4)) ) f = f**(1.0_pReal/3.0_pReal) ho = ax(1:3) * f @@ -916,7 +927,7 @@ function ax2cu(ax) result(cu) real(pReal), intent(in), dimension(4) :: ax real(pReal), dimension(3) :: cu - cu = ho2cu(ax2ho(ax)) + cu = ho2cu(ax2ho(ax)) end function ax2cu @@ -929,7 +940,7 @@ pure function ro2qu(ro) result(qu) real(pReal), intent(in), dimension(4) :: ro real(pReal), dimension(4) :: qu - + qu = ax2qu(ro2ax(ro)) end function ro2qu @@ -957,7 +968,7 @@ pure function ro2eu(ro) result(eu) real(pReal), intent(in), dimension(4) :: ro real(pReal), dimension(3) :: eu - + eu = om2eu(ro2om(ro)) end function ro2eu @@ -971,14 +982,14 @@ pure function ro2ax(ro) result(ax) real(pReal), intent(in), dimension(4) :: ro real(pReal), dimension(4) :: ax - + real(pReal) :: ta, angle - + ta = ro(4) - + if (.not. IEEE_is_finite(ta)) then ax = [ ro(1), ro(2), ro(3), PI ] - elseif (dEq0(ta)) then + elseif (dEq0(ta)) then ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ] else angle = 2.0_pReal*atan(ta) @@ -997,9 +1008,9 @@ pure function ro2ho(ro) result(ho) real(pReal), intent(in), dimension(4) :: ro real(pReal), dimension(3) :: ho - + real(pReal) :: f - + if (dEq0(norm2(ro(1:3)))) then ho = [ 0.0_pReal, 0.0_pReal, 0.0_pReal ] else @@ -1074,26 +1085,26 @@ pure function ho2ax(ho) result(ax) real(pReal), intent(in), dimension(3) :: ho real(pReal), dimension(4) :: ax - + integer :: i real(pReal) :: hmag_squared, s, hm real(pReal), parameter, dimension(16) :: & - tfit = [ 1.0000000000018852_pReal, -0.5000000002194847_pReal, & - -0.024999992127593126_pReal, -0.003928701544781374_pReal, & - -0.0008152701535450438_pReal, -0.0002009500426119712_pReal, & - -0.00002397986776071756_pReal, -0.00008202868926605841_pReal, & - +0.00012448715042090092_pReal, -0.0001749114214822577_pReal, & - +0.0001703481934140054_pReal, -0.00012062065004116828_pReal, & - +0.000059719705868660826_pReal, -0.00001980756723965647_pReal, & + tfit = [ 1.0000000000018852_pReal, -0.5000000002194847_pReal, & + -0.024999992127593126_pReal, -0.003928701544781374_pReal, & + -0.0008152701535450438_pReal, -0.0002009500426119712_pReal, & + -0.00002397986776071756_pReal, -0.00008202868926605841_pReal, & + +0.00012448715042090092_pReal, -0.0001749114214822577_pReal, & + +0.0001703481934140054_pReal, -0.00012062065004116828_pReal, & + +0.000059719705868660826_pReal, -0.00001980756723965647_pReal, & +0.000003953714684212874_pReal, -0.00000036555001439719544_pReal ] - + ! normalize h and store the magnitude hmag_squared = sum(ho**2.0_pReal) if (dEq0(hmag_squared)) then ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ] else hm = hmag_squared - + ! convert the magnitude to the rotation angle s = tfit(1) + tfit(2) * hmag_squared do i=3,16 @@ -1120,16 +1131,56 @@ pure function ho2ro(ho) result(ro) end function ho2ro -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief convert homochoric to cubochoric -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------- pure function ho2cu(ho) result(cu) - real(pReal), intent(in), dimension(3) :: ho - real(pReal), dimension(3) :: cu + real(pReal), intent(in), dimension(3) :: ho + 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 @@ -1204,25 +1255,93 @@ pure function cu2ro(cu) result(ro) end function cu2ro -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------- !> @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) - real(pReal), intent(in), dimension(3) :: cu - real(pReal), dimension(3) :: ho + real(pReal), intent(in), dimension(3) :: cu + 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 +!-------------------------------------------------------------------------- +!> @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 !-------------------------------------------------------------------------------------------------- subroutine unitTest - + type(rotation) :: R real(pReal), dimension(4) :: qu, ax, ro real(pReal), dimension(3) :: x, eu, ho, v3 @@ -1233,7 +1352,7 @@ subroutine unitTest integer :: i do i=1,10 - + msg = '' #if defined(__GFORTRAN__) && __GNUC__<9 @@ -1307,15 +1426,15 @@ subroutine unitTest #endif call R%fromMatrix(om) - + call random_number(v3) if(all(dNeq(R%rotVector(R%rotVector(v3),active=.true.),v3,1.0e-12_pReal))) & msg = trim(msg)//'rotVector,' - + call random_number(t33) if(all(dNeq(R%rotTensor2(R%rotTensor2(t33),active=.true.),t33,1.0e-12_pReal))) & msg = trim(msg)//'rotTensor2,' - + call random_number(t3333) if(all(dNeq(R%rotTensor4(R%rotTensor4(t3333),active=.true.),t3333,1.0e-12_pReal))) & msg = trim(msg)//'rotTensor4,'