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

Merge branch 'development' into YAML-compatible-debug

This commit is contained in:
Martin Diehl 2020-07-03 16:19:38 +02:00
parent 5a423103a0 e1bbaac0d7
commit 76f0c5fc5e
34 changed files with 3398 additions and 1100 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
VERSION
View File

@ -1 +1 @@
v2.0.3-2717-g52aacf37 v2.0.3-2783-g1a8feab2

75
processing/post/perceptualUniformColorMap.py
View File

@ -1,75 +0,0 @@
#!/usr/bin/env python3
import os
import sys
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
#Borland, D., & Taylor, R. M. (2007). Rainbow Color Map (Still) Considered Harmful. Computer Graphics and Applications, IEEE, 27(2), 14--17.
#Moreland, K. (2009). Diverging Color Maps for Scientific Visualization. In Proc. 5th Int. Symp. Visual Computing (pp. 92--103).
outtypes = ['paraview','gmsh','raw','GOM']
extensions = ['.json','.msh','.txt','.legend']
colormodels = ['RGB','HSL','XYZ','CIELAB','MSH']
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
Produces perceptually linear diverging and sequential colormaps in formats suitable for visualization software
or simply as a list of interpolated colors.
""", version = scriptID)
parser.add_option('-N','--steps', dest='N', type='int', nargs=1, metavar='int',
help='number of interpolation steps [%default]')
parser.add_option('-t','--trim', dest='trim', type='float', nargs=2, metavar='float float',
help='relative trim of colormap range [%default]')
parser.add_option('-l','--left', dest='left', type='float', nargs=3, metavar='float float float',
help='left color [%default]')
parser.add_option('-r','--right', dest='right', type='float', nargs=3, metavar='float float float',
help='right color [%default]')
parser.add_option('-c','--colormodel', dest='colormodel', metavar='string',
help='colormodel: '+', '.join(colormodels)+' [%default]')
parser.add_option('-p','--predefined', dest='predefined', metavar='string',
help='predefined colormap')
parser.add_option('-f','--format', dest='format', metavar='string',
help='output format: '+', '.join(outtypes)+' [%default]')
parser.set_defaults(colormodel = 'RGB')
parser.set_defaults(predefined = None)
parser.set_defaults(basename = None)
parser.set_defaults(format = 'paraview')
parser.set_defaults(N = 10)
parser.set_defaults(trim = (-1.0,1.0))
parser.set_defaults(left = (1.0,1.0,1.0))
parser.set_defaults(right = (0.0,0.0,0.0))
(options,filename) = parser.parse_args()
if options.format not in outtypes:
parser.error('invalid format: "{}" (choices: {}).'.format(options.format,', '.join(outtypes)))
if options.N < 2:
parser.error('too few steps (need at least 2).')
if options.trim[0] < -1.0 or \
options.trim[1] > 1.0 or \
options.trim[0] >= options.trim[1]:
parser.error('invalid trim range (-1 +1).')
name = options.format if filename == [] \
else filename[0]
output = sys.stdout if filename == [] \
else open(os.path.basename(filename[0])+extensions[outtypes.index(options.format)],'w')
colorLeft = damask.Color(options.colormodel.upper(), list(options.left))
colorRight = damask.Color(options.colormodel.upper(), list(options.right))
colormap = damask.Colormap(colorLeft, colorRight, predefined=options.predefined)
output.write(colormap.export(name,options.format,options.N,list(options.trim)))
output.close()

2
python/damask/__init__.py
View File

@ -10,7 +10,7 @@ with open(_Path(__file__).parent/_Path('VERSION')) as _f:
from ._environment import Environment # noqa from ._environment import Environment # noqa
from ._table import Table # noqa from ._table import Table # noqa
from ._vtk import VTK # noqa from ._vtk import VTK # noqa
from ._colormaps import Colormap, Color # noqa from ._colormap import Colormap # noqa
from ._rotation import Rotation # noqa from ._rotation import Rotation # noqa
from ._lattice import Symmetry, Lattice# noqa from ._lattice import Symmetry, Lattice# noqa
from ._orientation import Orientation # noqa from ._orientation import Orientation # noqa

612
python/damask/_colormap.py Normal file
View File

@ -0,0 +1,612 @@
import json
import functools
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib import cm
import damask
from . import Table
_eps = 216./24389.
_kappa = 24389./27.
_ref_white = np.array([.95047, 1.00000, 1.08883]) # Observer = 2, Illuminant = D65
# ToDo (if needed)
# - support alpha channel (paraview/ASCII/input)
# - support NaN color (paraview)
class Colormap(mpl.colors.ListedColormap):
def __add__(self,other):
"""Concatenate colormaps."""
return Colormap(np.vstack((self.colors,other.colors)),
f'{self.name}+{other.name}')
def __iadd__(self,other):
"""Concatenate colormaps."""
return self.__add__(other)
def __invert__(self):
"""Return inverted colormap."""
return self.reversed()
@staticmethod
def from_range(low,high,name='DAMASK colormap',N=256,model='rgb'):
"""
Create a perceptually uniform colormap between given (inclusive) bounds.
Colors are internally stored as R(ed) G(green) B(lue) values.
The colormap can be used in matplotlib/seaborn or exported to
file for external use.
Parameters
----------
low : numpy.ndarray of shape (3)
Color definition for minimum value.
high : numpy.ndarray of shape (3)
Color definition for maximum value.
N : integer, optional
The number of color quantization levels. Defaults to 256.
name : str, optional
The name of the colormap. Defaults to `DAMASK colormap`.
model : {'rgb', 'hsv', 'hsl', 'xyz', 'lab', 'msh'}
Colormodel used for input color definitions. Defaults to `rgb`.
The available color models are:
- 'rgb': R(ed) G(green) B(lue).
- 'hsv': H(ue) S(aturation) V(alue).
- 'hsl': H(ue) S(aturation) L(uminance).
- 'xyz': CIE Xyz.
- 'lab': CIE Lab.
- 'msh': Msh (for perceptual uniform interpolation).
"""
low_high = np.vstack((low,high))
if model.lower() == 'rgb':
if np.any(low_high<0) or np.any(low_high>1):
raise ValueError(f'RGB color {low} | {high} are out of range.')
low_,high_ = map(Colormap._rgb2msh,low_high)
elif model.lower() == 'hsv':
if np.any(low_high<0) or np.any(low_high[:,1:3]>1) or np.any(low_high[:,0]>360):
raise ValueError(f'HSV color {low} | {high} are out of range.')
low_,high_ = map(Colormap._hsv2msh,low_high)
elif model.lower() == 'hsl':
if np.any(low_high<0) or np.any(low_high[:,1:3]>1) or np.any(low_high[:,0]>360):
raise ValueError(f'HSL color {low} | {high} are out of range.')
low_,high_ = map(Colormap._hsl2msh,low_high)
elif model.lower() == 'xyz':
low_,high_ = map(Colormap._xyz2msh,low_high)
elif model.lower() == 'lab':
if np.any(low_high[:,0]<0):
raise ValueError(f'CIE Lab color {low} | {high} are out of range.')
low_,high_ = map(Colormap._lab2msh,low_high)
elif model.lower() == 'msh':
low_,high_ = low_high[0],low_high[1]
else:
raise ValueError(f'Invalid color model: {model}.')
msh = map(functools.partial(Colormap._interpolate_msh,low=low_,high=high_),np.linspace(0,1,N))
rgb = np.array(list(map(Colormap._msh2rgb,msh)))
return Colormap(rgb,name=name)
@staticmethod
def from_predefined(name,N=256):
"""
Select from a set of predefined colormaps.
Predefined colormaps include native matplotlib colormaps
and common DAMASK colormaps.
Parameters
----------
name : str
The name of the colormap.
N : int, optional
The number of color quantization levels. Defaults to 256.
This parameter is not used for matplotlib colormaps
that are of type `ListedColormap`.
"""
# matplotlib presets
for cat in Colormap._predefined_mpl:
for n in cat[1]:
if n == name:
colormap = cm.__dict__[name]
if isinstance(colormap,mpl.colors.LinearSegmentedColormap):
return Colormap(np.array(list(map(colormap,np.linspace(0,1,N)))),name=name)
else:
return Colormap(np.array(colormap.colors),name=name)
# DAMASK presets
definition = Colormap._predefined_DAMASK[name]
return Colormap.from_range(definition['low'],definition['high'],name,N)
@staticmethod
def list_predefined():
"""
List predefined colormaps by category.
References
----------
.. [1] DAMASK colormap theory
https://www.kennethmoreland.com/color-maps/ColorMapsExpanded.pdf
.. [2] DAMASK colormaps first use
https://doi.org/10.1016/j.ijplas.2012.09.012
.. [3] Matplotlib colormaps overview
https://matplotlib.org/tutorials/colors/colormaps.html
"""
print('DAMASK colormaps')
print(' '+', '.join(Colormap._predefined_DAMASK.keys()))
for cat in Colormap._predefined_mpl:
print(f'{cat[0]}')
print(' '+', '.join(cat[1]))
def show(self,aspect=10,vertical=False):
"""Show colormap as matplotlib figure."""
fig = plt.figure(figsize=(5/aspect,5) if vertical else (5,5/aspect))
ax1 = fig.add_axes([0, 0, 1, 1])
ax1.set_axis_off()
ax1.imshow(np.linspace(1 if vertical else 0,
0 if vertical else 1,
self.N).reshape((-1,1) if vertical else (1,-1)),
aspect='auto', cmap=self, interpolation='nearest')
plt.show()
def reversed(self,name=None):
"""
Make a reversed instance of the colormap.
Parameters
----------
name : str, optional
The name for the reversed colormap.
A name of None will be replaced by the name of the parent colormap + "_r".
Returns
-------
damask.Colormap
The reversed colormap.
"""
rev = super(Colormap,self).reversed(name)
return Colormap(np.array(rev.colors),rev.name[:-4] if rev.name.endswith('_r_r') else rev.name)
def to_file(self,fname=None,format='ParaView'):
"""
Export colormap to file for use in external programs.
Parameters
----------
fname : file, str, or pathlib.Path, optional.
Filename to store results. If not given, the filename will
consist of the name of the colormap and an extension that
depends on the file format.
format : {'ParaView', 'ASCII', 'GOM', 'gmsh'}, optional
File format, defaults to 'ParaView'. Available formats are:
- ParaView: JSON file, extension '.json'.
- ASCII: Plain text file, extension '.txt'.
- GOM: Aramis GOM (DIC), extension '.legend'.
- Gmsh: Gmsh FEM mesh-generator, extension '.msh'.
"""
if fname is not None:
try:
f = open(fname,'w')
except TypeError:
f = fname
else:
f = None
if format.lower() == 'paraview':
Colormap._export_paraview(self,f)
elif format.lower() == 'ascii':
Colormap._export_ASCII(self,f)
elif format.lower() == 'gom':
Colormap._export_GOM(self,f)
elif format.lower() == 'gmsh':
Colormap._export_gmsh(self,f)
else:
raise ValueError('Unknown output format: {format}.')
@staticmethod
def _export_paraview(colormap,fhandle=None):
"""Write colormap to JSON file for Paraview."""
colors = []
for i,c in enumerate(np.round(colormap.colors,6).tolist()):
colors+=[i]+c
out = [{
'Creator':f'damask.Colormap v{damask.version}',
'ColorSpace':'RGB',
'Name':colormap.name,
'DefaultMap':True,
'RGBPoints':colors
}]
if fhandle is None:
with open(colormap.name.replace(' ','_')+'.json', 'w') as f:
json.dump(out, f,indent=4)
else:
json.dump(out,fhandle,indent=4)
@staticmethod
def _export_ASCII(colormap,fhandle=None):
"""Write colormap to ASCII table."""
labels = {'RGBA':4} if colormap.colors.shape[1] == 4 else {'RGB': 3}
t = Table(colormap.colors,labels,f'Creator: damask.Colormap v{damask.version}')
if fhandle is None:
with open(colormap.name.replace(' ','_')+'.txt', 'w') as f:
t.to_ASCII(f,True)
else:
t.to_ASCII(fhandle,True)
@staticmethod
def _export_GOM(colormap,fhandle=None):
"""Write colormap to GOM Aramis compatible format."""
# ToDo: test in GOM
GOM_str = f'1 1 {colormap.name.replace(" ","_")} 9 {colormap.name.replace(" ","_")} ' \
+ '0 1 0 3 0 0 -1 9 \\ 0 0 0 255 255 255 0 0 255 ' \
+ f'30 NO_UNIT 1 1 64 64 64 255 1 0 0 0 0 0 0 3 0 {len(colormap.colors)}' \
+ ' '.join([f' 0 {c[0]} {c[1]} {c[2]} 255 1' for c in reversed((colormap.colors*255).astype(int))]) \
+ '\n'
if fhandle is None:
with open(colormap.name.replace(' ','_')+'.legend', 'w') as f:
f.write(GOM_str)
else:
fhandle.write(GOM_str)
@staticmethod
def _export_gmsh(colormap,fhandle=None):
"""Write colormap to Gmsh compatible format."""
# ToDo: test in gmsh
gmsh_str = 'View.ColorTable = {\n' \
+'\n'.join([f'{c[0]},{c[1]},{c[2]},' for c in colormap.colors[:,:3]*255]) \
+'\n}\n'
if fhandle is None:
with open(colormap.name.replace(' ','_')+'.msh', 'w') as f:
f.write(gmsh_str)
else:
fhandle.write(gmsh_str)
@staticmethod
def _interpolate_msh(frac,low,high):
"""
Interpolate in Msh color space.
This interpolation gives a perceptually uniform colormap.
References
----------
https://www.kennethmoreland.com/color-maps/ColorMapsExpanded.pdf
https://www.kennethmoreland.com/color-maps/diverging_map.py
"""
def rad_diff(a,b):
return abs(a[2]-b[2])
def adjust_hue(msh_sat, msh_unsat):
"""If saturation of one of the two colors is much less than the other, hue of the less."""
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
lo = np.array(low)
hi = np.array(high)
if (lo[1] > 0.05 and hi[1] > 0.05 and rad_diff(lo,hi) > np.pi/3.0):
M_mid = max(lo[0],hi[0],88.0)
if frac < 0.5:
hi = np.array([M_mid,0.0,0.0])
frac *= 2.0
else:
lo = np.array([M_mid,0.0,0.0])
frac = 2.0*frac - 1.0
if lo[1] < 0.05 and hi[1] > 0.05:
lo[2] = adjust_hue(hi,lo)
elif lo[1] > 0.05 and hi[1] < 0.05:
hi[2] = adjust_hue(lo,hi)
return (1.0 - frac) * lo + frac * hi
_predefined_mpl= [('Perceptually Uniform Sequential', [
'viridis', 'plasma', 'inferno', 'magma', 'cividis']),
('Sequential', [
'Greys', 'Purples', 'Blues', 'Greens', 'Oranges', 'Reds',
'YlOrBr', 'YlOrRd', 'OrRd', 'PuRd', 'RdPu', 'BuPu',
'GnBu', 'PuBu', 'YlGnBu', 'PuBuGn', 'BuGn', 'YlGn']),
('Sequential (2)', [
'binary', 'gist_yarg', 'gist_gray', 'gray', 'bone', 'pink',
'spring', 'summer', 'autumn', 'winter', 'cool', 'Wistia',
'hot', 'afmhot', 'gist_heat', 'copper']),
('Diverging', [
'PiYG', 'PRGn', 'BrBG', 'PuOr', 'RdGy', 'RdBu',
'RdYlBu', 'RdYlGn', 'Spectral', 'coolwarm', 'bwr', 'seismic']),
('Cyclic', ['twilight', 'twilight_shifted', 'hsv']),
('Qualitative', [
'Pastel1', 'Pastel2', 'Paired', 'Accent',
'Dark2', 'Set1', 'Set2', 'Set3',
'tab10', 'tab20', 'tab20b', 'tab20c']),
('Miscellaneous', [
'flag', 'prism', 'ocean', 'gist_earth', 'terrain', 'gist_stern',
'gnuplot', 'gnuplot2', 'CMRmap', 'cubehelix', 'brg',
'gist_rainbow', 'rainbow', 'jet', 'nipy_spectral', 'gist_ncar'])]
_predefined_DAMASK = {'orientation': {'low': [0.933334,0.878432,0.878431],
'high': [0.250980,0.007843,0.000000]},
'strain': {'low': [0.941177,0.941177,0.870588],
'high': [0.266667,0.266667,0.000000]},
'stress': {'low': [0.878432,0.874511,0.949019],
'high': [0.000002,0.000000,0.286275]}}
@staticmethod
def _hsv2rgb(hsv):
"""
H(ue) S(aturation) V(alue) to R(red) G(reen) B(lue).
References
----------
https://www.rapidtables.com/convert/color/hsv-to-rgb.html
"""
sextant = np.clip(int(hsv[0]/60.),0,5)
c = hsv[1]*hsv[2]
x = c*(1.0 - abs((hsv[0]/60.)%2 - 1.))
return np.array([
[c, x, 0],
[x, c, 0],
[0, c, x],
[0, x, c],
[x, 0, c],
[c, 0, x],
])[sextant] + hsv[2] - c
@staticmethod
def _rgb2hsv(rgb):
"""
R(ed) G(reen) B(lue) to H(ue) S(aturation) V(alue).
References
----------
https://www.rapidtables.com/convert/color/rgb-to-hsv.html
"""
C_max = rgb.max()
C_min = rgb.min()
Delta = C_max - C_min
v = C_max
s = 0. if np.isclose(C_max,0.) else Delta/C_max
if np.isclose(Delta,0.):
h = 0.
elif rgb.argmax() == 0:
h = (rgb[1]-rgb[2])/Delta%6
elif rgb.argmax() == 1:
h = (rgb[2]-rgb[0])/Delta + 2.
elif rgb.argmax() == 2:
h = (rgb[0]-rgb[1])/Delta + 4.
h = np.clip(h,0.,6.) * 60.
return np.array([h,s,v])
@staticmethod
def _hsl2rgb(hsl):
"""
H(ue) S(aturation) L(uminance) to R(red) G(reen) B(lue).
References
----------
https://www.rapidtables.com/convert/color/hsl-to-rgb.html
"""
sextant = np.clip(int(hsl[0]/60.),0,5)
c = (1.0 - abs(2.0 * hsl[2] - 1.))*hsl[1]
x = c*(1.0 - abs((hsl[0]/60.)%2 - 1.))
m = hsl[2] - 0.5*c
return np.array([
[c+m, x+m, m],
[x+m, c+m, m],
[m, c+m, x+m],
[m, x+m, c+m],
[x+m, m, c+m],
[c+m, m, x+m],
])[sextant]
@staticmethod
def _rgb2hsl(rgb):
"""
R(ed) G(reen) B(lue) to H(ue) S(aturation) L(uminance).
References
----------
https://www.rapidtables.com/convert/color/rgb-to-hsl.html
"""
C_max = rgb.max()
C_min = rgb.min()
Delta = C_max - C_min
l = np.clip((C_max + C_min)*.5,0.,1.) # noqa
s = 0. if np.isclose(C_max,C_min) else Delta/(1.-np.abs(2*l-1.))
if np.isclose(Delta,0.):
h = 0.
elif rgb.argmax() == 0:
h = (rgb[1]-rgb[2])/Delta%6
elif rgb.argmax() == 1:
h = (rgb[2]-rgb[0])/Delta + 2.
elif rgb.argmax() == 2:
h = (rgb[0]-rgb[1])/Delta + 4.
h = np.clip(h,0.,6.) * 60.
return np.array([h,s,l])
@staticmethod
def _xyz2rgb(xyz):
"""
CIE Xyz to R(ed) G(reen) B(lue).
References
----------
http://www.ryanjuckett.com/programming/rgb-color-space-conversion
"""
rgb_lin = np.dot(np.array([
[ 3.240969942,-1.537383178,-0.498610760],
[-0.969243636, 1.875967502, 0.041555057],
[ 0.055630080,-0.203976959, 1.056971514]
]),xyz)
with np.errstate(invalid='ignore'):
rgb = np.where(rgb_lin>0.0031308,rgb_lin**(1.0/2.4)*1.0555-0.0555,rgb_lin*12.92)
return np.clip(rgb,0.,1.)
@staticmethod
def _rgb2xyz(rgb):
"""
R(ed) G(reen) B(lue) to CIE Xyz.
References
----------
http://www.ryanjuckett.com/programming/rgb-color-space-conversion
"""
rgb_lin = np.where(rgb>0.04045,((rgb+0.0555)/1.0555)**2.4,rgb/12.92)
return np.dot(np.array([
[0.412390799,0.357584339,0.180480788],
[0.212639006,0.715168679,0.072192315],
[0.019330819,0.119194780,0.950532152]
]),rgb_lin)
@staticmethod
def _lab2xyz(lab,ref_white=None):
"""
CIE Lab to CIE Xyz.
References
----------
http://www.brucelindbloom.com/index.html?Eqn_Lab_to_XYZ.html
"""
f_x = (lab[0]+16.)/116. + lab[1]/500.
f_z = (lab[0]+16.)/116. - lab[2]/200.
return np.array([
f_x**3. if f_x**3. > _eps else (116.*f_x-16.)/_kappa,
((lab[0]+16.)/116.)**3 if lab[0]>_kappa*_eps else lab[0]/_kappa,
f_z**3. if f_z**3. > _eps else (116.*f_z-16.)/_kappa
])*(ref_white if ref_white is not None else _ref_white)
@staticmethod
def _xyz2lab(xyz,ref_white=None):
"""
CIE Xyz to CIE Lab.
References
----------
http://www.brucelindbloom.com/index.html?Eqn_Lab_to_XYZ.html
"""
ref_white = ref_white if ref_white is not None else _ref_white
f = np.where(xyz/ref_white > _eps,(xyz/ref_white)**(1./3.),(_kappa*xyz/ref_white+16.)/116.)
return np.array([
116.0 * f[1] - 16.0,
500.0 * (f[0] - f[1]),
200.0 * (f[1] - f[2])
])
@staticmethod
def _lab2msh(lab):
"""
CIE Lab to Msh.
References
----------
https://www.kennethmoreland.com/color-maps/ColorMapsExpanded.pdf
https://www.kennethmoreland.com/color-maps/diverging_map.py
"""
M = np.linalg.norm(lab)
return np.array([
M,
np.arccos(lab[0]/M) if M>1e-8 else 0.,
np.arctan2(lab[2],lab[1]) if M>1e-8 else 0.,
])
@staticmethod
def _msh2lab(msh):
"""
Msh to CIE Lab.
References
----------
https://www.kennethmoreland.com/color-maps/ColorMapsExpanded.pdf
https://www.kennethmoreland.com/color-maps/diverging_map.py
"""
return np.array([
msh[0] * np.cos(msh[1]),
msh[0] * np.sin(msh[1]) * np.cos(msh[2]),
msh[0] * np.sin(msh[1]) * np.sin(msh[2])
])
@staticmethod
def _lab2rgb(lab):
return Colormap._xyz2rgb(Colormap._lab2xyz(lab))
@staticmethod
def _rgb2lab(rgb):
return Colormap._xyz2lab(Colormap._rgb2xyz(rgb))
@staticmethod
def _msh2rgb(msh):
return Colormap._lab2rgb(Colormap._msh2lab(msh))
@staticmethod
def _rgb2msh(rgb):
return Colormap._lab2msh(Colormap._rgb2lab(rgb))
@staticmethod
def _hsv2msh(hsv):
return Colormap._rgb2msh(Colormap._hsv2rgb(hsv))
@staticmethod
def _hsl2msh(hsl):
return Colormap._rgb2msh(Colormap._hsl2rgb(hsl))
@staticmethod
def _xyz2msh(xyz):
return Colormap._lab2msh(Colormap._xyz2lab(xyz))

541
python/damask/_colormaps.py
View File

@ -1,541 +0,0 @@
import numpy as np
from . import util
class Color:
"""Color representation in and conversion between different color-spaces."""
__slots__ = [
'model',
'color',
'__dict__',
]
def __init__(self,
model = 'RGB',
color = np.zeros(3)):
"""
Create a Color object.
Parameters
----------
model : string
color model
color : numpy.ndarray
vector representing the color according to the selected model
"""
self.__transforms__ = \
{'HSV': {'index': 0, 'next': self._HSV2HSL},
'HSL': {'index': 1, 'next': self._HSL2RGB, 'prev': self._HSL2HSV},
'RGB': {'index': 2, 'next': self._RGB2XYZ, 'prev': self._RGB2HSL},
'XYZ': {'index': 3, 'next': self._XYZ2CIELAB, 'prev': self._XYZ2RGB},
'CIELAB': {'index': 4, 'next': self._CIELAB2MSH, 'prev': self._CIELAB2XYZ},
'MSH': {'index': 5, 'prev': self._MSH2CIELAB},
}
model = model.upper()
if model not in list(self.__transforms__.keys()): model = 'RGB'
if model == 'RGB' and max(color) > 1.0: # are we RGB255 ?
for i in range(3):
color[i] /= 255.0 # rescale to RGB
if model == 'HSL': # are we HSL ?
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] < 0.0: color[0] += 1.0 # rewind to proper range
self.model = model
self.color = np.array(color,'d')
def __repr__(self):
"""Color model and values."""
return 'Model: %s Color: %s'%(self.model,str(self.color))
def __str__(self):
"""Color model and values."""
return self.__repr__()
def convert_to(self,toModel = 'RGB'):
"""
Change the color model permanently.
Parameters
----------
toModel : string
color model
"""
toModel = toModel.upper()
if toModel not in list(self.__transforms__.keys()): return
sourcePos = self.__transforms__[self.model]['index']
targetPos = self.__transforms__[toModel]['index']
while sourcePos < targetPos:
self.__transforms__[self.model]['next']()
sourcePos += 1
while sourcePos > targetPos:
self.__transforms__[self.model]['prev']()
sourcePos -= 1
return self
def express_as(self,asModel = 'RGB'):
"""
Return the color in a different model.
Parameters
----------
asModel : string
color model
"""
return self.__class__(self.model,self.color).convert_to(asModel)
def _HSV2HSL(self):
"""
Convert H(ue) S(aturation) V(alue or brightness) to H(ue) S(aturation) L(uminance).
All values are in the range [0,1]
http://codeitdown.com/hsl-hsb-hsv-color
"""
if self.model != 'HSV': return
converted = Color('HSL',np.array([
self.color[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.)),
0.5*self.color[2]*(2.-self.color[1]),
]))
self.model = converted.model
self.color = converted.color
def _HSL2HSV(self):
"""
Convert H(ue) S(aturation) L(uminance) to H(ue) S(aturation) V(alue or brightness).
All values are in the range [0,1]
http://codeitdown.com/hsl-hsb-hsv-color
"""
if self.model != 'HSL': return
h = self.color[0]
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
converted = Color('HSV',np.array([h,s,b]))
self.model = converted.model
self.color = converted.color
def _HSL2RGB(self):
"""
Convert H(ue) S(aturation) L(uminance) to R(red) G(reen) B(lue).
All values are in the range [0,1]
from http://en.wikipedia.org/wiki/HSL_and_HSV
"""
if self.model != 'HSL': return
sextant = self.color[0]*6.0
c = (1.0 - abs(2.0 * self.color[2] - 1.0))*self.color[1]
x = c*(1.0 - abs(sextant%2 - 1.0))
m = self.color[2] - 0.5*c
converted = Color('RGB',np.array([
[c+m, x+m, m],
[x+m, c+m, m],
[m, c+m, x+m],
[m, x+m, c+m],
[x+m, m, c+m],
[c+m, m, x+m],
][int(sextant)]))
self.model = converted.model
self.color = converted.color
def _RGB2HSL(self):
"""
Convert R(ed) G(reen) B(lue) to H(ue) S(aturation) L(uminance).
All values are in the range [0,1]
from http://130.113.54.154/~monger/hsl-rgb.html
"""
if self.model != 'RGB': return
HSL = np.zeros(3)
maxcolor = self.color.max()
mincolor = self.color.min()
HSL[2] = (maxcolor + mincolor)/2.0
if(mincolor == maxcolor):
HSL[0] = 0.0
HSL[1] = 0.0
else:
if (HSL[2]<0.5):
HSL[1] = (maxcolor - mincolor)/(maxcolor + mincolor)
else:
HSL[1] = (maxcolor - mincolor)/(2.0 - maxcolor - mincolor)
if (maxcolor == self.color[0]):
HSL[0] = 0.0 + (self.color[1] - self.color[2])/(maxcolor - mincolor)
elif (maxcolor == self.color[1]):
HSL[0] = 2.0 + (self.color[2] - self.color[0])/(maxcolor - mincolor)
elif (maxcolor == self.color[2]):
HSL[0] = 4.0 + (self.color[0] - self.color[1])/(maxcolor - mincolor)
HSL[0] = HSL[0]*60.0 # scaling to 360 might be dangerous for small values
if (HSL[0] < 0.0):
HSL[0] = HSL[0] + 360.0
HSL[1:] = np.clip(HSL[1:],0.0,1.0)
converted = Color('HSL', HSL)
self.model = converted.model
self.color = converted.color
def _RGB2XYZ(self):
"""
Convert R(ed) G(reen) B(lue) to CIE XYZ.
All values are in the range [0,1]
from http://www.cs.rit.edu/~ncs/color/t_convert.html
"""
if self.model != 'RGB': return
XYZ = np.zeros(3)
RGB_lin = np.zeros(3)
convert = np.array([[0.412453,0.357580,0.180423],
[0.212671,0.715160,0.072169],
[0.019334,0.119193,0.950227]])
for i in range(3):
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
XYZ = np.dot(convert,RGB_lin)
XYZ = np.clip(XYZ,0.0,None)
converted = Color('XYZ', XYZ)
self.model = converted.model
self.color = converted.color
def _XYZ2RGB(self):
"""
Convert CIE XYZ to R(ed) G(reen) B(lue).
All values are in the range [0,1]
from http://www.cs.rit.edu/~ncs/color/t_convert.html
"""
if self.model != 'XYZ': return
convert = np.array([[ 3.240479,-1.537150,-0.498535],
[-0.969256, 1.875992, 0.041556],
[ 0.055648,-0.204043, 1.057311]])
RGB_lin = np.dot(convert,self.color)
RGB = np.zeros(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
else: RGB[i] = RGB_lin[i] *12.92
RGB = np.clip(RGB,0.0,1.0)
maxVal = max(RGB) # clipping colors according to the display gamut
if (maxVal > 1.0): RGB /= maxVal
converted = Color('RGB', RGB)
self.model = converted.model
self.color = converted.color
def _CIELAB2XYZ(self):
"""
Convert CIE Lab to CIE XYZ.
All values are in the range [0,1]
from http://www.easyrgb.com/index.php?X=MATH&H=07#text7
"""
if self.model != 'CIELAB': return
ref_white = np.array([.95047, 1.00000, 1.08883]) # Observer = 2, Illuminant = D65
XYZ = np.zeros(3)
XYZ[1] = (self.color[0] + 16.0 ) / 116.0
XYZ[0] = XYZ[1] + self.color[1]/ 500.0
XYZ[2] = XYZ[1] - self.color[2]/ 200.0
for i in range(len(XYZ)):
if (XYZ[i] > 6./29. ): XYZ[i] = XYZ[i]**3.
else: XYZ[i] = 108./841. * (XYZ[i] - 4./29.)
converted = Color('XYZ', XYZ*ref_white)
self.model = converted.model
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]) # 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):
"""
Convert CIE Lab to Msh colorspace.
from http://www.cs.unm.edu/~kmorel/documents/ColorMaps/DivergingColorMapWorkshop.xls
"""
if self.model != 'CIELAB': return
Msh = np.zeros(3)
Msh[0] = np.sqrt(np.dot(self.color,self.color))
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)
self.model = converted.model
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)
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
):
"""
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):
"""Left and right value of colormap."""
return 'Left: %s Right: %s'%(self.left,self.right)
def invert(self):
"""Switch left/minimum with right/maximum."""
(self.left, self.right) = (self.right, self.left)
return self
def show_predefined(self):
"""Show the labels of the predefined colormaps."""
print('\n'.join(self.__predefined__.keys()))
def color(self,fraction = 0.5):
def interpolate_Msh(lo, hi, frac):
def rad_diff(a,b):
return abs(a[2]-b[2])
def adjust_hue(Msh_sat, Msh_unsat):
"""If saturation of one of the two colors is too less than the other, hue of the less."""
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[:])
Msh2 = np.array(hi[:])
if (Msh1[1] > 0.05 and Msh2[1] > 0.05 and rad_diff(Msh1,Msh2) > np.pi/3.0):
M_mid = max(Msh1[0],Msh2[0],88.0)
if frac < 0.5:
Msh2 = np.array([M_mid,0.0,0.0])
frac *= 2.0
else:
Msh1 = np.array([M_mid,0.0,0.0])
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
return Color('MSH',Msh)
def interpolate_linear(lo, hi, frac):
"""Linear interpolation between lo and hi color at given fraction; output in model of lo color."""
interpolation = (1.0 - frac) * np.array(lo.color[:]) \
+ frac * np.array(hi.express_as(lo.model).color[:])
return Color(lo.model,interpolation)
if self.interpolate == 'perceptualuniform':
return interpolate_Msh(self.left.express_as('MSH').color,
self.right.express_as('MSH').color,fraction)
elif self.interpolate == 'linear':
return interpolate_linear(self.left,self.right,fraction)
else:
raise NameError('unknown color interpolation method')
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.
Arguments: name, format, steps, crop.
Format is one of (paraview, gmsh, gom, raw, list).
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 = [f'[\n {{\n "ColorSpace": "RGB", "Name": "{name}", "DefaultMap": true,\n "RGBPoints" : ['] \
+ [f' {i:4d},{color[0]:8.6f},{color[1]:8.6f},{color[2]:8.6f}{"," if i+1<len(colors) else ""}' \
for i,color in enumerate(colors)] \
+ [' ]\n }\n]']
elif format == 'gmsh':
colormap = ['View.ColorTable = {'] \
+ [',\n'.join([','.join([str(x*255.0) for x in color]) for color in colors])] \
+ ['}']
elif format == 'gom':
colormap = [ f'1 1 {name}'
+ f' 9 {name}'
+ ' 0 1 0 3 0 0 -1 9 \\ 0 0 0 255 255 255 0 0 255 '
+ f'30 NO_UNIT 1 1 64 64 64 255 1 0 0 0 0 0 0 3 0 {len(colors)}'
+ ' '.join([f' 0 {util.srepr((255*np.array(c)).astype(int)," ")} 255 1' for c 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

32
python/damask/_environment.py
View File

@ -1,19 +1,32 @@
import os import os
from pathlib import Path from pathlib import Path
import tkinter
class Environment: class Environment:
# ToDo: Probably, we don't need a class (just a module with a few functions)
def __init__(self): def __init__(self):
"""Read and provide values of DAMASK configuration.""" """Do Nothing."""
pass
@property
def screen_size(self):
width = 1024
height = 768
try: try:
tk = tkinter.Tk() import wx
self.screen_width = tk.winfo_screenwidth() _ = wx.App(False) # noqa
self.screen_height = tk.winfo_screenheight() width, height = wx.GetDisplaySize()
tk.destroy() except ImportError:
except tkinter.TclError: try:
self.screen_width = 1024 import tkinter
self.screen_height = 768 tk = tkinter.Tk()
width = tk.winfo_screenwidth()
height = tk.winfo_screenheight()
tk.destroy()
except Exception as e:
pass
return (width,height)
@property @property
def options(self): def options(self):
@ -26,6 +39,7 @@ class Environment:
return options return options
@property @property
def root_dir(self): def root_dir(self):
"""Return DAMASK root path.""" """Return DAMASK root path."""

261
python/damask/_lattice.py
View File

@ -1,11 +1,9 @@
import numpy as np import numpy as np
from . import Rotation
class Symmetry: class Symmetry:
""" """
Symmetry operations for lattice systems. Symmetry-related operations for crystal systems.
References References
---------- ----------
@ -13,34 +11,34 @@ class Symmetry:
""" """
lattices = [None,'orthorhombic','tetragonal','hexagonal','cubic',] crystal_systems = [None,'orthorhombic','tetragonal','hexagonal','cubic']
def __init__(self, symmetry = None): def __init__(self, system = None):
""" """
Symmetry Definition. Symmetry Definition.
Parameters Parameters
---------- ----------
symmetry : str, optional system : {None,'orthorhombic','tetragonal','hexagonal','cubic'}, optional
label of the crystal system Name of the crystal system. Defaults to 'None'.
""" """
if symmetry is not None and symmetry.lower() not in Symmetry.lattices: if system is not None and system.lower() not in self.crystal_systems:
raise KeyError(f'Symmetry/crystal system "{symmetry}" is unknown') raise KeyError(f'Crystal system "{system}" is unknown')
self.lattice = symmetry.lower() if isinstance(symmetry,str) else symmetry self.system = system.lower() if isinstance(system,str) else system
def __copy__(self): def __copy__(self):
"""Copy.""" """Copy."""
return self.__class__(self.lattice) return self.__class__(self.system)
copy = __copy__ copy = __copy__
def __repr__(self): def __repr__(self):
"""Readable string.""" """Readable string."""
return f'{self.lattice}' return f'{self.system}'
def __eq__(self, other): def __eq__(self, other):
@ -53,7 +51,7 @@ class Symmetry:
Symmetry to check for equality. Symmetry to check for equality.
""" """
return self.lattice == other.lattice return self.system == other.system
def __neq__(self, other): def __neq__(self, other):
""" """
@ -77,14 +75,16 @@ class Symmetry:
Symmetry to check for for order. Symmetry to check for for order.
""" """
myOrder = Symmetry.lattices.index(self.lattice) myOrder = self.crystal_systems.index(self.system)
otherOrder = Symmetry.lattices.index(other.lattice) otherOrder = self.crystal_systems.index(other.system)
return (myOrder > otherOrder) - (myOrder < otherOrder) return (myOrder > otherOrder) - (myOrder < otherOrder)
def symmetryOperations(self,members=[]):
"""List (or single element) of symmetry operations as rotations.""" @property
if self.lattice == 'cubic': def symmetry_operations(self):
symQuats = [ """Symmetry operations as quaternions."""
if self.system == 'cubic':
sym_quats = [
[ 1.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, 1.0, 0.0, 0.0 ],
[ 0.0, 0.0, 1.0, 0.0 ], [ 0.0, 0.0, 1.0, 0.0 ],
@ -110,8 +110,8 @@ class Symmetry:
[-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 ],
[-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': elif self.system == 'hexagonal':
symQuats = [ sym_quats = [
[ 1.0, 0.0, 0.0, 0.0 ], [ 1.0, 0.0, 0.0, 0.0 ],
[-0.5*np.sqrt(3), 0.0, 0.0, -0.5 ], [-0.5*np.sqrt(3), 0.0, 0.0, -0.5 ],
[ 0.5, 0.0, 0.0, 0.5*np.sqrt(3) ], [ 0.5, 0.0, 0.0, 0.5*np.sqrt(3) ],
@ -125,8 +125,8 @@ class Symmetry:
[ 0.0, -0.5, -0.5*np.sqrt(3), 0.0 ], [ 0.0, -0.5, -0.5*np.sqrt(3), 0.0 ],
[ 0.0, 0.5*np.sqrt(3), 0.5, 0.0 ], [ 0.0, 0.5*np.sqrt(3), 0.5, 0.0 ],
] ]
elif self.lattice == 'tetragonal': elif self.system == 'tetragonal':
symQuats = [ sym_quats = [
[ 1.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, 1.0, 0.0, 0.0 ],
[ 0.0, 0.0, 1.0, 0.0 ], [ 0.0, 0.0, 1.0, 0.0 ],
@ -136,64 +136,54 @@ class Symmetry:
[ 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.0, 0.5*np.sqrt(2) ], [-0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
] ]
elif self.lattice == 'orthorhombic': elif self.system == 'orthorhombic':
symQuats = [ sym_quats = [
[ 1.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,1.0,0.0,0.0 ],
[ 0.0,0.0,1.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,1.0 ],
] ]
else: else:
symQuats = [ sym_quats = [
[ 1.0,0.0,0.0,0.0 ], [ 1.0,0.0,0.0,0.0 ],
] ]
return np.array(sym_quats)
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): def in_FZ(self,rho):
""" """
Check whether given Rodrigues-Frank vector falls into fundamental zone of own symmetry. Check whether given Rodrigues-Frank vector falls into fundamental zone.
Fundamental zone in Rodrigues space is point symmetric around origin. Fundamental zone in Rodrigues space is point symmetric around origin.
""" """
if (len(rodrigues) != 3): if(rho.shape[-1] != 3):
raise ValueError('Input is not a Rodrigues-Frank vector.\n') raise ValueError('Input is not a Rodrigues-Frank vector field.')
if np.any(rodrigues == np.inf): return False rho_abs = np.abs(rho)
Rabs = abs(rodrigues) with np.errstate(invalid='ignore'):
# using '*'/prod for 'and'
if self.lattice == 'cubic': if self.system == 'cubic':
return np.sqrt(2.0)-1.0 >= Rabs[0] \ return np.where(np.prod(np.sqrt(2)-1. >= rho_abs,axis=-1) * \
and np.sqrt(2.0)-1.0 >= Rabs[1] \ (1. >= np.sum(rho_abs,axis=-1)),True,False)
and np.sqrt(2.0)-1.0 >= Rabs[2] \ elif self.system == 'hexagonal':
and 1.0 >= Rabs[0] + Rabs[1] + Rabs[2] return np.where(np.prod(1. >= rho_abs,axis=-1) * \
elif self.lattice == 'hexagonal': (2. >= np.sqrt(3)*rho_abs[...,0] + rho_abs[...,1]) * \
return 1.0 >= Rabs[0] and 1.0 >= Rabs[1] and 1.0 >= Rabs[2] \ (2. >= np.sqrt(3)*rho_abs[...,1] + rho_abs[...,0]) * \
and 2.0 >= np.sqrt(3)*Rabs[0] + Rabs[1] \ (2. >= np.sqrt(3) + rho_abs[...,2]),True,False)
and 2.0 >= np.sqrt(3)*Rabs[1] + Rabs[0] \ elif self.system == 'tetragonal':
and 2.0 >= np.sqrt(3) + Rabs[2] return np.where(np.prod(1. >= rho_abs[...,:2],axis=-1) * \
elif self.lattice == 'tetragonal': (np.sqrt(2) >= rho_abs[...,0] + rho_abs[...,1]) * \
return 1.0 >= Rabs[0] and 1.0 >= Rabs[1] \ (np.sqrt(2) >= rho_abs[...,2] + 1.),True,False)
and np.sqrt(2.0) >= Rabs[0] + Rabs[1] \ elif self.system == 'orthorhombic':
and np.sqrt(2.0) >= Rabs[2] + 1.0 return np.where(np.prod(1. >= rho_abs,axis=-1),True,False)
elif self.lattice == 'orthorhombic': else:
return 1.0 >= Rabs[0] and 1.0 >= Rabs[1] and 1.0 >= Rabs[2] return np.where(np.all(np.isfinite(rho_abs),axis=-1),True,False)
else:
return True
def inDisorientationSST(self,rodrigues): def in_disorientation_SST(self,rho):
""" """
Check whether given Rodrigues-Frank vector (of misorientation) falls into standard stereographic triangle of own symmetry. Check whether given Rodrigues-Frank vector (of misorientation) falls into standard stereographic triangle.
References References
---------- ----------
@ -201,27 +191,33 @@ class Symmetry:
https://doi.org/10.1107/S0108767391006864 https://doi.org/10.1107/S0108767391006864
""" """
if (len(rodrigues) != 3): if(rho.shape[-1] != 3):
raise ValueError('Input is not a Rodrigues-Frank vector.\n') raise ValueError('Input is not a Rodrigues-Frank vector field.')
R = rodrigues
epsilon = 0.0 with np.errstate(invalid='ignore'):
if self.lattice == 'cubic': # using '*' for 'and'
return R[0] >= R[1]+epsilon and R[1] >= R[2]+epsilon and R[2] >= epsilon if self.system == 'cubic':
elif self.lattice == 'hexagonal': return np.where((rho[...,0] >= rho[...,1]) * \
return R[0] >= np.sqrt(3)*(R[1]-epsilon) and R[1] >= epsilon and R[2] >= epsilon (rho[...,1] >= rho[...,2]) * \
elif self.lattice == 'tetragonal': (rho[...,2] >= 0),True,False)
return R[0] >= R[1]-epsilon and R[1] >= epsilon and R[2] >= epsilon elif self.system == 'hexagonal':
elif self.lattice == 'orthorhombic': return np.where((rho[...,0] >= rho[...,1]*np.sqrt(3)) * \
return R[0] >= epsilon and R[1] >= epsilon and R[2] >= epsilon (rho[...,1] >= 0) * \
else: (rho[...,2] >= 0),True,False)
return True elif self.system == 'tetragonal':
return np.where((rho[...,0] >= rho[...,1]) * \
(rho[...,1] >= 0) * \
(rho[...,2] >= 0),True,False)
elif self.system == 'orthorhombic':
return np.where((rho[...,0] >= 0) * \
(rho[...,1] >= 0) * \
(rho[...,2] >= 0),True,False)
else:
return np.ones_like(rho[...,0],dtype=bool)
def inSST(self, #ToDo: IPF color in separate function
vector, def in_SST(self,vector,proper=False,color=False):
proper = False,
color = False):
""" """
Check whether given vector falls into standard stereographic triangle of own symmetry. Check whether given vector falls into standard stereographic triangle of own symmetry.
@ -244,7 +240,10 @@ class Symmetry:
... } ... }
""" """
if self.lattice == 'cubic': if(vector.shape[-1] != 3):
raise ValueError('Input is not a 3D vector field.')
if self.system == 'cubic':
basis = {'improper':np.array([ [-1. , 0. , 1. ], basis = {'improper':np.array([ [-1. , 0. , 1. ],
[ np.sqrt(2.) , -np.sqrt(2.) , 0. ], [ np.sqrt(2.) , -np.sqrt(2.) , 0. ],
[ 0. , np.sqrt(3.) , 0. ] ]), [ 0. , np.sqrt(3.) , 0. ] ]),
@ -252,7 +251,7 @@ class Symmetry:
[-np.sqrt(2.) , np.sqrt(2.) , 0. ], [-np.sqrt(2.) , np.sqrt(2.) , 0. ],
[ np.sqrt(3.) , 0. , 0. ] ]), [ np.sqrt(3.) , 0. , 0. ] ]),
} }
elif self.lattice == 'hexagonal': elif self.system == 'hexagonal':
basis = {'improper':np.array([ [ 0. , 0. , 1. ], basis = {'improper':np.array([ [ 0. , 0. , 1. ],
[ 1. , -np.sqrt(3.) , 0. ], [ 1. , -np.sqrt(3.) , 0. ],
[ 0. , 2. , 0. ] ]), [ 0. , 2. , 0. ] ]),
@ -260,7 +259,7 @@ class Symmetry:
[-1. , np.sqrt(3.) , 0. ], [-1. , np.sqrt(3.) , 0. ],
[ np.sqrt(3.) , -1. , 0. ] ]), [ np.sqrt(3.) , -1. , 0. ] ]),
} }
elif self.lattice == 'tetragonal': elif self.system == 'tetragonal':
basis = {'improper':np.array([ [ 0. , 0. , 1. ], basis = {'improper':np.array([ [ 0. , 0. , 1. ],
[ 1. , -1. , 0. ], [ 1. , -1. , 0. ],
[ 0. , np.sqrt(2.) , 0. ] ]), [ 0. , np.sqrt(2.) , 0. ] ]),
@ -268,7 +267,7 @@ class Symmetry:
[-1. , 1. , 0. ], [-1. , 1. , 0. ],
[ np.sqrt(2.) , 0. , 0. ] ]), [ np.sqrt(2.) , 0. , 0. ] ]),
} }
elif self.lattice == 'orthorhombic': elif self.system == 'orthorhombic':
basis = {'improper':np.array([ [ 0., 0., 1.], basis = {'improper':np.array([ [ 0., 0., 1.],
[ 1., 0., 0.], [ 1., 0., 0.],
[ 0., 1., 0.] ]), [ 0., 1., 0.] ]),
@ -278,43 +277,41 @@ class Symmetry:
} }
else: # direct exit for unspecified symmetry else: # direct exit for unspecified symmetry
if color: if color:
return (True,np.zeros(3,'d')) return (np.ones_like(vector[...,0],bool),np.zeros_like(vector))
else: else:
return True return np.ones_like(vector[...,0],bool)
v = np.array(vector,dtype=float)
if proper: # check both improper ... b_i = np.broadcast_to(basis['improper'],vector.shape+(3,))
theComponents = np.around(np.dot(basis['improper'],v),12) if proper:
inSST = np.all(theComponents >= 0.0) b_p = np.broadcast_to(basis['proper'], vector.shape+(3,))
if not inSST: # ... and proper SST improper = np.all(np.around(np.einsum('...ji,...i',b_i,vector),12)>=0.0,axis=-1,keepdims=True)
theComponents = np.around(np.dot(basis['proper'],v),12) theComponents = np.where(np.broadcast_to(improper,vector.shape),
inSST = np.all(theComponents >= 0.0) np.around(np.einsum('...ji,...i',b_i,vector),12),
np.around(np.einsum('...ji,...i',b_p,vector),12))
else: else:
v[2] = abs(v[2]) # z component projects identical vector_ = np.block([vector[...,0:2],np.abs(vector[...,2:3])]) # z component projects identical
theComponents = np.around(np.dot(basis['improper'],v),12) # for positive and negative values theComponents = np.around(np.einsum('...ji,...i',b_i,vector_),12)
inSST = np.all(theComponents >= 0.0)
in_SST = np.all(theComponents >= 0.0,axis=-1)
if color: # have to return color array if color: # have to return color array
if inSST: with np.errstate(invalid='ignore',divide='ignore'):
rgb = np.power(theComponents/np.linalg.norm(theComponents),0.5) # smoothen color ramps rgb = (theComponents/np.linalg.norm(theComponents,axis=-1,keepdims=True))**0.5 # smoothen color ramps
rgb = np.minimum(np.ones(3,dtype=float),rgb) # limit to maximum intensity rgb = np.minimum(1.,rgb) # limit to maximum intensity
rgb /= max(rgb) # normalize to (HS)V = 1 rgb /= np.max(rgb,axis=-1,keepdims=True) # normalize to (HS)V = 1
else: rgb[np.broadcast_to(~in_SST.reshape(vector[...,0].shape+(1,)),vector.shape)] = 0.0
rgb = np.zeros(3,dtype=float) return (in_SST,rgb)
return (inSST,rgb)
else: else:
return inSST return in_SST
# code derived from https://github.com/ezag/pyeuclid
# suggested reading: http://web.mit.edu/2.998/www/QuaternionReport1.pdf
# ****************************************************************************************** # ******************************************************************************************
class Lattice: class Lattice: # ToDo: Make a subclass of Symmetry!
""" """
Lattice system. Bravais lattice.
Currently, this contains only a mapping from Bravais lattice to symmetry This contains only a mapping from Bravais lattice to symmetry
and orientation relationships. It could include twin and slip systems. and orientation relationships. It could include twin and slip systems.
References References
@ -324,15 +321,15 @@ class Lattice:
""" """
lattices = { lattices = {
'triclinic':{'symmetry':None}, 'triclinic':{'system':None},
'bct':{'symmetry':'tetragonal'}, 'bct': {'system':'tetragonal'},
'hex':{'symmetry':'hexagonal'}, 'hex': {'system':'hexagonal'},
'fcc':{'symmetry':'cubic','c/a':1.0}, 'fcc': {'system':'cubic','c/a':1.0},
'bcc':{'symmetry':'cubic','c/a':1.0}, 'bcc': {'system':'cubic','c/a':1.0},
} }
def __init__(self, lattice): def __init__(self,lattice,c_over_a=None):
""" """
New lattice of given type. New lattice of given type.
@ -343,18 +340,23 @@ class Lattice:
""" """
self.lattice = lattice self.lattice = lattice
self.symmetry = Symmetry(self.lattices[lattice]['symmetry']) self.symmetry = Symmetry(self.lattices[lattice]['system'])
# transition to subclass
self.system = self.symmetry.system
self.in_SST = self.symmetry.in_SST
self.in_FZ = self.symmetry.in_FZ
self.in_disorientation_SST = self.symmetry.in_disorientation_SST
def __repr__(self): def __repr__(self):
"""Report basic lattice information.""" """Report basic lattice information."""
return f'Bravais lattice {self.lattice} ({self.symmetry} symmetry)' return f'Bravais lattice {self.lattice} ({self.symmetry} crystal system)'
# Kurdjomov--Sachs orientation relationship for fcc <-> bcc transformation # Kurdjomov--Sachs orientation relationship for fcc <-> bcc transformation
# from S. Morito et al., Journal of Alloys and Compounds 577:s587-s592, 2013 # 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 # also see K. Kitahara et al., Acta Materialia 54:1279-1288, 2006
KS = {'mapping':{'fcc':0,'bcc':1}, _KS = {'mapping':{'fcc':0,'bcc':1},
'planes': np.array([ '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]],
@ -408,7 +410,7 @@ class Lattice:
# Greninger--Troiano orientation relationship for fcc <-> bcc transformation # Greninger--Troiano orientation relationship for fcc <-> bcc transformation
# from Y. He et al., Journal of Applied Crystallography 39:72-81, 2006 # from Y. He et al., Journal of Applied Crystallography 39:72-81, 2006
GT = {'mapping':{'fcc':0,'bcc':1}, _GT = {'mapping':{'fcc':0,'bcc':1},
'planes': np.array([ 'planes': np.array([
[[ 1, 1, 1],[ 1, 0, 1]], [[ 1, 1, 1],[ 1, 0, 1]],
[[ 1, 1, 1],[ 1, 1, 0]], [[ 1, 1, 1],[ 1, 1, 0]],
@ -462,7 +464,7 @@ class Lattice:
# Greninger--Troiano' orientation relationship for fcc <-> bcc transformation # Greninger--Troiano' orientation relationship for fcc <-> bcc transformation
# from Y. He et al., Journal of Applied Crystallography 39:72-81, 2006 # from Y. He et al., Journal of Applied Crystallography 39:72-81, 2006
GTprime = {'mapping':{'fcc':0,'bcc':1}, _GTprime = {'mapping':{'fcc':0,'bcc':1},
'planes': np.array([ 'planes': np.array([
[[ 7, 17, 17],[ 12, 5, 17]], [[ 7, 17, 17],[ 12, 5, 17]],
[[ 17, 7, 17],[ 17, 12, 5]], [[ 17, 7, 17],[ 17, 12, 5]],
@ -516,7 +518,7 @@ class Lattice:
# Nishiyama--Wassermann orientation relationship for fcc <-> bcc transformation # Nishiyama--Wassermann orientation relationship for fcc <-> bcc transformation
# from H. Kitahara et al., Materials Characterization 54:378-386, 2005 # from H. Kitahara et al., Materials Characterization 54:378-386, 2005
NW = {'mapping':{'fcc':0,'bcc':1}, _NW = {'mapping':{'fcc':0,'bcc':1},
'planes': np.array([ '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]],
@ -546,7 +548,7 @@ class Lattice:
# Pitsch orientation relationship for fcc <-> bcc transformation # Pitsch orientation relationship for fcc <-> bcc transformation
# from Y. He et al., Acta Materialia 53:1179-1190, 2005 # from Y. He et al., Acta Materialia 53:1179-1190, 2005
Pitsch = {'mapping':{'fcc':0,'bcc':1}, _Pitsch = {'mapping':{'fcc':0,'bcc':1},
'planes': np.array([ 'planes': np.array([
[[ 0, 1, 0],[ -1, 0, 1]], [[ 0, 1, 0],[ -1, 0, 1]],
[[ 0, 0, 1],[ 1, -1, 0]], [[ 0, 0, 1],[ 1, -1, 0]],
@ -576,7 +578,7 @@ class Lattice:
# Bain orientation relationship for fcc <-> bcc transformation # Bain orientation relationship for fcc <-> bcc transformation
# from Y. He et al., Journal of Applied Crystallography 39:72-81, 2006 # from Y. He et al., Journal of Applied Crystallography 39:72-81, 2006
Bain = {'mapping':{'fcc':0,'bcc':1}, _Bain = {'mapping':{'fcc':0,'bcc':1},
'planes': np.array([ 'planes': np.array([
[[ 1, 0, 0],[ 1, 0, 0]], [[ 1, 0, 0],[ 1, 0, 0]],
[[ 0, 1, 0],[ 0, 1, 0]], [[ 0, 1, 0],[ 0, 1, 0]],
@ -586,7 +588,8 @@ class Lattice:
[[ 0, 0, 1],[ 1, 0, 1]], [[ 0, 0, 1],[ 1, 0, 1]],
[[ 1, 0, 0],[ 1, 1, 0]]],dtype='float')} [[ 1, 0, 0],[ 1, 1, 0]]],dtype='float')}
def relationOperations(self,model):
def relation_operations(self,model):
""" """
Crystallographic orientation relationships for phase transformations. Crystallographic orientation relationships for phase transformations.
@ -608,8 +611,8 @@ class Lattice:
https://doi.org/10.1016/j.actamat.2004.11.021 https://doi.org/10.1016/j.actamat.2004.11.021
""" """
models={'KS':self.KS, 'GT':self.GT, 'GT_prime':self.GTprime, models={'KS':self._KS, 'GT':self._GT, 'GT_prime':self._GTprime,
'NW':self.NW, 'Pitsch': self.Pitsch, 'Bain':self.Bain} 'NW':self._NW, 'Pitsch': self._Pitsch, 'Bain':self._Bain}
try: try:
relationship = models[model] relationship = models[model]
except KeyError : except KeyError :
@ -635,6 +638,8 @@ class Lattice:
otherDir = miller[otherDir_id]/ np.linalg.norm(miller[otherDir_id]) otherDir = miller[otherDir_id]/ np.linalg.norm(miller[otherDir_id])
otherMatrix = np.array([otherDir,np.cross(otherPlane,otherDir),otherPlane]) otherMatrix = np.array([otherDir,np.cross(otherPlane,otherDir),otherPlane])
r['rotations'].append(Rotation.from_matrix(np.dot(otherMatrix.T,myMatrix))) r['rotations'].append(np.dot(otherMatrix.T,myMatrix))
r['rotations'] = np.array(r['rotations'])
return r return r

143
python/damask/_orientation.py
View File

@ -3,7 +3,7 @@ import numpy as np
from . import Lattice from . import Lattice
from . import Rotation from . import Rotation
class Orientation: class Orientation: # ToDo: make subclass of lattice and Rotation?
""" """
Crystallographic orientation. Crystallographic orientation.
@ -39,9 +39,12 @@ class Orientation:
else: else:
self.rotation = Rotation.from_quaternion(rotation) # assume quaternion self.rotation = Rotation.from_quaternion(rotation) # assume quaternion
if self.rotation.quaternion.shape != (4,): def __getitem__(self,item):
raise NotImplementedError('Support for multiple rotations missing') """Iterate over leading/leftmost dimension of Orientation array."""
return self.__class__(self.rotation[item],self.lattice)
# ToDo: Discuss vectorization/calling signature
def disorientation(self, def disorientation(self,
other, other,
SST = True, SST = True,
@ -58,8 +61,8 @@ class Orientation:
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.equivalent if SST else self.equivalent[0] #ToDo: This is just me! # take all or only first sym operation
otherSymEqs = other.equivalentOrientations() otherSymEqs = other.equivalent
for i,sA in enumerate(mySymEqs): for i,sA in enumerate(mySymEqs):
aInv = sA.rotation.inversed() aInv = sA.rotation.inversed()
@ -68,8 +71,8 @@ class Orientation:
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.as_Rodrigues(vector=True)) \ breaker = self.lattice.in_FZ(r.as_Rodrigues(vector=True)) \
and (not SST or other.lattice.symmetry.inDisorientationSST(r.as_Rodrigues(vector=True))) and (not SST or other.lattice.in_disorientation_SST(r.as_Rodrigues(vector=True)))
if breaker: break if breaker: break
if breaker: break if breaker: break
if breaker: break if breaker: break
@ -77,79 +80,123 @@ class Orientation:
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,
# self-->other: True, self<--other: False # self-->other: True, self<--other: False
def inFZ(self):
return self.lattice.symmetry.inFZ(self.rotation.as_Rodrigues(vector=True)) @property
def in_FZ(self):
"""Check if orientations fall into Fundamental Zone."""
return self.lattice.in_FZ(self.rotation.as_Rodrigues(vector=True))
def equivalentOrientations(self,members=[]): @property
"""List of orientations which are symmetrically equivalent.""" def equivalent(self):
try: """
iter(members) # asking for (even empty) list of members? Orientations which are symmetrically equivalent.
except TypeError:
return self.__class__(self.lattice.symmetry.symmetryOperations(members)*self.rotation,self.lattice) # no, return rotation object
else:
return [self.__class__(q*self.rotation,self.lattice) \
for q in self.lattice.symmetry.symmetryOperations(members)] # yes, return list of rotations
def relatedOrientations(self,model): One dimension (length according to number of symmetrically equivalent orientations)
"""List of orientations related by the given orientation relationship.""" is added to the left of the Rotation array.
r = self.lattice.relationOperations(model)
return [self.__class__(o*self.rotation,r['lattice']) for o in r['rotations']] """
o = self.lattice.symmetry.symmetry_operations
o = o.reshape(o.shape[:1]+(1,)*len(self.rotation.shape)+(4,))
o = Rotation(np.broadcast_to(o,o.shape[:1]+self.rotation.quaternion.shape))
s = np.broadcast_to(self.rotation.quaternion,o.shape[:1]+self.rotation.quaternion.shape)
return self.__class__(o@Rotation(s),self.lattice)
def related(self,model):
"""
Orientations related by the given orientation relationship.
One dimension (length according to number of related orientations)
is added to the left of the Rotation array.
"""
o = Rotation.from_matrix(self.lattice.relation_operations(model)['rotations']).as_quaternion()
o = o.reshape(o.shape[:1]+(1,)*len(self.rotation.shape)+(4,))
o = Rotation(np.broadcast_to(o,o.shape[:1]+self.rotation.quaternion.shape))
s = np.broadcast_to(self.rotation.quaternion,o.shape[:1]+self.rotation.quaternion.shape)
return self.__class__(o@Rotation(s),self.lattice.relation_operations(model)['lattice'])
@property
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(): eq = self.equivalent
if self.lattice.symmetry.inFZ(me.rotation.as_Rodrigues(vector=True)): break in_FZ = eq.in_FZ
return self.__class__(me.rotation,self.lattice) # remove duplicates (occur for highly symmetric orientations)
found = np.zeros_like(in_FZ[0],dtype=bool)
q = self.rotation.quaternion[0]
for s in range(in_FZ.shape[0]):
#something fishy... why does q needs to be initialized?
q = np.where(np.expand_dims(np.logical_and(in_FZ[s],~found),-1),eq.rotation.quaternion[s],q)
found = np.logical_or(in_FZ[s],found)
return self.__class__(q,self.lattice)
def inversePole(self, def inverse_pole(self,axis,proper=False,SST=True):
axis,
proper = False,
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:
for i,o in enumerate(self.equivalentOrientations()): # test all symmetric equivalent quaternions eq = self.equivalent
pole = o.rotation*axis # align crystal direction to axis pole = eq.rotation @ np.broadcast_to(axis/np.linalg.norm(axis),eq.rotation.shape+(3,))
if self.lattice.symmetry.inSST(pole,proper): break # found SST version in_SST = self.lattice.in_SST(pole,proper=proper)
# remove duplicates (occur for highly symmetric orientations)
found = np.zeros_like(in_SST[0],dtype=bool)
p = pole[0]
for s in range(in_SST.shape[0]):
p = np.where(np.expand_dims(np.logical_and(in_SST[s],~found),-1),pole[s],p)
found = np.logical_or(in_SST[s],found)
return p
else: else:
pole = self.rotation*axis # align crystal direction to axis return self.rotation @ np.broadcast_to(axis/np.linalg.norm(axis),self.rotation.shape+(3,))
return (pole,i if SST else 0)
def IPFcolor(self,axis):
def IPF_color(self,axis): #ToDo axis or direction?
"""TSL color of inverse pole figure for given axis.""" """TSL color of inverse pole figure for given axis."""
color = np.zeros(3,'d') eq = self.equivalent
pole = eq.rotation @ np.broadcast_to(axis/np.linalg.norm(axis),eq.rotation.shape+(3,))
in_SST, color = self.lattice.in_SST(pole,color=True)
for o in self.equivalentOrientations(): # remove duplicates (occur for highly symmetric orientations)
pole = o.rotation*axis # align crystal direction to axis found = np.zeros_like(in_SST[0],dtype=bool)
inSST,color = self.lattice.symmetry.inSST(pole,color=True) c = color[0]
if inSST: break for s in range(in_SST.shape[0]):
c = np.where(np.expand_dims(np.logical_and(in_SST[s],~found),-1),color[s],c)
found = np.logical_or(in_SST[s],found)
return color return c
# ToDo: Discuss vectorization/calling signature
@staticmethod @staticmethod
def fromAverage(orientations, def from_average(orientations,
weights = []): weights = []):
"""Create orientation from average of list of orientations.""" """Create orientation from average of list of orientations."""
# further read: Orientation distribution analysis in deformed grains
# https://doi.org/10.1107/S0021889801003077
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.equivalent[
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.from_average(closest,weights),ref.lattice)
# ToDo: Discuss vectorization/calling signature
def average(self,other): def average(self,other):
"""Calculate the average rotation.""" """Calculate the average rotation."""
return Orientation.fromAverage([self,other]) return Orientation.from_average([self,other])

45
python/damask/_result.py
View File

@ -1,4 +1,4 @@
import multiprocessing import multiprocessing as mp
import re import re
import inspect import inspect
import glob import glob
@ -11,7 +11,9 @@ from functools import partial
import h5py import h5py
import numpy as np import numpy as np
from numpy.lib import recfunctions as rfn
import damask
from . import VTK from . import VTK
from . import Table from . import Table
from . import Rotation from . import Rotation
@ -20,7 +22,6 @@ from . import Environment
from . import grid_filters from . import grid_filters
from . import mechanics from . import mechanics
from . import util from . import util
from . import version
class Result: class Result:
@ -413,17 +414,18 @@ class Result:
for i in self.iterate('increments'): for i in self.iterate('increments'):
message += f'\n{i} ({self.times[self.increments.index(i)]}s)\n' message += f'\n{i} ({self.times[self.increments.index(i)]}s)\n'
for o,p in zip(['constituents','materialpoints'],['con_physics','mat_physics']): for o,p in zip(['constituents','materialpoints'],['con_physics','mat_physics']):
message += f' {o[:-1]}\n'
for oo in self.iterate(o): for oo in self.iterate(o):
message += f' {oo}\n' message += f' {oo}\n'
for pp in self.iterate(p): for pp in self.iterate(p):
message += f' {pp}\n' message += f' {pp}\n'
group = '/'.join([i,o[:-1],oo,pp]) # o[:-1]: plural/singular issue group = '/'.join([i,o[:-1],oo,pp]) # o[:-1]: plural/singular issue
for d in f[group].keys(): for d in f[group].keys():
try: try:
dataset = f['/'.join([group,d])] dataset = f['/'.join([group,d])]
unit = f" / {dataset.attrs['Unit'].decode()}" if 'Unit' in dataset.attrs else '' unit = f" / {dataset.attrs['Unit'].decode()}" if 'Unit' in dataset.attrs else ''
description = dataset.attrs['Description'].decode() description = dataset.attrs['Description'].decode()
message += f' {d}{unit}: {description}\n' message += f' {d}{unit}: {description}\n'
except KeyError: except KeyError:
pass pass
return message return message
@ -703,7 +705,6 @@ class Result:
label,p = 'intermediate',1 label,p = 'intermediate',1
elif eigenvalue == 'min': elif eigenvalue == 'min':
label,p = 'minimum',0 label,p = 'minimum',0
print('p',eigenvalue)
return { return {
'data': mechanics.eigenvectors(T_sym['data'])[:,p], 'data': mechanics.eigenvectors(T_sym['data'])[:,p],
'label': f"v_{eigenvalue}({T_sym['label']})", 'label': f"v_{eigenvalue}({T_sym['label']})",
@ -731,29 +732,23 @@ class Result:
@staticmethod @staticmethod
def _add_IPFcolor(q,l): def _add_IPF_color(q,l):
d = np.array(l) m = util.scale_to_coprime(np.array(l))
d_unit = d/np.linalg.norm(d)
m = util.scale_to_coprime(d)
colors = np.empty((len(q['data']),3),np.uint8)
lattice = q['meta']['Lattice'] o = Orientation(Rotation(rfn.structured_to_unstructured(q['data'])),
lattice = q['meta']['Lattice'])
for i,qu in enumerate(q['data']):
o = Orientation(np.array([qu['w'],qu['x'],qu['y'],qu['z']]),lattice).reduced()
colors[i] = np.uint8(o.IPFcolor(d_unit)*255)
return { return {
'data': colors, 'data': np.uint8(o.IPF_color(l)*255),
'label': 'IPFcolor_[{} {} {}]'.format(*m), 'label': 'IPFcolor_[{} {} {}]'.format(*m),
'meta' : { 'meta' : {
'Unit': 'RGB (8bit)', 'Unit': '8-bit RGB',
'Lattice': lattice, 'Lattice': q['meta']['Lattice'],
'Description': 'Inverse Pole Figure (IPF) colors along sample direction [{} {} {}]'.format(*m), 'Description': 'Inverse Pole Figure (IPF) colors along sample direction [{} {} {}]'.format(*m),
'Creator': inspect.stack()[0][3][1:] 'Creator': inspect.stack()[0][3][1:]
} }
} }
def add_IPFcolor(self,q,l): def add_IPF_color(self,q,l):
""" """
Add RGB color tuple of inverse pole figure (IPF) color. Add RGB color tuple of inverse pole figure (IPF) color.
@ -765,7 +760,7 @@ class Result:
Lab frame direction for inverse pole figure. Lab frame direction for inverse pole figure.
""" """
self._add_generic_pointwise(self._add_IPFcolor,{'q':q},{'l':l}) self._add_generic_pointwise(self._add_IPF_color,{'q':q},{'l':l})
@staticmethod @staticmethod
@ -1066,8 +1061,8 @@ class Result:
""" """
num_threads = Environment().options['DAMASK_NUM_THREADS'] num_threads = Environment().options['DAMASK_NUM_THREADS']
pool = multiprocessing.Pool(int(num_threads) if num_threads is not None else None) pool = mp.Pool(int(num_threads) if num_threads is not None else None)
lock = multiprocessing.Manager().Lock() lock = mp.Manager().Lock()
groups = self.groups_with_datasets(datasets.values()) groups = self.groups_with_datasets(datasets.values())
default_arg = partial(self._job,func=func,datasets=datasets,args=args,lock=lock) default_arg = partial(self._job,func=func,datasets=datasets,args=args,lock=lock)
@ -1090,7 +1085,7 @@ class Result:
for l,v in result[1]['meta'].items(): for l,v in result[1]['meta'].items():
dataset.attrs[l]=v.encode() dataset.attrs[l]=v.encode()
creator = f"damask.Result.{dataset.attrs['Creator'].decode()} v{version}" creator = f"damask.Result.{dataset.attrs['Creator'].decode()} v{damask.version}"
dataset.attrs['Creator'] = creator.encode() dataset.attrs['Creator'] = creator.encode()
except (OSError,RuntimeError) as err: except (OSError,RuntimeError) as err:
@ -1222,7 +1217,7 @@ class Result:
elif mode.lower()=='point': elif mode.lower()=='point':
v = VTK.from_polyData(self.cell_coordinates()) v = VTK.from_polyData(self.cell_coordinates())
N_digits = int(np.floor(np.log10(int(self.increments[-1][3:]))))+1 N_digits = int(np.floor(np.log10(max(1,int(self.increments[-1][3:])))))+1
for inc in util.show_progress(self.iterate('increments'),len(self.selection['increments'])): for inc in util.show_progress(self.iterate('increments'),len(self.selection['increments'])):

46
python/damask/_rotation.py
View File

@ -53,6 +53,8 @@ class Rotation:
Use .from_quaternion to perform a sanity check. Use .from_quaternion to perform a sanity check.
""" """
if quaternion.shape[-1] != 4:
raise ValueError('Not a quaternion')
self.quaternion = quaternion.copy() self.quaternion = quaternion.copy()
@ -61,6 +63,7 @@ class Rotation:
return self.quaternion.shape[:-1] return self.quaternion.shape[:-1]
# ToDo: Check difference __copy__ vs __deepcopy__
def __copy__(self): def __copy__(self):
"""Copy.""" """Copy."""
return self.__class__(self.quaternion) return self.__class__(self.quaternion)
@ -71,7 +74,7 @@ class Rotation:
def __repr__(self): def __repr__(self):
"""Orientation displayed as unit quaternion, rotation matrix, and Bunge-Euler angles.""" """Orientation displayed as unit quaternion, rotation matrix, and Bunge-Euler angles."""
if self.quaternion.shape != (4,): if self.quaternion.shape != (4,):
raise NotImplementedError('Support for multiple rotations missing') return 'Quaternions:\n'+str(self.quaternion) # ToDo: could be nicer ...
return '\n'.join([ return '\n'.join([
'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)), 'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)),
'Matrix:\n{}'.format(np.round(self.as_matrix(),8)), 'Matrix:\n{}'.format(np.round(self.as_matrix(),8)),
@ -79,6 +82,19 @@ class Rotation:
]) ])
def __getitem__(self,item):
"""Iterate over leading/leftmost dimension of Rotation array."""
if self.shape == (): return self.copy()
if isinstance(item,tuple) and len(item) >= len(self):
raise IndexError('Too many indices')
return self.__class__(self.quaternion[item])
def __len__(self):
"""Length of leading/leftmost dimension of Rotation array."""
return 0 if self.shape == () else self.shape[0]
def __matmul__(self, other): def __matmul__(self, other):
""" """
Rotation of vector, second or fourth order tensor, or rotation object. Rotation of vector, second or fourth order tensor, or rotation object.
@ -88,6 +104,11 @@ class Rotation:
other : numpy.ndarray or Rotation other : numpy.ndarray or Rotation
Vector, second or fourth order tensor, or rotation object that is rotated. Vector, second or fourth order tensor, or rotation object that is rotated.
Returns
-------
other_rot : numpy.ndarray or Rotation
Rotated vector, second or fourth order tensor, or rotation object.
""" """
if isinstance(other, Rotation): if isinstance(other, Rotation):
q_m = self.quaternion[...,0:1] q_m = self.quaternion[...,0:1]
@ -155,7 +176,7 @@ class Rotation:
def broadcast_to(self,shape): def broadcast_to(self,shape):
if isinstance(shape,int): shape = (shape,) if isinstance(shape,(int,np.integer)): shape = (shape,)
if self.shape == (): if self.shape == ():
q = np.broadcast_to(self.quaternion,shape+(4,)) q = np.broadcast_to(self.quaternion,shape+(4,))
else: else:
@ -178,7 +199,7 @@ class Rotation:
""" """
if self.quaternion.shape != (4,) or other.quaternion.shape != (4,): if self.quaternion.shape != (4,) or other.quaternion.shape != (4,):
raise NotImplementedError('Support for multiple rotations missing') raise NotImplementedError('Support for multiple rotations missing')
return Rotation.fromAverage([self,other]) return Rotation.from_average([self,other])
################################################################################################ ################################################################################################
@ -273,7 +294,11 @@ class Rotation:
""" """
ro = Rotation._qu2ro(self.quaternion) ro = Rotation._qu2ro(self.quaternion)
return ro[...,:3]*ro[...,3] if vector else ro if vector:
with np.errstate(invalid='ignore'):
return ro[...,:3]*ro[...,3:4]
else:
return ro
def as_homochoric(self): def as_homochoric(self):
""" """
@ -299,6 +324,7 @@ class Rotation:
""" """
return Rotation._qu2cu(self.quaternion) return Rotation._qu2cu(self.quaternion)
@property
def M(self): # ToDo not sure about the name: as_M or M? we do not have a from_M def M(self): # ToDo not sure about the name: as_M or M? we do not have a from_M
""" """
Intermediate representation supporting quaternion averaging. Intermediate representation supporting quaternion averaging.
@ -555,7 +581,7 @@ class Rotation:
@staticmethod @staticmethod
def fromAverage(rotations,weights = None): def from_average(rotations,weights = None):
""" """
Average rotation. Average rotation.
@ -580,8 +606,8 @@ class Rotation:
weights = np.ones(N,dtype='i') weights = np.ones(N,dtype='i')
for i,(r,n) in enumerate(zip(rotations,weights)): for i,(r,n) in enumerate(zip(rotations,weights)):
M = r.M() * n if i == 0 \ M = r.M * n if i == 0 \
else M + r.M() * n # noqa add (multiples) of this rotation to average noqa else M + r.M * n # noqa add (multiples) of this rotation to average noqa
eig, vec = np.linalg.eig(M/N) eig, vec = np.linalg.eig(M/N)
return Rotation.from_quaternion(np.real(vec.T[eig.argmax()]),accept_homomorph = True) return Rotation.from_quaternion(np.real(vec.T[eig.argmax()]),accept_homomorph = True)
@ -593,7 +619,7 @@ class Rotation:
elif hasattr(shape, '__iter__'): elif hasattr(shape, '__iter__'):
r = np.random.random(tuple(shape)+(3,)) r = np.random.random(tuple(shape)+(3,))
else: else:
r = np.random.random((shape,3)) r = np.random.rand(shape,3)
A = np.sqrt(r[...,2]) A = np.sqrt(r[...,2])
B = np.sqrt(1.0-r[...,2]) B = np.sqrt(1.0-r[...,2])
@ -604,11 +630,9 @@ class Rotation:
return Rotation(q.reshape(r.shape[:-1]+(4,)) if shape is not None else q)._standardize() return Rotation(q.reshape(r.shape[:-1]+(4,)) if shape is not None else q)._standardize()
# for compatibility (old names do not follow convention) # for compatibility (old names do not follow convention)
asM = M
fromQuaternion = from_quaternion
fromEulers = from_Eulers fromEulers = from_Eulers
fromQuaternion = from_quaternion
asAxisAngle = as_axis_angle asAxisAngle = as_axis_angle
__mul__ = __matmul__ __mul__ = __matmul__

5
python/damask/_table.py
View File

@ -19,11 +19,12 @@ class Table:
Data. Column labels from a pandas.DataFrame will be replaced. Data. Column labels from a pandas.DataFrame will be replaced.
shapes : dict with str:tuple pairs shapes : dict with str:tuple pairs
Shapes of the columns. Example 'F':(3,3) for a deformation gradient. Shapes of the columns. Example 'F':(3,3) for a deformation gradient.
comments : iterable of str, optional comments : str or iterable of str, optional
Additional, human-readable information. Additional, human-readable information.
""" """
self.comments = [] if comments is None else [c for c in comments] comments_ = [comments] if isinstance(comments,str) else comments
self.comments = [] if comments_ is None else [c for c in comments_]
self.data = pd.DataFrame(data=data) self.data = pd.DataFrame(data=data)
self.shapes = { k:(v,) if isinstance(v,(np.int,int)) else v for k,v in shapes.items() } self.shapes = { k:(v,) if isinstance(v,(np.int,int)) else v for k,v in shapes.items() }
self._label_condensed() self._label_condensed()

32
python/damask/_vtk.py
View File

@ -1,3 +1,4 @@
import multiprocessing as mp
from pathlib import Path from pathlib import Path
import pandas as pd import pandas as pd
@ -6,9 +7,9 @@ import vtk
from vtk.util.numpy_support import numpy_to_vtk as np_to_vtk from vtk.util.numpy_support import numpy_to_vtk as np_to_vtk
from vtk.util.numpy_support import numpy_to_vtkIdTypeArray as np_to_vtkIdTypeArray from vtk.util.numpy_support import numpy_to_vtkIdTypeArray as np_to_vtkIdTypeArray
import damask
from . import Table from . import Table
from . import Environment from . import Environment
from . import version
class VTK: class VTK:
@ -159,8 +160,11 @@ class VTK:
return VTK(geom) return VTK(geom)
@staticmethod
def write(self,fname): def _write(writer):
"""Wrapper for parallel writing."""
writer.Write()
def write(self,fname,parallel=True):
""" """
Write to file. Write to file.
@ -168,6 +172,8 @@ class VTK:
---------- ----------
fname : str or pathlib.Path fname : str or pathlib.Path
Filename for writing. Filename for writing.
parallel : boolean, optional
Write data in parallel background process. Defaults to True.
""" """
if isinstance(self.geom,vtk.vtkRectilinearGrid): if isinstance(self.geom,vtk.vtkRectilinearGrid):
@ -185,8 +191,11 @@ class VTK:
writer.SetCompressorTypeToZLib() writer.SetCompressorTypeToZLib()
writer.SetDataModeToBinary() writer.SetDataModeToBinary()
writer.SetInputData(self.geom) writer.SetInputData(self.geom)
if parallel:
writer.Write() mp_writer = mp.Process(target=self._write,args=(writer,))
mp_writer.start()
else:
writer.Write()
# Check https://blog.kitware.com/ghost-and-blanking-visibility-changes/ for missing data # Check https://blog.kitware.com/ghost-and-blanking-visibility-changes/ for missing data
@ -197,14 +206,21 @@ class VTK:
N_cells = self.geom.GetNumberOfCells() N_cells = self.geom.GetNumberOfCells()
if isinstance(data,np.ndarray): if isinstance(data,np.ndarray):
d = np_to_vtk(num_array=data.reshape(data.shape[0],-1),deep=True)
if label is None: if label is None:
raise ValueError('No label defined for numpy.ndarray') raise ValueError('No label defined for numpy.ndarray')
if data.dtype in [np.float64, np.float128]: # avoid large files
d = np_to_vtk(num_array=data.astype(np.float32).reshape(data.shape[0],-1),deep=True)
else:
d = np_to_vtk(num_array=data.reshape(data.shape[0],-1),deep=True)
d.SetName(label) d.SetName(label)
if data.shape[0] == N_cells: if data.shape[0] == N_cells:
self.geom.GetCellData().AddArray(d) self.geom.GetCellData().AddArray(d)
elif data.shape[0] == N_points: elif data.shape[0] == N_points:
self.geom.GetPointData().AddArray(d) self.geom.GetPointData().AddArray(d)
else:
raise ValueError(f'Invalid shape {data.shape[0]}')
elif isinstance(data,pd.DataFrame): elif isinstance(data,pd.DataFrame):
raise NotImplementedError('pd.DataFrame') raise NotImplementedError('pd.DataFrame')
elif isinstance(data,Table): elif isinstance(data,Table):
@ -216,7 +232,7 @@ class VTK:
def __repr__(self): def __repr__(self):
"""ASCII representation of the VTK data.""" """ASCII representation of the VTK data."""
writer = vtk.vtkDataSetWriter() writer = vtk.vtkDataSetWriter()
writer.SetHeader(f'# DAMASK.VTK v{version}') writer.SetHeader(f'# damask.VTK v{damask.version}')
writer.WriteToOutputStringOn() writer.WriteToOutputStringOn()
writer.SetInputData(self.geom) writer.SetInputData(self.geom)
writer.Write() writer.Write()
@ -242,7 +258,7 @@ class VTK:
ren.AddActor(actor) ren.AddActor(actor)
ren.SetBackground(0.2,0.2,0.2) ren.SetBackground(0.2,0.2,0.2)
window.SetSize(Environment().screen_width,Environment().screen_height) window.SetSize(Environment().screen_size[0],Environment().screen_size[1])
iren = vtk.vtkRenderWindowInteractor() iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(window) iren.SetRenderWindow(window)

95
python/tests/conftest.py
View File

@ -1,7 +1,13 @@
from pathlib import Path from pathlib import Path
import numpy as np
import pytest import pytest
# Use to monkeypatch damask.version (for comparsion to reference files that contain version information)
def pytest_configure():
pytest.dummy_version = '99.99.99-9999-pytest'
def pytest_addoption(parser): def pytest_addoption(parser):
parser.addoption("--update", parser.addoption("--update",
action="store_true", action="store_true",
@ -16,3 +22,92 @@ def update(request):
def reference_dir_base(): def reference_dir_base():
"""Directory containing reference results.""" """Directory containing reference results."""
return Path(__file__).parent/'reference' return Path(__file__).parent/'reference'
@pytest.fixture
def set_of_quaternions():
"""A set of n random rotations."""
def random_quaternions(N):
r = np.random.rand(N,3)
A = np.sqrt(r[:,2])
B = np.sqrt(1.0-r[:,2])
qu = np.column_stack([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])
qu[:,0]*=np.sign(qu[:,0])
return qu
n = 1100
scatter=1.e-2
specials = np.array([
[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,-1.0, 0.0, 0.0],
[0.0, 0.0,-1.0, 0.0],
[0.0, 0.0, 0.0,-1.0],
#----------------------
[1.0, 1.0, 0.0, 0.0],
[1.0, 0.0, 1.0, 0.0],
[1.0, 0.0, 0.0, 1.0],
[0.0, 1.0, 1.0, 0.0],
[0.0, 1.0, 0.0, 1.0],
[0.0, 0.0, 1.0, 1.0],
#----------------------
[1.0,-1.0, 0.0, 0.0],
[1.0, 0.0,-1.0, 0.0],
[1.0, 0.0, 0.0,-1.0],
[0.0, 1.0,-1.0, 0.0],
[0.0, 1.0, 0.0,-1.0],
[0.0, 0.0, 1.0,-1.0],
#----------------------
[0.0, 1.0,-1.0, 0.0],
[0.0, 1.0, 0.0,-1.0],
[0.0, 0.0, 1.0,-1.0],
#----------------------
[0.0,-1.0,-1.0, 0.0],
[0.0,-1.0, 0.0,-1.0],
[0.0, 0.0,-1.0,-1.0],
#----------------------
[1.0, 1.0, 1.0, 0.0],
[1.0, 1.0, 0.0, 1.0],
[1.0, 0.0, 1.0, 1.0],
[1.0,-1.0, 1.0, 0.0],
[1.0,-1.0, 0.0, 1.0],
[1.0, 0.0,-1.0, 1.0],
[1.0, 1.0,-1.0, 0.0],
[1.0, 1.0, 0.0,-1.0],
[1.0, 0.0, 1.0,-1.0],
[1.0,-1.0,-1.0, 0.0],
[1.0,-1.0, 0.0,-1.0],
[1.0, 0.0,-1.0,-1.0],
#----------------------
[0.0, 1.0, 1.0, 1.0],
[0.0, 1.0,-1.0, 1.0],
[0.0, 1.0, 1.0,-1.0],
[0.0,-1.0, 1.0, 1.0],
[0.0,-1.0,-1.0, 1.0],
[0.0,-1.0, 1.0,-1.0],
[0.0,-1.0,-1.0,-1.0],
#----------------------
[1.0, 1.0, 1.0, 1.0],
[1.0,-1.0, 1.0, 1.0],
[1.0, 1.0,-1.0, 1.0],
[1.0, 1.0, 1.0,-1.0],
[1.0,-1.0,-1.0, 1.0],
[1.0,-1.0, 1.0,-1.0],
[1.0, 1.0,-1.0,-1.0],
[1.0,-1.0,-1.0,-1.0],
])
specials /= np.linalg.norm(specials,axis=1).reshape(-1,1)
specials_scatter = specials + np.broadcast_to(np.random.rand(4)*scatter,specials.shape)
specials_scatter /= np.linalg.norm(specials_scatter,axis=1).reshape(-1,1)
specials_scatter[specials_scatter[:,0]<0]*=-1
return np.array([s for s in specials] + \
[s for s in specials_scatter] + \
[s for s in random_quaternions(n-2*len(specials))])

1290
python/tests/reference/Colormap/binary.json Normal file
View File

File diff suppressed because it is too large Load Diff

1
python/tests/reference/Colormap/binary.legend Normal file
View File

File diff suppressed because one or more lines are too long

258
python/tests/reference/Colormap/binary.msh Normal file
View File

@ -0,0 +1,258 @@
View.ColorTable = {
255.0,255.0,255.0,
254.0,254.0,254.0,
253.0,253.0,253.0,
252.0,252.0,252.0,
251.0,251.0,251.0,
250.0,250.0,250.0,
249.0,249.0,249.0,
248.0,248.0,248.0,
247.0,247.0,247.0,
246.0,246.0,246.0,
245.0,245.0,245.0,
244.0,244.0,244.0,
243.0,243.0,243.0,
242.0,242.0,242.0,
241.0,241.0,241.0,
240.0,240.0,240.0,
239.0,239.0,239.0,
238.0,238.0,238.0,
237.0,237.0,237.0,
236.0,236.0,236.0,
235.00000000000003,235.00000000000003,235.00000000000003,
234.0,234.0,234.0,
233.0,233.0,233.0,
232.0,232.0,232.0,
231.0,231.0,231.0,
230.0,230.0,230.0,
229.0,229.0,229.0,
228.0,228.0,228.0,
227.0,227.0,227.0,
226.0,226.0,226.0,
225.0,225.0,225.0,
224.0,224.0,224.0,
223.0,223.0,223.0,
222.0,222.0,222.0,
221.0,221.0,221.0,
220.0,220.0,220.0,
219.00000000000003,219.00000000000003,219.00000000000003,
218.00000000000003,218.00000000000003,218.00000000000003,
217.0,217.0,217.0,
216.0,216.0,216.0,
215.0,215.0,215.0,
214.0,214.0,214.0,
213.0,213.0,213.0,
211.99999999999997,211.99999999999997,211.99999999999997,
211.0,211.0,211.0,
210.0,210.0,210.0,
209.0,209.0,209.0,
208.0,208.0,208.0,
207.0,207.0,207.0,
206.0,206.0,206.0,
205.0,205.0,205.0,
204.0,204.0,204.0,
203.00000000000003,203.00000000000003,203.00000000000003,
202.00000000000003,202.00000000000003,202.00000000000003,
201.0,201.0,201.0,
200.0,200.0,200.0,
199.0,199.0,199.0,
198.0,198.0,198.0,
197.0,197.0,197.0,
195.99999999999997,195.99999999999997,195.99999999999997,
195.0,195.0,195.0,
194.0,194.0,194.0,
193.0,193.0,193.0,
192.0,192.0,192.0,
191.0,191.0,191.0,
190.0,190.0,190.0,
189.0,189.0,189.0,
188.0,188.0,188.0,
187.00000000000003,187.00000000000003,187.00000000000003,
186.00000000000003,186.00000000000003,186.00000000000003,
185.0,185.0,185.0,
184.0,184.0,184.0,
183.0,183.0,183.0,
182.0,182.0,182.0,
181.0,181.0,181.0,
179.99999999999997,179.99999999999997,179.99999999999997,
179.0,179.0,179.0,
178.0,178.0,178.0,
177.0,177.0,177.0,
176.0,176.0,176.0,
175.0,175.0,175.0,
174.0,174.0,174.0,
173.0,173.0,173.0,
172.0,172.0,172.0,
171.00000000000003,171.00000000000003,171.00000000000003,
170.00000000000003,170.00000000000003,170.00000000000003,
169.0,169.0,169.0,
168.0,168.0,168.0,
167.0,167.0,167.0,
166.0,166.0,166.0,
165.0,165.0,165.0,
163.99999999999997,163.99999999999997,163.99999999999997,
163.0,163.0,163.0,
162.0,162.0,162.0,
161.0,161.0,161.0,
160.0,160.0,160.0,
159.0,159.0,159.0,
158.0,158.0,158.0,
157.0,157.0,157.0,
156.0,156.0,156.0,
155.00000000000003,155.00000000000003,155.00000000000003,
154.00000000000003,154.00000000000003,154.00000000000003,
153.0,153.0,153.0,
152.0,152.0,152.0,
151.0,151.0,151.0,
150.0,150.0,150.0,
149.0,149.0,149.0,
147.99999999999997,147.99999999999997,147.99999999999997,
147.0,147.0,147.0,
146.0,146.0,146.0,
145.0,145.0,145.0,
144.0,144.0,144.0,
143.0,143.0,143.0,
142.0,142.0,142.0,
141.0,141.0,141.0,
140.0,140.0,140.0,
139.00000000000003,139.00000000000003,139.00000000000003,
138.00000000000003,138.00000000000003,138.00000000000003,
137.0,137.0,137.0,
136.0,136.0,136.0,
135.0,135.0,135.0,
134.0,134.0,134.0,
133.0,133.0,133.0,
131.99999999999997,131.99999999999997,131.99999999999997,
131.0,131.0,131.0,
130.0,130.0,130.0,
129.0,129.0,129.0,
128.0,128.0,128.0,
127.0,127.0,127.0,
126.0,126.0,126.0,
125.00000000000001,125.00000000000001,125.00000000000001,
124.00000000000001,124.00000000000001,124.00000000000001,
123.00000000000001,123.00000000000001,123.00000000000001,
121.99999999999999,121.99999999999999,121.99999999999999,
121.0,121.0,121.0,
120.0,120.0,120.0,
119.0,119.0,119.0,
118.0,118.0,118.0,
117.00000000000001,117.00000000000001,117.00000000000001,
116.00000000000001,116.00000000000001,116.00000000000001,
114.99999999999999,114.99999999999999,114.99999999999999,
113.99999999999999,113.99999999999999,113.99999999999999,
113.0,113.0,113.0,
112.0,112.0,112.0,
111.0,111.0,111.0,
110.0,110.0,110.0,
109.00000000000001,109.00000000000001,109.00000000000001,
108.00000000000001,108.00000000000001,108.00000000000001,
107.00000000000001,107.00000000000001,107.00000000000001,
105.99999999999999,105.99999999999999,105.99999999999999,
105.0,105.0,105.0,
104.0,104.0,104.0,
103.0,103.0,103.0,
102.0,102.0,102.0,
101.00000000000001,101.00000000000001,101.00000000000001,
100.00000000000001,100.00000000000001,100.00000000000001,
98.99999999999999,98.99999999999999,98.99999999999999,
97.99999999999999,97.99999999999999,97.99999999999999,
97.0,97.0,97.0,
96.0,96.0,96.0,
95.0,95.0,95.0,
94.0,94.0,94.0,
93.00000000000001,93.00000000000001,93.00000000000001,
92.00000000000001,92.00000000000001,92.00000000000001,
91.00000000000001,91.00000000000001,91.00000000000001,
89.99999999999999,89.99999999999999,89.99999999999999,
89.0,89.0,89.0,
88.0,88.0,88.0,
87.0,87.0,87.0,
86.0,86.0,86.0,
85.00000000000001,85.00000000000001,85.00000000000001,
84.00000000000001,84.00000000000001,84.00000000000001,
82.99999999999999,82.99999999999999,82.99999999999999,
81.99999999999999,81.99999999999999,81.99999999999999,
81.0,81.0,81.0,
80.0,80.0,80.0,
79.0,79.0,79.0,
78.0,78.0,78.0,
77.00000000000001,77.00000000000001,77.00000000000001,
76.00000000000001,76.00000000000001,76.00000000000001,
75.00000000000001,75.00000000000001,75.00000000000001,
73.99999999999999,73.99999999999999,73.99999999999999,
73.0,73.0,73.0,
72.0,72.0,72.0,
71.0,71.0,71.0,
70.0,70.0,70.0,
69.00000000000001,69.00000000000001,69.00000000000001,
68.00000000000001,68.00000000000001,68.00000000000001,
66.99999999999999,66.99999999999999,66.99999999999999,
65.99999999999999,65.99999999999999,65.99999999999999,
65.0,65.0,65.0,
64.0,64.0,64.0,
63.0,63.0,63.0,
62.00000000000001,62.00000000000001,62.00000000000001,
61.00000000000001,61.00000000000001,61.00000000000001,
60.000000000000014,60.000000000000014,60.000000000000014,
59.000000000000014,59.000000000000014,59.000000000000014,
57.99999999999999,57.99999999999999,57.99999999999999,
56.99999999999999,56.99999999999999,56.99999999999999,
56.0,56.0,56.0,
55.0,55.0,55.0,
54.00000000000001,54.00000000000001,54.00000000000001,
53.00000000000001,53.00000000000001,53.00000000000001,
52.000000000000014,52.000000000000014,52.000000000000014,
50.999999999999986,50.999999999999986,50.999999999999986,
49.99999999999999,49.99999999999999,49.99999999999999,
48.99999999999999,48.99999999999999,48.99999999999999,
48.0,48.0,48.0,
47.0,47.0,47.0,
46.00000000000001,46.00000000000001,46.00000000000001,
45.00000000000001,45.00000000000001,45.00000000000001,
44.000000000000014,44.000000000000014,44.000000000000014,
43.000000000000014,43.000000000000014,43.000000000000014,
41.99999999999999,41.99999999999999,41.99999999999999,
40.99999999999999,40.99999999999999,40.99999999999999,
40.0,40.0,40.0,
39.0,39.0,39.0,
38.00000000000001,38.00000000000001,38.00000000000001,
37.00000000000001,37.00000000000001,37.00000000000001,
36.000000000000014,36.000000000000014,36.000000000000014,
34.999999999999986,34.999999999999986,34.999999999999986,
33.99999999999999,33.99999999999999,33.99999999999999,
32.99999999999999,32.99999999999999,32.99999999999999,
32.0,32.0,32.0,
31.000000000000004,31.000000000000004,31.000000000000004,
30.000000000000007,30.000000000000007,30.000000000000007,
29.00000000000001,29.00000000000001,29.00000000000001,
28.000000000000014,28.000000000000014,28.000000000000014,
27.000000000000018,27.000000000000018,27.000000000000018,
25.999999999999993,25.999999999999993,25.999999999999993,
24.999999999999996,24.999999999999996,24.999999999999996,
24.0,24.0,24.0,
23.000000000000004,23.000000000000004,23.000000000000004,
22.000000000000007,22.000000000000007,22.000000000000007,
21.00000000000001,21.00000000000001,21.00000000000001,
20.000000000000014,20.000000000000014,20.000000000000014,
18.99999999999999,18.99999999999999,18.99999999999999,
17.999999999999993,17.999999999999993,17.999999999999993,
16.999999999999996,16.999999999999996,16.999999999999996,
16.0,16.0,16.0,
15.000000000000004,15.000000000000004,15.000000000000004,
14.000000000000007,14.000000000000007,14.000000000000007,
13.00000000000001,13.00000000000001,13.00000000000001,
12.000000000000014,12.000000000000014,12.000000000000014,
11.000000000000018,11.000000000000018,11.000000000000018,
9.999999999999993,9.999999999999993,9.999999999999993,
8.999999999999996,8.999999999999996,8.999999999999996,
8.0,8.0,8.0,
7.0000000000000036,7.0000000000000036,7.0000000000000036,
6.000000000000007,6.000000000000007,6.000000000000007,
5.000000000000011,5.000000000000011,5.000000000000011,
4.000000000000014,4.000000000000014,4.000000000000014,
2.9999999999999893,2.9999999999999893,2.9999999999999893,
1.999999999999993,1.999999999999993,1.999999999999993,
0.9999999999999964,0.9999999999999964,0.9999999999999964,
0.0,0.0,0.0,
}

258
python/tests/reference/Colormap/binary.txt Normal file
View File

@ -0,0 +1,258 @@
# Creator: damask.Colormap v99.99.99-9999-pytest
1_RGBA 2_RGBA 3_RGBA 4_RGBA
1.0 1.0 1.0 1.0
0.996078431372549 0.996078431372549 0.996078431372549 1.0
0.9921568627450981 0.9921568627450981 0.9921568627450981 1.0
0.9882352941176471 0.9882352941176471 0.9882352941176471 1.0
0.9843137254901961 0.9843137254901961 0.9843137254901961 1.0
0.9803921568627451 0.9803921568627451 0.9803921568627451 1.0
0.9764705882352941 0.9764705882352941 0.9764705882352941 1.0
0.9725490196078431 0.9725490196078431 0.9725490196078431 1.0
0.9686274509803922 0.9686274509803922 0.9686274509803922 1.0
0.9647058823529412 0.9647058823529412 0.9647058823529412 1.0
0.9607843137254902 0.9607843137254902 0.9607843137254902 1.0
0.9568627450980393 0.9568627450980393 0.9568627450980393 1.0
0.9529411764705882 0.9529411764705882 0.9529411764705882 1.0
0.9490196078431372 0.9490196078431372 0.9490196078431372 1.0
0.9450980392156862 0.9450980392156862 0.9450980392156862 1.0
0.9411764705882353 0.9411764705882353 0.9411764705882353 1.0
0.9372549019607843 0.9372549019607843 0.9372549019607843 1.0
0.9333333333333333 0.9333333333333333 0.9333333333333333 1.0
0.9294117647058824 0.9294117647058824 0.9294117647058824 1.0
0.9254901960784314 0.9254901960784314 0.9254901960784314 1.0
0.9215686274509804 0.9215686274509804 0.9215686274509804 1.0
0.9176470588235294 0.9176470588235294 0.9176470588235294 1.0
0.9137254901960784 0.9137254901960784 0.9137254901960784 1.0
0.9098039215686274 0.9098039215686274 0.9098039215686274 1.0
0.9058823529411765 0.9058823529411765 0.9058823529411765 1.0
0.9019607843137255 0.9019607843137255 0.9019607843137255 1.0
0.8980392156862745 0.8980392156862745 0.8980392156862745 1.0
0.8941176470588236 0.8941176470588236 0.8941176470588236 1.0
0.8901960784313725 0.8901960784313725 0.8901960784313725 1.0
0.8862745098039215 0.8862745098039215 0.8862745098039215 1.0
0.8823529411764706 0.8823529411764706 0.8823529411764706 1.0
0.8784313725490196 0.8784313725490196 0.8784313725490196 1.0
0.8745098039215686 0.8745098039215686 0.8745098039215686 1.0
0.8705882352941177 0.8705882352941177 0.8705882352941177 1.0
0.8666666666666667 0.8666666666666667 0.8666666666666667 1.0
0.8627450980392157 0.8627450980392157 0.8627450980392157 1.0
0.8588235294117648 0.8588235294117648 0.8588235294117648 1.0
0.8549019607843138 0.8549019607843138 0.8549019607843138 1.0
0.8509803921568627 0.8509803921568627 0.8509803921568627 1.0
0.8470588235294118 0.8470588235294118 0.8470588235294118 1.0
0.8431372549019608 0.8431372549019608 0.8431372549019608 1.0
0.8392156862745098 0.8392156862745098 0.8392156862745098 1.0
0.8352941176470589 0.8352941176470589 0.8352941176470589 1.0
0.8313725490196078 0.8313725490196078 0.8313725490196078 1.0
0.8274509803921568 0.8274509803921568 0.8274509803921568 1.0
0.8235294117647058 0.8235294117647058 0.8235294117647058 1.0
0.8196078431372549 0.8196078431372549 0.8196078431372549 1.0
0.8156862745098039 0.8156862745098039 0.8156862745098039 1.0
0.8117647058823529 0.8117647058823529 0.8117647058823529 1.0
0.807843137254902 0.807843137254902 0.807843137254902 1.0
0.803921568627451 0.803921568627451 0.803921568627451 1.0
0.8 0.8 0.8 1.0
0.7960784313725491 0.7960784313725491 0.7960784313725491 1.0
0.7921568627450981 0.7921568627450981 0.7921568627450981 1.0
0.788235294117647 0.788235294117647 0.788235294117647 1.0
0.7843137254901961 0.7843137254901961 0.7843137254901961 1.0
0.7803921568627451 0.7803921568627451 0.7803921568627451 1.0
0.7764705882352941 0.7764705882352941 0.7764705882352941 1.0
0.7725490196078432 0.7725490196078432 0.7725490196078432 1.0
0.7686274509803921 0.7686274509803921 0.7686274509803921 1.0
0.7647058823529411 0.7647058823529411 0.7647058823529411 1.0
0.7607843137254902 0.7607843137254902 0.7607843137254902 1.0
0.7568627450980392 0.7568627450980392 0.7568627450980392 1.0
0.7529411764705882 0.7529411764705882 0.7529411764705882 1.0
0.7490196078431373 0.7490196078431373 0.7490196078431373 1.0
0.7450980392156863 0.7450980392156863 0.7450980392156863 1.0
0.7411764705882353 0.7411764705882353 0.7411764705882353 1.0
0.7372549019607844 0.7372549019607844 0.7372549019607844 1.0
0.7333333333333334 0.7333333333333334 0.7333333333333334 1.0
0.7294117647058824 0.7294117647058824 0.7294117647058824 1.0
0.7254901960784313 0.7254901960784313 0.7254901960784313 1.0
0.7215686274509804 0.7215686274509804 0.7215686274509804 1.0
0.7176470588235294 0.7176470588235294 0.7176470588235294 1.0
0.7137254901960784 0.7137254901960784 0.7137254901960784 1.0
0.7098039215686275 0.7098039215686275 0.7098039215686275 1.0
0.7058823529411764 0.7058823529411764 0.7058823529411764 1.0
0.7019607843137254 0.7019607843137254 0.7019607843137254 1.0
0.6980392156862745 0.6980392156862745 0.6980392156862745 1.0
0.6941176470588235 0.6941176470588235 0.6941176470588235 1.0
0.6901960784313725 0.6901960784313725 0.6901960784313725 1.0
0.6862745098039216 0.6862745098039216 0.6862745098039216 1.0
0.6823529411764706 0.6823529411764706 0.6823529411764706 1.0
0.6784313725490196 0.6784313725490196 0.6784313725490196 1.0
0.6745098039215687 0.6745098039215687 0.6745098039215687 1.0
0.6705882352941177 0.6705882352941177 0.6705882352941177 1.0
0.6666666666666667 0.6666666666666667 0.6666666666666667 1.0
0.6627450980392157 0.6627450980392157 0.6627450980392157 1.0
0.6588235294117647 0.6588235294117647 0.6588235294117647 1.0
0.6549019607843137 0.6549019607843137 0.6549019607843137 1.0
0.6509803921568628 0.6509803921568628 0.6509803921568628 1.0
0.6470588235294118 0.6470588235294118 0.6470588235294118 1.0
0.6431372549019607 0.6431372549019607 0.6431372549019607 1.0
0.6392156862745098 0.6392156862745098 0.6392156862745098 1.0
0.6352941176470588 0.6352941176470588 0.6352941176470588 1.0
0.6313725490196078 0.6313725490196078 0.6313725490196078 1.0
0.6274509803921569 0.6274509803921569 0.6274509803921569 1.0
0.6235294117647059 0.6235294117647059 0.6235294117647059 1.0
0.6196078431372549 0.6196078431372549 0.6196078431372549 1.0
0.615686274509804 0.615686274509804 0.615686274509804 1.0
0.611764705882353 0.611764705882353 0.611764705882353 1.0
0.607843137254902 0.607843137254902 0.607843137254902 1.0
0.603921568627451 0.603921568627451 0.603921568627451 1.0
0.6 0.6 0.6 1.0
0.596078431372549 0.596078431372549 0.596078431372549 1.0
0.592156862745098 0.592156862745098 0.592156862745098 1.0
0.5882352941176471 0.5882352941176471 0.5882352941176471 1.0
0.5843137254901961 0.5843137254901961 0.5843137254901961 1.0
0.580392156862745 0.580392156862745 0.580392156862745 1.0
0.5764705882352941 0.5764705882352941 0.5764705882352941 1.0
0.5725490196078431 0.5725490196078431 0.5725490196078431 1.0
0.5686274509803921 0.5686274509803921 0.5686274509803921 1.0
0.5647058823529412 0.5647058823529412 0.5647058823529412 1.0
0.5607843137254902 0.5607843137254902 0.5607843137254902 1.0
0.5568627450980392 0.5568627450980392 0.5568627450980392 1.0
0.5529411764705883 0.5529411764705883 0.5529411764705883 1.0
0.5490196078431373 0.5490196078431373 0.5490196078431373 1.0
0.5450980392156863 0.5450980392156863 0.5450980392156863 1.0
0.5411764705882354 0.5411764705882354 0.5411764705882354 1.0
0.5372549019607843 0.5372549019607843 0.5372549019607843 1.0
0.5333333333333333 0.5333333333333333 0.5333333333333333 1.0
0.5294117647058824 0.5294117647058824 0.5294117647058824 1.0
0.5254901960784314 0.5254901960784314 0.5254901960784314 1.0
0.5215686274509804 0.5215686274509804 0.5215686274509804 1.0
0.5176470588235293 0.5176470588235293 0.5176470588235293 1.0
0.5137254901960784 0.5137254901960784 0.5137254901960784 1.0
0.5098039215686274 0.5098039215686274 0.5098039215686274 1.0
0.5058823529411764 0.5058823529411764 0.5058823529411764 1.0
0.5019607843137255 0.5019607843137255 0.5019607843137255 1.0
0.4980392156862745 0.4980392156862745 0.4980392156862745 1.0
0.49411764705882355 0.49411764705882355 0.49411764705882355 1.0
0.4901960784313726 0.4901960784313726 0.4901960784313726 1.0
0.4862745098039216 0.4862745098039216 0.4862745098039216 1.0
0.48235294117647065 0.48235294117647065 0.48235294117647065 1.0
0.4784313725490196 0.4784313725490196 0.4784313725490196 1.0
0.4745098039215686 0.4745098039215686 0.4745098039215686 1.0
0.47058823529411764 0.47058823529411764 0.47058823529411764 1.0
0.4666666666666667 0.4666666666666667 0.4666666666666667 1.0
0.4627450980392157 0.4627450980392157 0.4627450980392157 1.0
0.45882352941176474 0.45882352941176474 0.45882352941176474 1.0
0.4549019607843138 0.4549019607843138 0.4549019607843138 1.0
0.4509803921568627 0.4509803921568627 0.4509803921568627 1.0
0.44705882352941173 0.44705882352941173 0.44705882352941173 1.0
0.44313725490196076 0.44313725490196076 0.44313725490196076 1.0
0.4392156862745098 0.4392156862745098 0.4392156862745098 1.0
0.43529411764705883 0.43529411764705883 0.43529411764705883 1.0
0.43137254901960786 0.43137254901960786 0.43137254901960786 1.0
0.4274509803921569 0.4274509803921569 0.4274509803921569 1.0
0.42352941176470593 0.42352941176470593 0.42352941176470593 1.0
0.41960784313725497 0.41960784313725497 0.41960784313725497 1.0
0.4156862745098039 0.4156862745098039 0.4156862745098039 1.0
0.4117647058823529 0.4117647058823529 0.4117647058823529 1.0
0.40784313725490196 0.40784313725490196 0.40784313725490196 1.0
0.403921568627451 0.403921568627451 0.403921568627451 1.0
0.4 0.4 0.4 1.0
0.39607843137254906 0.39607843137254906 0.39607843137254906 1.0
0.3921568627450981 0.3921568627450981 0.3921568627450981 1.0
0.388235294117647 0.388235294117647 0.388235294117647 1.0
0.38431372549019605 0.38431372549019605 0.38431372549019605 1.0
0.3803921568627451 0.3803921568627451 0.3803921568627451 1.0
0.3764705882352941 0.3764705882352941 0.3764705882352941 1.0
0.37254901960784315 0.37254901960784315 0.37254901960784315 1.0
0.3686274509803922 0.3686274509803922 0.3686274509803922 1.0
0.3647058823529412 0.3647058823529412 0.3647058823529412 1.0
0.36078431372549025 0.36078431372549025 0.36078431372549025 1.0
0.3568627450980393 0.3568627450980393 0.3568627450980393 1.0
0.3529411764705882 0.3529411764705882 0.3529411764705882 1.0
0.34901960784313724 0.34901960784313724 0.34901960784313724 1.0
0.34509803921568627 0.34509803921568627 0.34509803921568627 1.0
0.3411764705882353 0.3411764705882353 0.3411764705882353 1.0
0.33725490196078434 0.33725490196078434 0.33725490196078434 1.0
0.33333333333333337 0.33333333333333337 0.33333333333333337 1.0
0.3294117647058824 0.3294117647058824 0.3294117647058824 1.0
0.3254901960784313 0.3254901960784313 0.3254901960784313 1.0
0.32156862745098036 0.32156862745098036 0.32156862745098036 1.0
0.3176470588235294 0.3176470588235294 0.3176470588235294 1.0
0.3137254901960784 0.3137254901960784 0.3137254901960784 1.0
0.30980392156862746 0.30980392156862746 0.30980392156862746 1.0
0.3058823529411765 0.3058823529411765 0.3058823529411765 1.0
0.3019607843137255 0.3019607843137255 0.3019607843137255 1.0
0.29803921568627456 0.29803921568627456 0.29803921568627456 1.0
0.2941176470588236 0.2941176470588236 0.2941176470588236 1.0
0.2901960784313725 0.2901960784313725 0.2901960784313725 1.0
0.28627450980392155 0.28627450980392155 0.28627450980392155 1.0
0.2823529411764706 0.2823529411764706 0.2823529411764706 1.0
0.2784313725490196 0.2784313725490196 0.2784313725490196 1.0
0.27450980392156865 0.27450980392156865 0.27450980392156865 1.0
0.2705882352941177 0.2705882352941177 0.2705882352941177 1.0
0.2666666666666667 0.2666666666666667 0.2666666666666667 1.0
0.26274509803921564 0.26274509803921564 0.26274509803921564 1.0
0.2588235294117647 0.2588235294117647 0.2588235294117647 1.0
0.2549019607843137 0.2549019607843137 0.2549019607843137 1.0
0.25098039215686274 0.25098039215686274 0.25098039215686274 1.0
0.24705882352941178 0.24705882352941178 0.24705882352941178 1.0
0.2431372549019608 0.2431372549019608 0.2431372549019608 1.0
0.23921568627450984 0.23921568627450984 0.23921568627450984 1.0
0.23529411764705888 0.23529411764705888 0.23529411764705888 1.0
0.2313725490196079 0.2313725490196079 0.2313725490196079 1.0
0.22745098039215683 0.22745098039215683 0.22745098039215683 1.0
0.22352941176470587 0.22352941176470587 0.22352941176470587 1.0
0.2196078431372549 0.2196078431372549 0.2196078431372549 1.0
0.21568627450980393 0.21568627450980393 0.21568627450980393 1.0
0.21176470588235297 0.21176470588235297 0.21176470588235297 1.0
0.207843137254902 0.207843137254902 0.207843137254902 1.0
0.20392156862745103 0.20392156862745103 0.20392156862745103 1.0
0.19999999999999996 0.19999999999999996 0.19999999999999996 1.0
0.196078431372549 0.196078431372549 0.196078431372549 1.0
0.19215686274509802 0.19215686274509802 0.19215686274509802 1.0
0.18823529411764706 0.18823529411764706 0.18823529411764706 1.0
0.1843137254901961 0.1843137254901961 0.1843137254901961 1.0
0.18039215686274512 0.18039215686274512 0.18039215686274512 1.0
0.17647058823529416 0.17647058823529416 0.17647058823529416 1.0
0.1725490196078432 0.1725490196078432 0.1725490196078432 1.0
0.16862745098039222 0.16862745098039222 0.16862745098039222 1.0
0.16470588235294115 0.16470588235294115 0.16470588235294115 1.0
0.16078431372549018 0.16078431372549018 0.16078431372549018 1.0
0.1568627450980392 0.1568627450980392 0.1568627450980392 1.0
0.15294117647058825 0.15294117647058825 0.15294117647058825 1.0
0.14901960784313728 0.14901960784313728 0.14901960784313728 1.0
0.14509803921568631 0.14509803921568631 0.14509803921568631 1.0
0.14117647058823535 0.14117647058823535 0.14117647058823535 1.0
0.13725490196078427 0.13725490196078427 0.13725490196078427 1.0
0.1333333333333333 0.1333333333333333 0.1333333333333333 1.0
0.12941176470588234 0.12941176470588234 0.12941176470588234 1.0
0.12549019607843137 0.12549019607843137 0.12549019607843137 1.0
0.1215686274509804 0.1215686274509804 0.1215686274509804 1.0
0.11764705882352944 0.11764705882352944 0.11764705882352944 1.0
0.11372549019607847 0.11372549019607847 0.11372549019607847 1.0
0.1098039215686275 0.1098039215686275 0.1098039215686275 1.0
0.10588235294117654 0.10588235294117654 0.10588235294117654 1.0
0.10196078431372546 0.10196078431372546 0.10196078431372546 1.0
0.0980392156862745 0.0980392156862745 0.0980392156862745 1.0
0.09411764705882353 0.09411764705882353 0.09411764705882353 1.0
0.09019607843137256 0.09019607843137256 0.09019607843137256 1.0
0.0862745098039216 0.0862745098039216 0.0862745098039216 1.0
0.08235294117647063 0.08235294117647063 0.08235294117647063 1.0
0.07843137254901966 0.07843137254901966 0.07843137254901966 1.0
0.07450980392156858 0.07450980392156858 0.07450980392156858 1.0
0.07058823529411762 0.07058823529411762 0.07058823529411762 1.0
0.06666666666666665 0.06666666666666665 0.06666666666666665 1.0
0.06274509803921569 0.06274509803921569 0.06274509803921569 1.0
0.05882352941176472 0.05882352941176472 0.05882352941176472 1.0
0.05490196078431375 0.05490196078431375 0.05490196078431375 1.0
0.050980392156862786 0.050980392156862786 0.050980392156862786 1.0
0.04705882352941182 0.04705882352941182 0.04705882352941182 1.0
0.04313725490196085 0.04313725490196085 0.04313725490196085 1.0
0.039215686274509776 0.039215686274509776 0.039215686274509776 1.0
0.03529411764705881 0.03529411764705881 0.03529411764705881 1.0
0.03137254901960784 0.03137254901960784 0.03137254901960784 1.0
0.027450980392156876 0.027450980392156876 0.027450980392156876 1.0
0.02352941176470591 0.02352941176470591 0.02352941176470591 1.0
0.019607843137254943 0.019607843137254943 0.019607843137254943 1.0
0.015686274509803977 0.015686274509803977 0.015686274509803977 1.0
0.0117647058823529 0.0117647058823529 0.0117647058823529 1.0
0.007843137254901933 0.007843137254901933 0.007843137254901933 1.0
0.0039215686274509665 0.0039215686274509665 0.0039215686274509665 1.0
0.0 0.0 0.0 1.0

6
python/tests/reference/VTK/polyData.vtp
View File

@ -3,11 +3,11 @@
<PolyData> <PolyData>
<Piece NumberOfPoints="10" NumberOfVerts="0" NumberOfLines="0" NumberOfStrips="0" NumberOfPolys="0"> <Piece NumberOfPoints="10" NumberOfVerts="0" NumberOfLines="0" NumberOfStrips="0" NumberOfPolys="0">
<PointData> <PointData>
<DataArray type="Float64" Name="coordinates" NumberOfComponents="3" format="binary" RangeMin="0.7453559924999299" RangeMax="2.449489742783178"> <DataArray type="Float32" Name="coordinates" NumberOfComponents="3" format="binary" RangeMin="0.7453560147132696" RangeMax="2.449489742783178">
AQAAAACAAADwAAAAagAAAA==eF5jYMAGPuyXOV4IRHvsIfQZe8u+xxZ9j1/sh/Eh9B37IjDj4f5QMLhqD6Gf2odB+Pth6iD0G3sFiLn7ofqg+j/CxOH6IfRX+xCouRYQ+6H0D/sCqH6YuRD6D1wd1B9QmsEBxgcAJsNfhw== AQAAAACAAAB4AAAAVgAAAA==eF5jYICBhv2WfY9tLfuS7Ypk3PeDaCDf7okF3/7Vq1bZrV6lZQ+k94HEgHL2QHovUM7+iUUfiG0LlQdhkH77Ipnj9iB5qFp7kBjQDiBmcADRANsaLXM=
<InformationKey name="L2_NORM_RANGE" location="vtkDataArray" length="2"> <InformationKey name="L2_NORM_RANGE" location="vtkDataArray" length="2">
<Value index="0"> <Value index="0">
0.7453559925 0.74535601471
</Value> </Value>
<Value index="1"> <Value index="1">
2.4494897428 2.4494897428

14
python/tests/reference/VTK/rectilinearGrid.vtr
View File

@ -3,27 +3,27 @@
<RectilinearGrid WholeExtent="0 5 0 6 0 7"> <RectilinearGrid WholeExtent="0 5 0 6 0 7">
<Piece Extent="0 5 0 6 0 7"> <Piece Extent="0 5 0 6 0 7">
<PointData> <PointData>
<DataArray type="Float64" Name="node" NumberOfComponents="3" format="binary" RangeMin="0" RangeMax="1.268857754044952"> <DataArray type="Float32" Name="node" NumberOfComponents="3" format="binary" RangeMin="0" RangeMax="1.268857765318962">
AQAAAACAAACAHwAAywMAAA==eF51mTGKVUEQRX88izByB8byY0MDlyG4CTfgFswMDWQigw8/Ejr40E7woGkYXqLJW4LDb8qmz71VDDjvWMGhijvd0KeTr8dXn/++f/x59rwIf3j6+untw1PS34S/udez8KgP97r+///w8bwIDx/f34SHD3nU4DXxIS/CVx/2N+GrTxWfUV18PC/C132xvwlf9zV51PDck/mQF+HrfNjfhK/z2cXn273+iI/nRXj4+P4mPHzI1zqSfZEX4eu+2N+Er/uanPXl9buXn+9n5n3lRTjzvvY34cx78Pgee7ye6eN5Ec6804ecefc+NfEhL8KZd+9TE58qPqO6+HhehDPv9CFn3v189mQ+5EU48+7ns4sPefhE7ujjeRHOvNOHnHmnz6gj2Rd5Ec6804ecefc+kbtLkvcLfCb3eb/AZ3Kf90uS9+njOfN+SfI+fch93ulTEx9y5p0+7Gfe6VPFZ1QXH8+Zd+6L/cw799XFZ3ju4uM58875sJ9553x28VlzN308Z94vSd6nD7nPO/d1iI/nzDv3xX7mnfs6Ep/Tafvx8eXnl+R95UU48772N+HMe/D4Hnu8nunjeRHOvNOHnHn3PjXxIS/CmXfvUxOfKj6juvh4XoQz7/QhZ979fPZkPuRFOPPu57OLD3n4RO7o43kRzrzTh5x5p8+oI9kXeRHOvNOHnHn3PnHO3iTvK+f5fkvO9xt8Jmfeg8f32GOcs9PHc57vN8k7fciZd+9TEx9ynu/0YT/Pd/pU8RnVxcdznu/cF/t5vnNfXXyG5y4+nvN853zYz/Od89nFZz1np4/nPN9vknf6kDPv9Bl1iI/nPN+5L/bzfOe+jsTndLr/Gdh+S95XXoQz72t/E868B4/vscfrmT6eF+HMO33ImXfvUxMf8iKcefc+NfGp4jOqi4/nRTjzTh9y5t3PZ0/mQ16EM+9+Prv4kIdP5I4+nhfhzDt9yJl3+ow6kn2RF+HMO33ImXfvE/fqTfK+ct7nt+Q+v8FncuY9eHyPPca9evp4zvv8JnmnDznz7n1q4kPO+zx92M/7PH2q+Izq4uM57/PcF/t5n+e+uvgMz118POd9nvNhP+/znM8uPpE7+njO+/wmeacPOfNOn1GH+HjO+zz3xX7e57mvI/GJ6pL3lfM9rkve1/4mnHkPHr+NPca72PTxnO9xXfK+9jfhzLv3qYkPOd/j6MP+Jpx5p8/6zX2R8z2O+2J/E868r//yPY7zIed7HOfD/iaceadP5I4+nvM9rkve6UPOvNNn1CE+nvM9jvtifxPOvAf/B5cwXsg= AQAAAACAAADADwAAFQMAAA==eF5Vl7GtFEEQBTeG8wgAhxxux8bEOIsYwDu/M0A6D0JofMAmjU0BF5+d7anXjzP+L6NV4l3pj8S29efL77/35ucO//nwS3zeiL99fTM2+3zPd/u6uTc/d3h67EY8PXB50jxpnjRPmifNk/Kcn4Gn+do18NiNeO3StvPfJk/ztUseuxGvXfI8Hg95mp87PD12I54euD5hu0IeuHbpRly74r9mb9+/7nQvru6T6b5uxHSfPH/PdniaqzseuxHTvT1pnjRPmifNk+ZJec7PsF3Ddg3bxY2Y7rZLnubqbrvkgemOZ7bD01zd8diNmO69K2xXyAPTvXeF7QrzzHa3vbtPpvtt7+7Xjbi73/b1/ex4muleHrsRd3c8aZ40T5onzZPmSXm2q512Dds1bBc34u6uXfI001275IG7e3mqXXma6V4euxF3d3aF7Qp54O7OrrBdYZ5t+/npo7oXV/fJdF83YrpPXt/Pjqe5uuOxGzHd25PmSfOkedI8aZ6U5/wM2zVs17Bd3Ijpbrvkaa7utksemO54Zjs8zdUdj92I6d67wnaFPDDde1fYrjDP9Vare7Heeft7f6n7ZHvn7e/9pe54YHvn1R0PXJ40T5onzZPmSfOkPFu91epuu4bt4kZs77z9vWuXPLC98+puu+RZb7W644HtnVd3PHDNCtsV8sD2zqt77wrzbNvn44e6F1f3yXRfN2K6T17fz46nubrjsRsx3duT5knzpHnSPGmelOf8DNs1bNewXdyI6W675Gmu7rZLHpjueGY7PM3VHY/diOneu8J2hTww3XtX2K4wz3yrD3Uv5p0/7J0/1H1yv/OHuuNp5p0/7J0/1B0PXJ40T5onzZPmSfOkPNv1VmvXsF3DdnEj7ndeu+Rp5p3XLnngfucPdcfTzDt/2Dt/qDseuGaF7Qp54H7n2RW2K8xT3xHdi5/67ui+bsR0nzx/zHZ4mqs7HrsR0709aZ40T5onzZPmSXnWb3YNPDDd8cB0t13yNFd32yUPTHc8sx2eZv0/btAdD0z33hW2K+SB6d67wnYV/wMyiXC6
<InformationKey name="L2_NORM_RANGE" location="vtkDataArray" length="2"> <InformationKey name="L2_NORM_RANGE" location="vtkDataArray" length="2">
<Value index="0"> <Value index="0">
0 0
</Value> </Value>
<Value index="1"> <Value index="1">
1.268857754 1.2688577653
</Value> </Value>
</InformationKey> </InformationKey>
</DataArray> </DataArray>
</PointData> </PointData>
<CellData> <CellData>
<DataArray type="Float64" Name="cell" NumberOfComponents="3" format="binary" RangeMin="0.10871961482881586" RangeMax="1.160792402743735"> <DataArray type="Float32" Name="cell" NumberOfComponents="3" format="binary" RangeMin="0.10871961651677305" RangeMax="1.1607924233069467">
AQAAAACAAACwEwAAZQIAAA==eF511r2KFEEYRmHjjY2N9Ao0lQ79yQy8jBVjc29gL0Ezhc1maDBfU0FB3NVgNyqQuged1leZ56sqBpo5HRx4+wS13nn989l6vjzfzm45u/vk9+/NcvL17cuHJx8Lf3D/cD4XfvPq9vmj68vCH19vLwpf/3pvbedT8crjlccrj1ce77vtXBavPF55vPJ45fG+3/hN8crjlccrj1d+vHOb7NyKV368cyte+XrUVZ901YtXftxVL175Ss/f96dX+9MPpedwew6353B7DrdnvXJ71iu3Z73pTa/cnvXK7VlvetMrt2e9cnse79wmO7fildvzeOdWvOlt3FUvXrk9j7vqE+9uWQ/46qL0HG7P4fYcbs/h9qxXbs965fasN73plduzXrk9601veuX2rFduz+Od22TnVrxyex7v3Io3vY276sUrt+dxV33i3f37/vYcbs/h9hxuz+H2rFduz3rl9pynPYfbs165PeuV27NeuT3rlduz3j//22TnVrxyew6353B7DrdnvXJ71iu353uHa8jZl9JzuD2H23O4PYfbs165PeuV27Pe9KZXbs965fasN73plduzXrk9j3duk51b8crtebxzK970Nu6qF6/cnsdd9Yl3tzzdLtbfSs/h9hxuz+H2HG7PeuX2rFduz3rTm165PeuV27Pe9KZXbs965fY83rlNdm7FK7fn8c6teNPbuKtevHJ7HnfVJ97d8uJwDdn/KD2H23O4PYfbc7g965Xbs165PetNb3rl9qxXbs9605teuT3rldvzeOc22bkVr9yexzu34k1v46568crtedzVf/4LDINsdg== AQAAAACAAADYCQAA8QEAAA==eF5N0CGOFlEQAOHVK0nmBnAFgnn/PLsSsZd4CSQIPDfYZN2iUbMePI4DEMR/BbgDNJ2ve1yZSir18P3jeD6O8eruxfj99s0Ff356Kh63v4o/jNsdP/xzb24+Xbg4XBwuDheHe3//s1wcLg4Xh4vT3fZ2k9NNTjc53eRsnuXibJ7l4mye5T4fq1ycr1a5OF+tk3uMb++u9TnY52Cfg30O9pmLfeZin7nxjYt95mKf2932dpN9bjfZ526e5WKfu3mWi33uV6tc7HO/Wif3GO+vry8+B/sc7HOwz8E+c7HPXOwzN75xsc9c7HO7295uss/tJvvczbNc7HM3z3Kxz/1qlYt97lfr5B7/f/kc7HOwz8E+B/vMxT5zsc/c+MbFPnOxz+1ue7vJPreb7HM3z3Kxz908y8U+96tVLva5X62Te4y7xy/1OdjnYJ+DfQ72mYt95mKfufGNi33mYp/b3fZ2k31uN9nnbp7lYp+7eZaLfe5Xq1zsc79aJ/cYjy9/1Odgn4N9DvY52Gcu9pmLfebGNy72mYt9bnfb2032ud1kn7t5lot97uZZLva5X61ysc/9ap3cY1y//qnPwT4H+xzsc7DPXOwzF/vMjW9c7DMX+9zutreb7HO7yT538ywX+9zNs1zsc79a5WKf+1XyX6K10A4=
<InformationKey name="L2_NORM_RANGE" location="vtkDataArray" length="2"> <InformationKey name="L2_NORM_RANGE" location="vtkDataArray" length="2">
<Value index="0"> <Value index="0">
0.10871961483 0.10871961652
</Value> </Value>
<Value index="1"> <Value index="1">
1.1607924027 1.1607924233
</Value> </Value>
</InformationKey> </InformationKey>
</DataArray> </DataArray>

151
python/tests/test_Colormap.py Normal file
View File

@ -0,0 +1,151 @@
import os
import filecmp
import time
import numpy as np
import pytest
import damask
from damask import Colormap
@pytest.fixture
def reference_dir(reference_dir_base):
"""Directory containing reference results."""
return reference_dir_base/'Colormap'
class TestColormap:
def test_conversion(self):
specials = np.array([[0.,0.,0.],
[1.,0.,0.],
[0.,1.,0.],
[0.,0.,1.],
[1.,1.,0.],
[0.,1.,1.],
[1.,0.,1.],
[1.,1.,1.]
])
rgbs = np.vstack((specials,np.random.rand(100,3)))
pass # class not integrated
for rgb in rgbs:
print('rgb',rgb)
# rgb2hsv2rgb
hsv = Colormap._rgb2hsv(rgb)
print('hsv',hsv)
assert np.allclose(Colormap._hsv2rgb(hsv),rgb)
# rgb2hsl2rgb
hsl = Colormap._rgb2hsl(rgb)
print('hsl',hsl)
assert np.allclose(Colormap._hsl2rgb(hsl),rgb)
# rgb2xyz2rgb
xyz = Colormap._rgb2xyz(rgb)
print('xyz',xyz)
assert np.allclose(Colormap._xyz2rgb(xyz),rgb,atol=1.e-6,rtol=0)
# xyz2lab2xyz
lab = Colormap._xyz2lab(xyz)
print('lab',lab)
assert np.allclose(Colormap._lab2xyz(lab),xyz)
# lab2msh2lab
msh = Colormap._lab2msh(lab)
print('msh',msh)
assert np.allclose(Colormap._msh2lab(msh),lab)
# lab2rgb2lab
assert np.allclose(Colormap._rgb2lab(Colormap._lab2rgb(lab)),lab,atol=1.e-6,rtol=0)
# rgb2msh2rgb
assert np.allclose(Colormap._msh2rgb(Colormap._rgb2msh(rgb)),rgb,atol=1.e-6,rtol=0)
# hsv2msh
assert np.allclose(Colormap._hsv2msh(hsv),msh,atol=1.e-6,rtol=0)
# hsl2msh
assert np.allclose(Colormap._hsv2msh(hsv),msh,atol=1.e-6,rtol=0)
# xyz2msh
assert np.allclose(Colormap._xyz2msh(xyz),msh,atol=1.e-6,rtol=0)
@pytest.mark.parametrize('format',['ASCII','paraview','GOM','Gmsh'])
@pytest.mark.parametrize('model',['rgb','hsv','hsl','xyz','lab','msh'])
def test_from_range(self,model,format,tmpdir):
N = np.random.randint(2,256)
c = Colormap.from_range(np.random.rand(3),np.random.rand(3),model=model,N=N)
c.to_file(tmpdir/'color_out',format=format)
@pytest.mark.parametrize('format',['ASCII','paraview','GOM','Gmsh'])
@pytest.mark.parametrize('name',['strain','gnuplot','Greys','PRGn','viridis'])
def test_from_predefined(self,name,format,tmpdir):
N = np.random.randint(2,256)
c = Colormap.from_predefined(name,N)
os.chdir(tmpdir)
c.to_file(format=format)
@pytest.mark.parametrize('format,name',[('ASCII','test.txt'),
('paraview','test.json'),
('GOM','test.legend'),
('Gmsh','test.msh')
])
def test_write_filehandle(self,format,name,tmpdir):
c = Colormap.from_predefined('Dark2')
fname = tmpdir/name
with open(fname,'w') as f:
c.to_file(f,format=format)
for i in range(10):
if fname.exists(): return
time.sleep(.5)
assert False
def test_write_invalid_format(self):
c = Colormap.from_predefined('Dark2')
with pytest.raises(ValueError):
c.to_file(format='invalid')
@pytest.mark.parametrize('model',['rgb','hsv','hsl','lab','invalid'])
def test_invalid_color(self,model):
with pytest.raises(ValueError):
c = Colormap.from_range(-2.+np.random.rand(3),np.random.rand(3),N=10,model=model) # noqa
def test_reversed(self):
c_1 = Colormap.from_predefined('stress')
c_2 = c_1.reversed()
assert (not np.allclose(c_1.colors,c_2.colors)) and \
np.allclose(c_1.colors,c_2.reversed().colors)
def test_invert(self):
c_1 = Colormap.from_predefined('strain')
c_2 = ~c_1
assert (not np.allclose(c_1.colors,c_2.colors)) and \
np.allclose(c_1.colors,(~c_2).colors)
def test_add(self):
c = Colormap.from_predefined('jet')
c += c
assert (np.allclose(c.colors[:len(c.colors)//2],c.colors[len(c.colors)//2:]))
def test_list(self):
Colormap.list_predefined()
@pytest.mark.parametrize('format,ext',[('ASCII','.txt'),
('paraview','.json'),
('GOM','.legend'),
('Gmsh','.msh')
])
def test_compare_reference(self,format,ext,tmpdir,reference_dir,update,monkeypatch):
monkeypatch.setattr(damask, 'version', pytest.dummy_version)
name = 'binary'
c = Colormap.from_predefined(name)
if update:
os.chdir(reference_dir)
c.to_file(format=format)
else:
os.chdir(tmpdir)
c.to_file(format=format)
time.sleep(.5)
assert filecmp.cmp(tmpdir/(name+ext),reference_dir/(name+ext))

140
python/tests/test_Lattice.py
View File

@ -3,38 +3,154 @@ import random
import pytest import pytest
import numpy as np import numpy as np
from damask import Rotation
from damask import Symmetry from damask import Symmetry
def in_FZ(system,rho):
"""Non-vectorized version of 'in_FZ'."""
rho_abs = abs(rho)
if system == 'cubic':
return np.sqrt(2.0)-1.0 >= rho_abs[0] \
and np.sqrt(2.0)-1.0 >= rho_abs[1] \
and np.sqrt(2.0)-1.0 >= rho_abs[2] \
and 1.0 >= rho_abs[0] + rho_abs[1] + rho_abs[2]
elif system == 'hexagonal':
return 1.0 >= rho_abs[0] and 1.0 >= rho_abs[1] and 1.0 >= rho_abs[2] \
and 2.0 >= np.sqrt(3)*rho_abs[0] + rho_abs[1] \
and 2.0 >= np.sqrt(3)*rho_abs[1] + rho_abs[0] \
and 2.0 >= np.sqrt(3) + rho_abs[2]
elif system == 'tetragonal':
return 1.0 >= rho_abs[0] and 1.0 >= rho_abs[1] \
and np.sqrt(2.0) >= rho_abs[0] + rho_abs[1] \
and np.sqrt(2.0) >= rho_abs[2] + 1.0
elif system == 'orthorhombic':
return 1.0 >= rho_abs[0] and 1.0 >= rho_abs[1] and 1.0 >= rho_abs[2]
else:
return np.all(np.isfinite(rho_abs))
def in_disorientation_SST(system,rho):
"""Non-vectorized version of 'in_Disorientation_SST'."""
epsilon = 0.0
if system == 'cubic':
return rho[0] >= rho[1]+epsilon and rho[1] >= rho[2]+epsilon and rho[2] >= epsilon
elif system == 'hexagonal':
return rho[0] >= np.sqrt(3)*(rho[1]-epsilon) and rho[1] >= epsilon and rho[2] >= epsilon
elif system == 'tetragonal':
return rho[0] >= rho[1]-epsilon and rho[1] >= epsilon and rho[2] >= epsilon
elif system == 'orthorhombic':
return rho[0] >= epsilon and rho[1] >= epsilon and rho[2] >= epsilon
else:
return True
def in_SST(system,vector,proper = False):
"""Non-vectorized version of 'in_SST'."""
if system == '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 system == '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 system == '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 system == '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:
return True
v = np.array(vector,dtype=float)
if proper:
theComponents = np.around(np.dot(basis['improper'],v),12)
inSST = np.all(theComponents >= 0.0)
if not inSST:
theComponents = np.around(np.dot(basis['proper'],v),12)
inSST = np.all(theComponents >= 0.0)
else:
v[2] = abs(v[2])
theComponents = np.around(np.dot(basis['improper'],v),12)
inSST = np.all(theComponents >= 0.0)
return inSST
@pytest.fixture
def set_of_rodrigues(set_of_quaternions):
return Rotation(set_of_quaternions).as_Rodrigues(vector=True)[:200]
class TestSymmetry: class TestSymmetry:
@pytest.mark.parametrize('system',Symmetry.crystal_systems)
def test_in_FZ_vectorize(self,set_of_rodrigues,system):
result = Symmetry(system).in_FZ(set_of_rodrigues.reshape(50,4,3)).reshape(200)
for i,r in enumerate(result):
assert r == in_FZ(system,set_of_rodrigues[i])
@pytest.mark.parametrize('system',Symmetry.crystal_systems)
def test_in_disorientation_SST_vectorize(self,set_of_rodrigues,system):
result = Symmetry(system).in_disorientation_SST(set_of_rodrigues.reshape(50,4,3)).reshape(200)
for i,r in enumerate(result):
assert r == in_disorientation_SST(system,set_of_rodrigues[i])
@pytest.mark.parametrize('proper',[True,False])
@pytest.mark.parametrize('system',Symmetry.crystal_systems)
def test_in_SST_vectorize(self,system,proper):
vecs = np.random.rand(20,4,3)
result = Symmetry(system).in_SST(vecs,proper).reshape(20*4)
for i,r in enumerate(result):
assert r == in_SST(system,vecs.reshape(20*4,3)[i],proper)
@pytest.mark.parametrize('invalid_symmetry',['fcc','bcc','hello']) @pytest.mark.parametrize('invalid_symmetry',['fcc','bcc','hello'])
def test_invalid_symmetry(self,invalid_symmetry): def test_invalid_symmetry(self,invalid_symmetry):
with pytest.raises(KeyError): with pytest.raises(KeyError):
s = Symmetry(invalid_symmetry) # noqa s = Symmetry(invalid_symmetry) # noqa
def test_equal(self): def test_equal(self):
symmetry = random.choice(Symmetry.lattices) symmetry = random.choice(Symmetry.crystal_systems)
print(symmetry) print(symmetry)
assert Symmetry(symmetry) == Symmetry(symmetry) assert Symmetry(symmetry) == Symmetry(symmetry)
def test_not_equal(self): def test_not_equal(self):
symmetries = random.sample(Symmetry.lattices,k=2) symmetries = random.sample(Symmetry.crystal_systems,k=2)
assert Symmetry(symmetries[0]) != Symmetry(symmetries[1]) assert Symmetry(symmetries[0]) != Symmetry(symmetries[1])
@pytest.mark.parametrize('lattice',Symmetry.lattices) @pytest.mark.parametrize('system',Symmetry.crystal_systems)
def test_inFZ(self,lattice): def test_in_FZ(self,system):
assert Symmetry(lattice).inFZ(np.zeros(3)) assert Symmetry(system).in_FZ(np.zeros(3))
@pytest.mark.parametrize('lattice',Symmetry.lattices) @pytest.mark.parametrize('system',Symmetry.crystal_systems)
def test_inDisorientationSST(self,lattice): def test_in_disorientation_SST(self,system):
assert Symmetry(lattice).inDisorientationSST(np.zeros(3)) assert Symmetry(system).in_disorientation_SST(np.zeros(3))
@pytest.mark.parametrize('lattice',Symmetry.lattices) @pytest.mark.parametrize('system',Symmetry.crystal_systems)
@pytest.mark.parametrize('proper',[True,False]) @pytest.mark.parametrize('proper',[True,False])
def test_inSST(self,lattice,proper): def test_in_SST(self,system,proper):
assert Symmetry(lattice).inSST(np.zeros(3),proper) assert Symmetry(system).in_SST(np.zeros(3),proper)
@pytest.mark.parametrize('function',['inFZ','inDisorientationSST']) @pytest.mark.parametrize('function',['in_FZ','in_disorientation_SST','in_SST'])
def test_invalid_argument(self,function): def test_invalid_argument(self,function):
s = Symmetry() # noqa s = Symmetry() # noqa
with pytest.raises(ValueError): with pytest.raises(ValueError):

94
python/tests/test_Orientation.py
View File

@ -11,6 +11,25 @@ from damask import Lattice
n = 1000 n = 1000
def IPF_color(orientation,direction):
"""TSL color of inverse pole figure for given axis (non-vectorized)."""
for o in orientation.equivalent:
pole = o.rotation@direction
inSST,color = orientation.lattice.in_SST(pole,color=True)
if inSST: break
return color
def inverse_pole(orientation,axis,proper=False,SST=True):
if SST:
for eq in orientation.equivalent:
pole = eq.rotation @ axis/np.linalg.norm(axis)
if orientation.lattice.in_SST(pole,proper=proper):
return pole
else:
return orientation.rotation @ axis/np.linalg.norm(axis)
@pytest.fixture @pytest.fixture
def reference_dir(reference_dir_base): def reference_dir(reference_dir_base):
"""Directory containing reference results.""" """Directory containing reference results."""
@ -19,6 +38,31 @@ def reference_dir(reference_dir_base):
class TestOrientation: class TestOrientation:
@pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])
@pytest.mark.parametrize('lattice',['fcc','bcc'])
def test_relationship_vectorize(self,set_of_quaternions,lattice,model):
result = Orientation(set_of_quaternions[:200].reshape(50,4,4),lattice).related(model)
ref_qu = result.rotation.quaternion.reshape(-1,200,4)
for i in range(200):
single = Orientation(set_of_quaternions[i],lattice).related(model).rotation.quaternion
assert np.allclose(ref_qu[:,i,:],single)
@pytest.mark.parametrize('lattice',Lattice.lattices)
def test_IPF_vectorize(self,set_of_quaternions,lattice):
direction = np.random.random(3)*2.0-1
oris = Orientation(Rotation(set_of_quaternions),lattice)[:200]
for i,color in enumerate(oris.IPF_color(direction)):
assert np.allclose(color,IPF_color(oris[i],direction))
@pytest.mark.parametrize('SST',[False,True])
@pytest.mark.parametrize('proper',[True,False])
@pytest.mark.parametrize('lattice',Lattice.lattices)
def test_inverse_pole_vectorize(self,set_of_quaternions,lattice,SST,proper):
axis = np.random.random(3)*2.0-1
oris = Orientation(Rotation(set_of_quaternions),lattice)[:200]
for i,pole in enumerate(oris.inverse_pole(axis,SST=SST)):
assert np.allclose(pole,inverse_pole(oris[i],axis,SST=SST))
@pytest.mark.parametrize('color',[{'label':'red', 'RGB':[1,0,0],'direction':[0,0,1]}, @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':'green','RGB':[0,1,0],'direction':[0,1,1]},
{'label':'blue', 'RGB':[0,0,1],'direction':[1,1,1]}]) {'label':'blue', 'RGB':[0,0,1],'direction':[1,1,1]}])
@ -26,35 +70,63 @@ class TestOrientation:
def test_IPF_cubic(self,color,lattice): def test_IPF_cubic(self,color,lattice):
cube = damask.Orientation(damask.Rotation(),lattice) cube = damask.Orientation(damask.Rotation(),lattice)
for direction in set(permutations(np.array(color['direction']))): for direction in set(permutations(np.array(color['direction']))):
assert np.allclose(cube.IPFcolor(np.array(direction)),np.array(color['RGB'])) assert np.allclose(cube.IPF_color(np.array(direction)),np.array(color['RGB']))
@pytest.mark.parametrize('lattice',Lattice.lattices) @pytest.mark.parametrize('lattice',Lattice.lattices)
def test_IPF(self,lattice): def test_IPF_equivalent(self,set_of_quaternions,lattice):
direction = np.random.random(3)*2.0-1 direction = np.random.random(3)*2.0-1
for rot in [Rotation.from_random() for r in range(n//100)]: for ori in Orientation(Rotation(set_of_quaternions),lattice)[:200]:
R = damask.Orientation(rot,lattice) color = ori.IPF_color(direction)
color = R.IPFcolor(direction) for equivalent in ori.equivalent:
for equivalent in R.equivalentOrientations(): assert np.allclose(color,equivalent.IPF_color(direction))
assert np.allclose(color,R.IPFcolor(direction))
@pytest.mark.parametrize('lattice',Lattice.lattices)
def test_reduced(self,set_of_quaternions,lattice):
oris = Orientation(Rotation(set_of_quaternions),lattice)
reduced = oris.reduced
assert np.all(reduced.in_FZ) and oris.rotation.shape == reduced.rotation.shape
@pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch']) @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])
@pytest.mark.parametrize('lattice',['fcc','bcc']) @pytest.mark.parametrize('lattice',['fcc','bcc'])
def test_relationship_forward_backward(self,model,lattice): def test_relationship_forward_backward(self,model,lattice):
ori = Orientation(Rotation.from_random(),lattice) ori = Orientation(Rotation.from_random(),lattice)
for i,r in enumerate(ori.relatedOrientations(model)): for i,r in enumerate(ori.related(model)):
ori2 = r.relatedOrientations(model)[i] ori2 = r.related(model)[i]
misorientation = ori.rotation.misorientation(ori2.rotation) misorientation = ori.rotation.misorientation(ori2.rotation)
assert misorientation.asAxisAngle(degrees=True)[3]<1.0e-5 assert misorientation.as_axis_angle(degrees=True)[3]<1.0e-5
@pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch']) @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])
@pytest.mark.parametrize('lattice',['fcc','bcc']) @pytest.mark.parametrize('lattice',['fcc','bcc'])
def test_relationship_reference(self,update,reference_dir,model,lattice): def test_relationship_reference(self,update,reference_dir,model,lattice):
reference = os.path.join(reference_dir,f'{lattice}_{model}.txt') reference = os.path.join(reference_dir,f'{lattice}_{model}.txt')
ori = Orientation(Rotation(),lattice) ori = Orientation(Rotation(),lattice)
eu = np.array([o.rotation.as_Eulers(degrees=True) for o in ori.relatedOrientations(model)]) eu = np.array([o.rotation.as_Eulers(degrees=True) for o in ori.related(model)])
if update: if update:
coords = np.array([(1,i+1) for i,x in enumerate(eu)]) coords = np.array([(1,i+1) for i,x in enumerate(eu)])
table = damask.Table(eu,{'Eulers':(3,)}) table = damask.Table(eu,{'Eulers':(3,)})
table.add('pos',coords) table.add('pos',coords)
table.to_ASCII(reference) table.to_ASCII(reference)
assert np.allclose(eu,damask.Table.from_ASCII(reference).get('Eulers')) assert np.allclose(eu,damask.Table.from_ASCII(reference).get('Eulers'))
@pytest.mark.parametrize('lattice',Lattice.lattices)
def test_disorientation360(self,lattice):
R_1 = Orientation(Rotation(),lattice)
R_2 = Orientation(damask.Rotation.from_Eulers([360,0,0],degrees=True),lattice)
assert np.allclose(R_1.disorientation(R_2).as_matrix(),np.eye(3))
@pytest.mark.parametrize('lattice',Lattice.lattices)
@pytest.mark.parametrize('angle',[10,20,30,40])
def test_average(self,angle,lattice):
R_1 = Orientation(Rotation.from_axis_angle([0,0,1,10],degrees=True),lattice)
R_2 = Orientation(Rotation.from_axis_angle([0,0,1,angle],degrees=True),lattice)
avg_angle = R_1.average(R_2).rotation.as_axis_angle(degrees=True,pair=True)[1]
assert np.isclose(avg_angle,10+(angle-10)/2.)
@pytest.mark.parametrize('lattice',Lattice.lattices)
def test_from_average(self,lattice):
R_1 = Orientation(Rotation.from_random(),lattice)
eqs = [r for r in R_1.equivalent]
R_2 = damask.Orientation.from_average(eqs)
assert np.allclose(R_1.rotation.quaternion,R_2.rotation.quaternion)

10
python/tests/test_Result.py
View File

@ -153,16 +153,16 @@ class TestResult:
assert np.allclose(in_memory,in_file) assert np.allclose(in_memory,in_file)
@pytest.mark.parametrize('d',[[1,0,0],[0,1,0],[0,0,1]]) @pytest.mark.parametrize('d',[[1,0,0],[0,1,0],[0,0,1]])
def test_add_IPFcolor(self,default,d): def test_add_IPF_color(self,default,d):
default.add_IPFcolor('orientation',d) default.add_IPF_color('orientation',d)
loc = {'orientation': default.get_dataset_location('orientation'), loc = {'orientation': default.get_dataset_location('orientation'),
'color': default.get_dataset_location('IPFcolor_[{} {} {}]'.format(*d))} 'color': default.get_dataset_location('IPFcolor_[{} {} {}]'.format(*d))}
qu = default.read_dataset(loc['orientation']).view(np.double).reshape(-1,4) qu = default.read_dataset(loc['orientation']).view(np.double).reshape(-1,4)
crystal_structure = default.get_crystal_structure() crystal_structure = default.get_crystal_structure()
in_memory = np.empty((qu.shape[0],3),np.uint8) in_memory = np.empty((qu.shape[0],3),np.uint8)
for i,q in enumerate(qu): for i,q in enumerate(qu):
o = damask.Orientation(q,crystal_structure).reduced() o = damask.Orientation(q,crystal_structure).reduced
in_memory[i] = np.uint8(o.IPFcolor(np.array(d))*255) in_memory[i] = np.uint8(o.IPF_color(np.array(d))*255)
in_file = default.read_dataset(loc['color']) in_file = default.read_dataset(loc['color'])
assert np.allclose(in_memory,in_file) assert np.allclose(in_memory,in_file)
@ -319,4 +319,4 @@ class TestResult:
def test_XDMF(self,tmp_path,single_phase): def test_XDMF(self,tmp_path,single_phase):
os.chdir(tmp_path) os.chdir(tmp_path)
single_phase.write_XDMF single_phase.write_XDMF()

187
python/tests/test_Rotation.py
View File

@ -6,94 +6,21 @@ import numpy as np
from damask import Rotation from damask import Rotation
from damask import _rotation from damask import _rotation
n = 1100 n = 1100
atol=1.e-4 atol=1.e-4
scatter=1.e-2
@pytest.fixture
def default():
"""A set of n rotations (corner cases and random)."""
specials = np.array([
[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,-1.0, 0.0, 0.0],
[0.0, 0.0,-1.0, 0.0],
[0.0, 0.0, 0.0,-1.0],
#----------------------
[1.0, 1.0, 0.0, 0.0],
[1.0, 0.0, 1.0, 0.0],
[1.0, 0.0, 0.0, 1.0],
[0.0, 1.0, 1.0, 0.0],
[0.0, 1.0, 0.0, 1.0],
[0.0, 0.0, 1.0, 1.0],
#----------------------
[1.0,-1.0, 0.0, 0.0],
[1.0, 0.0,-1.0, 0.0],
[1.0, 0.0, 0.0,-1.0],
[0.0, 1.0,-1.0, 0.0],
[0.0, 1.0, 0.0,-1.0],
[0.0, 0.0, 1.0,-1.0],
#----------------------
[0.0, 1.0,-1.0, 0.0],
[0.0, 1.0, 0.0,-1.0],
[0.0, 0.0, 1.0,-1.0],
#----------------------
[0.0,-1.0,-1.0, 0.0],
[0.0,-1.0, 0.0,-1.0],
[0.0, 0.0,-1.0,-1.0],
#----------------------
[1.0, 1.0, 1.0, 0.0],
[1.0, 1.0, 0.0, 1.0],
[1.0, 0.0, 1.0, 1.0],
[1.0,-1.0, 1.0, 0.0],
[1.0,-1.0, 0.0, 1.0],
[1.0, 0.0,-1.0, 1.0],
[1.0, 1.0,-1.0, 0.0],
[1.0, 1.0, 0.0,-1.0],
[1.0, 0.0, 1.0,-1.0],
[1.0,-1.0,-1.0, 0.0],
[1.0,-1.0, 0.0,-1.0],
[1.0, 0.0,-1.0,-1.0],
#----------------------
[0.0, 1.0, 1.0, 1.0],
[0.0, 1.0,-1.0, 1.0],
[0.0, 1.0, 1.0,-1.0],
[0.0,-1.0, 1.0, 1.0],
[0.0,-1.0,-1.0, 1.0],
[0.0,-1.0, 1.0,-1.0],
[0.0,-1.0,-1.0,-1.0],
#----------------------
[1.0, 1.0, 1.0, 1.0],
[1.0,-1.0, 1.0, 1.0],
[1.0, 1.0,-1.0, 1.0],
[1.0, 1.0, 1.0,-1.0],
[1.0,-1.0,-1.0, 1.0],
[1.0,-1.0, 1.0,-1.0],
[1.0, 1.0,-1.0,-1.0],
[1.0,-1.0,-1.0,-1.0],
])
specials /= np.linalg.norm(specials,axis=1).reshape(-1,1)
specials_scatter = specials + np.broadcast_to(np.random.rand(4)*scatter,specials.shape)
specials_scatter /= np.linalg.norm(specials_scatter,axis=1).reshape(-1,1)
specials_scatter[specials_scatter[:,0]<0]*=-1
return [Rotation.from_quaternion(s) for s in specials] + \
[Rotation.from_quaternion(s) for s in specials_scatter] + \
[Rotation.from_random() for _ in range(n-len(specials)-len(specials_scatter))]
@pytest.fixture @pytest.fixture
def reference_dir(reference_dir_base): def reference_dir(reference_dir_base):
"""Directory containing reference results.""" """Directory containing reference results."""
return os.path.join(reference_dir_base,'Rotation') return os.path.join(reference_dir_base,'Rotation')
@pytest.fixture
def set_of_rotations(set_of_quaternions):
return [Rotation.from_quaternion(s) for s in set_of_quaternions]
#################################################################################################### ####################################################################################################
# Code below available according to the following conditions on https://github.com/MarDiehl/3Drotations # Code below available according to the following conditions
#################################################################################################### ####################################################################################################
# Copyright (c) 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH # Copyright (c) 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH
# Copyright (c) 2013-2014, Marc De Graef/Carnegie Mellon University # Copyright (c) 2013-2014, Marc De Graef/Carnegie Mellon University
@ -567,9 +494,9 @@ class TestRotation:
(Rotation._qu2ro,Rotation._ro2qu), (Rotation._qu2ro,Rotation._ro2qu),
(Rotation._qu2ho,Rotation._ho2qu), (Rotation._qu2ho,Rotation._ho2qu),
(Rotation._qu2cu,Rotation._cu2qu)]) (Rotation._qu2cu,Rotation._cu2qu)])
def test_quaternion_internal(self,default,forward,backward): def test_quaternion_internal(self,set_of_rotations,forward,backward):
"""Ensure invariance of conversion from quaternion and back.""" """Ensure invariance of conversion from quaternion and back."""
for rot in default: for rot in set_of_rotations:
m = rot.as_quaternion() m = rot.as_quaternion()
o = backward(forward(m)) o = backward(forward(m))
ok = np.allclose(m,o,atol=atol) ok = np.allclose(m,o,atol=atol)
@ -584,9 +511,9 @@ class TestRotation:
(Rotation._om2ro,Rotation._ro2om), (Rotation._om2ro,Rotation._ro2om),
(Rotation._om2ho,Rotation._ho2om), (Rotation._om2ho,Rotation._ho2om),
(Rotation._om2cu,Rotation._cu2om)]) (Rotation._om2cu,Rotation._cu2om)])
def test_matrix_internal(self,default,forward,backward): def test_matrix_internal(self,set_of_rotations,forward,backward):
"""Ensure invariance of conversion from rotation matrix and back.""" """Ensure invariance of conversion from rotation matrix and back."""
for rot in default: for rot in set_of_rotations:
m = rot.as_matrix() m = rot.as_matrix()
o = backward(forward(m)) o = backward(forward(m))
ok = np.allclose(m,o,atol=atol) ok = np.allclose(m,o,atol=atol)
@ -599,9 +526,9 @@ class TestRotation:
(Rotation._eu2ro,Rotation._ro2eu), (Rotation._eu2ro,Rotation._ro2eu),
(Rotation._eu2ho,Rotation._ho2eu), (Rotation._eu2ho,Rotation._ho2eu),
(Rotation._eu2cu,Rotation._cu2eu)]) (Rotation._eu2cu,Rotation._cu2eu)])
def test_Eulers_internal(self,default,forward,backward): def test_Eulers_internal(self,set_of_rotations,forward,backward):
"""Ensure invariance of conversion from Euler angles and back.""" """Ensure invariance of conversion from Euler angles and back."""
for rot in default: for rot in set_of_rotations:
m = rot.as_Eulers() m = rot.as_Eulers()
o = backward(forward(m)) o = backward(forward(m))
u = np.array([np.pi*2,np.pi,np.pi*2]) u = np.array([np.pi*2,np.pi,np.pi*2])
@ -619,9 +546,9 @@ class TestRotation:
(Rotation._ax2ro,Rotation._ro2ax), (Rotation._ax2ro,Rotation._ro2ax),
(Rotation._ax2ho,Rotation._ho2ax), (Rotation._ax2ho,Rotation._ho2ax),
(Rotation._ax2cu,Rotation._cu2ax)]) (Rotation._ax2cu,Rotation._cu2ax)])
def test_axis_angle_internal(self,default,forward,backward): def test_axis_angle_internal(self,set_of_rotations,forward,backward):
"""Ensure invariance of conversion from axis angle angles pair and back.""" """Ensure invariance of conversion from axis angle angles pair and back."""
for rot in default: for rot in set_of_rotations:
m = rot.as_axis_angle() m = rot.as_axis_angle()
o = backward(forward(m)) o = backward(forward(m))
ok = np.allclose(m,o,atol=atol) ok = np.allclose(m,o,atol=atol)
@ -636,10 +563,10 @@ class TestRotation:
(Rotation._ro2ax,Rotation._ax2ro), (Rotation._ro2ax,Rotation._ax2ro),
(Rotation._ro2ho,Rotation._ho2ro), (Rotation._ro2ho,Rotation._ho2ro),
(Rotation._ro2cu,Rotation._cu2ro)]) (Rotation._ro2cu,Rotation._cu2ro)])
def test_Rodrigues_internal(self,default,forward,backward): def test_Rodrigues_internal(self,set_of_rotations,forward,backward):
"""Ensure invariance of conversion from Rodrigues-Frank vector and back.""" """Ensure invariance of conversion from Rodrigues-Frank vector and back."""
cutoff = np.tan(np.pi*.5*(1.-1e-4)) cutoff = np.tan(np.pi*.5*(1.-1e-4))
for rot in default: for rot in set_of_rotations:
m = rot.as_Rodrigues() m = rot.as_Rodrigues()
o = backward(forward(m)) o = backward(forward(m))
ok = np.allclose(np.clip(m,None,cutoff),np.clip(o,None,cutoff),atol=atol) ok = np.allclose(np.clip(m,None,cutoff),np.clip(o,None,cutoff),atol=atol)
@ -653,9 +580,9 @@ class TestRotation:
(Rotation._ho2ax,Rotation._ax2ho), (Rotation._ho2ax,Rotation._ax2ho),
(Rotation._ho2ro,Rotation._ro2ho), (Rotation._ho2ro,Rotation._ro2ho),
(Rotation._ho2cu,Rotation._cu2ho)]) (Rotation._ho2cu,Rotation._cu2ho)])
def test_homochoric_internal(self,default,forward,backward): def test_homochoric_internal(self,set_of_rotations,forward,backward):
"""Ensure invariance of conversion from homochoric vector and back.""" """Ensure invariance of conversion from homochoric vector and back."""
for rot in default: for rot in set_of_rotations:
m = rot.as_homochoric() m = rot.as_homochoric()
o = backward(forward(m)) o = backward(forward(m))
ok = np.allclose(m,o,atol=atol) ok = np.allclose(m,o,atol=atol)
@ -668,9 +595,9 @@ class TestRotation:
(Rotation._cu2ax,Rotation._ax2cu), (Rotation._cu2ax,Rotation._ax2cu),
(Rotation._cu2ro,Rotation._ro2cu), (Rotation._cu2ro,Rotation._ro2cu),
(Rotation._cu2ho,Rotation._ho2cu)]) (Rotation._cu2ho,Rotation._ho2cu)])
def test_cubochoric_internal(self,default,forward,backward): def test_cubochoric_internal(self,set_of_rotations,forward,backward):
"""Ensure invariance of conversion from cubochoric vector and back.""" """Ensure invariance of conversion from cubochoric vector and back."""
for rot in default: for rot in set_of_rotations:
m = rot.as_cubochoric() m = rot.as_cubochoric()
o = backward(forward(m)) o = backward(forward(m))
ok = np.allclose(m,o,atol=atol) ok = np.allclose(m,o,atol=atol)
@ -684,9 +611,9 @@ class TestRotation:
(Rotation._qu2ax,qu2ax), (Rotation._qu2ax,qu2ax),
(Rotation._qu2ro,qu2ro), (Rotation._qu2ro,qu2ro),
(Rotation._qu2ho,qu2ho)]) (Rotation._qu2ho,qu2ho)])
def test_quaternion_vectorization(self,default,vectorized,single): def test_quaternion_vectorization(self,set_of_quaternions,vectorized,single):
"""Check vectorized implementation for quaternion against single point calculation.""" """Check vectorized implementation for quaternion against single point calculation."""
qu = np.array([rot.as_quaternion() for rot in default]) qu = np.array(set_of_quaternions)
vectorized(qu.reshape(qu.shape[0]//2,-1,4)) vectorized(qu.reshape(qu.shape[0]//2,-1,4))
co = vectorized(qu) co = vectorized(qu)
for q,c in zip(qu,co): for q,c in zip(qu,co):
@ -697,9 +624,9 @@ class TestRotation:
@pytest.mark.parametrize('vectorized, single',[(Rotation._om2qu,om2qu), @pytest.mark.parametrize('vectorized, single',[(Rotation._om2qu,om2qu),
(Rotation._om2eu,om2eu), (Rotation._om2eu,om2eu),
(Rotation._om2ax,om2ax)]) (Rotation._om2ax,om2ax)])
def test_matrix_vectorization(self,default,vectorized,single): def test_matrix_vectorization(self,set_of_rotations,vectorized,single):
"""Check vectorized implementation for rotation matrix against single point calculation.""" """Check vectorized implementation for rotation matrix against single point calculation."""
om = np.array([rot.as_matrix() for rot in default]) om = np.array([rot.as_matrix() for rot in set_of_rotations])
vectorized(om.reshape(om.shape[0]//2,-1,3,3)) vectorized(om.reshape(om.shape[0]//2,-1,3,3))
co = vectorized(om) co = vectorized(om)
for o,c in zip(om,co): for o,c in zip(om,co):
@ -710,9 +637,9 @@ class TestRotation:
(Rotation._eu2om,eu2om), (Rotation._eu2om,eu2om),
(Rotation._eu2ax,eu2ax), (Rotation._eu2ax,eu2ax),
(Rotation._eu2ro,eu2ro)]) (Rotation._eu2ro,eu2ro)])
def test_Eulers_vectorization(self,default,vectorized,single): def test_Eulers_vectorization(self,set_of_rotations,vectorized,single):
"""Check vectorized implementation for Euler angles against single point calculation.""" """Check vectorized implementation for Euler angles against single point calculation."""
eu = np.array([rot.as_Eulers() for rot in default]) eu = np.array([rot.as_Eulers() for rot in set_of_rotations])
vectorized(eu.reshape(eu.shape[0]//2,-1,3)) vectorized(eu.reshape(eu.shape[0]//2,-1,3))
co = vectorized(eu) co = vectorized(eu)
for e,c in zip(eu,co): for e,c in zip(eu,co):
@ -723,9 +650,9 @@ class TestRotation:
(Rotation._ax2om,ax2om), (Rotation._ax2om,ax2om),
(Rotation._ax2ro,ax2ro), (Rotation._ax2ro,ax2ro),
(Rotation._ax2ho,ax2ho)]) (Rotation._ax2ho,ax2ho)])
def test_axis_angle_vectorization(self,default,vectorized,single): def test_axis_angle_vectorization(self,set_of_rotations,vectorized,single):
"""Check vectorized implementation for axis angle pair against single point calculation.""" """Check vectorized implementation for axis angle pair against single point calculation."""
ax = np.array([rot.as_axis_angle() for rot in default]) ax = np.array([rot.as_axis_angle() for rot in set_of_rotations])
vectorized(ax.reshape(ax.shape[0]//2,-1,4)) vectorized(ax.reshape(ax.shape[0]//2,-1,4))
co = vectorized(ax) co = vectorized(ax)
for a,c in zip(ax,co): for a,c in zip(ax,co):
@ -735,9 +662,9 @@ class TestRotation:
@pytest.mark.parametrize('vectorized, single',[(Rotation._ro2ax,ro2ax), @pytest.mark.parametrize('vectorized, single',[(Rotation._ro2ax,ro2ax),
(Rotation._ro2ho,ro2ho)]) (Rotation._ro2ho,ro2ho)])
def test_Rodrigues_vectorization(self,default,vectorized,single): def test_Rodrigues_vectorization(self,set_of_rotations,vectorized,single):
"""Check vectorized implementation for Rodrigues-Frank vector against single point calculation.""" """Check vectorized implementation for Rodrigues-Frank vector against single point calculation."""
ro = np.array([rot.as_Rodrigues() for rot in default]) ro = np.array([rot.as_Rodrigues() for rot in set_of_rotations])
vectorized(ro.reshape(ro.shape[0]//2,-1,4)) vectorized(ro.reshape(ro.shape[0]//2,-1,4))
co = vectorized(ro) co = vectorized(ro)
for r,c in zip(ro,co): for r,c in zip(ro,co):
@ -746,9 +673,9 @@ class TestRotation:
@pytest.mark.parametrize('vectorized, single',[(Rotation._ho2ax,ho2ax), @pytest.mark.parametrize('vectorized, single',[(Rotation._ho2ax,ho2ax),
(Rotation._ho2cu,ho2cu)]) (Rotation._ho2cu,ho2cu)])
def test_homochoric_vectorization(self,default,vectorized,single): def test_homochoric_vectorization(self,set_of_rotations,vectorized,single):
"""Check vectorized implementation for homochoric vector against single point calculation.""" """Check vectorized implementation for homochoric vector against single point calculation."""
ho = np.array([rot.as_homochoric() for rot in default]) ho = np.array([rot.as_homochoric() for rot in set_of_rotations])
vectorized(ho.reshape(ho.shape[0]//2,-1,3)) vectorized(ho.reshape(ho.shape[0]//2,-1,3))
co = vectorized(ho) co = vectorized(ho)
for h,c in zip(ho,co): for h,c in zip(ho,co):
@ -756,9 +683,9 @@ class TestRotation:
assert np.allclose(single(h),c) and np.allclose(single(h),vectorized(h)) assert np.allclose(single(h),c) and np.allclose(single(h),vectorized(h))
@pytest.mark.parametrize('vectorized, single',[(Rotation._cu2ho,cu2ho)]) @pytest.mark.parametrize('vectorized, single',[(Rotation._cu2ho,cu2ho)])
def test_cubochoric_vectorization(self,default,vectorized,single): def test_cubochoric_vectorization(self,set_of_rotations,vectorized,single):
"""Check vectorized implementation for cubochoric vector against single point calculation.""" """Check vectorized implementation for cubochoric vector against single point calculation."""
cu = np.array([rot.as_cubochoric() for rot in default]) cu = np.array([rot.as_cubochoric() for rot in set_of_rotations])
vectorized(cu.reshape(cu.shape[0]//2,-1,3)) vectorized(cu.reshape(cu.shape[0]//2,-1,3))
co = vectorized(cu) co = vectorized(cu)
for u,c in zip(cu,co): for u,c in zip(cu,co):
@ -766,8 +693,8 @@ class TestRotation:
assert np.allclose(single(u),c) and np.allclose(single(u),vectorized(u)) assert np.allclose(single(u),c) and np.allclose(single(u),vectorized(u))
@pytest.mark.parametrize('degrees',[True,False]) @pytest.mark.parametrize('degrees',[True,False])
def test_Eulers(self,default,degrees): def test_Eulers(self,set_of_rotations,degrees):
for rot in default: for rot in set_of_rotations:
m = rot.as_quaternion() m = rot.as_quaternion()
o = Rotation.from_Eulers(rot.as_Eulers(degrees),degrees).as_quaternion() o = Rotation.from_Eulers(rot.as_Eulers(degrees),degrees).as_quaternion()
ok = np.allclose(m,o,atol=atol) ok = np.allclose(m,o,atol=atol)
@ -779,9 +706,9 @@ class TestRotation:
@pytest.mark.parametrize('P',[1,-1]) @pytest.mark.parametrize('P',[1,-1])
@pytest.mark.parametrize('normalise',[True,False]) @pytest.mark.parametrize('normalise',[True,False])
@pytest.mark.parametrize('degrees',[True,False]) @pytest.mark.parametrize('degrees',[True,False])
def test_axis_angle(self,default,degrees,normalise,P): def test_axis_angle(self,set_of_rotations,degrees,normalise,P):
c = np.array([P*-1,P*-1,P*-1,1.]) c = np.array([P*-1,P*-1,P*-1,1.])
for rot in default: for rot in set_of_rotations:
m = rot.as_Eulers() m = rot.as_Eulers()
o = Rotation.from_axis_angle(rot.as_axis_angle(degrees)*c,degrees,normalise,P).as_Eulers() o = Rotation.from_axis_angle(rot.as_axis_angle(degrees)*c,degrees,normalise,P).as_Eulers()
u = np.array([np.pi*2,np.pi,np.pi*2]) u = np.array([np.pi*2,np.pi,np.pi*2])
@ -793,8 +720,8 @@ class TestRotation:
print(m,o,rot.as_quaternion()) print(m,o,rot.as_quaternion())
assert ok and (np.zeros(3)-1.e-9 <= o).all() and (o <= np.array([np.pi*2.,np.pi,np.pi*2.])+1.e-9).all() assert ok and (np.zeros(3)-1.e-9 <= o).all() and (o <= np.array([np.pi*2.,np.pi,np.pi*2.])+1.e-9).all()
def test_matrix(self,default): def test_matrix(self,set_of_rotations):
for rot in default: for rot in set_of_rotations:
m = rot.as_axis_angle() m = rot.as_axis_angle()
o = Rotation.from_axis_angle(rot.as_axis_angle()).as_axis_angle() o = Rotation.from_axis_angle(rot.as_axis_angle()).as_axis_angle()
ok = np.allclose(m,o,atol=atol) ok = np.allclose(m,o,atol=atol)
@ -805,9 +732,9 @@ class TestRotation:
@pytest.mark.parametrize('P',[1,-1]) @pytest.mark.parametrize('P',[1,-1])
@pytest.mark.parametrize('normalise',[True,False]) @pytest.mark.parametrize('normalise',[True,False])
def test_Rodrigues(self,default,normalise,P): def test_Rodrigues(self,set_of_rotations,normalise,P):
c = np.array([P*-1,P*-1,P*-1,1.]) c = np.array([P*-1,P*-1,P*-1,1.])
for rot in default: for rot in set_of_rotations:
m = rot.as_matrix() m = rot.as_matrix()
o = Rotation.from_Rodrigues(rot.as_Rodrigues()*c,normalise,P).as_matrix() o = Rotation.from_Rodrigues(rot.as_Rodrigues()*c,normalise,P).as_matrix()
ok = np.allclose(m,o,atol=atol) ok = np.allclose(m,o,atol=atol)
@ -815,9 +742,9 @@ class TestRotation:
assert ok and np.isclose(np.linalg.det(o),1.0) assert ok and np.isclose(np.linalg.det(o),1.0)
@pytest.mark.parametrize('P',[1,-1]) @pytest.mark.parametrize('P',[1,-1])
def test_homochoric(self,default,P): def test_homochoric(self,set_of_rotations,P):
cutoff = np.tan(np.pi*.5*(1.-1e-4)) cutoff = np.tan(np.pi*.5*(1.-1e-4))
for rot in default: for rot in set_of_rotations:
m = rot.as_Rodrigues() m = rot.as_Rodrigues()
o = Rotation.from_homochoric(rot.as_homochoric()*P*-1,P).as_Rodrigues() o = Rotation.from_homochoric(rot.as_homochoric()*P*-1,P).as_Rodrigues()
ok = np.allclose(np.clip(m,None,cutoff),np.clip(o,None,cutoff),atol=atol) ok = np.allclose(np.clip(m,None,cutoff),np.clip(o,None,cutoff),atol=atol)
@ -826,8 +753,8 @@ class TestRotation:
assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) assert ok and np.isclose(np.linalg.norm(o[:3]),1.0)
@pytest.mark.parametrize('P',[1,-1]) @pytest.mark.parametrize('P',[1,-1])
def test_cubochoric(self,default,P): def test_cubochoric(self,set_of_rotations,P):
for rot in default: for rot in set_of_rotations:
m = rot.as_homochoric() m = rot.as_homochoric()
o = Rotation.from_cubochoric(rot.as_cubochoric()*P*-1,P).as_homochoric() o = Rotation.from_cubochoric(rot.as_cubochoric()*P*-1,P).as_homochoric()
ok = np.allclose(m,o,atol=atol) ok = np.allclose(m,o,atol=atol)
@ -836,9 +763,9 @@ class TestRotation:
@pytest.mark.parametrize('P',[1,-1]) @pytest.mark.parametrize('P',[1,-1])
@pytest.mark.parametrize('accept_homomorph',[True,False]) @pytest.mark.parametrize('accept_homomorph',[True,False])
def test_quaternion(self,default,P,accept_homomorph): def test_quaternion(self,set_of_rotations,P,accept_homomorph):
c = np.array([1,P*-1,P*-1,P*-1]) * (-1 if accept_homomorph else 1) c = np.array([1,P*-1,P*-1,P*-1]) * (-1 if accept_homomorph else 1)
for rot in default: for rot in set_of_rotations:
m = rot.as_cubochoric() m = rot.as_cubochoric()
o = Rotation.from_quaternion(rot.as_quaternion()*c,accept_homomorph,P).as_cubochoric() o = Rotation.from_quaternion(rot.as_quaternion()*c,accept_homomorph,P).as_cubochoric()
ok = np.allclose(m,o,atol=atol) ok = np.allclose(m,o,atol=atol)
@ -848,8 +775,8 @@ class TestRotation:
assert ok and o.max() < np.pi**(2./3.)*0.5+1.e-9 assert ok and o.max() < np.pi**(2./3.)*0.5+1.e-9
@pytest.mark.parametrize('reciprocal',[True,False]) @pytest.mark.parametrize('reciprocal',[True,False])
def test_basis(self,default,reciprocal): def test_basis(self,set_of_rotations,reciprocal):
for rot in default: for rot in set_of_rotations:
om = rot.as_matrix() + 0.1*np.eye(3) om = rot.as_matrix() + 0.1*np.eye(3)
rot = Rotation.from_basis(om,False,reciprocal=reciprocal) rot = Rotation.from_basis(om,False,reciprocal=reciprocal)
assert np.isclose(np.linalg.det(rot.as_matrix()),1.0) assert np.isclose(np.linalg.det(rot.as_matrix()),1.0)
@ -923,8 +850,8 @@ class TestRotation:
@pytest.mark.parametrize('data',[np.random.rand(5,3), @pytest.mark.parametrize('data',[np.random.rand(5,3),
np.random.rand(5,3,3), np.random.rand(5,3,3),
np.random.rand(5,3,3,3,3)]) np.random.rand(5,3,3,3,3)])
def test_rotate_vectorization(self,default,data): def test_rotate_vectorization(self,set_of_rotations,data):
for rot in default: for rot in set_of_rotations:
v = rot.broadcast_to((5,)) @ data v = rot.broadcast_to((5,)) @ data
for i in range(data.shape[0]): for i in range(data.shape[0]):
print(i-data[i]) print(i-data[i])
@ -978,3 +905,15 @@ class TestRotation:
def test_misorientation(self): def test_misorientation(self):
R = Rotation.from_random() R = Rotation.from_random()
assert np.allclose(R.misorientation(R).as_matrix(),np.eye(3)) assert np.allclose(R.misorientation(R).as_matrix(),np.eye(3))
def test_misorientation360(self):
R_1 = Rotation()
R_2 = Rotation.from_Eulers([360,0,0],degrees=True)
assert np.allclose(R_1.misorientation(R_2).as_matrix(),np.eye(3))
@pytest.mark.parametrize('angle',[10,20,30,40,50,60,70,80,90,100,120])
def test_average(self,angle):
R_1 = Rotation.from_axis_angle([0,0,1,10],degrees=True)
R_2 = Rotation.from_axis_angle([0,0,1,angle],degrees=True)
avg_angle = R_1.average(R_2).as_axis_angle(degrees=True,pair=True)[1]
assert np.isclose(avg_angle,10+(angle-10)/2.)

30
python/tests/test_VTK.py
View File

@ -1,3 +1,7 @@
import os
import filecmp
import time
import pytest import pytest
import numpy as np import numpy as np
@ -17,7 +21,7 @@ class TestVTK:
origin = np.random.random(3) origin = np.random.random(3)
v = VTK.from_rectilinearGrid(grid,size,origin) v = VTK.from_rectilinearGrid(grid,size,origin)
string = v.__repr__() string = v.__repr__()
v.write(tmp_path/'rectilinearGrid') v.write(tmp_path/'rectilinearGrid',False)
vtr = VTK.from_file(tmp_path/'rectilinearGrid.vtr') vtr = VTK.from_file(tmp_path/'rectilinearGrid.vtr')
with open(tmp_path/'rectilinearGrid.vtk','w') as f: with open(tmp_path/'rectilinearGrid.vtk','w') as f:
f.write(string) f.write(string)
@ -25,10 +29,10 @@ class TestVTK:
assert(string == vtr.__repr__() == vtk.__repr__()) assert(string == vtr.__repr__() == vtk.__repr__())
def test_polyData(self,tmp_path): def test_polyData(self,tmp_path):
points = np.random.rand(3,100) points = np.random.rand(100,3)
v = VTK.from_polyData(points) v = VTK.from_polyData(points)
string = v.__repr__() string = v.__repr__()
v.write(tmp_path/'polyData') v.write(tmp_path/'polyData',False)
vtp = VTK.from_file(tmp_path/'polyData.vtp') vtp = VTK.from_file(tmp_path/'polyData.vtp')
with open(tmp_path/'polyData.vtk','w') as f: with open(tmp_path/'polyData.vtk','w') as f:
f.write(string) f.write(string)
@ -47,14 +51,30 @@ class TestVTK:
connectivity = np.random.choice(np.arange(n),n,False).reshape(-1,n) connectivity = np.random.choice(np.arange(n),n,False).reshape(-1,n)
v = VTK.from_unstructuredGrid(nodes,connectivity,cell_type) v = VTK.from_unstructuredGrid(nodes,connectivity,cell_type)
string = v.__repr__() string = v.__repr__()
v.write(tmp_path/'unstructuredGrid') v.write(tmp_path/'unstructuredGrid',False)
vtu = VTK.from_file(tmp_path/'unstructuredGrid.vtu') vtu = VTK.from_file(tmp_path/'unstructuredGrid.vtu')
with open(tmp_path/'unstructuredGrid.vtk','w') as f: with open(tmp_path/'unstructuredGrid.vtk','w') as f:
f.write(string) f.write(string)
vtk = VTK.from_file(tmp_path/'unstructuredGrid.vtk','unstructuredgrid') vtk = VTK.from_file(tmp_path/'unstructuredGrid.vtk','unstructuredgrid')
assert(string == vtu.__repr__() == vtk.__repr__()) assert(string == vtu.__repr__() == vtk.__repr__())
@pytest.mark.parametrize('name,dataset_type',[('this_file_does_not_exist.vtk',None),
def test_parallel_out(self,tmp_path):
points = np.random.rand(102,3)
v = VTK.from_polyData(points)
fname_s = tmp_path/'single.vtp'
fname_p = tmp_path/'parallel.vtp'
v.write(fname_s,False)
v.write(fname_p,True)
for i in range(10):
if os.path.isfile(fname_p) and filecmp.cmp(fname_s,fname_p):
assert(True)
return
time.sleep(.5)
assert(False)
@pytest.mark.parametrize('name,dataset_type',[('this_file_does_not_exist.vtk', None),
('this_file_does_not_exist.vtk','vtk'), ('this_file_does_not_exist.vtk','vtk'),
('this_file_does_not_exist.vtx', None)]) ('this_file_does_not_exist.vtx', None)])
def test_invalid_dataset_type(self,dataset_type,name): def test_invalid_dataset_type(self,dataset_type,name):

12
src/DAMASK_interface.f90
View File

@ -6,7 +6,7 @@
!> @brief Interfacing between the PETSc-based solvers and the material subroutines provided !> @brief Interfacing between the PETSc-based solvers and the material subroutines provided
!! by DAMASK !! by DAMASK
!> @details Interfacing between the PETSc-based solvers and the material subroutines provided !> @details Interfacing between the PETSc-based solvers and the material subroutines provided
!> by DAMASK. Interpretating the command line arguments to get load case, geometry file, !> by DAMASK. Interpreting the command line arguments to get load case, geometry file,
!> and working directory. !> and working directory.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
#define PETSC_MAJOR 3 #define PETSC_MAJOR 3
@ -20,11 +20,11 @@ module DAMASK_interface
use prec use prec
use system_routines use system_routines
implicit none implicit none
private private
logical, volatile, public, protected :: & logical, volatile, public, protected :: &
SIGTERM, & !< termination signal SIGTERM, & !< termination signal
SIGUSR1, & !< 1. user-defined signal SIGUSR1, & !< 1. user-defined signal
SIGUSR2 !< 2. user-defined signal SIGUSR2 !< 2. user-defined signal
integer, public, protected :: & integer, public, protected :: &
@ -159,14 +159,14 @@ subroutine DAMASK_interface_init
call MPI_Type_size(MPI_INTEGER,typeSize,err) call MPI_Type_size(MPI_INTEGER,typeSize,err)
if (err /= 0) call quit(1) if (err /= 0) call quit(1)
if (typeSize*8 /= bit_size(0)) then if (typeSize*8 /= bit_size(0)) then
write(6,'(a)') ' Mismatch between MPI and DAMASK integer' write(6,'(a)') ' Mismatch between MPI and DAMASK integer'
call quit(1) call quit(1)
endif endif
call MPI_Type_size(MPI_DOUBLE,typeSize,err) call MPI_Type_size(MPI_DOUBLE,typeSize,err)
if (err /= 0) call quit(1) if (err /= 0) call quit(1)
if (typeSize*8 /= storage_size(0.0_pReal)) then if (typeSize*8 /= storage_size(0.0_pReal)) then
write(6,'(a)') ' Mismatch between MPI and DAMASK real' write(6,'(a)') ' Mismatch between MPI and DAMASK real'
call quit(1) call quit(1)
endif endif
@ -340,7 +340,7 @@ end function getGeometryFile
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief relative path of loadcase from command line arguments !> @brief relative path of load case from command line arguments
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function getLoadCaseFile(loadCaseParameter) function getLoadCaseFile(loadCaseParameter)

6
src/constitutive.f90
View File

@ -400,7 +400,7 @@ subroutine constitutive_init
constitutive_source_maxSizeDotState = 0 constitutive_source_maxSizeDotState = 0
PhaseLoop2:do ph = 1,material_Nphase PhaseLoop2:do ph = 1,material_Nphase
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! partition and inititalize state ! partition and initialize state
plasticState(ph)%partionedState0 = plasticState(ph)%state0 plasticState(ph)%partionedState0 = plasticState(ph)%state0
plasticState(ph)%state = plasticState(ph)%partionedState0 plasticState(ph)%state = plasticState(ph)%partionedState0
forall(s = 1:phase_Nsources(ph)) forall(s = 1:phase_Nsources(ph))
@ -475,7 +475,7 @@ end subroutine constitutive_dependentState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient !> @brief contains the constitutive equation for calculating the velocity gradient
! ToDo: Discuss wheter it makes sense if crystallite handles the configuration conversion, i.e. ! ToDo: Discuss whether it makes sense if crystallite handles the configuration conversion, i.e.
! Mp in, dLp_dMp out ! Mp in, dLp_dMp out
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
@ -686,7 +686,7 @@ end subroutine constitutive_SandItsTangents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to !> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
!> the elastic and intermeidate deformation gradients using Hookes law !> the elastic and intermediate deformation gradients using Hooke's law
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
Fe, Fi, ipc, ip, el) Fe, Fi, ipc, ip, el)

2
src/constitutive_plastic_nonlocal.f90
View File

@ -37,7 +37,7 @@ submodule(constitutive) plastic_nonlocal
! BEGIN DEPRECATED ! BEGIN DEPRECATED
integer, dimension(:,:,:), allocatable :: & integer, dimension(:,:,:), allocatable :: &
iRhoU, & !< state indices for unblocked density iRhoU, & !< state indices for unblocked density
iV, & !< state indices for dislcation velocities iV, & !< state indices for dislocation velocities
iD !< state indices for stable dipole height iD !< state indices for stable dipole height
!END DEPRECATED !END DEPRECATED

2
src/kinematics_thermal_expansion.f90
View File

@ -100,7 +100,7 @@ end function kinematics_thermal_expansion_initialStrain
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient !> @brief constitutive equation for calculating the velocity gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el) subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el)

76
src/list.f90
View File

@ -1,18 +1,18 @@
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief linked list !> @brief Linked list
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module list module list
use prec use prec
use IO use IO
implicit none implicit none
private private
type, private :: tPartitionedString type, private :: tPartitionedString
character(len=:), allocatable :: val character(len=:), allocatable :: val
integer, dimension(:), allocatable :: pos integer, dimension(:), allocatable :: pos
end type tPartitionedString end type tPartitionedString
type, public :: tPartitionedStringList type, public :: tPartitionedStringList
type(tPartitionedString) :: string type(tPartitionedString) :: string
type(tPartitionedStringList), pointer :: next => null() type(tPartitionedStringList), pointer :: next => null()
@ -20,31 +20,31 @@ module list
procedure :: add => add procedure :: add => add
procedure :: show => show procedure :: show => show
procedure :: free => free procedure :: free => free
! currently, a finalize is needed for all shapes of tPartitionedStringList. ! currently, a finalize is needed for all shapes of tPartitionedStringList.
! with Fortran 2015, we can define one recursive elemental function ! with Fortran 2015, we can define one recursive elemental function
! https://software.intel.com/en-us/forums/intel-visual-fortran-compiler-for-windows/topic/543326 ! https://software.intel.com/en-us/forums/intel-visual-fortran-compiler-for-windows/topic/543326
final :: finalize, & final :: finalize, &
finalizeArray finalizeArray
procedure :: keyExists => keyExists procedure :: keyExists => keyExists
procedure :: countKeys => countKeys procedure :: countKeys => countKeys
procedure :: getFloat => getFloat procedure :: getFloat => getFloat
procedure :: getInt => getInt procedure :: getInt => getInt
procedure :: getString => getString procedure :: getString => getString
procedure :: getFloats => getFloats procedure :: getFloats => getFloats
procedure :: getInts => getInts procedure :: getInts => getInts
procedure :: getStrings => getStrings procedure :: getStrings => getStrings
end type tPartitionedStringList end type tPartitionedStringList
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief add element !> @brief add element
!> @details Adds a string together with the start/end position of chunks in this string. The new !> @details Adds a string together with the start/end position of chunks in this string. The new
!! element is added at the end of the list. Empty strings are not added. All strings are converted !! element is added at the end of the list. Empty strings are not added. All strings are converted
!! to lower case. The data is not stored in the new element but in the current. !! to lower case. The data is not stored in the new element but in the current.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -177,7 +177,7 @@ end function countKeys
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief gets float value of for a given key from a linked list !> @brief gets float value of for a given key from a linked list
!> @details gets the last value if the key occurs more than once. If key is not found exits with !> @details gets the last value if the key occurs more than once. If key is not found exits with
!! error unless default is given !! error unless default is given
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
real(pReal) function getFloat(this,key,defaultVal) real(pReal) function getFloat(this,key,defaultVal)
@ -187,11 +187,11 @@ real(pReal) function getFloat(this,key,defaultVal)
real(pReal), intent(in), optional :: defaultVal real(pReal), intent(in), optional :: defaultVal
type(tPartitionedStringList), pointer :: item type(tPartitionedStringList), pointer :: item
logical :: found logical :: found
getFloat = huge(1.0) ! suppress warning about unitialized value getFloat = huge(1.0) ! suppress warning about unitialized value
found = present(defaultVal) found = present(defaultVal)
if (found) getFloat = defaultVal if (found) getFloat = defaultVal
item => this item => this
do while (associated(item%next)) do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
@ -209,7 +209,7 @@ end function getFloat
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief gets integer value of for a given key from a linked list !> @brief gets integer value of for a given key from a linked list
!> @details gets the last value if the key occurs more than once. If key is not found exits with !> @details gets the last value if the key occurs more than once. If key is not found exits with
!! error unless default is given !! error unless default is given
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
integer function getInt(this,key,defaultVal) integer function getInt(this,key,defaultVal)
@ -223,7 +223,7 @@ integer function getInt(this,key,defaultVal)
getInt = huge(1) ! suppress warning about unitialized value getInt = huge(1) ! suppress warning about unitialized value
found = present(defaultVal) found = present(defaultVal)
if (found) getInt = defaultVal if (found) getInt = defaultVal
item => this item => this
do while (associated(item%next)) do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
@ -241,12 +241,12 @@ end function getInt
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief gets string value of for a given key from a linked list !> @brief gets string value of for a given key from a linked list
!> @details gets the last value if the key occurs more than once. If key is not found exits with !> @details gets the last value if the key occurs more than once. If key is not found exits with
!! error unless default is given. If raw is true, the the complete string is returned, otherwise !! error unless default is given. If raw is true, the the complete string is returned, otherwise
!! the individual chunks are returned !! the individual chunks are returned
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
character(len=pStringLen) function getString(this,key,defaultVal,raw) character(len=pStringLen) function getString(this,key,defaultVal,raw)
class(tPartitionedStringList), target, intent(in) :: this class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key character(len=*), intent(in) :: key
character(len=*), intent(in), optional :: defaultVal character(len=*), intent(in), optional :: defaultVal
@ -259,19 +259,19 @@ character(len=pStringLen) function getString(this,key,defaultVal,raw)
else else
whole = .false. whole = .false.
endif endif
found = present(defaultVal) found = present(defaultVal)
if (found) then if (found) then
if (len_trim(defaultVal) > len(getString)) call IO_error(0,ext_msg='getString') if (len_trim(defaultVal) > len(getString)) call IO_error(0,ext_msg='getString')
getString = trim(defaultVal) getString = trim(defaultVal)
endif endif
item => this item => this
do while (associated(item%next)) do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true. found = .true.
if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key) if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key)
if (whole) then if (whole) then
getString = trim(item%string%val(item%string%pos(4):)) ! raw string starting a second chunk getString = trim(item%string%val(item%string%pos(4):)) ! raw string starting a second chunk
else else
@ -280,7 +280,7 @@ character(len=pStringLen) function getString(this,key,defaultVal,raw)
endif endif
item => item%next item => item%next
enddo enddo
if (.not. found) call IO_error(140,ext_msg=key) if (.not. found) call IO_error(140,ext_msg=key)
end function getString end function getString
@ -292,7 +292,7 @@ end function getString
!! values from the last occurrence. If key is not found exits with error unless default is given. !! values from the last occurrence. If key is not found exits with error unless default is given.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function getFloats(this,key,defaultVal,requiredSize) function getFloats(this,key,defaultVal,requiredSize)
real(pReal), dimension(:), allocatable :: getFloats real(pReal), dimension(:), allocatable :: getFloats
class(tPartitionedStringList), target, intent(in) :: this class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key character(len=*), intent(in) :: key
@ -302,12 +302,12 @@ function getFloats(this,key,defaultVal,requiredSize)
integer :: i integer :: i
logical :: found, & logical :: found, &
cumulative cumulative
cumulative = (key(1:1) == '(' .and. key(len_trim(key):len_trim(key)) == ')') cumulative = (key(1:1) == '(' .and. key(len_trim(key):len_trim(key)) == ')')
found = .false. found = .false.
allocate(getFloats(0)) allocate(getFloats(0))
item => this item => this
do while (associated(item%next)) do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
@ -320,7 +320,7 @@ function getFloats(this,key,defaultVal,requiredSize)
endif endif
item => item%next item => item%next
enddo enddo
if (.not. found) then if (.not. found) then
if (present(defaultVal)) then; getFloats = defaultVal; else; call IO_error(140,ext_msg=key); endif if (present(defaultVal)) then; getFloats = defaultVal; else; call IO_error(140,ext_msg=key); endif
endif endif
@ -337,7 +337,7 @@ end function getFloats
!! values from the last occurrence. If key is not found exits with error unless default is given. !! values from the last occurrence. If key is not found exits with error unless default is given.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function getInts(this,key,defaultVal,requiredSize) function getInts(this,key,defaultVal,requiredSize)
integer, dimension(:), allocatable :: getInts integer, dimension(:), allocatable :: getInts
class(tPartitionedStringList), target, intent(in) :: this class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key character(len=*), intent(in) :: key
@ -347,12 +347,12 @@ function getInts(this,key,defaultVal,requiredSize)
integer :: i integer :: i
logical :: found, & logical :: found, &
cumulative cumulative
cumulative = (key(1:1) == '(' .and. key(len_trim(key):len_trim(key)) == ')') cumulative = (key(1:1) == '(' .and. key(len_trim(key):len_trim(key)) == ')')
found = .false. found = .false.
allocate(getInts(0)) allocate(getInts(0))
item => this item => this
do while (associated(item%next)) do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
@ -365,7 +365,7 @@ function getInts(this,key,defaultVal,requiredSize)
endif endif
item => item%next item => item%next
enddo enddo
if (.not. found) then if (.not. found) then
if (present(defaultVal)) then; getInts = defaultVal; else; call IO_error(140,ext_msg=key); endif if (present(defaultVal)) then; getInts = defaultVal; else; call IO_error(140,ext_msg=key); endif
endif endif
@ -383,7 +383,7 @@ end function getInts
!! If raw is true, the the complete string is returned, otherwise the individual chunks are returned !! If raw is true, the the complete string is returned, otherwise the individual chunks are returned
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function getStrings(this,key,defaultVal,raw) function getStrings(this,key,defaultVal,raw)
character(len=pStringLen),dimension(:), allocatable :: getStrings character(len=pStringLen),dimension(:), allocatable :: getStrings
class(tPartitionedStringList),target, intent(in) :: this class(tPartitionedStringList),target, intent(in) :: this
character(len=*), intent(in) :: key character(len=*), intent(in) :: key
@ -395,7 +395,7 @@ function getStrings(this,key,defaultVal,raw)
logical :: found, & logical :: found, &
whole, & whole, &
cumulative cumulative
cumulative = (key(1:1) == '(' .and. key(len_trim(key):len_trim(key)) == ')') cumulative = (key(1:1) == '(' .and. key(len_trim(key):len_trim(key)) == ')')
if (present(raw)) then if (present(raw)) then
whole = raw whole = raw
@ -403,14 +403,14 @@ function getStrings(this,key,defaultVal,raw)
whole = .false. whole = .false.
endif endif
found = .false. found = .false.
item => this item => this
do while (associated(item%next)) do while (associated(item%next))
if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
found = .true. found = .true.
if (allocated(getStrings) .and. .not. cumulative) deallocate(getStrings) if (allocated(getStrings) .and. .not. cumulative) deallocate(getStrings)
if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key) if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key)
notAllocated: if (.not. allocated(getStrings)) then notAllocated: if (.not. allocated(getStrings)) then
if (whole) then if (whole) then
str = item%string%val(item%string%pos(4):) str = item%string%val(item%string%pos(4):)
@ -437,7 +437,7 @@ function getStrings(this,key,defaultVal,raw)
endif endif
item => item%next item => item%next
enddo enddo
if (.not. found) then if (.not. found) then
if (present(defaultVal)) then if (present(defaultVal)) then
if (len(defaultVal) > len(getStrings)) call IO_error(0,ext_msg='getStrings') if (len(defaultVal) > len(getStrings)) call IO_error(0,ext_msg='getStrings')

