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

whitespace cleaning and other polishing

This commit is contained in:
Martin Diehl 2020-03-14 21:53:48 +01:00
parent 38b755740b
commit b4679fabfc
6 changed files with 535 additions and 525 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

844
python/damask/colormaps.py
View File

@ -1,6 +1,6 @@
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__ = [
@ -13,544 +13,538 @@ class Color():
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 convert_to(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 express_as(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).convert_to(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):
converted = Color('HSL', HSL) HSL[i+1] = min(HSL[i+1],1.0)
self.model = converted.model HSL[i+1] = max(HSL[i+1],0.0)
self.color = converted.color
converted = Color('HSL', HSL)
self.model = converted.model
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)
self.model = converted.model
self.color = converted.color
converted = Color('XYZ', XYZ)
self.model = converted.model
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)
self.model = converted.model
self.color = converted.color
converted = Color('RGB', RGB)
self.model = converted.model
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): def _XYZ2CIELAB(self):
""" """
Convert CIE XYZ to CIE Lab. Convert CIE XYZ to CIE Lab.
All values are in the range [0,1] All values are in the range [0,1]
from http://en.wikipedia.org/wiki/Lab_color_space, from http://en.wikipedia.org/wiki/Lab_color_space,
http://www.cs.rit.edu/~ncs/color/t_convert.html http://www.cs.rit.edu/~ncs/color/t_convert.html
""" """
if self.model != 'XYZ': return if self.model != 'XYZ': 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 = self.color/ref_white XYZ = self.color/ref_white
for i in range(len(XYZ)): for i in range(len(XYZ)):
if (XYZ[i] > 216./24389 ): XYZ[i] = XYZ[i]**(1.0/3.0) if (XYZ[i] > 216./24389 ): XYZ[i] = XYZ[i]**(1.0/3.0)
else: XYZ[i] = (841./108. * XYZ[i]) + 16.0/116.0 else: XYZ[i] = (841./108. * XYZ[i]) + 16.0/116.0
converted = Color('CIELAB', np.array([ 116.0 * XYZ[1] - 16.0, converted = Color('CIELAB', np.array([ 116.0 * XYZ[1] - 16.0,
500.0 * (XYZ[0] - XYZ[1]), 500.0 * (XYZ[0] - XYZ[1]),
200.0 * (XYZ[1] - XYZ[2]) ])) 200.0 * (XYZ[1] - XYZ[2]) ]))
self.model = converted.model self.model = converted.model
self.color = converted.color self.color = converted.color
def _CIELAB2MSH(self): def _CIELAB2MSH(self):
""" """
Convert CIE Lab to Msh colorspace. Convert CIE Lab to Msh colorspace.
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 != 'CIELAB': return
Msh = np.zeros(3,'d') Msh = np.zeros(3,'d')
Msh[0] = np.sqrt(np.dot(self.color,self.color)) Msh[0] = np.sqrt(np.dot(self.color,self.color))
if (Msh[0] > 0.001): if (Msh[0] > 0.001):
Msh[1] = np.arccos(self.color[0]/Msh[0]) Msh[1] = np.arccos(self.color[0]/Msh[0])
if (self.color[1] != 0.0): if (self.color[1] != 0.0):
Msh[2] = np.arctan2(self.color[2],self.color[1]) Msh[2] = np.arctan2(self.color[2],self.color[1])
converted = Color('MSH', Msh) converted = Color('MSH', Msh)
self.model = converted.model self.model = converted.model
self.color = converted.color self.color = converted.color
def _MSH2CIELAB(self): 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:
"""Perceptually uniform diverging or sequential colormap."""
__slots__ = [
'left',
'right',
'interpolate',
]
__predefined__ = {
'gray': {'left': Color('HSL',[0,1,1]),
'right': Color('HSL',[0,0,0.15]),
'interpolate': 'perceptualuniform'},
'grey': {'left': Color('HSL',[0,1,1]),
'right': Color('HSL',[0,0,0.15]),
'interpolate': 'perceptualuniform'},
'red': {'left': Color('HSL',[0,1,0.14]),
'right': Color('HSL',[0,0.35,0.91]),
'interpolate': 'perceptualuniform'},
'green': {'left': Color('HSL',[0.33333,1,0.14]),
'right': Color('HSL',[0.33333,0.35,0.91]),
'interpolate': 'perceptualuniform'},
'blue': {'left': Color('HSL',[0.66,1,0.14]),
'right': Color('HSL',[0.66,0.35,0.91]),
'interpolate': 'perceptualuniform'},
'seaweed': {'left': Color('HSL',[0.78,1.0,0.1]),
'right': Color('HSL',[0.40000,0.1,0.9]),
'interpolate': 'perceptualuniform'},
'bluebrown': {'left': Color('HSL',[0.65,0.53,0.49]),
'right': Color('HSL',[0.11,0.75,0.38]),
'interpolate': 'perceptualuniform'},
'redgreen': {'left': Color('HSL',[0.97,0.96,0.36]),
'right': Color('HSL',[0.33333,1.0,0.14]),
'interpolate': 'perceptualuniform'},
'bluered': {'left': Color('HSL',[0.65,0.53,0.49]),
'right': Color('HSL',[0.97,0.96,0.36]),
'interpolate': 'perceptualuniform'},
'blueredrainbow':{'left': Color('HSL',[2.0/3.0,1,0.5]),
'right': Color('HSL',[0,1,0.5]),
'interpolate': 'linear' },
'orientation': {'left': Color('RGB',[0.933334,0.878432,0.878431]),
'right': Color('RGB',[0.250980,0.007843,0.000000]),
'interpolate': 'perceptualuniform'},
'strain': {'left': Color('RGB',[0.941177,0.941177,0.870588]),
'right': Color('RGB',[0.266667,0.266667,0.000000]),
'interpolate': 'perceptualuniform'},
'stress': {'left': Color('RGB',[0.878432,0.874511,0.949019]),
'right': Color('RGB',[0.000002,0.000000,0.286275]),
'interpolate': 'perceptualuniform'},
}
# ------------------------------------------------------------------
def __init__(self,
left = Color('RGB',[1,1,1]),
right = Color('RGB',[0,0,0]),
interpolate = 'perceptualuniform',
predefined = None
):
""" """
Convert Msh colorspace to CIE Lab. Create a Colormap object.
Parameters
----------
left : Color
left color (minimum value)
right : Color
right color (maximum value)
interpolate : str
interpolation scheme (either 'perceptualuniform' or 'linear')
predefined : bool
ignore other arguments and use predefined definition
with s,h in radians
from http://www.cs.unm.edu/~kmorel/documents/ColorMaps/DivergingColorMapWorkshop.xls
""" """
if self.model != 'MSH': return if predefined is not None:
left = self.__predefined__[predefined.lower()]['left']
right= self.__predefined__[predefined.lower()]['right']
interpolate = self.__predefined__[predefined.lower()]['interpolate']
Lab = np.zeros(3,'d') if left.__class__.__name__ != 'Color':
Lab[0] = self.color[0] * np.cos(self.color[1]) left = Color()
Lab[1] = self.color[0] * np.sin(self.color[1]) * np.cos(self.color[2]) if right.__class__.__name__ != 'Color':
Lab[2] = self.color[0] * np.sin(self.color[1]) * np.sin(self.color[2]) right = Color()
converted = Color('CIELAB', Lab) self.left = left
self.model = converted.model self.right = right
self.color = converted.color self.interpolate = interpolate
class Colormap(): # ------------------------------------------------------------------
"""Perceptually uniform diverging or sequential colormap.""" def __repr__(self):
"""Left and right value of colormap."""
__slots__ = [ return 'Left: %s Right: %s'%(self.left,self.right)
'left',
'right',
'interpolate',
]
__predefined__ = {
'gray': {'left': Color('HSL',[0,1,1]),
'right': Color('HSL',[0,0,0.15]),
'interpolate': 'perceptualuniform'},
'grey': {'left': Color('HSL',[0,1,1]),
'right': Color('HSL',[0,0,0.15]),
'interpolate': 'perceptualuniform'},
'red': {'left': Color('HSL',[0,1,0.14]),
'right': Color('HSL',[0,0.35,0.91]),
'interpolate': 'perceptualuniform'},
'green': {'left': Color('HSL',[0.33333,1,0.14]),
'right': Color('HSL',[0.33333,0.35,0.91]),
'interpolate': 'perceptualuniform'},
'blue': {'left': Color('HSL',[0.66,1,0.14]),
'right': Color('HSL',[0.66,0.35,0.91]),
'interpolate': 'perceptualuniform'},
'seaweed': {'left': Color('HSL',[0.78,1.0,0.1]),
'right': Color('HSL',[0.40000,0.1,0.9]),
'interpolate': 'perceptualuniform'},
'bluebrown': {'left': Color('HSL',[0.65,0.53,0.49]),
'right': Color('HSL',[0.11,0.75,0.38]),
'interpolate': 'perceptualuniform'},
'redgreen': {'left': Color('HSL',[0.97,0.96,0.36]),
'right': Color('HSL',[0.33333,1.0,0.14]),
'interpolate': 'perceptualuniform'},
'bluered': {'left': Color('HSL',[0.65,0.53,0.49]),
'right': Color('HSL',[0.97,0.96,0.36]),
'interpolate': 'perceptualuniform'},
'blueredrainbow':{'left': Color('HSL',[2.0/3.0,1,0.5]),
'right': Color('HSL',[0,1,0.5]),
'interpolate': 'linear' },
'orientation': {'left': Color('RGB',[0.933334,0.878432,0.878431]),
'right': Color('RGB',[0.250980,0.007843,0.000000]),
'interpolate': 'perceptualuniform'},
'strain': {'left': Color('RGB',[0.941177,0.941177,0.870588]),
'right': Color('RGB',[0.266667,0.266667,0.000000]),
'interpolate': 'perceptualuniform'},
'stress': {'left': Color('RGB',[0.878432,0.874511,0.949019]),
'right': Color('RGB',[0.000002,0.000000,0.286275]),
'interpolate': 'perceptualuniform'},
}
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def __init__(self, def invert(self):
left = Color('RGB',[1,1,1]), """Switch left/minimum with right/maximum."""
right = Color('RGB',[0,0,0]), (self.left, self.right) = (self.right, self.left)
interpolate = 'perceptualuniform', return self
predefined = None
):
"""
Create a Colormap object.
Parameters
----------
left : Color
left color (minimum value)
right : Color
right color (maximum value)
interpolate : str
interpolation scheme (either 'perceptualuniform' or 'linear')
predefined : bool
ignore other arguments and use predefined definition
"""
if predefined is not None:
left = self.__predefined__[predefined.lower()]['left']
right= self.__predefined__[predefined.lower()]['right']
interpolate = self.__predefined__[predefined.lower()]['interpolate']
if left.__class__.__name__ != 'Color':
left = Color()
if right.__class__.__name__ != 'Color':
right = Color()
self.left = left
self.right = right
self.interpolate = interpolate
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def __repr__(self): def show_predefined(self):
"""Left and right value of colormap.""" """Show the labels of the predefined colormaps."""
return 'Left: %s Right: %s'%(self.left,self.right) print('\n'.join(self.__predefined__.keys()))
# ------------------------------------------------------------------
def color(self,fraction = 0.5):
# ------------------------------------------------------------------ def interpolate_Msh(lo, hi, frac):
def invert(self):
"""Switch left/minimum with right/maximum."""
(self.left, self.right) = (self.right, self.left)
return self
def rad_diff(a,b):
return abs(a[2]-b[2])
# if saturation of one of the two colors is too less than the other, hue of the less
def adjust_hue(Msh_sat, Msh_unsat):
if Msh_sat[0] >= Msh_unsat[0]:
return Msh_sat[2]
else:
hSpin = Msh_sat[1]/np.sin(Msh_sat[1])*np.sqrt(Msh_unsat[0]**2.0-Msh_sat[0]**2)/Msh_sat[0]
if Msh_sat[2] < - np.pi/3.0: hSpin *= -1.0
return Msh_sat[2] + hSpin
# ------------------------------------------------------------------ Msh1 = np.array(lo[:])
def show_predefined(self): Msh2 = np.array(hi[:])
"""Show the labels of the predefined colormaps."""
print('\n'.join(self.__predefined__.keys()))
# ------------------------------------------------------------------ if (Msh1[1] > 0.05 and Msh2[1] > 0.05 and rad_diff(Msh1,Msh2) > np.pi/3.0):
def color(self,fraction = 0.5): M_mid = max(Msh1[0],Msh2[0],88.0)
if frac < 0.5:
Msh2 = np.array([M_mid,0.0,0.0],'d')
frac *= 2.0
else:
Msh1 = np.array([M_mid,0.0,0.0],'d')
frac = 2.0*frac - 1.0
if Msh1[1] < 0.05 and Msh2[1] > 0.05: Msh1[2] = adjust_hue(Msh2,Msh1)
elif Msh1[1] > 0.05 and Msh2[1] < 0.05: Msh2[2] = adjust_hue(Msh1,Msh2)
Msh = (1.0 - frac) * Msh1 + frac * Msh2
def interpolate_Msh(lo, hi, frac): return Color('MSH',Msh)
def rad_diff(a,b): def interpolate_linear(lo, hi, frac):
return abs(a[2]-b[2]) """Linear interpolation between lo and hi color at given fraction; output in model of lo color."""
# if saturation of one of the two colors is too less than the other, hue of the less interpolation = (1.0 - frac) * np.array(lo.color[:]) \
def adjust_hue(Msh_sat, Msh_unsat): + frac * np.array(hi.express_as(lo.model).color[:])
if Msh_sat[0] >= Msh_unsat[0]:
return Msh_sat[2]
else:
hSpin = Msh_sat[1]/np.sin(Msh_sat[1])*np.sqrt(Msh_unsat[0]**2.0-Msh_sat[0]**2)/Msh_sat[0]
if Msh_sat[2] < - np.pi/3.0: hSpin *= -1.0
return Msh_sat[2] + hSpin
Msh1 = np.array(lo[:]) return Color(lo.model,interpolation)
Msh2 = np.array(hi[:])
if (Msh1[1] > 0.05 and Msh2[1] > 0.05 and rad_diff(Msh1,Msh2) > np.pi/3.0): if self.interpolate == 'perceptualuniform':
M_mid = max(Msh1[0],Msh2[0],88.0) return interpolate_Msh(self.left.express_as('MSH').color,
if frac < 0.5: self.right.express_as('MSH').color,fraction)
Msh2 = np.array([M_mid,0.0,0.0],'d') elif self.interpolate == 'linear':
frac *= 2.0 return interpolate_linear(self.left,
else: self.right,fraction)
Msh1 = np.array([M_mid,0.0,0.0],'d') else:
frac = 2.0*frac - 1.0 raise NameError('unknown color interpolation method')
if Msh1[1] < 0.05 and Msh2[1] > 0.05: Msh1[2] = adjust_hue(Msh2,Msh1)
elif Msh1[1] > 0.05 and Msh2[1] < 0.05: Msh2[2] = adjust_hue(Msh1,Msh2)
Msh = (1.0 - frac) * Msh1 + frac * Msh2
return Color('MSH',Msh) # ------------------------------------------------------------------
def export(self,name = 'uniformPerceptualColorMap',\
format = 'paraview',\
steps = 2,\
crop = [-1.0,1.0],
model = 'RGB'):
"""
[RGB] colormap for use in paraview or gmsh, or as raw string, or array.
def interpolate_linear(lo, hi, frac): Arguments: name, format, steps, crop.
"""Linear interpolation between lo and hi color at given fraction; output in model of lo color.""" Format is one of (paraview, gmsh, raw, list).
interpolation = (1.0 - frac) * np.array(lo.color[:]) \ Crop selects a (sub)range in [-1.0,1.0].
+ frac * np.array(hi.express_as(lo.model).color[:]) Generates sequential map if one limiting color is either white or black,
diverging map otherwise.
"""
format = format.lower() # consistent comparison basis
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]).express_as(model).color for i in range(steps)]
if format == 'paraview':
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],) \
for i,color in enumerate(colors[:-1])] \
+ [' {:4d},{:8.6f},{:8.6f},{:8.6f} '.format(len(colors),colors[-1][0],colors[-1][1],colors[-1][2],)] \
+ [' ]\n }\n]']
return Color(lo.model,interpolation) elif format == 'gmsh':
colormap = ['View.ColorTable = {'] \
+ [',\n'.join(['{%s}'%(','.join([str(x*255.0) for x in color])) for color in colors])] \
+ ['}']
if self.interpolate == 'perceptualuniform': elif format == 'gom':
return interpolate_Msh(self.left.express_as('MSH').color, colormap = ['1 1 ' + str(name)
self.right.express_as('MSH').color,fraction) + ' 9 ' + str(name)
elif self.interpolate == 'linear': + ' 0 1 0 3 0 0 -1 9 \\ 0 0 0 255 255 255 0 0 255 '
return interpolate_linear(self.left, + '30 NO_UNIT 1 1 64 64 64 255 1 0 0 0 0 0 0 3 0 ' + str(len(colors))
self.right,fraction) + ' '.join([' 0 %s 255 1'%(' '.join([str(int(x*255.0)) for x in color])) for color in reversed(colors)])]
else:
raise NameError('unknown color interpolation method')
# ------------------------------------------------------------------ elif format == 'raw':
def export(self,name = 'uniformPerceptualColorMap',\ colormap = ['\t'.join(map(str,color)) for color in colors]
format = 'paraview',\
steps = 2,\
crop = [-1.0,1.0],
model = 'RGB'):
"""
[RGB] colormap for use in paraview or gmsh, or as raw string, or array.
Arguments: name, format, steps, crop. elif format == 'list':
Format is one of (paraview, gmsh, raw, list). colormap = colors
Crop selects a (sub)range in [-1.0,1.0].
Generates sequential map if one limiting color is either white or black,
diverging map otherwise.
"""
format = format.lower() # consistent comparison basis
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]).express_as(model).color for i in range(steps)]
if format == 'paraview':
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],) \
for i,color in enumerate(colors[:-1])] \
+ [' {:4d},{:8.6f},{:8.6f},{:8.6f} '.format(len(colors),colors[-1][0],colors[-1][1],colors[-1][2],)] \
+ [' ]\n }\n]']
elif format == 'gmsh': else:
colormap = ['View.ColorTable = {'] \ raise NameError('unknown color export format')
+ [',\n'.join(['{%s}'%(','.join([str(x*255.0) for x in color])) for color in colors])] \
+ ['}']
elif format == 'gom': return '\n'.join(colormap) + '\n' if type(colormap[0]) is str else colormap
colormap = ['1 1 ' + str(name)
+ ' 9 ' + str(name)
+ ' 0 1 0 3 0 0 -1 9 \\ 0 0 0 255 255 255 0 0 255 '
+ '30 NO_UNIT 1 1 64 64 64 255 1 0 0 0 0 0 0 3 0 ' + str(len(colors))
+ ' '.join([' 0 %s 255 1'%(' '.join([str(int(x*255.0)) for x in color])) for color in reversed(colors)])]
elif format == 'raw':
colormap = ['\t'.join(map(str,color)) for color in colors]
elif format == 'list':
colormap = colors
else:
raise NameError('unknown color export format')
return '\n'.join(colormap) + '\n' if type(colormap[0]) is str else colormap

