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

Merge branch 'python-style' into 'development' 

Python style

See merge request damask/DAMASK!132
This commit is contained in:
Karo 2020-03-06 00:42:48 +01:00
parent c06ef00fe2 7a0c20b6de
commit 31136100c6
11 changed files with 2220 additions and 2157 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

100
python/damask/Lambert.py
View File

@ -52,37 +52,37 @@ def CubeToBall(cube):
""" """
if np.abs(np.max(cube))>np.pi**(2./3.) * 0.5: if np.abs(np.max(cube))>np.pi**(2./3.) * 0.5:
raise ValueError raise ValueError
# transform to the sphere grid via the curved square, and intercept the zero point # 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-300): if np.allclose(cube,0.0,rtol=0.0,atol=1.0e-300):
ball = np.zeros(3) ball = np.zeros(3)
else: else:
# get pyramide and scale by grid parameter ratio # get pyramide and scale by grid parameter ratio
p = get_order(cube) p = get_order(cube)
XYZ = cube[p] * sc XYZ = cube[p] * sc
# intercept all the points along the z-axis # intercept all the points along the z-axis
if np.allclose(XYZ[0:2],0.0,rtol=0.0,atol=1.0e-300): if np.allclose(XYZ[0:2],0.0,rtol=0.0,atol=1.0e-300):
ball = np.array([0.0, 0.0, np.sqrt(6.0/np.pi) * XYZ[2]]) ball = np.array([0.0, 0.0, np.sqrt(6.0/np.pi) * XYZ[2]])
else: else:
order = [1,0] if np.abs(XYZ[1]) <= np.abs(XYZ[0]) else [0,1] 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]] q = np.pi/12.0 * XYZ[order[0]]/XYZ[order[1]]
c = np.cos(q) c = np.cos(q)
s = np.sin(q) s = np.sin(q)
q = R1*2.0**0.25/beta * XYZ[order[1]] / np.sqrt(np.sqrt(2.0)-c) 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 T = np.array([ (np.sqrt(2.0)*c - 1.0), np.sqrt(2.0) * s]) * q
# transform to sphere grid (inverse Lambert) # 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 # note that there is no need to worry about dividing by zero, since XYZ[2] can not become zero
c = np.sum(T**2) c = np.sum(T**2)
s = c * np.pi/24.0 /XYZ[2]**2 s = c * np.pi/24.0 /XYZ[2]**2
c = c * np.sqrt(np.pi/24.0)/XYZ[2] c = c * np.sqrt(np.pi/24.0)/XYZ[2]
q = np.sqrt( 1.0 - s ) 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 ]) 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 # reverse the coordinates back to the regular order according to the original pyramid number
ball = ball[p] ball = ball[p]
return ball return ball
@ -104,39 +104,39 @@ def BallToCube(ball):
""" """
rs = np.linalg.norm(ball) rs = np.linalg.norm(ball)
if rs > R1: if rs > R1:
raise ValueError raise ValueError
if np.allclose(ball,0.0,rtol=0.0,atol=1.0e-300): if np.allclose(ball,0.0,rtol=0.0,atol=1.0e-300):
cube = np.zeros(3) cube = np.zeros(3)
else: else:
p = get_order(ball) p = get_order(ball)
xyz3 = ball[p] xyz3 = ball[p]
# inverse M_3 # inverse M_3
xyz2 = xyz3[0:2] * np.sqrt( 2.0*rs/(rs+np.abs(xyz3[2])) ) xyz2 = xyz3[0:2] * np.sqrt( 2.0*rs/(rs+np.abs(xyz3[2])) )
# inverse M_2 # inverse M_2
qxy = np.sum(xyz2**2) qxy = np.sum(xyz2**2)
if np.isclose(qxy,0.0,rtol=0.0,atol=1.0e-300): if np.isclose(qxy,0.0,rtol=0.0,atol=1.0e-300):
Tinv = np.zeros(2) Tinv = np.zeros(2)
else: else:
q2 = qxy + np.max(np.abs(xyz2))**2 q2 = qxy + np.max(np.abs(xyz2))**2
sq2 = np.sqrt(q2) sq2 = np.sqrt(q2)
q = (beta/np.sqrt(2.0)/R1) * np.sqrt(q2*qxy/(q2-np.max(np.abs(xyz2))*sq2)) 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) 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 \ 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]) np.array([np.arccos(tt)/np.pi*12.0,1.0])
Tinv = q * np.where(xyz2<0.0,-Tinv,Tinv) Tinv = q * np.where(xyz2<0.0,-Tinv,Tinv)
# inverse M_1 # 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 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
# reverse the coordinates back to the regular order according to the original pyramid number cube = cube[p]
cube = cube[p]
return cube return cube
def get_order(xyz): def get_order(xyz):
""" """
Get order of the coordinates. Get order of the coordinates.
@ -157,10 +157,10 @@ def get_order(xyz):
""" """
if (abs(xyz[0])<= xyz[2]) and (abs(xyz[1])<= xyz[2]) or \ if (abs(xyz[0])<= xyz[2]) and (abs(xyz[1])<= xyz[2]) or \
(abs(xyz[0])<=-xyz[2]) and (abs(xyz[1])<=-xyz[2]): (abs(xyz[0])<=-xyz[2]) and (abs(xyz[1])<=-xyz[2]):
return [0,1,2] return [0,1,2]
elif (abs(xyz[2])<= xyz[0]) and (abs(xyz[1])<= xyz[0]) or \ elif (abs(xyz[2])<= xyz[0]) and (abs(xyz[1])<= xyz[0]) or \
(abs(xyz[2])<=-xyz[0]) and (abs(xyz[1])<=-xyz[0]): (abs(xyz[2])<=-xyz[0]) and (abs(xyz[1])<=-xyz[0]):
return [1,2,0] return [1,2,0]
elif (abs(xyz[0])<= xyz[1]) and (abs(xyz[2])<= xyz[1]) or \ elif (abs(xyz[0])<= xyz[1]) and (abs(xyz[2])<= xyz[1]) or \
(abs(xyz[0])<=-xyz[1]) and (abs(xyz[2])<=-xyz[1]): (abs(xyz[0])<=-xyz[1]) and (abs(xyz[2])<=-xyz[1]):
return [2,0,1] return [2,0,1]

4
python/damask/__init__.py
View File

@ -13,7 +13,9 @@ from .asciitable import ASCIItable # noqa
from .config import Material # noqa from .config import Material # noqa
from .colormaps import Colormap, Color # noqa from .colormaps import Colormap, Color # noqa
from .orientation import Symmetry, Lattice, Rotation, Orientation # noqa from .rotation import Rotation # noqa
from .lattice import Symmetry, Lattice# noqa
from .orientation import Orientation # noqa
from .result import Result # noqa from .result import Result # noqa
from .result import Result as DADF5 # noqa from .result import Result as DADF5 # noqa

550
python/damask/colormaps.py
View File