2
src/source_damage_isoDuctile.f90
View File

@ -1,7 +1,7 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH !> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH !> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine incoprorating isotropic ductile damage source mechanism !> @brief material subroutine incorporating isotropic ductile damage source mechanism
!> @details to be done !> @details to be done
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module source_damage_isoDuctile module source_damage_isoDuctile

56
src/system_routines.f90
View File

@ -1,14 +1,14 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief provides wrappers to C routines !> @brief Wrappers to C routines for system operations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module system_routines module system_routines
use, intrinsic :: ISO_C_Binding use, intrinsic :: ISO_C_Binding
use prec use prec
implicit none implicit none
public :: & public :: &
signalterm_C, & signalterm_C, &
signalusr1_C, & signalusr1_C, &
@ -17,74 +17,74 @@ module system_routines
getCWD, & getCWD, &
getHostName, & getHostName, &
setCWD setCWD
interface interface
function isDirectory_C(path) bind(C) function isDirectory_C(path) bind(C)
use, intrinsic :: ISO_C_Binding, only: & use, intrinsic :: ISO_C_Binding, only: &
C_INT, & C_INT, &
C_CHAR C_CHAR
use prec use prec
integer(C_INT) :: isDirectory_C integer(C_INT) :: isDirectory_C
character(kind=C_CHAR), dimension(pPathLen), intent(in) :: path ! C string is an array character(kind=C_CHAR), dimension(pPathLen), intent(in) :: path ! C string is an array
end function isDirectory_C end function isDirectory_C
subroutine getCurrentWorkDir_C(path, stat) bind(C) subroutine getCurrentWorkDir_C(path, stat) bind(C)
use, intrinsic :: ISO_C_Binding, only: & use, intrinsic :: ISO_C_Binding, only: &
C_INT, & C_INT, &
C_CHAR C_CHAR
use prec use prec
character(kind=C_CHAR), dimension(pPathLen), intent(out) :: path ! C string is an array character(kind=C_CHAR), dimension(pPathLen), intent(out) :: path ! C string is an array
integer(C_INT), intent(out) :: stat integer(C_INT), intent(out) :: stat
end subroutine getCurrentWorkDir_C end subroutine getCurrentWorkDir_C
subroutine getHostName_C(str, stat) bind(C) subroutine getHostName_C(str, stat) bind(C)
use, intrinsic :: ISO_C_Binding, only: & use, intrinsic :: ISO_C_Binding, only: &
C_INT, & C_INT, &
C_CHAR C_CHAR
use prec use prec
character(kind=C_CHAR), dimension(pStringLen), intent(out) :: str ! C string is an array character(kind=C_CHAR), dimension(pStringLen), intent(out) :: str ! C string is an array
integer(C_INT), intent(out) :: stat integer(C_INT), intent(out) :: stat
end subroutine getHostName_C end subroutine getHostName_C
function chdir_C(path) bind(C) function chdir_C(path) bind(C)
use, intrinsic :: ISO_C_Binding, only: & use, intrinsic :: ISO_C_Binding, only: &
C_INT, & C_INT, &
C_CHAR C_CHAR
use prec use prec
integer(C_INT) :: chdir_C integer(C_INT) :: chdir_C
character(kind=C_CHAR), dimension(pPathLen), intent(in) :: path ! C string is an array character(kind=C_CHAR), dimension(pPathLen), intent(in) :: path ! C string is an array
end function chdir_C end function chdir_C
subroutine signalterm_C(handler) bind(C) subroutine signalterm_C(handler) bind(C)
use, intrinsic :: ISO_C_Binding, only: & use, intrinsic :: ISO_C_Binding, only: &
C_FUNPTR C_FUNPTR
type(C_FUNPTR), intent(in), value :: handler type(C_FUNPTR), intent(in), value :: handler
end subroutine signalterm_C end subroutine signalterm_C
subroutine signalusr1_C(handler) bind(C) subroutine signalusr1_C(handler) bind(C)
use, intrinsic :: ISO_C_Binding, only: & use, intrinsic :: ISO_C_Binding, only: &
C_FUNPTR C_FUNPTR
type(C_FUNPTR), intent(in), value :: handler type(C_FUNPTR), intent(in), value :: handler
end subroutine signalusr1_C end subroutine signalusr1_C
subroutine signalusr2_C(handler) bind(C) subroutine signalusr2_C(handler) bind(C)
use, intrinsic :: ISO_C_Binding, only: & use, intrinsic :: ISO_C_Binding, only: &
C_FUNPTR C_FUNPTR
type(C_FUNPTR), intent(in), value :: handler type(C_FUNPTR), intent(in), value :: handler
end subroutine signalusr2_C end subroutine signalusr2_C
end interface end interface
contains contains
@ -96,8 +96,8 @@ logical function isDirectory(path)
character(len=*), intent(in) :: path character(len=*), intent(in) :: path
character(kind=C_CHAR), dimension(pPathLen) :: strFixedLength ! C string as array character(kind=C_CHAR), dimension(pPathLen) :: strFixedLength ! C string as array
integer :: i integer :: i
strFixedLength = repeat(C_NULL_CHAR,len(strFixedLength)) strFixedLength = repeat(C_NULL_CHAR,len(strFixedLength))
do i=1,len(path) ! copy array components do i=1,len(path) ! copy array components
strFixedLength(i)=path(i:i) strFixedLength(i)=path(i:i)
@ -116,9 +116,9 @@ function getCWD()
character(len=:), allocatable :: getCWD character(len=:), allocatable :: getCWD
integer(C_INT) :: stat integer(C_INT) :: stat
integer :: i integer :: i
call getCurrentWorkDir_C(charArray,stat) call getCurrentWorkDir_C(charArray,stat)
if (stat /= 0_C_INT) then if (stat /= 0_C_INT) then
getCWD = 'Error occured when getting currend working directory' getCWD = 'Error occured when getting currend working directory'
else else
@ -145,7 +145,7 @@ function getHostName()
character(len=:), allocatable :: getHostName character(len=:), allocatable :: getHostName
integer(C_INT) :: stat integer(C_INT) :: stat
integer :: i integer :: i
call getHostName_C(charArray,stat) call getHostName_C(charArray,stat)
if (stat /= 0_C_INT) then if (stat /= 0_C_INT) then
@ -173,7 +173,7 @@ logical function setCWD(path)
character(len=*), intent(in) :: path character(len=*), intent(in) :: path
character(kind=C_CHAR), dimension(pPathLen) :: strFixedLength ! C string is an array character(kind=C_CHAR), dimension(pPathLen) :: strFixedLength ! C string is an array
integer :: i integer :: i
strFixedLength = repeat(C_NULL_CHAR,len(strFixedLength)) strFixedLength = repeat(C_NULL_CHAR,len(strFixedLength))
do i=1,len(path) ! copy array components do i=1,len(path) ! copy array components
strFixedLength(i)=path(i:i) strFixedLength(i)=path(i:i)