58
python/damask/geom.py
View File

@ -36,7 +36,7 @@ class Geom():
self.set_origin(origin) self.set_origin(origin)
self.set_homogenization(homogenization) self.set_homogenization(homogenization)
self.set_comments(comments) self.set_comments(comments)
def __repr__(self): def __repr__(self):
"""Basic information on geometry definition.""" """Basic information on geometry definition."""
@ -49,6 +49,7 @@ class Geom():
'max microstructure: {}'.format(np.nanmax(self.microstructure)), 'max microstructure: {}'.format(np.nanmax(self.microstructure)),
]) ])
def update(self,microstructure=None,size=None,origin=None,rescale=False): def update(self,microstructure=None,size=None,origin=None,rescale=False):
""" """
Updates microstructure and size. Updates microstructure and size.
@ -70,7 +71,7 @@ class Geom():
origin_old = self.get_origin() origin_old = self.get_origin()
unique_old = len(np.unique(self.microstructure)) unique_old = len(np.unique(self.microstructure))
max_old = np.nanmax(self.microstructure) max_old = np.nanmax(self.microstructure)
if size is not None and rescale: if size is not None and rescale:
raise ValueError('Either set size explicitly or rescale automatically') raise ValueError('Either set size explicitly or rescale automatically')
@ -108,9 +109,10 @@ class Geom():
if max_old != np.nanmax(self.microstructure): if max_old != np.nanmax(self.microstructure):
message[-1] = util.delete(message[-1]) message[-1] = util.delete(message[-1])
message.append(util.emph('max microstructure: {}'.format(np.nanmax(self.microstructure)))) message.append(util.emph('max microstructure: {}'.format(np.nanmax(self.microstructure))))
return util.return_message(message) return util.return_message(message)
def set_comments(self,comments): def set_comments(self,comments):
""" """
Replaces all existing comments. Replaces all existing comments.
@ -123,7 +125,8 @@ class Geom():
""" """
self.comments = [] self.comments = []
self.add_comments(comments) self.add_comments(comments)
def add_comments(self,comments): def add_comments(self,comments):
""" """
Appends comments to existing comments. Appends comments to existing comments.
@ -136,6 +139,7 @@ class Geom():
""" """
self.comments += [str(c) for c in comments] if isinstance(comments,list) else [str(comments)] self.comments += [str(c) for c in comments] if isinstance(comments,list) else [str(comments)]
def set_microstructure(self,microstructure): def set_microstructure(self,microstructure):
""" """
Replaces the existing microstructure representation. Replaces the existing microstructure representation.
@ -154,6 +158,7 @@ class Geom():
else: else:
self.microstructure = np.copy(microstructure) self.microstructure = np.copy(microstructure)
def set_size(self,size): def set_size(self,size):
""" """
Replaces the existing size information. Replaces the existing size information.
@ -173,6 +178,7 @@ class Geom():
else: else:
self.size = np.array(size) self.size = np.array(size)
def set_origin(self,origin): def set_origin(self,origin):
""" """
Replaces the existing origin information. Replaces the existing origin information.
@ -189,6 +195,7 @@ class Geom():
else: else:
self.origin = np.array(origin) self.origin = np.array(origin)
def set_homogenization(self,homogenization): def set_homogenization(self,homogenization):
""" """
Replaces the existing homogenization index. Replaces the existing homogenization index.
@ -205,34 +212,42 @@ class Geom():
else: else:
self.homogenization = homogenization self.homogenization = homogenization
@property @property
def grid(self): def grid(self):
return self.get_grid() return self.get_grid()
def get_microstructure(self): def get_microstructure(self):
"""Return the microstructure representation.""" """Return the microstructure representation."""
return np.copy(self.microstructure) return np.copy(self.microstructure)
def get_size(self): def get_size(self):
"""Return the physical size in meter.""" """Return the physical size in meter."""
return np.copy(self.size) return np.copy(self.size)
def get_origin(self): def get_origin(self):
"""Return the origin in meter.""" """Return the origin in meter."""
return np.copy(self.origin) return np.copy(self.origin)
def get_grid(self): def get_grid(self):
"""Return the grid discretization.""" """Return the grid discretization."""
return np.array(self.microstructure.shape) return np.array(self.microstructure.shape)
def get_homogenization(self): def get_homogenization(self):
"""Return the homogenization index.""" """Return the homogenization index."""
return self.homogenization return self.homogenization
def get_comments(self): def get_comments(self):
"""Return the comments.""" """Return the comments."""
return self.comments[:] return self.comments[:]
def get_header(self): def get_header(self):
"""Return the full header (grid, size, origin, homogenization, comments).""" """Return the full header (grid, size, origin, homogenization, comments)."""
header = ['{} header'.format(len(self.comments)+4)] + self.comments header = ['{} header'.format(len(self.comments)+4)] + self.comments
@ -241,7 +256,8 @@ class Geom():
header.append('origin x {} y {} z {}'.format(*self.get_origin())) header.append('origin x {} y {} z {}'.format(*self.get_origin()))
header.append('homogenization {}'.format(self.get_homogenization())) header.append('homogenization {}'.format(self.get_homogenization()))
return header return header
@staticmethod @staticmethod
def from_file(fname): def from_file(fname):
""" """
@ -266,7 +282,7 @@ class Geom():
if not keyword.startswith('head') or header_length < 3: if not keyword.startswith('head') or header_length < 3:
raise TypeError('Header length information missing or invalid') raise TypeError('Header length information missing or invalid')
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 items else '' key = items[0] if items else ''
@ -295,14 +311,14 @@ class Geom():
else: items = list(map(float,items)) else: items = list(map(float,items))
microstructure[i:i+len(items)] = items microstructure[i:i+len(items)] = items
i += len(items) i += len(items)
if i != grid.prod(): if i != grid.prod():
raise TypeError('Invalid file: expected {} entries, found {}'.format(grid.prod(),i)) raise TypeError('Invalid file: expected {} entries, found {}'.format(grid.prod(),i))
microstructure = microstructure.reshape(grid,order='F') microstructure = microstructure.reshape(grid,order='F')
if not np.any(np.mod(microstructure.flatten(),1) != 0.0): # no float present if not np.any(np.mod(microstructure.flatten(),1) != 0.0): # no float present
microstructure = microstructure.astype('int') microstructure = microstructure.astype('int')
return Geom(microstructure.reshape(grid),size,origin,homogenization,comments) return Geom(microstructure.reshape(grid),size,origin,homogenization,comments)
@ -320,7 +336,7 @@ class Geom():
""" """
header = self.get_header() header = self.get_header()
grid = self.get_grid() grid = self.get_grid()
if pack is None: if pack is None:
plain = grid.prod()/np.unique(self.microstructure).size < 250 plain = grid.prod()/np.unique(self.microstructure).size < 250
else: else:
@ -371,7 +387,7 @@ class Geom():
elif compressType == 'of': elif compressType == 'of':
f.write('{} of {}\n'.format(reps,former)) f.write('{} of {}\n'.format(reps,former))
def to_vtk(self,fname=None): def to_vtk(self,fname=None):
""" """
Generates vtk file. Generates vtk file.
@ -391,7 +407,7 @@ class Geom():
np.linspace(0,size[1],grid[1]) + origin[1], np.linspace(0,size[1],grid[1]) + origin[1],
np.linspace(0,size[2],grid[2]) + origin[2] np.linspace(0,size[2],grid[2]) + origin[2]
] ]
rGrid = vtk.vtkRectilinearGrid() rGrid = vtk.vtkRectilinearGrid()
coordArray = [vtk.vtkDoubleArray(),vtk.vtkDoubleArray(),vtk.vtkDoubleArray()] coordArray = [vtk.vtkDoubleArray(),vtk.vtkDoubleArray(),vtk.vtkDoubleArray()]
@ -403,7 +419,7 @@ class Geom():
rGrid.SetXCoordinates(coordArray[0]) rGrid.SetXCoordinates(coordArray[0])
rGrid.SetYCoordinates(coordArray[1]) rGrid.SetYCoordinates(coordArray[1])
rGrid.SetZCoordinates(coordArray[2]) rGrid.SetZCoordinates(coordArray[2])
ms = numpy_support.numpy_to_vtk(num_array=self.microstructure.flatten(order='F'), ms = numpy_support.numpy_to_vtk(num_array=self.microstructure.flatten(order='F'),
array_type=vtk.VTK_INT if self.microstructure.dtype == int else vtk.VTK_FLOAT) array_type=vtk.VTK_INT if self.microstructure.dtype == int else vtk.VTK_FLOAT)
ms.SetName('microstructure') ms.SetName('microstructure')
@ -418,7 +434,7 @@ class Geom():
writer = vtk.vtkXMLRectilinearGridWriter() writer = vtk.vtkXMLRectilinearGridWriter()
writer.SetCompressorTypeToZLib() writer.SetCompressorTypeToZLib()
writer.SetDataModeToBinary() writer.SetDataModeToBinary()
ext = os.path.splitext(fname)[1] ext = os.path.splitext(fname)[1]
if ext == '': if ext == '':
name = fname + '.' + writer.GetDefaultFileExtension() name = fname + '.' + writer.GetDefaultFileExtension()
@ -427,13 +443,13 @@ class Geom():
else: else:
raise ValueError("unknown extension {}".format(ext)) raise ValueError("unknown extension {}".format(ext))
writer.SetFileName(name) writer.SetFileName(name)
writer.SetInputData(rGrid) writer.SetInputData(rGrid)
writer.Write() writer.Write()
if fname is None: return writer.GetOutputString() if fname is None: return writer.GetOutputString()
def show(self): def show(self):
"""Show raw content (as in file).""" """Show raw content (as in file)."""
f=StringIO() f=StringIO()
@ -469,9 +485,9 @@ class Geom():
ms = np.concatenate([ms,ms[:,limits[0]:limits[1]:-1,:]],1) ms = np.concatenate([ms,ms[:,limits[0]:limits[1]:-1,:]],1)
if 'x' in directions: if 'x' in directions:
ms = np.concatenate([ms,ms[limits[0]:limits[1]:-1,:,:]],0) ms = np.concatenate([ms,ms[limits[0]:limits[1]:-1,:,:]],0)
#self.add_comments('geom.py:mirror v{}'.format(version)
return self.update(ms,rescale=True) return self.update(ms,rescale=True)
#self.add_comments('tbd')
def scale(self,grid): def scale(self,grid):
@ -484,6 +500,7 @@ class Geom():
new grid dimension new grid dimension
""" """
#self.add_comments('geom.py:scale v{}'.format(version)
return self.update( return self.update(
ndimage.interpolation.zoom( ndimage.interpolation.zoom(
self.microstructure, self.microstructure,
@ -494,7 +511,6 @@ class Geom():
prefilter=False prefilter=False
) )
) )
#self.add_comments('tbd')
def clean(self,stencil=3): def clean(self,stencil=3):
@ -511,13 +527,13 @@ class Geom():
unique, inverse = np.unique(arr, return_inverse=True) unique, inverse = np.unique(arr, return_inverse=True)
return unique[np.argmax(np.bincount(inverse))] return unique[np.argmax(np.bincount(inverse))]
#self.add_comments('geom.py:clean v{}'.format(version)
return self.update(ndimage.filters.generic_filter( return self.update(ndimage.filters.generic_filter(
self.microstructure, self.microstructure,
mostFrequent, mostFrequent,
size=(stencil,)*3 size=(stencil,)*3
).astype(self.microstructure.dtype) ).astype(self.microstructure.dtype)
) )
#self.add_comments('tbd')
def renumber(self): def renumber(self):
@ -526,5 +542,5 @@ class Geom():
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)
#self.add_comments('geom.py:renumber v{}'.format(version)
return self.update(renumbered) return self.update(renumbered)
#self.add_comments('tbd')