@ -1,359 +1,355 @@
import numpy as np import numpy as np
class Color(): class Color():
"""Color representation in and conversion between different color-spaces.""" """Color representation in and conversion between different color-spaces."""
__slots__ = [ __slots__ = [
'model', 'model',
'color', 'color',
'__dict__', '__dict__',
] ]
# ------------------------------------------------------------------ def __init__(self,
def __init__(self, model = 'RGB',
model = 'RGB', color = np.zeros(3,'d')):
color = np.zeros(3,'d')): """
""" Create a Color object.
Create a Color object.
Parameters Parameters
---------- ----------
model : string model : string
color model color model
color : numpy.ndarray color : numpy.ndarray
vector representing the color according to the selected model vector representing the color according to the selected model
""" """
self.__transforms__ = \ self.__transforms__ = \
{'HSV': {'index': 0, 'next': self._HSV2HSL}, {'HSV': {'index': 0, 'next': self._HSV2HSL},
'HSL': {'index': 1, 'next': self._HSL2RGB, 'prev': self._HSL2HSV}, 'HSL': {'index': 1, 'next': self._HSL2RGB, 'prev': self._HSL2HSV},
'RGB': {'index': 2, 'next': self._RGB2XYZ, 'prev': self._RGB2HSL}, 'RGB': {'index': 2, 'next': self._RGB2XYZ, 'prev': self._RGB2HSL},
'XYZ': {'index': 3, 'next': self._XYZ2CIELAB, 'prev': self._XYZ2RGB}, 'XYZ': {'index': 3, 'next': self._XYZ2CIELAB, 'prev': self._XYZ2RGB},
'CIELAB': {'index': 4, 'next': self._CIELAB2MSH, 'prev': self._CIELAB2XYZ}, 'CIELAB': {'index': 4, 'next': self._CIELAB2MSH, 'prev': self._CIELAB2XYZ},
'MSH': {'index': 5, 'prev': self._MSH2CIELAB}, 'MSH': {'index': 5, 'prev': self._MSH2CIELAB},
} }
model = model.upper() model = model.upper()
if model not in list(self.__transforms__.keys()): model = 'RGB' if model not in list(self.__transforms__.keys()): model = 'RGB'
if model == 'RGB' and max(color) > 1.0: # are we RGB255 ? if model == 'RGB' and max(color) > 1.0: # are we RGB255 ?
for i in range(3): for i in range(3):
color[i] /= 255.0 # rescale to RGB color[i] /= 255.0 # rescale to RGB
if model == 'HSL': # are we HSL ? if model == 'HSL': # are we HSL ?
if abs(color[0]) > 1.0: color[0] /= 360.0 # with angular hue? if abs(color[0]) > 1.0: color[0] /= 360.0 # with angular hue?
while color[0] >= 1.0: color[0] -= 1.0 # rewind to proper range while color[0] >= 1.0: color[0] -= 1.0 # rewind to proper range
while color[0] < 0.0: color[0] += 1.0 # rewind to proper range while color[0] < 0.0: color[0] += 1.0 # rewind to proper range
self.model = model self.model = model
self.color = np.array(color,'d') self.color = np.array(color,'d')
# ------------------------------------------------------------------ def __repr__(self):
def __repr__(self): """Color model and values."""
"""Color model and values.""" return 'Model: %s Color: %s'%(self.model,str(self.color))
return 'Model: %s Color: %s'%(self.model,str(self.color))
# ------------------------------------------------------------------ def __str__(self):
def __str__(self): """Color model and values."""
"""Color model and values.""" return self.__repr__()
return self.__repr__()
# ------------------------------------------------------------------ def convert_to(self,toModel = 'RGB'):
def convertTo(self,toModel = 'RGB'): """
""" Change the color model permanently.
Change the color model permanently.
Parameters Parameters
---------- ----------
toModel : string toModel : string
color model color model
""" """
toModel = toModel.upper() toModel = toModel.upper()
if toModel not in list(self.__transforms__.keys()): return if toModel not in list(self.__transforms__.keys()): return
sourcePos = self.__transforms__[self.model]['index'] sourcePos = self.__transforms__[self.model]['index']
targetPos = self.__transforms__[toModel]['index'] targetPos = self.__transforms__[toModel]['index']
while sourcePos < targetPos: while sourcePos < targetPos:
self.__transforms__[self.model]['next']() self.__transforms__[self.model]['next']()
sourcePos += 1 sourcePos += 1
while sourcePos > targetPos: while sourcePos > targetPos:
self.__transforms__[self.model]['prev']() self.__transforms__[self.model]['prev']()
sourcePos -= 1 sourcePos -= 1
return self return self
# ------------------------------------------------------------------ def express_as(self,asModel = 'RGB'):
def expressAs(self,asModel = 'RGB'): """
""" Return the color in a different model.
Return the color in a different model.
Parameters Parameters
---------- ----------
asModel : string asModel : string
color model color model
""" """
return self.__class__(self.model,self.color).convertTo(asModel) return self.__class__(self.model,self.color).convert_to(asModel)
def _HSV2HSL(self): def _HSV2HSL(self):
""" """
Convert H(ue) S(aturation) V(alue or brightness) to H(ue) S(aturation) L(uminance). Convert H(ue) S(aturation) V(alue or brightness) to H(ue) S(aturation) L(uminance).
All values are in the range [0,1] All values are in the range [0,1]
http://codeitdown.com/hsl-hsb-hsv-color http://codeitdown.com/hsl-hsb-hsv-color
""" """
if self.model != 'HSV': return if self.model != 'HSV': return
converted = Color('HSL',np.array([ converted = Color('HSL',np.array([
self.color[0], self.color[0],
1. if self.color[2] == 0.0 or (self.color[1] == 0.0 and self.color[2] == 1.0) \ 1. if self.color[2] == 0.0 or (self.color[1] == 0.0 and self.color[2] == 1.0) \
else self.color[1]*self.color[2]/(1.-abs(self.color[2]*(2.-self.color[1])-1.)), else self.color[1]*self.color[2]/(1.-abs(self.color[2]*(2.-self.color[1])-1.)),
0.5*self.color[2]*(2.-self.color[1]), 0.5*self.color[2]*(2.-self.color[1]),
])) ]))
self.model = converted.model self.model = converted.model
self.color = converted.color self.color = converted.color
def _HSL2HSV(self): def _HSL2HSV(self):
""" """
Convert H(ue) S(aturation) L(uminance) to H(ue) S(aturation) V(alue or brightness). Convert H(ue) S(aturation) L(uminance) to H(ue) S(aturation) V(alue or brightness).
All values are in the range [0,1] All values are in the range [0,1]
http://codeitdown.com/hsl-hsb-hsv-color http://codeitdown.com/hsl-hsb-hsv-color
""" """
if self.model != 'HSL': return if self.model != 'HSL': return
h = self.color[0] h = self.color[0]
b = self.color[2]+0.5*(self.color[1]*(1.-abs(2*self.color[2]-1))) b = self.color[2]+0.5*(self.color[1]*(1.-abs(2*self.color[2]-1)))
s = 1.0 if b == 0.0 else 2.*(b-self.color[2])/b s = 1.0 if b == 0.0 else 2.*(b-self.color[2])/b
converted = Color('HSV',np.array([h,s,b])) converted = Color('HSV',np.array([h,s,b]))
self.model = converted.model self.model = converted.model
self.color = converted.color self.color = converted.color
def _HSL2RGB(self): def _HSL2RGB(self):
""" """
Convert H(ue) S(aturation) L(uminance) to R(red) G(reen) B(lue). Convert H(ue) S(aturation) L(uminance) to R(red) G(reen) B(lue).
All values are in the range [0,1] All values are in the range [0,1]
from http://en.wikipedia.org/wiki/HSL_and_HSV from http://en.wikipedia.org/wiki/HSL_and_HSV
""" """
if self.model != 'HSL': return if self.model != 'HSL': return
sextant = self.color[0]*6.0 sextant = self.color[0]*6.0
c = (1.0 - abs(2.0 * self.color[2] - 1.0))*self.color[1] c = (1.0 - abs(2.0 * self.color[2] - 1.0))*self.color[1]
x = c*(1.0 - abs(sextant%2 - 1.0)) x = c*(1.0 - abs(sextant%2 - 1.0))
m = self.color[2] - 0.5*c m = self.color[2] - 0.5*c
converted = Color('RGB',np.array([ converted = Color('RGB',np.array([
[c+m, x+m, m], [c+m, x+m, m],
[x+m, c+m, m], [x+m, c+m, m],
[m, c+m, x+m], [m, c+m, x+m],
[m, x+m, c+m], [m, x+m, c+m],
[x+m, m, c+m], [x+m, m, c+m],
[c+m, m, x+m], [c+m, m, x+m],
][int(sextant)],'d')) ][int(sextant)],'d'))
self.model = converted.model self.model = converted.model
self.color = converted.color self.color = converted.color
def _RGB2HSL(self): def _RGB2HSL(self):
""" """
Convert R(ed) G(reen) B(lue) to H(ue) S(aturation) L(uminance). Convert R(ed) G(reen) B(lue) to H(ue) S(aturation) L(uminance).
All values are in the range [0,1] All values are in the range [0,1]
from http://130.113.54.154/~monger/hsl-rgb.html from http://130.113.54.154/~monger/hsl-rgb.html
""" """
if self.model != 'RGB': return if self.model != 'RGB': return
HSL = np.zeros(3,'d') HSL = np.zeros(3,'d')
maxcolor = self.color.max() maxcolor = self.color.max()
mincolor = self.color.min() mincolor = self.color.min()
HSL[2] = (maxcolor + mincolor)/2.0 HSL[2] = (maxcolor + mincolor)/2.0
if(mincolor == maxcolor): if(mincolor == maxcolor):
HSL[0] = 0.0 HSL[0] = 0.0
HSL[1] = 0.0 HSL[1] = 0.0
else:
if (HSL[2]<0.5):
HSL[1] = (maxcolor - mincolor)/(maxcolor + mincolor)
else: else:
HSL[1] = (maxcolor - mincolor)/(2.0 - maxcolor - mincolor) if (HSL[2]<0.5):
if (maxcolor == self.color[0]): HSL[1] = (maxcolor - mincolor)/(maxcolor + mincolor)
HSL[0] = 0.0 + (self.color[1] - self.color[2])/(maxcolor - mincolor) else:
elif (maxcolor == self.color[1]): HSL[1] = (maxcolor - mincolor)/(2.0 - maxcolor - mincolor)
HSL[0] = 2.0 + (self.color[2] - self.color[0])/(maxcolor - mincolor) if (maxcolor == self.color[0]):
elif (maxcolor == self.color[2]): HSL[0] = 0.0 + (self.color[1] - self.color[2])/(maxcolor - mincolor)
HSL[0] = 4.0 + (self.color[0] - self.color[1])/(maxcolor - mincolor) elif (maxcolor == self.color[1]):
HSL[0] = HSL[0]*60.0 # scaling to 360 might be dangerous for small values HSL[0] = 2.0 + (self.color[2] - self.color[0])/(maxcolor - mincolor)
if (HSL[0] < 0.0): elif (maxcolor == self.color[2]):
HSL[0] = HSL[0] + 360.0 HSL[0] = 4.0 + (self.color[0] - self.color[1])/(maxcolor - mincolor)
for i in range(2): HSL[0] = HSL[0]*60.0 # scaling to 360 might be dangerous for small values
HSL[i+1] = min(HSL[i+1],1.0) if (HSL[0] < 0.0):
HSL[i+1] = max(HSL[i+1],0.0) HSL[0] = HSL[0] + 360.0
for i in range(2):
HSL[i+1] = min(HSL[i+1],1.0)
HSL[i+1] = max(HSL[i+1],0.0)
converted = Color('HSL', HSL) converted = Color('HSL', HSL)
self.model = converted.model self.model = converted.model
self.color = converted.color self.color = converted.color
def _RGB2XYZ(self): def _RGB2XYZ(self):
""" """
Convert R(ed) G(reen) B(lue) to CIE XYZ. Convert R(ed) G(reen) B(lue) to CIE XYZ.
All values are in the range [0,1] All values are in the range [0,1]
from http://www.cs.rit.edu/~ncs/color/t_convert.html from http://www.cs.rit.edu/~ncs/color/t_convert.html
""" """
if self.model != 'RGB': return if self.model != 'RGB': return
XYZ = np.zeros(3,'d') XYZ = np.zeros(3,'d')
RGB_lin = np.zeros(3,'d') RGB_lin = np.zeros(3,'d')
convert = np.array([[0.412453,0.357580,0.180423], convert = np.array([[0.412453,0.357580,0.180423],
[0.212671,0.715160,0.072169], [0.212671,0.715160,0.072169],
[0.019334,0.119193,0.950227]]) [0.019334,0.119193,0.950227]])
for i in range(3): for i in range(3):
if (self.color[i] > 0.04045): RGB_lin[i] = ((self.color[i]+0.0555)/1.0555)**2.4 if (self.color[i] > 0.04045): RGB_lin[i] = ((self.color[i]+0.0555)/1.0555)**2.4
else: RGB_lin[i] = self.color[i] /12.92 else: RGB_lin[i] = self.color[i] /12.92
XYZ = np.dot(convert,RGB_lin) XYZ = np.dot(convert,RGB_lin)
for i in range(3): for i in range(3):
XYZ[i] = max(XYZ[i],0.0) XYZ[i] = max(XYZ[i],0.0)
converted = Color('XYZ', XYZ) converted = Color('XYZ', XYZ)
self.model = converted.model self.model = converted.model
self.color = converted.color self.color = converted.color
def _XYZ2RGB(self): def _XYZ2RGB(self):
""" """
Convert CIE XYZ to R(ed) G(reen) B(lue). Convert CIE XYZ to R(ed) G(reen) B(lue).
All values are in the range [0,1] All values are in the range [0,1]
from http://www.cs.rit.edu/~ncs/color/t_convert.html from http://www.cs.rit.edu/~ncs/color/t_convert.html
""" """
if self.model != 'XYZ': if self.model != 'XYZ':
return return
convert = np.array([[ 3.240479,-1.537150,-0.498535], convert = np.array([[ 3.240479,-1.537150,-0.498535],
[-0.969256, 1.875992, 0.041556], [-0.969256, 1.875992, 0.041556],
[ 0.055648,-0.204043, 1.057311]]) [ 0.055648,-0.204043, 1.057311]])
RGB_lin = np.dot(convert,self.color) RGB_lin = np.dot(convert,self.color)
RGB = np.zeros(3,'d') RGB = np.zeros(3,'d')
for i in range(3): for i in range(3):
if (RGB_lin[i] > 0.0031308): RGB[i] = ((RGB_lin[i])**(1.0/2.4))*1.0555-0.0555 if (RGB_lin[i] > 0.0031308): RGB[i] = ((RGB_lin[i])**(1.0/2.4))*1.0555-0.0555
else: RGB[i] = RGB_lin[i] *12.92 else: RGB[i] = RGB_lin[i] *12.92
for i in range(3): for i in range(3):
RGB[i] = min(RGB[i],1.0) RGB[i] = min(RGB[i],1.0)
RGB[i] = max(RGB[i],0.0) RGB[i] = max(RGB[i],0.0)
maxVal = max(RGB) # clipping colors according to the display gamut maxVal = max(RGB) # clipping colors according to the display gamut
if (maxVal > 1.0): RGB /= maxVal if (maxVal > 1.0): RGB /= maxVal
converted = Color('RGB', RGB) converted = Color('RGB', RGB)
self.model = converted.model self.model = converted.model
self.color = converted.color self.color = converted.color
def _CIELAB2XYZ(self): def _CIELAB2XYZ(self):
""" """
Convert CIE Lab to CIE XYZ. Convert CIE Lab to CIE XYZ.
All values are in the range [0,1] All values are in the range [0,1]
from http://www.easyrgb.com/index.php?X=MATH&H=07#text7 from http://www.easyrgb.com/index.php?X=MATH&H=07#text7
""" """
if self.model != 'CIELAB': return if self.model != 'CIELAB': return
ref_white = np.array([.95047, 1.00000, 1.08883],'d') # Observer = 2, Illuminant = D65 ref_white = np.array([.95047, 1.00000, 1.08883],'d') # Observer = 2, Illuminant = D65
XYZ = np.zeros(3,'d') XYZ = np.zeros(3,'d')
XYZ[1] = (self.color[0] + 16.0 ) / 116.0 XYZ[1] = (self.color[0] + 16.0 ) / 116.0
XYZ[0] = XYZ[1] + self.color[1]/ 500.0 XYZ[0] = XYZ[1] + self.color[1]/ 500.0
XYZ[2] = XYZ[1] - self.color[2]/ 200.0 XYZ[2] = XYZ[1] - self.color[2]/ 200.0
for i in range(len(XYZ)): for i in range(len(XYZ)):
if (XYZ[i] > 6./29. ): XYZ[i] = XYZ[i]**3. if (XYZ[i] > 6./29. ): XYZ[i] = XYZ[i]**3.
else: XYZ[i] = 108./841. * (XYZ[i] - 4./29.) else: XYZ[i] = 108./841. * (XYZ[i] - 4./29.)
converted = Color('XYZ', XYZ*ref_white) converted = Color('XYZ', XYZ*ref_white)
self.model = converted.model self.model = converted.model
self.color = converted.color self.color = converted.color
def _XYZ2CIELAB(self):
"""
Convert CIE XYZ to CIE Lab.
All values are in the range [0,1]
from http://en.wikipedia.org/wiki/Lab_color_space,
http://www.cs.rit.edu/~ncs/color/t_convert.html
"""
if self.model != 'XYZ': return
ref_white = np.array([.95047, 1.00000, 1.08883],'d') # Observer = 2, Illuminant = D65
XYZ = self.color/ref_white
for i in range(len(XYZ)):
if (XYZ[i] > 216./24389 ): XYZ[i] = XYZ[i]**(1.0/3.0)
else: XYZ[i] = (841./108. * XYZ[i]) + 16.0/116.0
converted = Color('CIELAB', np.array([ 116.0 * XYZ[1] - 16.0,
500.0 * (XYZ[0] - XYZ[1]),
200.0 * (XYZ[1] - XYZ[2]) ]))
self.model = converted.model
self.color = converted.color
def _CIELAB2MSH(self): def _XYZ2CIELAB(self):
""" """
Convert CIE Lab to Msh colorspace. Convert CIE XYZ to CIE Lab.
from http://www.cs.unm.edu/~kmorel/documents/ColorMaps/DivergingColorMapWorkshop.xls All values are in the range [0,1]
""" from http://en.wikipedia.org/wiki/Lab_color_space,
if self.model != 'CIELAB': return http://www.cs.rit.edu/~ncs/color/t_convert.html
"""
if self.model != 'XYZ': return
Msh = np.zeros(3,'d') ref_white = np.array([.95047, 1.00000, 1.08883],'d') # Observer = 2, Illuminant = D65
Msh[0] = np.sqrt(np.dot(self.color,self.color)) XYZ = self.color/ref_white
if (Msh[0] > 0.001):
Msh[1] = np.arccos(self.color[0]/Msh[0])
if (self.color[1] != 0.0):
Msh[2] = np.arctan2(self.color[2],self.color[1])
converted = Color('MSH', Msh) for i in range(len(XYZ)):
self.model = converted.model if (XYZ[i] > 216./24389 ): XYZ[i] = XYZ[i]**(1.0/3.0)
self.color = converted.color else: XYZ[i] = (841./108. * XYZ[i]) + 16.0/116.0
converted = Color('CIELAB', np.array([ 116.0 * XYZ[1] - 16.0,
500.0 * (XYZ[0] - XYZ[1]),
200.0 * (XYZ[1] - XYZ[2]) ]))
self.model = converted.model
self.color = converted.color
def _MSH2CIELAB(self): def _CIELAB2MSH(self):
""" """
Convert Msh colorspace to CIE Lab. Convert CIE Lab to Msh colorspace.
with s,h in radians from http://www.cs.unm.edu/~kmorel/documents/ColorMaps/DivergingColorMapWorkshop.xls
from http://www.cs.unm.edu/~kmorel/documents/ColorMaps/DivergingColorMapWorkshop.xls """
""" if self.model != 'CIELAB': return
if self.model != 'MSH': return
Lab = np.zeros(3,'d') Msh = np.zeros(3,'d')
Lab[0] = self.color[0] * np.cos(self.color[1]) Msh[0] = np.sqrt(np.dot(self.color,self.color))
Lab[1] = self.color[0] * np.sin(self.color[1]) * np.cos(self.color[2]) if (Msh[0] > 0.001):
Lab[2] = self.color[0] * np.sin(self.color[1]) * np.sin(self.color[2]) Msh[1] = np.arccos(self.color[0]/Msh[0])
if (self.color[1] != 0.0):
Msh[2] = np.arctan2(self.color[2],self.color[1])
converted = Color('CIELAB', Lab) converted = Color('MSH', Msh)
self.model = converted.model self.model = converted.model
self.color = converted.color self.color = converted.color
def _MSH2CIELAB(self):
"""
Convert Msh colorspace to CIE Lab.
with s,h in radians
from http://www.cs.unm.edu/~kmorel/documents/ColorMaps/DivergingColorMapWorkshop.xls
"""
if self.model != 'MSH': return
Lab = np.zeros(3,'d')
Lab[0] = self.color[0] * np.cos(self.color[1])
Lab[1] = self.color[0] * np.sin(self.color[1]) * np.cos(self.color[2])
Lab[2] = self.color[0] * np.sin(self.color[1]) * np.sin(self.color[2])
converted = Color('CIELAB', Lab)
self.model = converted.model
self.color = converted.color
class Colormap(): class Colormap():
@ -498,13 +494,13 @@ class Colormap():
def interpolate_linear(lo, hi, frac): def interpolate_linear(lo, hi, frac):
"""Linear interpolation between lo and hi color at given fraction; output in model of lo color.""" """Linear interpolation between lo and hi color at given fraction; output in model of lo color."""
interpolation = (1.0 - frac) * np.array(lo.color[:]) \ interpolation = (1.0 - frac) * np.array(lo.color[:]) \
+ frac * np.array(hi.expressAs(lo.model).color[:]) + frac * np.array(hi.express_as(lo.model).color[:])
return Color(lo.model,interpolation) return Color(lo.model,interpolation)
if self.interpolate == 'perceptualuniform': if self.interpolate == 'perceptualuniform':
return interpolate_Msh(self.left.expressAs('MSH').color, return interpolate_Msh(self.left.express_as('MSH').color,
self.right.expressAs('MSH').color,fraction) self.right.express_as('MSH').color,fraction)
elif self.interpolate == 'linear': elif self.interpolate == 'linear':
return interpolate_linear(self.left, return interpolate_linear(self.left,
self.right,fraction) self.right,fraction)
@ -528,7 +524,7 @@ class Colormap():
""" """
format = format.lower() # consistent comparison basis format = format.lower() # consistent comparison basis
frac = 0.5*(np.array(crop) + 1.0) # rescale crop range to fractions frac = 0.5*(np.array(crop) + 1.0) # rescale crop range to fractions
colors = [self.color(float(i)/(steps-1)*(frac[1]-frac[0])+frac[0]).expressAs(model).color for i in range(steps)] colors = [self.color(float(i)/(steps-1)*(frac[1]-frac[0])+frac[0]).express_as(model).color for i in range(steps)]
if format == 'paraview': if format == 'paraview':
colormap = ['[\n {{\n "ColorSpace": "RGB", "Name": "{}", "DefaultMap": true,\n "RGBPoints" : ['.format(name)] \ colormap = ['[\n {{\n "ColorSpace": "RGB", "Name": "{}", "DefaultMap": true,\n "RGBPoints" : ['.format(name)] \
+ [' {:4d},{:8.6f},{:8.6f},{:8.6f},'.format(i,color[0],color[1],color[2],) \ + [' {:4d},{:8.6f},{:8.6f},{:8.6f},'.format(i,color[0],color[1],color[2],) \

6
python/damask/geom.py
View File

@ -80,7 +80,7 @@ class Geom():
if size is not None: if size is not None:
self.set_size(size) self.set_size(size)
elif rescale: elif rescale:
self.set_size(self.get_grid()/grid_old*self.size) self.set_size(self.get_grid()/grid_old*self.size)
message = ['grid a b c: {}'.format(' x '.join(map(str,grid_old)))] message = ['grid a b c: {}'.format(' x '.join(map(str,grid_old)))]
if np.any(grid_old != self.get_grid()): if np.any(grid_old != self.get_grid()):
@ -269,7 +269,7 @@ class Geom():
comments = [] comments = []
for i,line in enumerate(content[:header_length]): for i,line in enumerate(content[:header_length]):
items = line.lower().strip().split() items = line.lower().strip().split()
key = items[0] if len(items) > 0 else '' key = items[0] if items else ''
if key == 'grid': if key == 'grid':
grid = np.array([ int(dict(zip(items[1::2],items[2::2]))[i]) for i in ['a','b','c']]) grid = np.array([ int(dict(zip(items[1::2],items[2::2]))[i]) for i in ['a','b','c']])
elif key == 'size': elif key == 'size':
@ -524,7 +524,7 @@ class Geom():
"""Renumber sorted microstructure indices to 1,...,N.""" """Renumber sorted microstructure indices to 1,...,N."""
renumbered = np.empty(self.get_grid(),dtype=self.microstructure.dtype) renumbered = np.empty(self.get_grid(),dtype=self.microstructure.dtype)
for i, oldID in enumerate(np.unique(self.microstructure)): for i, oldID in enumerate(np.unique(self.microstructure)):
renumbered = np.where(self.microstructure == oldID, i+1, renumbered) renumbered = np.where(self.microstructure == oldID, i+1, renumbered)
return self.update(renumbered) return self.update(renumbered)
#self.add_comments('tbd') #self.add_comments('tbd')

24
python/damask/grid_filters.py
View File

@ -1,7 +1,7 @@
from scipy import spatial from scipy import spatial
import numpy as np import numpy as np
def __ks(size,grid,first_order=False): def _ks(size,grid,first_order=False):
""" """
Get wave numbers operator. Get wave numbers operator.
@ -34,7 +34,7 @@ def curl(size,field):
""" """
n = np.prod(field.shape[3:]) n = np.prod(field.shape[3:])
k_s = __ks(size,field.shape[:3],True) k_s = _ks(size,field.shape[:3],True)
e = np.zeros((3, 3, 3)) e = np.zeros((3, 3, 3))
e[0, 1, 2] = e[1, 2, 0] = e[2, 0, 1] = +1.0 # Levi-Civita symbol e[0, 1, 2] = e[1, 2, 0] = e[2, 0, 1] = +1.0 # Levi-Civita symbol
@ -58,7 +58,7 @@ def divergence(size,field):
""" """
n = np.prod(field.shape[3:]) n = np.prod(field.shape[3:])
k_s = __ks(size,field.shape[:3],True) k_s = _ks(size,field.shape[:3],True)
field_fourier = np.fft.rfftn(field,axes=(0,1,2)) field_fourier = np.fft.rfftn(field,axes=(0,1,2))
divergence = (np.einsum('ijkl,ijkl ->ijk', k_s,field_fourier)*2.0j*np.pi if n == 3 else # vector, 3 -> 1 divergence = (np.einsum('ijkl,ijkl ->ijk', k_s,field_fourier)*2.0j*np.pi if n == 3 else # vector, 3 -> 1
@ -78,7 +78,7 @@ def gradient(size,field):
""" """
n = np.prod(field.shape[3:]) n = np.prod(field.shape[3:])
k_s = __ks(size,field.shape[:3],True) k_s = _ks(size,field.shape[:3],True)
field_fourier = np.fft.rfftn(field,axes=(0,1,2)) field_fourier = np.fft.rfftn(field,axes=(0,1,2))
gradient = (np.einsum('ijkl,ijkm->ijkm', field_fourier,k_s)*2.0j*np.pi if n == 1 else # scalar, 1 -> 3 gradient = (np.einsum('ijkl,ijkm->ijkm', field_fourier,k_s)*2.0j*np.pi if n == 1 else # scalar, 1 -> 3
@ -110,6 +110,7 @@ def cell_coord0(grid,size,origin=np.zeros(3)):
return np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3) return np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3)
def cell_displacement_fluct(size,F): def cell_displacement_fluct(size,F):
""" """
Cell center displacement field from fluctuation part of the deformation gradient field. Cell center displacement field from fluctuation part of the deformation gradient field.
@ -124,7 +125,7 @@ def cell_displacement_fluct(size,F):
""" """
integrator = 0.5j*size/np.pi integrator = 0.5j*size/np.pi
k_s = __ks(size,F.shape[:3],False) k_s = _ks(size,F.shape[:3],False)
k_s_squared = np.einsum('...l,...l',k_s,k_s) k_s_squared = np.einsum('...l,...l',k_s,k_s)
k_s_squared[0,0,0] = 1.0 k_s_squared[0,0,0] = 1.0
@ -136,6 +137,7 @@ def cell_displacement_fluct(size,F):
return np.fft.irfftn(displacement,axes=(0,1,2),s=F.shape[:3]) return np.fft.irfftn(displacement,axes=(0,1,2),s=F.shape[:3])
def cell_displacement_avg(size,F): def cell_displacement_avg(size,F):
""" """
Cell center displacement field from average part of the deformation gradient field. Cell center displacement field from average part of the deformation gradient field.
@ -151,6 +153,7 @@ def cell_displacement_avg(size,F):
F_avg = np.average(F,axis=(0,1,2)) F_avg = np.average(F,axis=(0,1,2))
return np.einsum('ml,ijkl->ijkm',F_avg-np.eye(3),cell_coord0(F.shape[:3][::-1],size)) return np.einsum('ml,ijkl->ijkm',F_avg-np.eye(3),cell_coord0(F.shape[:3][::-1],size))
def cell_displacement(size,F): def cell_displacement(size,F):
""" """
Cell center displacement field from deformation gradient field. Cell center displacement field from deformation gradient field.
@ -165,6 +168,7 @@ def cell_displacement(size,F):
""" """
return cell_displacement_avg(size,F) + cell_displacement_fluct(size,F) return cell_displacement_avg(size,F) + cell_displacement_fluct(size,F)
def cell_coord(size,F,origin=np.zeros(3)): def cell_coord(size,F,origin=np.zeros(3)):
""" """
Cell center positions. Cell center positions.
@ -181,6 +185,7 @@ def cell_coord(size,F,origin=np.zeros(3)):
""" """
return cell_coord0(F.shape[:3][::-1],size,origin) + cell_displacement(size,F) return cell_coord0(F.shape[:3][::-1],size,origin) + cell_displacement(size,F)
def cell_coord0_gridSizeOrigin(coord0,ordered=True): def cell_coord0_gridSizeOrigin(coord0,ordered=True):
""" """
Return grid 'DNA', i.e. grid, size, and origin from array of cell positions. Return grid 'DNA', i.e. grid, size, and origin from array of cell positions.
@ -221,6 +226,7 @@ def cell_coord0_gridSizeOrigin(coord0,ordered=True):
return (grid,size,origin) return (grid,size,origin)
def coord0_check(coord0): def coord0_check(coord0):
""" """
Check whether coordinates lie on a regular grid. Check whether coordinates lie on a regular grid.
@ -234,7 +240,6 @@ def coord0_check(coord0):
cell_coord0_gridSizeOrigin(coord0,ordered=True) cell_coord0_gridSizeOrigin(coord0,ordered=True)
def node_coord0(grid,size,origin=np.zeros(3)): def node_coord0(grid,size,origin=np.zeros(3)):
""" """
Nodal positions (undeformed). Nodal positions (undeformed).
@ -256,6 +261,7 @@ def node_coord0(grid,size,origin=np.zeros(3)):
return np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3) return np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3)
def node_displacement_fluct(size,F): def node_displacement_fluct(size,F):
""" """
Nodal displacement field from fluctuation part of the deformation gradient field. Nodal displacement field from fluctuation part of the deformation gradient field.
@ -270,6 +276,7 @@ def node_displacement_fluct(size,F):
""" """
return cell_2_node(cell_displacement_fluct(size,F)) return cell_2_node(cell_displacement_fluct(size,F))
def node_displacement_avg(size,F): def node_displacement_avg(size,F):
""" """
Nodal displacement field from average part of the deformation gradient field. Nodal displacement field from average part of the deformation gradient field.
@ -285,6 +292,7 @@ def node_displacement_avg(size,F):
F_avg = np.average(F,axis=(0,1,2)) F_avg = np.average(F,axis=(0,1,2))
return np.einsum('ml,ijkl->ijkm',F_avg-np.eye(3),node_coord0(F.shape[:3][::-1],size)) return np.einsum('ml,ijkl->ijkm',F_avg-np.eye(3),node_coord0(F.shape[:3][::-1],size))
def node_displacement(size,F): def node_displacement(size,F):
""" """
Nodal displacement field from deformation gradient field. Nodal displacement field from deformation gradient field.
@ -299,6 +307,7 @@ def node_displacement(size,F):
""" """
return node_displacement_avg(size,F) + node_displacement_fluct(size,F) return node_displacement_avg(size,F) + node_displacement_fluct(size,F)
def node_coord(size,F,origin=np.zeros(3)): def node_coord(size,F,origin=np.zeros(3)):
""" """
Nodal positions. Nodal positions.
@ -315,6 +324,7 @@ def node_coord(size,F,origin=np.zeros(3)):
""" """
return node_coord0(F.shape[:3][::-1],size,origin) + node_displacement(size,F) return node_coord0(F.shape[:3][::-1],size,origin) + node_displacement(size,F)
def cell_2_node(cell_data): def cell_2_node(cell_data):
"""Interpolate periodic cell data to nodal data.""" """Interpolate periodic cell data to nodal data."""
n = ( cell_data + np.roll(cell_data,1,(0,1,2)) n = ( cell_data + np.roll(cell_data,1,(0,1,2))
@ -323,6 +333,7 @@ def cell_2_node(cell_data):
return np.pad(n,((0,1),(0,1),(0,1))+((0,0),)*len(cell_data.shape[3:]),mode='wrap') return np.pad(n,((0,1),(0,1),(0,1))+((0,0),)*len(cell_data.shape[3:]),mode='wrap')
def node_2_cell(node_data): def node_2_cell(node_data):
"""Interpolate periodic nodal data to cell data.""" """Interpolate periodic nodal data to cell data."""
c = ( node_data + np.roll(node_data,1,(0,1,2)) c = ( node_data + np.roll(node_data,1,(0,1,2))
@ -331,6 +342,7 @@ def node_2_cell(node_data):
return c[:-1,:-1,:-1] return c[:-1,:-1,:-1]
def node_coord0_gridSizeOrigin(coord0,ordered=False): def node_coord0_gridSizeOrigin(coord0,ordered=False):
""" """
Return grid 'DNA', i.e. grid, size, and origin from array of nodal positions. Return grid 'DNA', i.e. grid, size, and origin from array of nodal positions.

641
python/damask/lattice.py Normal file
View File

@ -0,0 +1,641 @@
import numpy as np
from .rotation import Rotation
P = -1
# ******************************************************************************************
class Symmetry:
"""
Symmetry operations for lattice systems.
References
----------
https://en.wikipedia.org/wiki/Crystal_system
"""
lattices = [None,'orthorhombic','tetragonal','hexagonal','cubic',]
def __init__(self, symmetry = None):
"""
Symmetry Definition.
Parameters
----------
symmetry : str, optional
label of the crystal system
"""
if symmetry is not None and symmetry.lower() not in Symmetry.lattices:
raise KeyError('Symmetry/crystal system "{}" is unknown'.format(symmetry))
self.lattice = symmetry.lower() if isinstance(symmetry,str) else symmetry
def __copy__(self):
"""Copy."""
return self.__class__(self.lattice)
copy = __copy__
def __repr__(self):
"""Readable string."""
return '{}'.format(self.lattice)
def __eq__(self, other):
"""
Equal to other.
Parameters
----------
other : Symmetry
Symmetry to check for equality.
"""
return self.lattice == other.lattice
def __neq__(self, other):
"""
Not Equal to other.
Parameters
----------
other : Symmetry
Symmetry to check for inequality.
"""
return not self.__eq__(other)
def __cmp__(self,other):
"""
Linear ordering.
Parameters
----------
other : Symmetry
Symmetry to check for for order.
"""
myOrder = Symmetry.lattices.index(self.lattice)
otherOrder = Symmetry.lattices.index(other.lattice)
return (myOrder > otherOrder) - (myOrder < otherOrder)
def symmetryOperations(self,members=[]):
"""List (or single element) of symmetry operations as rotations."""
if self.lattice == 'cubic':
symQuats = [
[ 1.0, 0.0, 0.0, 0.0 ],
[ 0.0, 1.0, 0.0, 0.0 ],
[ 0.0, 0.0, 1.0, 0.0 ],
[ 0.0, 0.0, 0.0, 1.0 ],
[ 0.0, 0.0, 0.5*np.sqrt(2), 0.5*np.sqrt(2) ],
[ 0.0, 0.0, 0.5*np.sqrt(2),-0.5*np.sqrt(2) ],
[ 0.0, 0.5*np.sqrt(2), 0.0, 0.5*np.sqrt(2) ],
[ 0.0, 0.5*np.sqrt(2), 0.0, -0.5*np.sqrt(2) ],
[ 0.0, 0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0 ],
[ 0.0, -0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0 ],
[ 0.5, 0.5, 0.5, 0.5 ],
[-0.5, 0.5, 0.5, 0.5 ],
[-0.5, 0.5, 0.5, -0.5 ],
[-0.5, 0.5, -0.5, 0.5 ],
[-0.5, -0.5, 0.5, 0.5 ],
[-0.5, -0.5, 0.5, -0.5 ],
[-0.5, -0.5, -0.5, 0.5 ],
[-0.5, 0.5, -0.5, -0.5 ],
[-0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
[ 0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
[-0.5*np.sqrt(2), 0.0, 0.5*np.sqrt(2), 0.0 ],
[-0.5*np.sqrt(2), 0.0, -0.5*np.sqrt(2), 0.0 ],
[-0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0, 0.0 ],
[-0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0, 0.0 ],
]
elif self.lattice == 'hexagonal':
symQuats = [
[ 1.0, 0.0, 0.0, 0.0 ],
[-0.5*np.sqrt(3), 0.0, 0.0, -0.5 ],
[ 0.5, 0.0, 0.0, 0.5*np.sqrt(3) ],
[ 0.0, 0.0, 0.0, 1.0 ],
[-0.5, 0.0, 0.0, 0.5*np.sqrt(3) ],
[-0.5*np.sqrt(3), 0.0, 0.0, 0.5 ],
[ 0.0, 1.0, 0.0, 0.0 ],
[ 0.0, -0.5*np.sqrt(3), 0.5, 0.0 ],
[ 0.0, 0.5, -0.5*np.sqrt(3), 0.0 ],
[ 0.0, 0.0, 1.0, 0.0 ],
[ 0.0, -0.5, -0.5*np.sqrt(3), 0.0 ],
[ 0.0, 0.5*np.sqrt(3), 0.5, 0.0 ],
]
elif self.lattice == 'tetragonal':
symQuats = [
[ 1.0, 0.0, 0.0, 0.0 ],
[ 0.0, 1.0, 0.0, 0.0 ],
[ 0.0, 0.0, 1.0, 0.0 ],
[ 0.0, 0.0, 0.0, 1.0 ],
[ 0.0, 0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0 ],
[ 0.0, -0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0 ],
[ 0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
[-0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
]
elif self.lattice == 'orthorhombic':
symQuats = [
[ 1.0,0.0,0.0,0.0 ],
[ 0.0,1.0,0.0,0.0 ],
[ 0.0,0.0,1.0,0.0 ],
[ 0.0,0.0,0.0,1.0 ],
]
else:
symQuats = [
[ 1.0,0.0,0.0,0.0 ],
]
symOps = list(map(Rotation,
np.array(symQuats)[np.atleast_1d(members) if members != [] else range(len(symQuats))]))
try:
iter(members) # asking for (even empty) list of members?
except TypeError:
return symOps[0] # no, return rotation object
else:
return symOps # yes, return list of rotations
def inFZ(self,rodrigues):
"""
Check whether given Rodriques-Frank vector falls into fundamental zone of own symmetry.
Fundamental zone in Rodrigues space is point symmetric around origin.
"""
if (len(rodrigues) != 3):
raise ValueError('Input is not a Rodriques-Frank vector.\n')
if np.any(rodrigues == np.inf): return False
Rabs = abs(rodrigues)
if self.lattice == 'cubic':
return np.sqrt(2.0)-1.0 >= Rabs[0] \
and np.sqrt(2.0)-1.0 >= Rabs[1] \
and np.sqrt(2.0)-1.0 >= Rabs[2] \
and 1.0 >= Rabs[0] + Rabs[1] + Rabs[2]
elif self.lattice == 'hexagonal':
return 1.0 >= Rabs[0] and 1.0 >= Rabs[1] and 1.0 >= Rabs[2] \
and 2.0 >= np.sqrt(3)*Rabs[0] + Rabs[1] \
and 2.0 >= np.sqrt(3)*Rabs[1] + Rabs[0] \
and 2.0 >= np.sqrt(3) + Rabs[2]
elif self.lattice == 'tetragonal':
return 1.0 >= Rabs[0] and 1.0 >= Rabs[1] \
and np.sqrt(2.0) >= Rabs[0] + Rabs[1] \
and np.sqrt(2.0) >= Rabs[2] + 1.0
elif self.lattice == 'orthorhombic':
return 1.0 >= Rabs[0] and 1.0 >= Rabs[1] and 1.0 >= Rabs[2]
else:
return True
def inDisorientationSST(self,rodrigues):
"""
Check whether given Rodriques-Frank vector (of misorientation) falls into standard stereographic triangle of own symmetry.
References
----------
A. Heinz and P. Neumann, Acta Crystallographica Section A 47:780-789, 1991
https://doi.org/10.1107/S0108767391006864
"""
if (len(rodrigues) != 3):
raise ValueError('Input is not a Rodriques-Frank vector.\n')
R = rodrigues
epsilon = 0.0
if self.lattice == 'cubic':
return R[0] >= R[1]+epsilon and R[1] >= R[2]+epsilon and R[2] >= epsilon
elif self.lattice == 'hexagonal':
return R[0] >= np.sqrt(3)*(R[1]-epsilon) and R[1] >= epsilon and R[2] >= epsilon
elif self.lattice == 'tetragonal':
return R[0] >= R[1]-epsilon and R[1] >= epsilon and R[2] >= epsilon
elif self.lattice == 'orthorhombic':
return R[0] >= epsilon and R[1] >= epsilon and R[2] >= epsilon
else:
return True
def inSST(self,
vector,
proper = False,
color = False):
"""
Check whether given vector falls into standard stereographic triangle of own symmetry.
proper considers only vectors with z >= 0, hence uses two neighboring SSTs.
Return inverse pole figure color if requested.
Bases are computed from
basis = {'cubic' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red
[1.,0.,1.]/np.sqrt(2.), # direction of green
[1.,1.,1.]/np.sqrt(3.)]).T), # direction of blue
'hexagonal' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red
[1.,0.,0.], # direction of green
[np.sqrt(3.),1.,0.]/np.sqrt(4.)]).T), # direction of blue
'tetragonal' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red
[1.,0.,0.], # direction of green
[1.,1.,0.]/np.sqrt(2.)]).T), # direction of blue
'orthorhombic' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red
[1.,0.,0.], # direction of green
[0.,1.,0.]]).T), # direction of blue
}
"""
if self.lattice == 'cubic':
basis = {'improper':np.array([ [-1. , 0. , 1. ],
[ np.sqrt(2.) , -np.sqrt(2.) , 0. ],
[ 0. , np.sqrt(3.) , 0. ] ]),
'proper':np.array([ [ 0. , -1. , 1. ],
[-np.sqrt(2.) , np.sqrt(2.) , 0. ],
[ np.sqrt(3.) , 0. , 0. ] ]),
}
elif self.lattice == 'hexagonal':
basis = {'improper':np.array([ [ 0. , 0. , 1. ],
[ 1. , -np.sqrt(3.) , 0. ],
[ 0. , 2. , 0. ] ]),
'proper':np.array([ [ 0. , 0. , 1. ],
[-1. , np.sqrt(3.) , 0. ],
[ np.sqrt(3.) , -1. , 0. ] ]),
}
elif self.lattice == 'tetragonal':
basis = {'improper':np.array([ [ 0. , 0. , 1. ],
[ 1. , -1. , 0. ],
[ 0. , np.sqrt(2.) , 0. ] ]),
'proper':np.array([ [ 0. , 0. , 1. ],
[-1. , 1. , 0. ],
[ np.sqrt(2.) , 0. , 0. ] ]),
}
elif self.lattice == 'orthorhombic':
basis = {'improper':np.array([ [ 0., 0., 1.],
[ 1., 0., 0.],
[ 0., 1., 0.] ]),
'proper':np.array([ [ 0., 0., 1.],
[-1., 0., 0.],
[ 0., 1., 0.] ]),
}
else: # direct exit for unspecified symmetry
if color:
return (True,np.zeros(3,'d'))
else:
return True
v = np.array(vector,dtype=float)
if proper: # check both improper ...
theComponents = np.around(np.dot(basis['improper'],v),12)
inSST = np.all(theComponents >= 0.0)
if not inSST: # ... and proper SST
theComponents = np.around(np.dot(basis['proper'],v),12)
inSST = np.all(theComponents >= 0.0)
else:
v[2] = abs(v[2]) # z component projects identical
theComponents = np.around(np.dot(basis['improper'],v),12) # for positive and negative values
inSST = np.all(theComponents >= 0.0)
if color: # have to return color array
if inSST:
rgb = np.power(theComponents/np.linalg.norm(theComponents),0.5) # smoothen color ramps
rgb = np.minimum(np.ones(3,dtype=float),rgb) # limit to maximum intensity
rgb /= max(rgb) # normalize to (HS)V = 1
else:
rgb = np.zeros(3,dtype=float)
return (inSST,rgb)
else:
return inSST
# code derived from https://github.com/ezag/pyeuclid
# suggested reading: http://web.mit.edu/2.998/www/QuaternionReport1.pdf
# ******************************************************************************************
class Lattice:
"""
Lattice system.
Currently, this contains only a mapping from Bravais lattice to symmetry
and orientation relationships. It could include twin and slip systems.
References
----------
https://en.wikipedia.org/wiki/Bravais_lattice
"""
lattices = {
'triclinic':{'symmetry':None},
'bct':{'symmetry':'tetragonal'},
'hex':{'symmetry':'hexagonal'},
'fcc':{'symmetry':'cubic','c/a':1.0},
'bcc':{'symmetry':'cubic','c/a':1.0},
}
def __init__(self, lattice):
"""
New lattice of given type.
Parameters
----------
lattice : str
Bravais lattice.
"""
self.lattice = lattice
self.symmetry = Symmetry(self.lattices[lattice]['symmetry'])
def __repr__(self):
"""Report basic lattice information."""
return 'Bravais lattice {} ({} symmetry)'.format(self.lattice,self.symmetry)
# Kurdjomov--Sachs orientation relationship for fcc <-> bcc transformation
# from S. Morito et al., Journal of Alloys and Compounds 577:s587-s592, 2013
# also see K. Kitahara et al., Acta Materialia 54:1279-1288, 2006
KS = {'mapping':{'fcc':0,'bcc':1},
'planes': np.array([
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, -1],[ 0, 1, 1]],
[[ 1, 1, -1],[ 0, 1, 1]],
[[ 1, 1, -1],[ 0, 1, 1]],
[[ 1, 1, -1],[ 0, 1, 1]],
[[ 1, 1, -1],[ 0, 1, 1]],
[[ 1, 1, -1],[ 0, 1, 1]]],dtype='float'),
'directions': np.array([
[[ -1, 0, 1],[ -1, -1, 1]],
[[ -1, 0, 1],[ -1, 1, -1]],
[[ 0, 1, -1],[ -1, -1, 1]],
[[ 0, 1, -1],[ -1, 1, -1]],
[[ 1, -1, 0],[ -1, -1, 1]],
[[ 1, -1, 0],[ -1, 1, -1]],
[[ 1, 0, -1],[ -1, -1, 1]],
[[ 1, 0, -1],[ -1, 1, -1]],
[[ -1, -1, 0],[ -1, -1, 1]],
[[ -1, -1, 0],[ -1, 1, -1]],
[[ 0, 1, 1],[ -1, -1, 1]],
[[ 0, 1, 1],[ -1, 1, -1]],
[[ 0, -1, 1],[ -1, -1, 1]],
[[ 0, -1, 1],[ -1, 1, -1]],
[[ -1, 0, -1],[ -1, -1, 1]],
[[ -1, 0, -1],[ -1, 1, -1]],
[[ 1, 1, 0],[ -1, -1, 1]],
[[ 1, 1, 0],[ -1, 1, -1]],
[[ -1, 1, 0],[ -1, -1, 1]],
[[ -1, 1, 0],[ -1, 1, -1]],
[[ 0, -1, -1],[ -1, -1, 1]],
[[ 0, -1, -1],[ -1, 1, -1]],
[[ 1, 0, 1],[ -1, -1, 1]],
[[ 1, 0, 1],[ -1, 1, -1]]],dtype='float')}
# Greninger--Troiano orientation relationship for fcc <-> bcc transformation
# from Y. He et al., Journal of Applied Crystallography 39:72-81, 2006
GT = {'mapping':{'fcc':0,'bcc':1},
'planes': np.array([
[[ 1, 1, 1],[ 1, 0, 1]],
[[ 1, 1, 1],[ 1, 1, 0]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ -1, -1, 1],[ -1, 0, 1]],
[[ -1, -1, 1],[ -1, -1, 0]],
[[ -1, -1, 1],[ 0, -1, 1]],
[[ -1, 1, 1],[ -1, 0, 1]],
[[ -1, 1, 1],[ -1, 1, 0]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 1, 0, 1]],
[[ 1, -1, 1],[ 1, -1, 0]],
[[ 1, -1, 1],[ 0, -1, 1]],
[[ 1, 1, 1],[ 1, 1, 0]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, 1],[ 1, 0, 1]],
[[ -1, -1, 1],[ -1, -1, 0]],
[[ -1, -1, 1],[ 0, -1, 1]],
[[ -1, -1, 1],[ -1, 0, 1]],
[[ -1, 1, 1],[ -1, 1, 0]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ -1, 0, 1]],
[[ 1, -1, 1],[ 1, -1, 0]],
[[ 1, -1, 1],[ 0, -1, 1]],
[[ 1, -1, 1],[ 1, 0, 1]]],dtype='float'),
'directions': np.array([
[[ -5,-12, 17],[-17, -7, 17]],
[[ 17, -5,-12],[ 17,-17, -7]],
[[-12, 17, -5],[ -7, 17,-17]],
[[ 5, 12, 17],[ 17, 7, 17]],
[[-17, 5,-12],[-17, 17, -7]],
[[ 12,-17, -5],[ 7,-17,-17]],
[[ -5, 12,-17],[-17, 7,-17]],
[[ 17, 5, 12],[ 17, 17, 7]],
[[-12,-17, 5],[ -7,-17, 17]],
[[ 5,-12,-17],[ 17, -7,-17]],
[[-17, -5, 12],[-17,-17, 7]],
[[ 12, 17, 5],[ 7, 17, 17]],
[[ -5, 17,-12],[-17, 17, -7]],
[[-12, -5, 17],[ -7,-17, 17]],
[[ 17,-12, -5],[ 17, -7,-17]],
[[ 5,-17,-12],[ 17,-17, -7]],
[[ 12, 5, 17],[ 7, 17, 17]],
[[-17, 12, -5],[-17, 7,-17]],
[[ -5,-17, 12],[-17,-17, 7]],
[[-12, 5,-17],[ -7, 17,-17]],
[[ 17, 12, 5],[ 17, 7, 17]],
[[ 5, 17, 12],[ 17, 17, 7]],
[[ 12, -5,-17],[ 7,-17,-17]],
[[-17,-12, 5],[-17,-7, 17]]],dtype='float')}
# Greninger--Troiano' orientation relationship for fcc <-> bcc transformation
# from Y. He et al., Journal of Applied Crystallography 39:72-81, 2006
GTprime = {'mapping':{'fcc':0,'bcc':1},
'planes': np.array([
[[ 7, 17, 17],[ 12, 5, 17]],
[[ 17, 7, 17],[ 17, 12, 5]],
[[ 17, 17, 7],[ 5, 17, 12]],
[[ -7,-17, 17],[-12, -5, 17]],
[[-17, -7, 17],[-17,-12, 5]],
[[-17,-17, 7],[ -5,-17, 12]],
[[ 7,-17,-17],[ 12, -5,-17]],
[[ 17, -7,-17],[ 17,-12, -5]],
[[ 17,-17, -7],[ 5,-17,-12]],
[[ -7, 17,-17],[-12, 5,-17]],
[[-17, 7,-17],[-17, 12, -5]],
[[-17, 17, -7],[ -5, 17,-12]],
[[ 7, 17, 17],[ 12, 17, 5]],
[[ 17, 7, 17],[ 5, 12, 17]],
[[ 17, 17, 7],[ 17, 5, 12]],
[[ -7,-17, 17],[-12,-17, 5]],
[[-17, -7, 17],[ -5,-12, 17]],
[[-17,-17, 7],[-17, -5, 12]],
[[ 7,-17,-17],[ 12,-17, -5]],
[[ 17, -7,-17],[ 5, -12,-17]],
[[ 17,-17, -7],[ 17, -5,-12]],
[[ -7, 17,-17],[-12, 17, -5]],
[[-17, 7,-17],[ -5, 12,-17]],
[[-17, 17, -7],[-17, 5,-12]]],dtype='float'),
'directions': np.array([
[[ 0, 1, -1],[ 1, 1, -1]],
[[ -1, 0, 1],[ -1, 1, 1]],
[[ 1, -1, 0],[ 1, -1, 1]],
[[ 0, -1, -1],[ -1, -1, -1]],
[[ 1, 0, 1],[ 1, -1, 1]],
[[ 1, -1, 0],[ 1, -1, -1]],
[[ 0, 1, -1],[ -1, 1, -1]],
[[ 1, 0, 1],[ 1, 1, 1]],
[[ -1, -1, 0],[ -1, -1, 1]],
[[ 0, -1, -1],[ 1, -1, -1]],
[[ -1, 0, 1],[ -1, -1, 1]],
[[ -1, -1, 0],[ -1, -1, -1]],
[[ 0, -1, 1],[ 1, -1, 1]],
[[ 1, 0, -1],[ 1, 1, -1]],
[[ -1, 1, 0],[ -1, 1, 1]],
[[ 0, 1, 1],[ -1, 1, 1]],
[[ -1, 0, -1],[ -1, -1, -1]],
[[ -1, 1, 0],[ -1, 1, -1]],
[[ 0, -1, 1],[ -1, -1, 1]],
[[ -1, 0, -1],[ -1, 1, -1]],
[[ 1, 1, 0],[ 1, 1, 1]],
[[ 0, 1, 1],[ 1, 1, 1]],
[[ 1, 0, -1],[ 1, -1, -1]],
[[ 1, 1, 0],[ 1, 1, -1]]],dtype='float')}
# Nishiyama--Wassermann orientation relationship for fcc <-> bcc transformation
# from H. Kitahara et al., Materials Characterization 54:378-386, 2005
NW = {'mapping':{'fcc':0,'bcc':1},
'planes': np.array([
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ -1, -1, 1],[ 0, 1, 1]],
[[ -1, -1, 1],[ 0, 1, 1]],
[[ -1, -1, 1],[ 0, 1, 1]]],dtype='float'),
'directions': np.array([
[[ 2, -1, -1],[ 0, -1, 1]],
[[ -1, 2, -1],[ 0, -1, 1]],
[[ -1, -1, 2],[ 0, -1, 1]],
[[ -2, -1, -1],[ 0, -1, 1]],
[[ 1, 2, -1],[ 0, -1, 1]],
[[ 1, -1, 2],[ 0, -1, 1]],
[[ 2, 1, -1],[ 0, -1, 1]],
[[ -1, -2, -1],[ 0, -1, 1]],
[[ -1, 1, 2],[ 0, -1, 1]],
[[ 2, -1, 1],[ 0, -1, 1]], #It is wrong in the paper, but matrix is correct
[[ -1, 2, 1],[ 0, -1, 1]],
[[ -1, -1, -2],[ 0, -1, 1]]],dtype='float')}
# Pitsch orientation relationship for fcc <-> bcc transformation
# from Y. He et al., Acta Materialia 53:1179-1190, 2005
Pitsch = {'mapping':{'fcc':0,'bcc':1},
'planes': np.array([
[[ 0, 1, 0],[ -1, 0, 1]],
[[ 0, 0, 1],[ 1, -1, 0]],
[[ 1, 0, 0],[ 0, 1, -1]],
[[ 1, 0, 0],[ 0, -1, -1]],
[[ 0, 1, 0],[ -1, 0, -1]],
[[ 0, 0, 1],[ -1, -1, 0]],
[[ 0, 1, 0],[ -1, 0, -1]],
[[ 0, 0, 1],[ -1, -1, 0]],
[[ 1, 0, 0],[ 0, -1, -1]],
[[ 1, 0, 0],[ 0, -1, 1]],
[[ 0, 1, 0],[ 1, 0, -1]],
[[ 0, 0, 1],[ -1, 1, 0]]],dtype='float'),
'directions': np.array([
[[ 1, 0, 1],[ 1, -1, 1]],
[[ 1, 1, 0],[ 1, 1, -1]],
[[ 0, 1, 1],[ -1, 1, 1]],
[[ 0, 1, -1],[ -1, 1, -1]],
[[ -1, 0, 1],[ -1, -1, 1]],
[[ 1, -1, 0],[ 1, -1, -1]],
[[ 1, 0, -1],[ 1, -1, -1]],
[[ -1, 1, 0],[ -1, 1, -1]],
[[ 0, -1, 1],[ -1, -1, 1]],
[[ 0, 1, 1],[ -1, 1, 1]],
[[ 1, 0, 1],[ 1, -1, 1]],
[[ 1, 1, 0],[ 1, 1, -1]]],dtype='float')}
# Bain orientation relationship for fcc <-> bcc transformation
# from Y. He et al., Journal of Applied Crystallography 39:72-81, 2006
Bain = {'mapping':{'fcc':0,'bcc':1},
'planes': np.array([
[[ 1, 0, 0],[ 1, 0, 0]],
[[ 0, 1, 0],[ 0, 1, 0]],
[[ 0, 0, 1],[ 0, 0, 1]]],dtype='float'),
'directions': np.array([
[[ 0, 1, 0],[ 0, 1, 1]],
[[ 0, 0, 1],[ 1, 0, 1]],
[[ 1, 0, 0],[ 1, 1, 0]]],dtype='float')}
def relationOperations(self,model):
"""
Crystallographic orientation relationships for phase transformations.
References
----------
S. Morito et al., Journal of Alloys and Compounds 577:s587-s592, 2013
https://doi.org/10.1016/j.jallcom.2012.02.004
K. Kitahara et al., Acta Materialia 54(5):1279-1288, 2006
https://doi.org/10.1016/j.actamat.2005.11.001
Y. He et al., Journal of Applied Crystallography 39:72-81, 2006
https://doi.org/10.1107/S0021889805038276
H. Kitahara et al., Materials Characterization 54(4-5):378-386, 2005
https://doi.org/10.1016/j.matchar.2004.12.015
Y. He et al., Acta Materialia 53(4):1179-1190, 2005
https://doi.org/10.1016/j.actamat.2004.11.021
"""
models={'KS':self.KS, 'GT':self.GT, 'GT_prime':self.GTprime,
'NW':self.NW, 'Pitsch': self.Pitsch, 'Bain':self.Bain}
try:
relationship = models[model]
except KeyError :
raise KeyError('Orientation relationship "{}" is unknown'.format(model))
if self.lattice not in relationship['mapping']:
raise ValueError('Relationship "{}" not supported for lattice "{}"'.format(model,self.lattice))
r = {'lattice':Lattice((set(relationship['mapping'])-{self.lattice}).pop()), # target lattice
'rotations':[] }
myPlane_id = relationship['mapping'][self.lattice]
otherPlane_id = (myPlane_id+1)%2
myDir_id = myPlane_id +2
otherDir_id = otherPlane_id +2
for miller in np.hstack((relationship['planes'],relationship['directions'])):
myPlane = miller[myPlane_id]/ np.linalg.norm(miller[myPlane_id])
myDir = miller[myDir_id]/ np.linalg.norm(miller[myDir_id])
myMatrix = np.array([myDir,np.cross(myPlane,myDir),myPlane])
otherPlane = miller[otherPlane_id]/ np.linalg.norm(miller[otherPlane_id])
otherDir = miller[otherDir_id]/ np.linalg.norm(miller[otherDir_id])
otherMatrix = np.array([otherDir,np.cross(otherPlane,otherDir),otherPlane])
r['rotations'].append(Rotation.fromMatrix(np.dot(otherMatrix.T,myMatrix)))
return r

1665
python/damask/orientation.py
View File

File diff suppressed because it is too large Load Diff

837
python/damask/rotation.py Normal file
View File

@ -0,0 +1,837 @@
import numpy as np
from . import Lambert
P = -1
def iszero(a):
return np.isclose(a,0.0,atol=1.0e-12,rtol=0.0)
class Rotation:
u"""
Orientation stored with functionality for conversion to different representations.
References
----------
D. Rowenhorst et al., Modelling and Simulation in Materials Science and Engineering 23:083501, 2015
https://doi.org/10.1088/0965-0393/23/8/083501
Conventions
-----------
Convention 1: Coordinate frames are right-handed.
Convention 2: A rotation angle ω is taken to be positive for a counterclockwise rotation
when viewing from the end point of the rotation axis towards the origin.
Convention 3: Rotations will be interpreted in the passive sense.
Convention 4: Euler angle triplets are implemented using the Bunge convention,
with the angular ranges as [0, 2π],[0, π],[0, 2π].
Convention 5: The rotation angle ω is limited to the interval [0, π].
Convention 6: the real part of a quaternion is positive, Re(q) > 0
Convention 7: P = -1 (as default).
Usage
-----
Vector "a" (defined in coordinate system "A") is passively rotated
resulting in new coordinates "b" when expressed in system "B".
b = Q * a
b = np.dot(Q.asMatrix(),a)
"""
__slots__ = ['quaternion']
def __init__(self,quaternion = np.array([1.0,0.0,0.0,0.0])):
"""
Initializes to identity unless specified.
Parameters
----------
quaternion : numpy.ndarray, optional
Unit quaternion that follows the conventions. Use .fromQuaternion to perform a sanity check.
"""
self.quaternion = quaternion.copy()
def __copy__(self):
"""Copy."""
return self.__class__(self.quaternion)
copy = __copy__
def __repr__(self):
"""Orientation displayed as unit quaternion, rotation matrix, and Bunge-Euler angles."""
return '\n'.join([
'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)),
'Matrix:\n{}'.format(self.asMatrix()),
'Bunge Eulers / deg: ({:3.2f}, {:3.2f}, {:3.2f})'.format(*self.asEulers(degrees=True)),
])
def __mul__(self, other):
"""
Multiplication.
Parameters
----------
other : numpy.ndarray or Rotation
Vector, second or fourth order tensor, or rotation object that is rotated.
Todo
----
Document details active/passive)
considere rotation of (3,3,3,3)-matrix
"""
if isinstance(other, Rotation): # rotate a rotation
self_q = self.quaternion[0]
self_p = self.quaternion[1:]
other_q = other.quaternion[0]
other_p = other.quaternion[1:]
R = self.__class__(np.append(self_q*other_q - np.dot(self_p,other_p),
self_q*other_p + other_q*self_p + P * np.cross(self_p,other_p)))
return R.standardize()
elif isinstance(other, (tuple,np.ndarray)):
if isinstance(other,tuple) or other.shape == (3,): # rotate a single (3)-vector or meshgrid
A = self.quaternion[0]**2.0 - np.dot(self.quaternion[1:],self.quaternion[1:])
B = 2.0 * ( self.quaternion[1]*other[0]
+ self.quaternion[2]*other[1]
+ self.quaternion[3]*other[2])
C = 2.0 * P*self.quaternion[0]
return np.array([
A*other[0] + B*self.quaternion[1] + C*(self.quaternion[2]*other[2] - self.quaternion[3]*other[1]),
A*other[1] + B*self.quaternion[2] + C*(self.quaternion[3]*other[0] - self.quaternion[1]*other[2]),
A*other[2] + B*self.quaternion[3] + C*(self.quaternion[1]*other[1] - self.quaternion[2]*other[0]),
])
elif other.shape == (3,3,): # rotate a single (3x3)-matrix
return np.dot(self.asMatrix(),np.dot(other,self.asMatrix().T))
elif other.shape == (3,3,3,3,):
raise NotImplementedError
else:
return NotImplemented
else:
return NotImplemented
def inverse(self):
"""In-place inverse rotation/backward rotation."""
self.quaternion[1:] *= -1
return self
def inversed(self):
"""Inverse rotation/backward rotation."""
return self.copy().inverse()
def standardize(self):
"""In-place quaternion representation with positive q."""
if self.quaternion[0] < 0.0: self.quaternion*=-1
return self
def standardized(self):
"""Quaternion representation with positive q."""
return self.copy().standardize()
def misorientation(self,other):
"""
Get Misorientation.
Parameters
----------
other : Rotation
Rotation to which the misorientation is computed.
"""
return other*self.inversed()
def average(self,other):
"""
Calculate the average rotation.
Parameters
----------
other : Rotation
Rotation from which the average is rotated.
"""
return Rotation.fromAverage([self,other])
################################################################################################
# convert to different orientation representations (numpy arrays)
def asQuaternion(self):
"""
Unit quaternion [q, p_1, p_2, p_3] unless quaternion == True: damask.quaternion object.
Parameters
----------
quaternion : bool, optional
return quaternion as DAMASK object.
"""
return self.quaternion
def asEulers(self,
degrees = False):
"""
Bunge-Euler angles: (φ_1, ϕ, φ_2).
Parameters
----------
degrees : bool, optional
return angles in degrees.
"""
eu = Rotation.qu2eu(self.quaternion)
if degrees: eu = np.degrees(eu)
return eu
def asAxisAngle(self,
degrees = False,
pair = False):
"""
Axis angle representation [n_1, n_2, n_3, ω] unless pair == True: ([n_1, n_2, n_3], ω).
Parameters
----------
degrees : bool, optional
return rotation angle in degrees.
pair : bool, optional
return tuple of axis and angle.
"""
ax = Rotation.qu2ax(self.quaternion)
if degrees: ax[3] = np.degrees(ax[3])
return (ax[:3],np.degrees(ax[3])) if pair else ax
def asMatrix(self):
"""Rotation matrix."""
return Rotation.qu2om(self.quaternion)
def asRodrigues(self,
vector = False):
"""
Rodrigues-Frank vector representation [n_1, n_2, n_3, tan(ω/2)] unless vector == True: [n_1, n_2, n_3] * tan(ω/2).
Parameters
----------
vector : bool, optional
return as actual Rodrigues--Frank vector, i.e. rotation axis scaled by tan(ω/2).
"""
ro = Rotation.qu2ro(self.quaternion)
return ro[:3]*ro[3] if vector else ro
def asHomochoric(self):
"""Homochoric vector: (h_1, h_2, h_3)."""
return Rotation.qu2ho(self.quaternion)
def asCubochoric(self):
"""Cubochoric vector: (c_1, c_2, c_3)."""
return Rotation.qu2cu(self.quaternion)
def asM(self):
"""
Intermediate representation supporting quaternion averaging.
References
----------
F. Landis Markley et al., Journal of Guidance, Control, and Dynamics 30(4):1193-1197, 2007
https://doi.org/10.2514/1.28949
"""
return np.outer(self.quaternion,self.quaternion)
################################################################################################
# static constructors. The input data needs to follow the convention, options allow to
# relax these convections
@staticmethod
def fromQuaternion(quaternion,
acceptHomomorph = False,
P = -1):
qu = quaternion if isinstance(quaternion,np.ndarray) and quaternion.dtype == np.dtype(float) \
else np.array(quaternion,dtype=float)
if P > 0: qu[1:4] *= -1 # convert from P=1 to P=-1
if qu[0] < 0.0:
if acceptHomomorph:
qu *= -1.
else:
raise ValueError('Quaternion has negative first component.\n{}'.format(qu[0]))
if not np.isclose(np.linalg.norm(qu), 1.0):
raise ValueError('Quaternion is not of unit length.\n{} {} {} {}'.format(*qu))
return Rotation(qu)
@staticmethod
def fromEulers(eulers,
degrees = False):
eu = eulers if isinstance(eulers, np.ndarray) and eulers.dtype == np.dtype(float) \
else np.array(eulers,dtype=float)
eu = np.radians(eu) if degrees else eu
if np.any(eu < 0.0) or np.any(eu > 2.0*np.pi) or eu[1] > np.pi:
raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π].\n{} {} {}.'.format(*eu))
return Rotation(Rotation.eu2qu(eu))
@staticmethod
def fromAxisAngle(angleAxis,
degrees = False,
normalise = False,
P = -1):
ax = angleAxis if isinstance(angleAxis, np.ndarray) and angleAxis.dtype == np.dtype(float) \
else np.array(angleAxis,dtype=float)
if P > 0: ax[0:3] *= -1 # convert from P=1 to P=-1
if degrees: ax[ 3] = np.radians(ax[3])
if normalise: ax[0:3] /= np.linalg.norm(ax[0:3])
if ax[3] < 0.0 or ax[3] > np.pi:
raise ValueError('Axis angle rotation angle outside of [0..π].\n'.format(ax[3]))
if not np.isclose(np.linalg.norm(ax[0:3]), 1.0):
raise ValueError('Axis angle rotation axis is not of unit length.\n{} {} {}'.format(*ax[0:3]))
return Rotation(Rotation.ax2qu(ax))
@staticmethod
def fromBasis(basis,
orthonormal = True,
reciprocal = False,
):
om = basis if isinstance(basis, np.ndarray) else np.array(basis).reshape((3,3))
if reciprocal:
om = np.linalg.inv(om.T/np.pi) # transform reciprocal basis set
orthonormal = False # contains stretch
if not orthonormal:
(U,S,Vh) = np.linalg.svd(om) # singular value decomposition
om = np.dot(U,Vh)
if not np.isclose(np.linalg.det(om),1.0):
raise ValueError('matrix is not a proper rotation.\n{}'.format(om))
if not np.isclose(np.dot(om[0],om[1]), 0.0) \
or not np.isclose(np.dot(om[1],om[2]), 0.0) \
or not np.isclose(np.dot(om[2],om[0]), 0.0):
raise ValueError('matrix is not orthogonal.\n{}'.format(om))
return Rotation(Rotation.om2qu(om))
@staticmethod
def fromMatrix(om,
):
return Rotation.fromBasis(om)
@staticmethod
def fromRodrigues(rodrigues,
normalise = False,
P = -1):
ro = rodrigues if isinstance(rodrigues, np.ndarray) and rodrigues.dtype == np.dtype(float) \
else np.array(rodrigues,dtype=float)
if P > 0: ro[0:3] *= -1 # convert from P=1 to P=-1
if normalise: ro[0:3] /= np.linalg.norm(ro[0:3])
if not np.isclose(np.linalg.norm(ro[0:3]), 1.0):
raise ValueError('Rodrigues rotation axis is not of unit length.\n{} {} {}'.format(*ro[0:3]))
if ro[3] < 0.0:
raise ValueError('Rodriques rotation angle not positive.\n'.format(ro[3]))
return Rotation(Rotation.ro2qu(ro))
@staticmethod
def fromHomochoric(homochoric,
P = -1):
ho = homochoric if isinstance(homochoric, np.ndarray) and homochoric.dtype == np.dtype(float) \
else np.array(homochoric,dtype=float)
if P > 0: ho *= -1 # convert from P=1 to P=-1
return Rotation(Rotation.ho2qu(ho))
@staticmethod
def fromCubochoric(cubochoric,
P = -1):
cu = cubochoric if isinstance(cubochoric, np.ndarray) and cubochoric.dtype == np.dtype(float) \
else np.array(cubochoric,dtype=float)
ho = Rotation.cu2ho(cu)
if P > 0: ho *= -1 # convert from P=1 to P=-1
return Rotation(Rotation.ho2qu(ho))
@staticmethod
def fromAverage(rotations,
weights = []):
"""
Average rotation.
References
----------
F. Landis Markley et al., Journal of Guidance, Control, and Dynamics 30(4):1193-1197, 2007
https://doi.org/10.2514/1.28949
Parameters
----------
rotations : list of Rotations
Rotations to average from
weights : list of floats, optional
Weights for each rotation used for averaging
"""
if not all(isinstance(item, Rotation) for item in rotations):
raise TypeError("Only instances of Rotation can be averaged.")
N = len(rotations)
if weights == [] or not weights:
weights = np.ones(N,dtype='i')
for i,(r,n) in enumerate(zip(rotations,weights)):
M = r.asM() * n if i == 0 \
else M + r.asM() * n # noqa add (multiples) of this rotation to average noqa
eig, vec = np.linalg.eig(M/N)
return Rotation.fromQuaternion(np.real(vec.T[eig.argmax()]),acceptHomomorph = True)
@staticmethod
def fromRandom():
r = np.random.random(3)
A = np.sqrt(r[2])
B = np.sqrt(1.0-r[2])
return Rotation(np.array([np.cos(2.0*np.pi*r[0])*A,
np.sin(2.0*np.pi*r[1])*B,
np.cos(2.0*np.pi*r[1])*B,
np.sin(2.0*np.pi*r[0])*A])).standardize()
####################################################################################################
# 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.
####################################################################################################
#---------- Quaternion ----------
@staticmethod
def qu2om(qu):
"""Quaternion to rotation matrix."""
qq = qu[0]**2-(qu[1]**2 + qu[2]**2 + qu[3]**2)
om = np.diag(qq + 2.0*np.array([qu[1],qu[2],qu[3]])**2)
om[1,0] = 2.0*(qu[2]*qu[1]+qu[0]*qu[3])
om[0,1] = 2.0*(qu[1]*qu[2]-qu[0]*qu[3])
om[2,1] = 2.0*(qu[3]*qu[2]+qu[0]*qu[1])
om[1,2] = 2.0*(qu[2]*qu[3]-qu[0]*qu[1])
om[0,2] = 2.0*(qu[1]*qu[3]+qu[0]*qu[2])
om[2,0] = 2.0*(qu[3]*qu[1]-qu[0]*qu[2])
return om if P > 0.0 else om.T
@staticmethod
def qu2eu(qu):
"""Quaternion to Bunge-Euler angles."""
q03 = qu[0]**2+qu[3]**2
q12 = qu[1]**2+qu[2]**2
chi = np.sqrt(q03*q12)
if iszero(chi):
eu = np.array([np.arctan2(-P*2.0*qu[0]*qu[3],qu[0]**2-qu[3]**2), 0.0, 0.0]) if iszero(q12) else \
np.array([np.arctan2(2.0*qu[1]*qu[2],qu[1]**2-qu[2]**2), np.pi, 0.0])
else:
eu = np.array([np.arctan2((-P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-P*qu[0]*qu[1]-qu[2]*qu[3])*chi ),
np.arctan2( 2.0*chi, q03-q12 ),
np.arctan2(( P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-P*qu[0]*qu[1]+qu[2]*qu[3])*chi )])
# reduce Euler angles to definition range, i.e a lower limit of 0.0
eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu)
return eu
@staticmethod
def qu2ax(qu):
"""
Quaternion to axis angle pair.
Modified version of the original formulation, should be numerically more stable
"""
if iszero(qu[1]**2+qu[2]**2+qu[3]**2): # set axis to [001] if the angle is 0/360
ax = [ 0.0, 0.0, 1.0, 0.0 ]
elif not iszero(qu[0]):
s = np.sign(qu[0])/np.sqrt(qu[1]**2+qu[2]**2+qu[3]**2)
omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0))
ax = [ qu[1]*s, qu[2]*s, qu[3]*s, omega ]
else:
ax = [ qu[1], qu[2], qu[3], np.pi]
return np.array(ax)
@staticmethod
def qu2ro(qu):
"""Quaternion to Rodriques-Frank vector."""
if iszero(qu[0]):
ro = [qu[1], qu[2], qu[3], np.inf]
else:
s = np.linalg.norm([qu[1],qu[2],qu[3]])
ro = [0.0,0.0,P,0.0] if iszero(s) else \
[ qu[1]/s, qu[2]/s, qu[3]/s, np.tan(np.arccos(np.clip(qu[0],-1.0,1.0)))]
return np.array(ro)
@staticmethod
def qu2ho(qu):
"""Quaternion to homochoric vector."""
omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0))
if iszero(omega):
ho = np.array([ 0.0, 0.0, 0.0 ])
else:
ho = np.array([qu[1], qu[2], qu[3]])
f = 0.75 * ( omega - np.sin(omega) )
ho = ho/np.linalg.norm(ho) * f**(1./3.)
return ho
@staticmethod
def qu2cu(qu):
"""Quaternion to cubochoric vector."""
return Rotation.ho2cu(Rotation.qu2ho(qu))
#---------- Rotation matrix ----------
@staticmethod
def om2qu(om):
"""
Rotation matrix to quaternion.
The original formulation (direct conversion) had (numerical?) issues
"""
return Rotation.eu2qu(Rotation.om2eu(om))
@staticmethod
def om2eu(om):
"""Rotation matrix to Bunge-Euler angles."""
if abs(om[2,2]) < 1.0:
zeta = 1.0/np.sqrt(1.0-om[2,2]**2)
eu = np.array([np.arctan2(om[2,0]*zeta,-om[2,1]*zeta),
np.arccos(om[2,2]),
np.arctan2(om[0,2]*zeta, om[1,2]*zeta)])
else:
eu = np.array([np.arctan2( om[0,1],om[0,0]), np.pi*0.5*(1-om[2,2]),0.0]) # following the paper, not the reference implementation
# reduce Euler angles to definition range, i.e a lower limit of 0.0
eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu)
return eu
@staticmethod
def om2ax(om):
"""Rotation matrix to axis angle pair."""
ax=np.empty(4)
# first get the rotation angle
t = 0.5*(om.trace() -1.0)
ax[3] = np.arccos(np.clip(t,-1.0,1.0))
if iszero(ax[3]):
ax = [ 0.0, 0.0, 1.0, 0.0]
else:
w,vr = np.linalg.eig(om)
# next, find the eigenvalue (1,0j)
i = np.where(np.isclose(w,1.0+0.0j))[0][0]
ax[0:3] = np.real(vr[0:3,i])
diagDelta = np.array([om[1,2]-om[2,1],om[2,0]-om[0,2],om[0,1]-om[1,0]])
ax[0:3] = np.where(iszero(diagDelta), ax[0:3],np.abs(ax[0:3])*np.sign(-P*diagDelta))
return np.array(ax)
@staticmethod
def om2ro(om):
"""Rotation matrix to Rodriques-Frank vector."""
return Rotation.eu2ro(Rotation.om2eu(om))
@staticmethod
def om2ho(om):
"""Rotation matrix to homochoric vector."""
return Rotation.ax2ho(Rotation.om2ax(om))
@staticmethod
def om2cu(om):
"""Rotation matrix to cubochoric vector."""
return Rotation.ho2cu(Rotation.om2ho(om))
#---------- Bunge-Euler angles ----------
@staticmethod
def eu2qu(eu):
"""Bunge-Euler angles to quaternion."""
ee = 0.5*eu
cPhi = np.cos(ee[1])
sPhi = np.sin(ee[1])
qu = np.array([ cPhi*np.cos(ee[0]+ee[2]),
-P*sPhi*np.cos(ee[0]-ee[2]),
-P*sPhi*np.sin(ee[0]-ee[2]),
-P*cPhi*np.sin(ee[0]+ee[2]) ])
if qu[0] < 0.0: qu*=-1
return qu
@staticmethod
def eu2om(eu):
"""Bunge-Euler angles to rotation matrix."""
c = np.cos(eu)
s = np.sin(eu)
om = np.array([[+c[0]*c[2]-s[0]*s[2]*c[1], +s[0]*c[2]+c[0]*s[2]*c[1], +s[2]*s[1]],
[-c[0]*s[2]-s[0]*c[2]*c[1], -s[0]*s[2]+c[0]*c[2]*c[1], +c[2]*s[1]],
[+s[0]*s[1], -c[0]*s[1], +c[1] ]])
om[np.where(iszero(om))] = 0.0
return om
@staticmethod
def eu2ax(eu):
"""Bunge-Euler angles to axis angle pair."""
t = np.tan(eu[1]*0.5)
sigma = 0.5*(eu[0]+eu[2])
delta = 0.5*(eu[0]-eu[2])
tau = np.linalg.norm([t,np.sin(sigma)])
alpha = np.pi if iszero(np.cos(sigma)) else \
2.0*np.arctan(tau/np.cos(sigma))
if iszero(alpha):
ax = np.array([ 0.0, 0.0, 1.0, 0.0 ])
else:
ax = -P/tau * np.array([ t*np.cos(delta), t*np.sin(delta), np.sin(sigma) ]) # passive axis angle pair so a minus sign in front
ax = np.append(ax,alpha)
if alpha < 0.0: ax *= -1.0 # ensure alpha is positive
return ax
@staticmethod
def eu2ro(eu):
"""Bunge-Euler angles to Rodriques-Frank vector."""
ro = Rotation.eu2ax(eu) # convert to axis angle pair representation
if ro[3] >= np.pi: # Differs from original implementation. check convention 5
ro[3] = np.inf
elif iszero(ro[3]):
ro = np.array([ 0.0, 0.0, P, 0.0 ])
else:
ro[3] = np.tan(ro[3]*0.5)
return ro
@staticmethod
def eu2ho(eu):
"""Bunge-Euler angles to homochoric vector."""
return Rotation.ax2ho(Rotation.eu2ax(eu))
@staticmethod
def eu2cu(eu):
"""Bunge-Euler angles to cubochoric vector."""
return Rotation.ho2cu(Rotation.eu2ho(eu))
#---------- Axis angle pair ----------
@staticmethod
def ax2qu(ax):
"""Axis angle pair to quaternion."""
if iszero(ax[3]):
qu = np.array([ 1.0, 0.0, 0.0, 0.0 ])
else:
c = np.cos(ax[3]*0.5)
s = np.sin(ax[3]*0.5)
qu = np.array([ c, ax[0]*s, ax[1]*s, ax[2]*s ])
return qu
@staticmethod
def ax2om(ax):
"""Axis angle pair to rotation matrix."""
c = np.cos(ax[3])
s = np.sin(ax[3])
omc = 1.0-c
om=np.diag(ax[0:3]**2*omc + c)
for idx in [[0,1,2],[1,2,0],[2,0,1]]:
q = omc*ax[idx[0]] * ax[idx[1]]
om[idx[0],idx[1]] = q + s*ax[idx[2]]
om[idx[1],idx[0]] = q - s*ax[idx[2]]
return om if P < 0.0 else om.T
@staticmethod
def ax2eu(ax):
"""Rotation matrix to Bunge Euler angles."""
return Rotation.om2eu(Rotation.ax2om(ax))
@staticmethod
def ax2ro(ax):
"""Axis angle pair to Rodriques-Frank vector."""
if iszero(ax[3]):
ro = [ 0.0, 0.0, P, 0.0 ]
else:
ro = [ax[0], ax[1], ax[2]]
# 180 degree case
ro += [np.inf] if np.isclose(ax[3],np.pi,atol=1.0e-15,rtol=0.0) else \
[np.tan(ax[3]*0.5)]
return np.array(ro)
@staticmethod
def ax2ho(ax):
"""Axis angle pair to homochoric vector."""
f = (0.75 * ( ax[3] - np.sin(ax[3]) ))**(1.0/3.0)
ho = ax[0:3] * f
return ho
@staticmethod
def ax2cu(ax):
"""Axis angle pair to cubochoric vector."""
return Rotation.ho2cu(Rotation.ax2ho(ax))
#---------- Rodrigues-Frank vector ----------
@staticmethod
def ro2qu(ro):
"""Rodriques-Frank vector to quaternion."""
return Rotation.ax2qu(Rotation.ro2ax(ro))
@staticmethod
def ro2om(ro):
"""Rodgrigues-Frank vector to rotation matrix."""
return Rotation.ax2om(Rotation.ro2ax(ro))
@staticmethod
def ro2eu(ro):
"""Rodriques-Frank vector to Bunge-Euler angles."""
return Rotation.om2eu(Rotation.ro2om(ro))
@staticmethod
def ro2ax(ro):
"""Rodriques-Frank vector to axis angle pair."""
ta = ro[3]
if iszero(ta):
ax = [ 0.0, 0.0, 1.0, 0.0 ]
elif not np.isfinite(ta):
ax = [ ro[0], ro[1], ro[2], np.pi ]
else:
angle = 2.0*np.arctan(ta)
ta = 1.0/np.linalg.norm(ro[0:3])
ax = [ ro[0]/ta, ro[1]/ta, ro[2]/ta, angle ]
return np.array(ax)
@staticmethod
def ro2ho(ro):
"""Rodriques-Frank vector to homochoric vector."""
if iszero(np.sum(ro[0:3]**2.0)):
ho = [ 0.0, 0.0, 0.0 ]
else:
f = 2.0*np.arctan(ro[3]) -np.sin(2.0*np.arctan(ro[3])) if np.isfinite(ro[3]) else np.pi
ho = ro[0:3] * (0.75*f)**(1.0/3.0)
return np.array(ho)
@staticmethod
def ro2cu(ro):
"""Rodriques-Frank vector to cubochoric vector."""
return Rotation.ho2cu(Rotation.ro2ho(ro))
#---------- Homochoric vector----------
@staticmethod
def ho2qu(ho):
"""Homochoric vector to quaternion."""
return Rotation.ax2qu(Rotation.ho2ax(ho))
@staticmethod
def ho2om(ho):
"""Homochoric vector to rotation matrix."""
return Rotation.ax2om(Rotation.ho2ax(ho))
@staticmethod
def ho2eu(ho):
"""Homochoric vector to Bunge-Euler angles."""
return Rotation.ax2eu(Rotation.ho2ax(ho))
@staticmethod
def ho2ax(ho):
"""Homochoric vector to axis angle pair."""
tfit = np.array([+1.0000000000018852, -0.5000000002194847,
-0.024999992127593126, -0.003928701544781374,
-0.0008152701535450438, -0.0002009500426119712,
-0.00002397986776071756, -0.00008202868926605841,
+0.00012448715042090092, -0.0001749114214822577,
+0.0001703481934140054, -0.00012062065004116828,
+0.000059719705868660826, -0.00001980756723965647,
+0.000003953714684212874, -0.00000036555001439719544])
# normalize h and store the magnitude
hmag_squared = np.sum(ho**2.)
if iszero(hmag_squared):
ax = np.array([ 0.0, 0.0, 1.0, 0.0 ])
else:
hm = hmag_squared
# convert the magnitude to the rotation angle
s = tfit[0] + tfit[1] * hmag_squared
for i in range(2,16):
hm *= hmag_squared
s += tfit[i] * hm
ax = np.append(ho/np.sqrt(hmag_squared),2.0*np.arccos(np.clip(s,-1.0,1.0)))
return ax
@staticmethod
def ho2ro(ho):
"""Axis angle pair to Rodriques-Frank vector."""
return Rotation.ax2ro(Rotation.ho2ax(ho))
@staticmethod
def ho2cu(ho):
"""Homochoric vector to cubochoric vector."""
return Lambert.BallToCube(ho)
#---------- Cubochoric ----------
@staticmethod
def cu2qu(cu):
"""Cubochoric vector to quaternion."""
return Rotation.ho2qu(Rotation.cu2ho(cu))
@staticmethod
def cu2om(cu):
"""Cubochoric vector to rotation matrix."""
return Rotation.ho2om(Rotation.cu2ho(cu))
@staticmethod
def cu2eu(cu):
"""Cubochoric vector to Bunge-Euler angles."""
return Rotation.ho2eu(Rotation.cu2ho(cu))
@staticmethod
def cu2ax(cu):
"""Cubochoric vector to axis angle pair."""
return Rotation.ho2ax(Rotation.cu2ho(cu))
@staticmethod
def cu2ro(cu):
"""Cubochoric vector to Rodriques-Frank vector."""
return Rotation.ho2ro(Rotation.cu2ho(cu))
@staticmethod
def cu2ho(cu):
"""Cubochoric vector to homochoric vector."""
return Lambert.CubeToBall(cu)

312
python/damask/util.py
View File

@ -40,162 +40,210 @@ class bcolors:
self.CROSSOUT = '' self.CROSSOUT = ''
# -----------------------------
def srepr(arg,glue = '\n'): def srepr(arg,glue = '\n'):
"""Joins arguments as individual lines.""" r"""
if (not hasattr(arg, "strip") and Join arguments as individual lines.
(hasattr(arg, "__getitem__") or
hasattr(arg, "__iter__"))): Parameters
return glue.join(str(x) for x in arg) ----------
return arg if isinstance(arg,str) else repr(arg) arg : iterable
Items to join.
glue : str, optional
Defaults to \n.
"""
if (not hasattr(arg, "strip") and
(hasattr(arg, "__getitem__") or
hasattr(arg, "__iter__"))):
return glue.join(str(x) for x in arg)
return arg if isinstance(arg,str) else repr(arg)
# -----------------------------
def croak(what, newline = True): def croak(what, newline = True):
"""Writes formated to stderr.""" """
if what is not None: Write formated to stderr.
sys.stderr.write(srepr(what,glue = '\n') + ('\n' if newline else ''))
sys.stderr.flush() Parameters
----------
what : str or iterable
Content to be displayed
newline : bool, optional
Separate items of what by newline. Defaults to True.
"""
if not what:
sys.stderr.write(srepr(what,glue = '\n') + ('\n' if newline else ''))
sys.stderr.flush()
# -----------------------------
def report(who = None, def report(who = None,
what = None): what = None):
"""Reports script and file name.""" """
croak( (emph(who)+': ' if who is not None else '') + (what if what is not None else '') + '\n' ) Reports script and file name.
DEPRECATED
"""
croak( (emph(who)+': ' if who is not None else '') + (what if what is not None else '') + '\n' )
# -----------------------------
def emph(what): def emph(what):
"""Formats string with emphasis.""" """Formats string with emphasis."""
return bcolors.BOLD+srepr(what)+bcolors.ENDC return bcolors.BOLD+srepr(what)+bcolors.ENDC
# -----------------------------
def deemph(what): def deemph(what):
"""Formats string with deemphasis.""" """Formats string with deemphasis."""
return bcolors.DIM+srepr(what)+bcolors.ENDC return bcolors.DIM+srepr(what)+bcolors.ENDC
# -----------------------------
def delete(what): def delete(what):
"""Formats string as deleted.""" """Formats string as deleted."""
return bcolors.DIM+srepr(what)+bcolors.ENDC return bcolors.DIM+srepr(what)+bcolors.ENDC
# -----------------------------
def strikeout(what): def strikeout(what):
"""Formats string as strikeout.""" """Formats string as strikeout."""
return bcolors.CROSSOUT+srepr(what)+bcolors.ENDC return bcolors.CROSSOUT+srepr(what)+bcolors.ENDC
# -----------------------------
def execute(cmd, def execute(cmd,
streamIn = None, streamIn = None,
wd = './', wd = './',
env = None): env = None):
"""Executes a command in given directory and returns stdout and stderr for optional stdin."""
initialPath = os.getcwd()
os.chdir(wd)
myEnv = os.environ if env is None else env
process = subprocess.Popen(shlex.split(cmd),
stdout = subprocess.PIPE,
stderr = subprocess.PIPE,
stdin = subprocess.PIPE,
env = myEnv)
out,error = [i for i in (process.communicate() if streamIn is None
else process.communicate(streamIn.read().encode('utf-8')))]
out = out.decode('utf-8').replace('\x08','')
error = error.decode('utf-8').replace('\x08','')
os.chdir(initialPath)
if process.returncode != 0: raise RuntimeError('{} failed with returncode {}'.format(cmd,process.returncode))
return out,error
# -----------------------------
class extendableOption(Option):
"""
Used for definition of new option parser action 'extend', which enables to take multiple option arguments.
Adopted from online tutorial http://docs.python.org/library/optparse.html
"""
ACTIONS = Option.ACTIONS + ("extend",)
STORE_ACTIONS = Option.STORE_ACTIONS + ("extend",)
TYPED_ACTIONS = Option.TYPED_ACTIONS + ("extend",)
ALWAYS_TYPED_ACTIONS = Option.ALWAYS_TYPED_ACTIONS + ("extend",)
def take_action(self, action, dest, opt, value, values, parser):
if action == "extend":
lvalue = value.split(",")
values.ensure_value(dest, []).extend(lvalue)
else:
Option.take_action(self, action, dest, opt, value, values, parser)
# Print iterations progress
# from https://gist.github.com/aubricus/f91fb55dc6ba5557fbab06119420dd6a
def progressBar(iteration, total, prefix='', bar_length=50):
"""
Call in a loop to create terminal progress bar.
@params:
iteration - Required : current iteration (Int)
total - Required : total iterations (Int)
prefix - Optional : prefix string (Str)
bar_length - Optional : character length of bar (Int)
"""
fraction = iteration / float(total)
if not hasattr(progressBar, "last_fraction"): # first call to function
progressBar.start_time = time.time()
progressBar.last_fraction = -1.0
remaining_time = ' n/a'
else:
if fraction <= progressBar.last_fraction or iteration == 0: # reset: called within a new loop
progressBar.start_time = time.time()
progressBar.last_fraction = -1.0
remaining_time = ' n/a'
else:
progressBar.last_fraction = fraction
remainder = (total - iteration) * (time.time()-progressBar.start_time)/iteration
remaining_time = '{: 3d}:'.format(int( remainder//3600)) + \
'{:02d}:'.format(int((remainder//60)%60)) + \
'{:02d}' .format(int( remainder %60))
filled_length = int(round(bar_length * fraction))
bar = '█' * filled_length + '░' * (bar_length - filled_length)
sys.stderr.write('\r{} {} {}'.format(prefix, bar, remaining_time)),
if iteration == total: sys.stderr.write('\n')
sys.stderr.flush()
def scale_to_coprime(v):
"""Scale vector to co-prime (relatively prime) integers."""
MAX_DENOMINATOR = 1000
def get_square_denominator(x):
"""Denominator of the square of a number."""
return Fraction(x ** 2).limit_denominator(MAX_DENOMINATOR).denominator
def lcm(a, b):
"""Least common multiple."""
return a * b // np.gcd(a, b)
denominators = [int(get_square_denominator(i)) for i in v]
s = reduce(lcm, denominators) ** 0.5
m = (np.array(v)*s).astype(np.int)
return m//reduce(np.gcd,m)
class return_message():
"""Object with formatted return message."""
def __init__(self,message):
""" """
Sets return message. Execute command.
Parameters Parameters
---------- ----------
message : str or list of str cmd : str
message for output to screen Command to be executed.
streanIn :, optional
Input (via pipe) for executed process.
wd : str, optional
Working directory of process. Defaults to ./ .
env :
Environment
""" """
self.message = message initialPath = os.getcwd()
os.chdir(wd)
myEnv = os.environ if env is None else env
process = subprocess.Popen(shlex.split(cmd),
stdout = subprocess.PIPE,
stderr = subprocess.PIPE,
stdin = subprocess.PIPE,
env = myEnv)
out,error = [i for i in (process.communicate() if streamIn is None
else process.communicate(streamIn.read().encode('utf-8')))]
out = out.decode('utf-8').replace('\x08','')
error = error.decode('utf-8').replace('\x08','')
os.chdir(initialPath)
if process.returncode != 0:
raise RuntimeError('{} failed with returncode {}'.format(cmd,process.returncode))
return out,error
class extendableOption(Option):
"""
Used for definition of new option parser action 'extend', which enables to take multiple option arguments.
Adopted from online tutorial http://docs.python.org/library/optparse.html
DEPRECATED
"""
ACTIONS = Option.ACTIONS + ("extend",)
STORE_ACTIONS = Option.STORE_ACTIONS + ("extend",)
TYPED_ACTIONS = Option.TYPED_ACTIONS + ("extend",)
ALWAYS_TYPED_ACTIONS = Option.ALWAYS_TYPED_ACTIONS + ("extend",)
def take_action(self, action, dest, opt, value, values, parser):
if action == "extend":
lvalue = value.split(",")
values.ensure_value(dest, []).extend(lvalue)
else:
Option.take_action(self, action, dest, opt, value, values, parser)
def progressBar(iteration, total, prefix='', bar_length=50):
"""
Call in a loop to create terminal progress bar.
From https://gist.github.com/aubricus/f91fb55dc6ba5557fbab06119420dd6a
Parameters
----------
iteration : int
Current iteration.
total : int
Total iterations.
prefix : str, optional
Prefix string.
bar_length : int, optional
Character length of bar. Defaults to 50.
"""
fraction = iteration / float(total)
if not hasattr(progressBar, "last_fraction"): # first call to function
progressBar.start_time = time.time()
progressBar.last_fraction = -1.0
remaining_time = ' n/a'
else:
if fraction <= progressBar.last_fraction or iteration == 0: # reset: called within a new loop
progressBar.start_time = time.time()
progressBar.last_fraction = -1.0
remaining_time = ' n/a'
else:
progressBar.last_fraction = fraction
remainder = (total - iteration) * (time.time()-progressBar.start_time)/iteration
remaining_time = '{: 3d}:'.format(int( remainder//3600)) + \
'{:02d}:'.format(int((remainder//60)%60)) + \
'{:02d}' .format(int( remainder %60))
filled_length = int(round(bar_length * fraction))
bar = '█' * filled_length + '░' * (bar_length - filled_length)
sys.stderr.write('\r{} {} {}'.format(prefix, bar, remaining_time)),
if iteration == total:
sys.stderr.write('\n')
sys.stderr.flush()
def scale_to_coprime(v):
"""Scale vector to co-prime (relatively prime) integers."""
MAX_DENOMINATOR = 1000
def get_square_denominator(x):
"""Denominator of the square of a number."""
return Fraction(x ** 2).limit_denominator(MAX_DENOMINATOR).denominator
def lcm(a, b):
"""Least common multiple."""
return a * b // np.gcd(a, b)
denominators = [int(get_square_denominator(i)) for i in v]
s = reduce(lcm, denominators) ** 0.5
m = (np.array(v)*s).astype(np.int)
return m//reduce(np.gcd,m)
class return_message():
"""Object with formatted return message."""
def __init__(self,message):
"""
Sets return message.
Parameters
----------
message : str or list of str
message for output to screen
"""
self.message = message
def __repr__(self):
"""Return message suitable for interactive shells."""
return srepr(self.message)
def __repr__(self):
"""Return message suitable for interactive shells."""
return srepr(self.message)

65
python/tests/test_Orientation.py Normal file
View File

@ -0,0 +1,65 @@
import os
from itertools import permutations
import pytest
import numpy as np
import damask
from damask import Rotation
from damask import Orientation
from damask import Lattice
n = 1000
@pytest.fixture
def default():
"""A set of n random rotations."""
return [Rotation.fromRandom() for r in range(n)]
@pytest.fixture
def reference_dir(reference_dir_base):
"""Directory containing reference results."""
return os.path.join(reference_dir_base,'Rotation')
class TestOrientation:
@pytest.mark.parametrize('color',[{'label':'red', 'RGB':[1,0,0],'direction':[0,0,1]},
{'label':'green','RGB':[0,1,0],'direction':[0,1,1]},
{'label':'blue', 'RGB':[0,0,1],'direction':[1,1,1]}])
@pytest.mark.parametrize('lattice',['fcc','bcc'])
def test_IPF_cubic(self,default,color,lattice):
cube = damask.Orientation(damask.Rotation(),lattice)
for direction in set(permutations(np.array(color['direction']))):
assert np.allclose(cube.IPFcolor(direction),np.array(color['RGB']))
@pytest.mark.parametrize('lattice',Lattice.lattices)
def test_IPF(self,lattice):
direction = np.random.random(3)*2.0-1
for rot in [Rotation.fromRandom() for r in range(n//100)]:
R = damask.Orientation(rot,lattice)
color = R.IPFcolor(direction)
for equivalent in R.equivalentOrientations():
assert np.allclose(color,R.IPFcolor(direction))
@pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])
@pytest.mark.parametrize('lattice',['fcc','bcc'])
def test_relationship_forward_backward(self,model,lattice):
ori = Orientation(Rotation.fromRandom(),lattice)
for i,r in enumerate(ori.relatedOrientations(model)):
ori2 = r.relatedOrientations(model)[i]
misorientation = ori.rotation.misorientation(ori2.rotation)
assert misorientation.asAxisAngle(degrees=True)[3]<1.0e-5
@pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])
@pytest.mark.parametrize('lattice',['fcc','bcc'])
def test_relationship_reference(self,update,reference_dir,model,lattice):
reference = os.path.join(reference_dir,'{}_{}.txt'.format(lattice,model))
ori = Orientation(Rotation(),lattice)
eu = np.array([o.rotation.asEulers(degrees=True) for o in ori.relatedOrientations(model)])
if update:
coords = np.array([(1,i+1) for i,x in enumerate(eu)])
table = damask.Table(eu,{'Eulers':(3,)})
table.add('pos',coords)
table.to_ASCII(reference)
assert np.allclose(eu,damask.Table.from_ASCII(reference).get('Eulers'))

45
python/tests/test_Rotation.py
View File

@ -1,13 +1,9 @@
import os import os
from itertools import permutations
import pytest import pytest
import numpy as np import numpy as np
import damask
from damask import Rotation from damask import Rotation
from damask import Orientation
from damask import Lattice
n = 1000 n = 1000
@ -58,44 +54,3 @@ class TestRotation:
for rot in default: for rot in default:
assert np.allclose(rot.asCubochoric(), assert np.allclose(rot.asCubochoric(),
Rotation.fromQuaternion(rot.asQuaternion()).asCubochoric()) Rotation.fromQuaternion(rot.asQuaternion()).asCubochoric())
@pytest.mark.parametrize('color',[{'label':'red', 'RGB':[1,0,0],'direction':[0,0,1]},
{'label':'green','RGB':[0,1,0],'direction':[0,1,1]},
{'label':'blue', 'RGB':[0,0,1],'direction':[1,1,1]}])
@pytest.mark.parametrize('lattice',['fcc','bcc'])
def test_IPF_cubic(self,default,color,lattice):
cube = damask.Orientation(damask.Rotation(),lattice)
for direction in set(permutations(np.array(color['direction']))):
assert np.allclose(cube.IPFcolor(direction),np.array(color['RGB']))
@pytest.mark.parametrize('lattice',Lattice.lattices)
def test_IPF(self,lattice):
direction = np.random.random(3)*2.0-1
for rot in [Rotation.fromRandom() for r in range(n//100)]:
R = damask.Orientation(rot,lattice)
color = R.IPFcolor(direction)
for equivalent in R.equivalentOrientations():
assert np.allclose(color,R.IPFcolor(direction))
@pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])
@pytest.mark.parametrize('lattice',['fcc','bcc'])
def test_relationship_forward_backward(self,model,lattice):
ori = Orientation(Rotation.fromRandom(),lattice)
for i,r in enumerate(ori.relatedOrientations(model)):
ori2 = r.relatedOrientations(model)[i]
misorientation = ori.rotation.misorientation(ori2.rotation)
assert misorientation.asAxisAngle(degrees=True)[3]<1.0e-5
@pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])
@pytest.mark.parametrize('lattice',['fcc','bcc'])
def test_relationship_reference(self,update,reference_dir,model,lattice):
reference = os.path.join(reference_dir,'{}_{}.txt'.format(lattice,model))
ori = Orientation(Rotation(),lattice)
eu = np.array([o.rotation.asEulers(degrees=True) for o in ori.relatedOrientations(model)])
if update:
coords = np.array([(1,i+1) for i,x in enumerate(eu)])
table = damask.Table(eu,{'Eulers':(3,)})
table.add('pos',coords)
table.to_ASCII(reference)
assert np.allclose(eu,damask.Table.from_ASCII(reference).get('Eulers'))