4
src/thermal_conduction.f90
View File

@ -73,7 +73,7 @@ end subroutine thermal_conduction_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns heat generation rate !> @brief return heat generation rate
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
@ -132,7 +132,7 @@ end subroutine thermal_conduction_getSourceAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns homogenized thermal conductivity in reference configuration !> @brief return homogenized thermal conductivity in reference configuration
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function thermal_conduction_getConductivity(ip,el) function thermal_conduction_getConductivity(ip,el)

8
src/thermal_isothermal.f90
View File

@ -8,14 +8,14 @@ module thermal_isothermal
implicit none implicit none
public public
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields, reads information from material configuration file !> @brief allocates fields, reads information from material configuration file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine thermal_isothermal_init subroutine thermal_isothermal_init
integer :: h,NofMyHomog integer :: h,NofMyHomog
write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_isothermal_label//' init -+>>>'; flush(6) write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_isothermal_label//' init -+>>>'; flush(6)
@ -28,7 +28,7 @@ subroutine thermal_isothermal_init
allocate(thermalState(h)%state0 (0,NofMyHomog)) allocate(thermalState(h)%state0 (0,NofMyHomog))
allocate(thermalState(h)%subState0(0,NofMyHomog)) allocate(thermalState(h)%subState0(0,NofMyHomog))
allocate(thermalState(h)%state (0,NofMyHomog)) allocate(thermalState(h)%state (0,NofMyHomog))
deallocate(temperature (h)%p) deallocate(temperature (h)%p)
allocate (temperature (h)%p(1), source=thermal_initialT(h)) allocate (temperature (h)%p(1), source=thermal_initialT(h))
deallocate(temperatureRate(h)%p) deallocate(temperatureRate(h)%p)