52
python/damask/mechanics.py
View File

@ -9,9 +9,9 @@ def Cauchy(P,F):
Parameters Parameters
---------- ----------
F : numpy.ndarray of shape (:,3,3) or (3,3) F : numpy.ndarray of shape (:,3,3) or (3,3)
Deformation gradient. Deformation gradient.
P : numpy.ndarray of shape (:,3,3) or (3,3) P : numpy.ndarray of shape (:,3,3) or (3,3)
1. Piola-Kirchhoff stress. First Piola-Kirchhoff stress.
""" """
if np.shape(F) == np.shape(P) == (3,3): if np.shape(F) == np.shape(P) == (3,3):
@ -28,7 +28,7 @@ def deviatoric_part(T):
Parameters Parameters
---------- ----------
T : numpy.ndarray of shape (:,3,3) or (3,3) T : numpy.ndarray of shape (:,3,3) or (3,3)
Tensor of which the deviatoric part is computed. Tensor of which the deviatoric part is computed.
""" """
return T - np.eye(3)*spherical_part(T) if np.shape(T) == (3,3) else \ return T - np.eye(3)*spherical_part(T) if np.shape(T) == (3,3) else \
@ -45,7 +45,7 @@ def eigenvalues(T_sym):
Parameters Parameters
---------- ----------
T_sym : numpy.ndarray of shape (:,3,3) or (3,3) T_sym : numpy.ndarray of shape (:,3,3) or (3,3)
Symmetric tensor of which the eigenvalues are computed. Symmetric tensor of which the eigenvalues are computed.
""" """
return np.linalg.eigvalsh(symmetric(T_sym)) return np.linalg.eigvalsh(symmetric(T_sym))
@ -60,9 +60,9 @@ def eigenvectors(T_sym,RHS=False):
Parameters Parameters
---------- ----------
T_sym : numpy.ndarray of shape (:,3,3) or (3,3) T_sym : numpy.ndarray of shape (:,3,3) or (3,3)
Symmetric tensor of which the eigenvectors are computed. Symmetric tensor of which the eigenvectors are computed.
RHS: bool, optional RHS: bool, optional
Enforce right-handed coordinate system. Default is False. Enforce right-handed coordinate system. Default is False.
""" """
(u,v) = np.linalg.eigh(symmetric(T_sym)) (u,v) = np.linalg.eigh(symmetric(T_sym))
@ -82,7 +82,7 @@ def left_stretch(T):
Parameters Parameters
---------- ----------
T : numpy.ndarray of shape (:,3,3) or (3,3) T : numpy.ndarray of shape (:,3,3) or (3,3)
Tensor of which the left stretch is computed. Tensor of which the left stretch is computed.
""" """
return __polar_decomposition(T,'V')[0] return __polar_decomposition(T,'V')[0]
@ -95,7 +95,7 @@ def maximum_shear(T_sym):
Parameters Parameters
---------- ----------
T_sym : numpy.ndarray of shape (:,3,3) or (3,3) T_sym : numpy.ndarray of shape (:,3,3) or (3,3)
Symmetric tensor of which the maximum shear is computed. Symmetric tensor of which the maximum shear is computed.
""" """
w = eigenvalues(T_sym) w = eigenvalues(T_sym)
@ -110,7 +110,7 @@ def Mises_strain(epsilon):
Parameters Parameters
---------- ----------
epsilon : numpy.ndarray of shape (:,3,3) or (3,3) epsilon : numpy.ndarray of shape (:,3,3) or (3,3)
Symmetric strain tensor of which the von Mises equivalent is computed. Symmetric strain tensor of which the von Mises equivalent is computed.
""" """
return __Mises(epsilon,2.0/3.0) return __Mises(epsilon,2.0/3.0)
@ -123,7 +123,7 @@ def Mises_stress(sigma):
Parameters Parameters
---------- ----------
sigma : numpy.ndarray of shape (:,3,3) or (3,3) sigma : numpy.ndarray of shape (:,3,3) or (3,3)
Symmetric stress tensor of which the von Mises equivalent is computed. Symmetric stress tensor of which the von Mises equivalent is computed.
""" """
return __Mises(sigma,3.0/2.0) return __Mises(sigma,3.0/2.0)
@ -136,9 +136,9 @@ def PK2(P,F):
Parameters Parameters
---------- ----------
P : numpy.ndarray of shape (:,3,3) or (3,3) P : numpy.ndarray of shape (:,3,3) or (3,3)
1. Piola-Kirchhoff stress. First Piola-Kirchhoff stress.
F : numpy.ndarray of shape (:,3,3) or (3,3) F : numpy.ndarray of shape (:,3,3) or (3,3)
Deformation gradient. Deformation gradient.
""" """
if np.shape(F) == np.shape(P) == (3,3): if np.shape(F) == np.shape(P) == (3,3):
@ -155,7 +155,7 @@ def right_stretch(T):
Parameters Parameters
---------- ----------
T : numpy.ndarray of shape (:,3,3) or (3,3) T : numpy.ndarray of shape (:,3,3) or (3,3)
Tensor of which the right stretch is computed. Tensor of which the right stretch is computed.
""" """
return __polar_decomposition(T,'U')[0] return __polar_decomposition(T,'U')[0]
@ -168,7 +168,7 @@ def rotational_part(T):
Parameters Parameters
---------- ----------
T : numpy.ndarray of shape (:,3,3) or (3,3) T : numpy.ndarray of shape (:,3,3) or (3,3)
Tensor of which the rotational part is computed. Tensor of which the rotational part is computed.
""" """
return __polar_decomposition(T,'R')[0] return __polar_decomposition(T,'R')[0]
@ -181,9 +181,9 @@ def spherical_part(T,tensor=False):
Parameters Parameters
---------- ----------
T : numpy.ndarray of shape (:,3,3) or (3,3) T : numpy.ndarray of shape (:,3,3) or (3,3)
Tensor of which the hydrostatic part is computed. Tensor of which the hydrostatic part is computed.
tensor : bool, optional tensor : bool, optional
Map spherical part onto identity tensor. Default is false Map spherical part onto identity tensor. Default is false
""" """
if T.shape == (3,3): if T.shape == (3,3):
@ -207,11 +207,11 @@ def strain_tensor(F,t,m):
Parameters Parameters
---------- ----------
F : numpy.ndarray of shape (:,3,3) or (3,3) F : numpy.ndarray of shape (:,3,3) or (3,3)
Deformation gradient. Deformation gradient.
t : {‘V’, ‘U’} t : {‘V’, ‘U’}
Type of the polar decomposition, ‘V’ for left stretch tensor and ‘U’ for right stretch tensor. Type of the polar decomposition, ‘V’ for left stretch tensor and ‘U’ for right stretch tensor.
m : float m : float
Order of the strain. Order of the strain.
""" """
F_ = F.reshape((1,3,3)) if F.shape == (3,3) else F F_ = F.reshape((1,3,3)) if F.shape == (3,3) else F
@ -242,7 +242,7 @@ def symmetric(T):
Parameters Parameters
---------- ----------
T : numpy.ndarray of shape (:,3,3) or (3,3) T : numpy.ndarray of shape (:,3,3) or (3,3)
Tensor of which the symmetrized values are computed. Tensor of which the symmetrized values are computed.
""" """
return (T+transpose(T))*0.5 return (T+transpose(T))*0.5
@ -255,7 +255,7 @@ def transpose(T):
Parameters Parameters
---------- ----------
T : numpy.ndarray of shape (:,3,3) or (3,3) T : numpy.ndarray of shape (:,3,3) or (3,3)
Tensor of which the transpose is computed. Tensor of which the transpose is computed.
""" """
return T.T if np.shape(T) == (3,3) else \ return T.T if np.shape(T) == (3,3) else \
@ -269,10 +269,10 @@ def __polar_decomposition(T,requested):
Parameters Parameters
---------- ----------
T : numpy.ndarray of shape (:,3,3) or (3,3) T : numpy.ndarray of shape (:,3,3) or (3,3)
Tensor of which the singular values are computed. Tensor of which the singular values are computed.
requested : iterable of str requested : iterable of str
Requested outputs: ‘R’ for the rotation tensor, Requested outputs: ‘R’ for the rotation tensor,
‘V’ for left stretch tensor and ‘U’ for right stretch tensor. ‘V’ for left stretch tensor and ‘U’ for right stretch tensor.
""" """
u, s, vh = np.linalg.svd(T) u, s, vh = np.linalg.svd(T)
@ -297,9 +297,9 @@ def __Mises(T_sym,s):
Parameters Parameters
---------- ----------
T_sym : numpy.ndarray of shape (:,3,3) or (3,3) T_sym : numpy.ndarray of shape (:,3,3) or (3,3)
Symmetric tensor of which the von Mises equivalent is computed. Symmetric tensor of which the von Mises equivalent is computed.
s : float s : float
Scaling factor (2/3 for strain, 3/2 for stress). Scaling factor (2/3 for strain, 3/2 for stress).
""" """
d = deviatoric_part(T_sym) d = deviatoric_part(T_sym)

54
python/damask/orientation.py
View File

@ -50,23 +50,23 @@ class Orientation:
Look into A. Heinz and P. Neumann 1991 for cases with differing sym.) Look into A. Heinz and P. Neumann 1991 for cases with differing sym.)
""" """
if self.lattice.symmetry != other.lattice.symmetry: if self.lattice.symmetry != other.lattice.symmetry:
raise NotImplementedError('disorientation between different symmetry classes not supported yet.') raise NotImplementedError('disorientation between different symmetry classes not supported yet.')
mySymEqs = self.equivalentOrientations() if SST else self.equivalentOrientations([0]) # take all or only first sym operation mySymEqs = self.equivalentOrientations() if SST else self.equivalentOrientations([0]) # take all or only first sym operation
otherSymEqs = other.equivalentOrientations() otherSymEqs = other.equivalentOrientations()
for i,sA in enumerate(mySymEqs): for i,sA in enumerate(mySymEqs):
aInv = sA.rotation.inversed() aInv = sA.rotation.inversed()
for j,sB in enumerate(otherSymEqs): for j,sB in enumerate(otherSymEqs):
b = sB.rotation b = sB.rotation
r = b*aInv r = b*aInv
for k in range(2): for k in range(2):
r.inverse() r.inverse()
breaker = self.lattice.symmetry.inFZ(r.asRodrigues(vector=True)) \ breaker = self.lattice.symmetry.inFZ(r.asRodrigues(vector=True)) \
and (not SST or other.lattice.symmetry.inDisorientationSST(r.asRodrigues(vector=True))) and (not SST or other.lattice.symmetry.inDisorientationSST(r.asRodrigues(vector=True)))
if breaker: break if breaker: break
if breaker: break
if breaker: break if breaker: break
if breaker: break
return (Orientation(r,self.lattice), i,j, k == 1) if symmetries else r # disorientation ... return (Orientation(r,self.lattice), i,j, k == 1) if symmetries else r # disorientation ...
# ... own sym, other sym, # ... own sym, other sym,
@ -78,11 +78,11 @@ class Orientation:
def equivalentOrientations(self,members=[]): def equivalentOrientations(self,members=[]):
"""List of orientations which are symmetrically equivalent.""" """List of orientations which are symmetrically equivalent."""
try: try:
iter(members) # asking for (even empty) list of members? iter(members) # asking for (even empty) list of members?
except TypeError: except TypeError:
return self.__class__(self.lattice.symmetry.symmetryOperations(members)*self.rotation,self.lattice) # no, return rotation object return self.__class__(self.lattice.symmetry.symmetryOperations(members)*self.rotation,self.lattice) # no, return rotation object
else: else:
return [self.__class__(q*self.rotation,self.lattice) \ return [self.__class__(q*self.rotation,self.lattice) \
for q in self.lattice.symmetry.symmetryOperations(members)] # yes, return list of rotations for q in self.lattice.symmetry.symmetryOperations(members)] # yes, return list of rotations
def relatedOrientations(self,model): def relatedOrientations(self,model):
@ -94,7 +94,7 @@ class Orientation:
def reduced(self): def reduced(self):
"""Transform orientation to fall into fundamental zone according to symmetry.""" """Transform orientation to fall into fundamental zone according to symmetry."""
for me in self.equivalentOrientations(): for me in self.equivalentOrientations():
if self.lattice.symmetry.inFZ(me.rotation.asRodrigues(vector=True)): break if self.lattice.symmetry.inFZ(me.rotation.asRodrigues(vector=True)): break
return self.__class__(me.rotation,self.lattice) return self.__class__(me.rotation,self.lattice)
@ -105,11 +105,11 @@ class Orientation:
SST = True): SST = True):
"""Axis rotated according to orientation (using crystal symmetry to ensure location falls into SST).""" """Axis rotated according to orientation (using crystal symmetry to ensure location falls into SST)."""
if SST: # pole requested to be within SST if SST: # pole requested to be within SST
for i,o in enumerate(self.equivalentOrientations()): # test all symmetric equivalent quaternions for i,o in enumerate(self.equivalentOrientations()): # test all symmetric equivalent quaternions
pole = o.rotation*axis # align crystal direction to axis pole = o.rotation*axis # align crystal direction to axis
if self.lattice.symmetry.inSST(pole,proper): break # found SST version if self.lattice.symmetry.inSST(pole,proper): break # found SST version
else: else:
pole = self.rotation*axis # align crystal direction to axis pole = self.rotation*axis # align crystal direction to axis
return (pole,i if SST else 0) return (pole,i if SST else 0)
@ -119,9 +119,9 @@ class Orientation:
color = np.zeros(3,'d') color = np.zeros(3,'d')
for o in self.equivalentOrientations(): for o in self.equivalentOrientations():
pole = o.rotation*axis # align crystal direction to axis pole = o.rotation*axis # align crystal direction to axis
inSST,color = self.lattice.symmetry.inSST(pole,color=True) inSST,color = self.lattice.symmetry.inSST(pole,color=True)
if inSST: break if inSST: break
return color return color
@ -131,15 +131,15 @@ class Orientation:
weights = []): weights = []):
"""Create orientation from average of list of orientations.""" """Create orientation from average of list of orientations."""
if not all(isinstance(item, Orientation) for item in orientations): if not all(isinstance(item, Orientation) for item in orientations):
raise TypeError("Only instances of Orientation can be averaged.") raise TypeError("Only instances of Orientation can be averaged.")
closest = [] closest = []
ref = orientations[0] ref = orientations[0]
for o in orientations: for o in orientations:
closest.append(o.equivalentOrientations( closest.append(o.equivalentOrientations(
ref.disorientation(o, ref.disorientation(o,
SST = False, # select (o[ther]'s) sym orientation SST = False, # select (o[ther]'s) sym orientation
symmetries = True)[2]).rotation) # with lowest misorientation symmetries = True)[2]).rotation) # with lowest misorientation
return Orientation(Rotation.fromAverage(closest,weights),ref.lattice) return Orientation(Rotation.fromAverage(closest,weights),ref.lattice)

4
python/damask/table.py
View File

@ -287,7 +287,7 @@ class Table:
Parameters Parameters
---------- ----------
other : Table other : Table
Table to append Table to append.
""" """
if self.shapes != other.shapes or not self.data.columns.equals(other.data.columns): if self.shapes != other.shapes or not self.data.columns.equals(other.data.columns):
@ -305,7 +305,7 @@ class Table:
Parameters Parameters
---------- ----------
other : Table other : Table
Table to join Table to join.
""" """
if set(self.shapes) & set(other.shapes) or self.data.shape[0] != other.data.shape[0]: if set(self.shapes) & set(other.shapes) or self.data.shape[0] != other.data.shape[0]:

48
python/damask/util.py
View File

@ -47,9 +47,9 @@ def srepr(arg,glue = '\n'):
Parameters Parameters
---------- ----------
arg : iterable arg : iterable
Items to join. Items to join.
glue : str, optional glue : str, optional
Defaults to \n. Defaults to \n.
""" """
if (not hasattr(arg, "strip") and if (not hasattr(arg, "strip") and
@ -66,9 +66,9 @@ def croak(what, newline = True):
Parameters Parameters
---------- ----------
what : str or iterable what : str or iterable
Content to be displayed Content to be displayed.
newline : bool, optional newline : bool, optional
Separate items of what by newline. Defaults to True. Separate items of what by newline. Defaults to True.
""" """
if not what: if not what:
@ -117,13 +117,13 @@ def execute(cmd,
Parameters Parameters
---------- ----------
cmd : str cmd : str
Command to be executed. Command to be executed.
streanIn :, optional streanIn :, optional
Input (via pipe) for executed process. Input (via pipe) for executed process.
wd : str, optional wd : str, optional
Working directory of process. Defaults to ./ . Working directory of process. Defaults to ./ .
env : env : dict, optional
Environment Environment for execution.
""" """
initialPath = os.getcwd() initialPath = os.getcwd()
@ -140,7 +140,7 @@ def execute(cmd,
error = error.decode('utf-8').replace('\x08','') error = error.decode('utf-8').replace('\x08','')
os.chdir(initialPath) os.chdir(initialPath)
if process.returncode != 0: if process.returncode != 0:
raise RuntimeError('{} failed with returncode {}'.format(cmd,process.returncode)) raise RuntimeError('{} failed with returncode {}'.format(cmd,process.returncode))
return out,error return out,error
@ -158,11 +158,11 @@ class extendableOption(Option):
ALWAYS_TYPED_ACTIONS = Option.ALWAYS_TYPED_ACTIONS + ("extend",) ALWAYS_TYPED_ACTIONS = Option.ALWAYS_TYPED_ACTIONS + ("extend",)
def take_action(self, action, dest, opt, value, values, parser): def take_action(self, action, dest, opt, value, values, parser):
if action == "extend": if action == "extend":
lvalue = value.split(",") lvalue = value.split(",")
values.ensure_value(dest, []).extend(lvalue) values.ensure_value(dest, []).extend(lvalue)
else: else:
Option.take_action(self, action, dest, opt, value, values, parser) Option.take_action(self, action, dest, opt, value, values, parser)
class _ProgressBar: class _ProgressBar:
@ -179,11 +179,11 @@ class _ProgressBar:
Parameters Parameters
---------- ----------
total : int total : int
Total # of iterations. Total # of iterations.
prefix : str prefix : str
Prefix string. Prefix string.
bar_length : int bar_length : int
Character length of bar. Character length of bar.
""" """
self.total = total self.total = total
@ -224,13 +224,13 @@ def show_progress(iterable,N_iter=None,prefix='',bar_length=50):
Parameters Parameters
---------- ----------
iterable : iterable/function with yield statement iterable : iterable/function with yield statement
Iterable (or function with yield statement) to be decorated. Iterable (or function with yield statement) to be decorated.
N_iter : int N_iter : int
Total # of iterations. Needed if number of iterations can not be obtained as len(iterable). Total # of iterations. Needed if number of iterations can not be obtained as len(iterable).
prefix : str, optional. prefix : str, optional.
Prefix string. Prefix string.
bar_length : int, optional bar_length : int, optional
Character length of bar. Defaults to 50. Character length of bar. Defaults to 50.
""" """
status = _ProgressBar(N_iter if N_iter else len(iterable),prefix,bar_length) status = _ProgressBar(N_iter if N_iter else len(iterable),prefix,bar_length)
@ -238,7 +238,7 @@ def show_progress(iterable,N_iter=None,prefix='',bar_length=50):
for i,item in enumerate(iterable): for i,item in enumerate(iterable):
yield item yield item
status.update(i) status.update(i)
def scale_to_coprime(v): def scale_to_coprime(v):
"""Scale vector to co-prime (relatively prime) integers.""" """Scale vector to co-prime (relatively prime) integers."""
@ -268,7 +268,7 @@ class return_message():
Parameters Parameters
---------- ----------
message : str or list of str message : str or list of str
message for output to screen message for output to screen
""" """
self.message = message self.message = message