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 nonSchmid-adaptation

This commit is contained in:
Philip Eisenlohr 2023-09-26 17:57:30 -04:00
parent 5f2b73f815 96ee2af279
commit f7e06660b1
12 changed files with 566 additions and 472 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
VERSION
View File

@ -1 +1 @@
3.0.0-alpha7-807-g104efaed5 3.0.0-alpha7-837-gc3d3ea658

60
python/damask/_configmaterial.py
View File

@ -106,15 +106,12 @@ class ConfigMaterial(Config):
Load DREAM.3D (HDF5) file. Load DREAM.3D (HDF5) file.
Data in DREAM.3D files can be stored per cell ('CellData') Data in DREAM.3D files can be stored per cell ('CellData')
and/or per grain ('Grain Data'). Per default, cell-wise data and/or per grain ('Grain Data'). Per default, i.e. if
is assumed. 'grain_data' is None, cell-wise data is assumed.
damask.Grid.load_DREAM3D allows to get the corresponding geometry
for the grid solver.
Parameters Parameters
---------- ----------
fname : str fname : str or pathlib.Path
Filename of the DREAM.3D (HDF5) file. Filename of the DREAM.3D (HDF5) file.
grain_data : str grain_data : str
Name of the group (folder) containing grain-wise data. Defaults Name of the group (folder) containing grain-wise data. Defaults
@ -140,36 +137,43 @@ class ConfigMaterial(Config):
and grain- or cell-wise data. Defaults to None, in which case and grain- or cell-wise data. Defaults to None, in which case
it is set as the path that contains _SIMPL_GEOMETRY/SPACING. it is set as the path that contains _SIMPL_GEOMETRY/SPACING.
Notes
-----
Homogenization and phase entries are emtpy and need to be defined separately.
Returns Returns
------- -------
loaded : damask.ConfigMaterial loaded : damask.ConfigMaterial
Material configuration from file. Material configuration from file.
Notes
-----
damask.Grid.load_DREAM3D gives the corresponding geometry for
the grid solver.
For cell-wise data, only unique combinations of
orientation and phase are considered.
Homogenization and phase entries are emtpy and need to be
defined separately.
""" """
b = util.DREAM3D_base_group(fname) if base_group is None else base_group with h5py.File(fname, 'r') as f:
c = util.DREAM3D_cell_data_group(fname) if cell_data is None else cell_data b = util.DREAM3D_base_group(f) if base_group is None else base_group
f = h5py.File(fname,'r') c = util.DREAM3D_cell_data_group(f) if cell_data is None else cell_data
if grain_data is None: if grain_data is None:
phase = f['/'.join([b,c,phases])][()].flatten() phase = f['/'.join([b,c,phases])][()].flatten()
O = Rotation.from_Euler_angles(f['/'.join([b,c,Euler_angles])]).as_quaternion().reshape(-1,4) # noqa O = Rotation.from_Euler_angles(f['/'.join([b,c,Euler_angles])]).as_quaternion().reshape(-1,4) # noqa
_,idx = np.unique(np.hstack([O,phase.reshape(-1,1)]),return_index=True,axis=0) _,idx = np.unique(np.hstack([O,phase.reshape(-1,1)]),return_index=True,axis=0)
idx = np.sort(idx) idx = np.sort(idx)
else: else:
phase = f['/'.join([b,grain_data,phases])][()] phase = f['/'.join([b,grain_data,phases])][()]
O = Rotation.from_Euler_angles(f['/'.join([b,grain_data,Euler_angles])]).as_quaternion() # noqa O = Rotation.from_Euler_angles(f['/'.join([b,grain_data,Euler_angles])]).as_quaternion() # noqa
idx = np.arange(phase.size) idx = np.arange(phase.size)
if cell_ensemble_data is not None and phase_names is not None: if cell_ensemble_data is not None and phase_names is not None:
try: try:
names = np.array([s.decode() for s in f['/'.join([b,cell_ensemble_data,phase_names])]]) names = np.array([s.decode() for s in f['/'.join([b,cell_ensemble_data,phase_names])]])
phase = names[phase] phase = names[phase]
except KeyError: except KeyError:
pass pass
base_config = ConfigMaterial({'phase':{k if isinstance(k,int) else str(k): None for k in np.unique(phase)}, base_config = ConfigMaterial({'phase':{k if isinstance(k,int) else str(k): None for k in np.unique(phase)},

44
python/damask/_grid.py
View File

@ -358,14 +358,14 @@ class Grid:
""" """
Load DREAM.3D (HDF5) file. Load DREAM.3D (HDF5) file.
Data in DREAM.3D files can be stored per cell ('CellData') and/or Data in DREAM.3D files can be stored per cell ('CellData')
per grain ('Grain Data'). Per default, cell-wise data is assumed. and/or per grain ('Grain Data'). Per default, i.e. if
'feature_IDs' is None, cell-wise data is assumed.
damask.ConfigMaterial.load_DREAM3D gives the corresponding material definition.
Parameters Parameters
---------- ----------
fname : str or or pathlib.Path fname : str or pathlib.Path
Filename of the DREAM.3D (HDF5) file. Filename of the DREAM.3D (HDF5) file.
feature_IDs : str, optional feature_IDs : str, optional
Name of the dataset containing the mapping between cells and Name of the dataset containing the mapping between cells and
@ -392,23 +392,31 @@ class Grid:
loaded : damask.Grid loaded : damask.Grid
Grid-based geometry from file. Grid-based geometry from file.
Notes
-----
damask.ConfigMaterial.load_DREAM3D gives the corresponding
material definition.
For cell-wise data, only unique combinations of
orientation and phase are considered.
""" """
b = util.DREAM3D_base_group(fname) if base_group is None else base_group with h5py.File(fname, 'r') as f:
c = util.DREAM3D_cell_data_group(fname) if cell_data is None else cell_data b = util.DREAM3D_base_group(f) if base_group is None else base_group
f = h5py.File(fname, 'r') c = util.DREAM3D_cell_data_group(f) if cell_data is None else cell_data
cells = f['/'.join([b,'_SIMPL_GEOMETRY','DIMENSIONS'])][()] cells = f['/'.join([b,'_SIMPL_GEOMETRY','DIMENSIONS'])][()]
size = f['/'.join([b,'_SIMPL_GEOMETRY','SPACING'])] * cells size = f['/'.join([b,'_SIMPL_GEOMETRY','SPACING'])] * cells
origin = f['/'.join([b,'_SIMPL_GEOMETRY','ORIGIN'])][()] origin = f['/'.join([b,'_SIMPL_GEOMETRY','ORIGIN'])][()]
if feature_IDs is None: if feature_IDs is None:
phase = f['/'.join([b,c,phases])][()].reshape(-1,1) phase = f['/'.join([b,c,phases])][()].reshape(-1,1)
O = Rotation.from_Euler_angles(f['/'.join([b,c,Euler_angles])]).as_quaternion().reshape(-1,4) # noqa O = Rotation.from_Euler_angles(f['/'.join([b,c,Euler_angles])]).as_quaternion().reshape(-1,4) # noqa
unique,unique_inverse = np.unique(np.hstack([O,phase]),return_inverse=True,axis=0) unique,unique_inverse = np.unique(np.hstack([O,phase]),return_inverse=True,axis=0)
ma = np.arange(cells.prod()) if len(unique) == cells.prod() else \ ma = np.arange(cells.prod()) if len(unique) == cells.prod() else \
np.arange(unique.size)[np.argsort(pd.unique(unique_inverse))][unique_inverse] np.arange(unique.size)[np.argsort(pd.unique(unique_inverse))][unique_inverse]
else: else:
ma = f['/'.join([b,c,feature_IDs])][()].flatten() ma = f['/'.join([b,c,feature_IDs])][()].flatten()
return Grid(material = ma.reshape(cells,order='F'), return Grid(material = ma.reshape(cells,order='F'),
size = size, size = size,

4
python/damask/_orientation.py
View File

@ -804,7 +804,7 @@ class Orientation(Rotation,Crystal):
blend += sym_ops.shape blend += sym_ops.shape
v = sym_ops.broadcast_to(shape) \ v = sym_ops.broadcast_to(shape) \
@ np.broadcast_to(v.reshape(util.shapeshifter(v.shape,shape+(3,))),shape+(3,)) @ np.broadcast_to(v.reshape(util.shapeshifter(v.shape,shape+(3,))),shape+(3,))
return ~(self.broadcast_to(blend))@ np.broadcast_to(v,blend+(3,)) return ~(self.broadcast_to(blend))@np.broadcast_to(v,blend+(3,))
def Schmid(self, *, def Schmid(self, *,
@ -833,7 +833,7 @@ class Orientation(Rotation,Crystal):
>>> import damask >>> import damask
>>> np.set_printoptions(3,suppress=True,floatmode='fixed') >>> np.set_printoptions(3,suppress=True,floatmode='fixed')
>>> O = damask.Orientation.from_Euler_angles(phi=[0,45,0],degrees=True,lattice='cF') >>> O = damask.Orientation.from_Euler_angles(phi=[0,45,0],degrees=True,lattice='cF')
>>> O.Schmid(N_slip=[1]) >>> O.Schmid(N_slip=[12])[0]
array([[ 0.000, 0.000, 0.000], array([[ 0.000, 0.000, 0.000],
[ 0.577, -0.000, 0.816], [ 0.577, -0.000, 0.816],
[ 0.000, 0.000, 0.000]]) [ 0.000, 0.000, 0.000]])

703
python/damask/_result.py
View File

@ -1,5 +1,3 @@
import multiprocessing as mp
from multiprocessing.synchronize import Lock
import re import re
import fnmatch import fnmatch
import os import os
@ -7,8 +5,8 @@ import copy
import datetime import datetime
import xml.etree.ElementTree as ET # noqa import xml.etree.ElementTree as ET # noqa
import xml.dom.minidom import xml.dom.minidom
import functools
from pathlib import Path from pathlib import Path
from functools import partial
from collections import defaultdict from collections import defaultdict
from collections.abc import Iterable from collections.abc import Iterable
from typing import Optional, Union, Callable, Any, Sequence, Literal, Dict, List, Tuple from typing import Optional, Union, Callable, Any, Sequence, Literal, Dict, List, Tuple
@ -601,17 +599,6 @@ class Result:
f['/geometry/T_c'].attrs['VTK_TYPE'].decode()) f['/geometry/T_c'].attrs['VTK_TYPE'].decode())
@staticmethod
def _add_absolute(x: Dict[str, Any]) -> Dict[str, Any]:
return {
'data': np.abs(x['data']),
'label': f'|{x["label"]}|',
'meta': {
'unit': x['meta']['unit'],
'description': f"absolute value of {x['label']} ({x['meta']['description']})",
'creator': 'add_absolute'
}
}
def add_absolute(self, x: str): def add_absolute(self, x: str):
""" """
Add absolute value. Add absolute value.
@ -622,28 +609,20 @@ class Result:
Name of scalar, vector, or tensor dataset to take absolute value of. Name of scalar, vector, or tensor dataset to take absolute value of.
""" """
self._add_generic_pointwise(self._add_absolute,{'x':x}) def absolute(x: Dict[str, Any]) -> Dict[str, Any]:
return {
'data': np.abs(x['data']),
'label': f'|{x["label"]}|',
'meta': {
'unit': x['meta']['unit'],
'description': f"absolute value of {x['label']} ({x['meta']['description']})",
'creator': 'add_absolute'
}
}
self._add_generic_pointwise(absolute,{'x':x})
@staticmethod
def _add_calculation(**kwargs) -> Dict[str, Any]:
formula = kwargs['formula']
for d in re.findall(r'#(.*?)#',formula):
formula = formula.replace(f'#{d}#',f"kwargs['{d}']['data']")
data = eval(formula)
if not hasattr(data,'shape') or data.shape[0] != kwargs[d]['data'].shape[0]:
raise ValueError('"{}" results in invalid shape'.format(kwargs['formula']))
return {
'data': data,
'label': kwargs['label'],
'meta': {
'unit': kwargs['unit'],
'description': f"{kwargs['description']} (formula: {kwargs['formula']})",
'creator': 'add_calculation'
}
}
def add_calculation(self, def add_calculation(self,
formula: str, formula: str,
name: str, name: str,
@ -692,24 +671,30 @@ class Result:
... 'Mises equivalent of the Cauchy stress') ... 'Mises equivalent of the Cauchy stress')
""" """
def calculation(**kwargs) -> Dict[str, Any]:
formula = kwargs['formula']
for d in re.findall(r'#(.*?)#',formula):
formula = formula.replace(f'#{d}#',f"kwargs['{d}']['data']")
data = eval(formula)
if not hasattr(data,'shape') or data.shape[0] != kwargs[d]['data'].shape[0]:
raise ValueError('"{}" results in invalid shape'.format(kwargs['formula']))
return {
'data': data,
'label': kwargs['label'],
'meta': {
'unit': kwargs['unit'],
'description': f"{kwargs['description']} (formula: {kwargs['formula']})",
'creator': 'add_calculation'
}
}
dataset_mapping = {d:d for d in set(re.findall(r'#(.*?)#',formula))} # datasets used in the formula dataset_mapping = {d:d for d in set(re.findall(r'#(.*?)#',formula))} # datasets used in the formula
args = {'formula':formula,'label':name,'unit':unit,'description':description} args = {'formula':formula,'label':name,'unit':unit,'description':description}
self._add_generic_pointwise(self._add_calculation,dataset_mapping,args) self._add_generic_pointwise(calculation,dataset_mapping,args)
@staticmethod
def _add_stress_Cauchy(P: Dict[str, Any], F: Dict[str, Any]) -> Dict[str, Any]:
return {
'data': mechanics.stress_Cauchy(P['data'],F['data']),
'label': 'sigma',
'meta': {
'unit': P['meta']['unit'],
'description': "Cauchy stress calculated "
f"from {P['label']} ({P['meta']['description']})"
f" and {F['label']} ({F['meta']['description']})",
'creator': 'add_stress_Cauchy'
}
}
def add_stress_Cauchy(self, def add_stress_Cauchy(self,
P: str = 'P', P: str = 'P',
F: str = 'F'): F: str = 'F'):
@ -726,20 +711,23 @@ class Result:
Defaults to 'F'. Defaults to 'F'.
""" """
self._add_generic_pointwise(self._add_stress_Cauchy,{'P':P,'F':F})
def stress_Cauchy(P: Dict[str, Any], F: Dict[str, Any]) -> Dict[str, Any]:
return {
'data': mechanics.stress_Cauchy(P['data'],F['data']),
'label': 'sigma',
'meta': {
'unit': P['meta']['unit'],
'description': "Cauchy stress calculated "
f"from {P['label']} ({P['meta']['description']})"
f" and {F['label']} ({F['meta']['description']})",
'creator': 'add_stress_Cauchy'
}
}
self._add_generic_pointwise(stress_Cauchy,{'P':P,'F':F})
@staticmethod
def _add_determinant(T: Dict[str, Any]) -> Dict[str, Any]:
return {
'data': np.linalg.det(T['data']),
'label': f"det({T['label']})",
'meta': {
'unit': T['meta']['unit'],
'description': f"determinant of tensor {T['label']} ({T['meta']['description']})",
'creator': 'add_determinant'
}
}
def add_determinant(self, T: str): def add_determinant(self, T: str):
""" """
Add the determinant of a tensor. Add the determinant of a tensor.
@ -758,20 +746,21 @@ class Result:
>>> r.add_determinant('F_p') >>> r.add_determinant('F_p')
""" """
self._add_generic_pointwise(self._add_determinant,{'T':T})
def determinant(T: Dict[str, Any]) -> Dict[str, Any]:
return {
'data': np.linalg.det(T['data']),
'label': f"det({T['label']})",
'meta': {
'unit': T['meta']['unit'],
'description': f"determinant of tensor {T['label']} ({T['meta']['description']})",
'creator': 'add_determinant'
}
}
self._add_generic_pointwise(determinant,{'T':T})
@staticmethod
def _add_deviator(T: Dict[str, Any]) -> Dict[str, Any]:
return {
'data': tensor.deviatoric(T['data']),
'label': f"s_{T['label']}",
'meta': {
'unit': T['meta']['unit'],
'description': f"deviator of tensor {T['label']} ({T['meta']['description']})",
'creator': 'add_deviator'
}
}
def add_deviator(self, T: str): def add_deviator(self, T: str):
""" """
Add the deviatoric part of a tensor. Add the deviatoric part of a tensor.
@ -790,29 +779,21 @@ class Result:
>>> r.add_deviator('sigma') >>> r.add_deviator('sigma')
""" """
self._add_generic_pointwise(self._add_deviator,{'T':T})
def deviator(T: Dict[str, Any]) -> Dict[str, Any]:
return {
'data': tensor.deviatoric(T['data']),
'label': f"s_{T['label']}",
'meta': {
'unit': T['meta']['unit'],
'description': f"deviator of tensor {T['label']} ({T['meta']['description']})",
'creator': 'add_deviator'
}
}
self._add_generic_pointwise(deviator,{'T':T})
@staticmethod
def _add_eigenvalue(T_sym: Dict[str, Any], eigenvalue: Literal['max, mid, min']) -> Dict[str, Any]:
if eigenvalue == 'max':
label,p = 'maximum',2
elif eigenvalue == 'mid':
label,p = 'intermediate',1
elif eigenvalue == 'min':
label,p = 'minimum',0
else:
raise ValueError(f'invalid eigenvalue: {eigenvalue}')
return {
'data': tensor.eigenvalues(T_sym['data'])[:,p],
'label': f"lambda_{eigenvalue}({T_sym['label']})",
'meta' : {
'unit': T_sym['meta']['unit'],
'description': f"{label} eigenvalue of {T_sym['label']} ({T_sym['meta']['description']})",
'creator': 'add_eigenvalue'
}
}
def add_eigenvalue(self, def add_eigenvalue(self,
T_sym: str, T_sym: str,
eigenvalue: Literal['max', 'mid', 'min'] = 'max'): eigenvalue: Literal['max', 'mid', 'min'] = 'max'):
@ -835,30 +816,30 @@ class Result:
>>> r.add_eigenvalue('sigma','min') >>> r.add_eigenvalue('sigma','min')
""" """
self._add_generic_pointwise(self._add_eigenvalue,{'T_sym':T_sym},{'eigenvalue':eigenvalue})
def eigenval(T_sym: Dict[str, Any], eigenvalue: Literal['max, mid, min']) -> Dict[str, Any]:
if eigenvalue == 'max':
label,p = 'maximum',2
elif eigenvalue == 'mid':
label,p = 'intermediate',1
elif eigenvalue == 'min':
label,p = 'minimum',0
else:
raise ValueError(f'invalid eigenvalue: {eigenvalue}')
return {
'data': tensor.eigenvalues(T_sym['data'])[:,p],
'label': f"lambda_{eigenvalue}({T_sym['label']})",
'meta' : {
'unit': T_sym['meta']['unit'],
'description': f"{label} eigenvalue of {T_sym['label']} ({T_sym['meta']['description']})",
'creator': 'add_eigenvalue'
}
}
self._add_generic_pointwise(eigenval,{'T_sym':T_sym},{'eigenvalue':eigenvalue})
@staticmethod
def _add_eigenvector(T_sym: Dict[str, Any], eigenvalue: Literal['max', 'mid', 'min']) -> Dict[str, Any]:
if eigenvalue == 'max':
label,p = 'maximum',2
elif eigenvalue == 'mid':
label,p = 'intermediate',1
elif eigenvalue == 'min':
label,p = 'minimum',0
else:
raise ValueError(f'invalid eigenvalue: {eigenvalue}')
return {
'data': tensor.eigenvectors(T_sym['data'])[:,p],
'label': f"v_{eigenvalue}({T_sym['label']})",
'meta' : {
'unit': '1',
'description': f"eigenvector corresponding to {label} eigenvalue"
f" of {T_sym['label']} ({T_sym['meta']['description']})",
'creator': 'add_eigenvector'
}
}
def add_eigenvector(self, def add_eigenvector(self,
T_sym: str, T_sym: str,
eigenvalue: Literal['max', 'mid', 'min'] = 'max'): eigenvalue: Literal['max', 'mid', 'min'] = 'max'):
@ -874,25 +855,31 @@ class Result:
Defaults to 'max'. Defaults to 'max'.
""" """
self._add_generic_pointwise(self._add_eigenvector,{'T_sym':T_sym},{'eigenvalue':eigenvalue})
def eigenvector(T_sym: Dict[str, Any], eigenvalue: Literal['max', 'mid', 'min']) -> Dict[str, Any]:
if eigenvalue == 'max':
label,p = 'maximum',2
elif eigenvalue == 'mid':
label,p = 'intermediate',1
elif eigenvalue == 'min':
label,p = 'minimum',0
else:
raise ValueError(f'invalid eigenvalue: {eigenvalue}')
return {
'data': tensor.eigenvectors(T_sym['data'])[:,p],
'label': f"v_{eigenvalue}({T_sym['label']})",
'meta' : {
'unit': '1',
'description': f"eigenvector corresponding to {label} eigenvalue"
f" of {T_sym['label']} ({T_sym['meta']['description']})",
'creator': 'add_eigenvector'
}
}
self._add_generic_pointwise(eigenvector,{'T_sym':T_sym},{'eigenvalue':eigenvalue})
@staticmethod
def _add_IPF_color(l: FloatSequence, q: Dict[str, Any]) -> Dict[str, Any]:
m = util.scale_to_coprime(np.array(l))
lattice = q['meta']['lattice']
o = Orientation(rotation = q['data'],lattice=lattice)
return {
'data': np.uint8(o.IPF_color(l)*255),
'label': 'IPFcolor_({} {} {})'.format(*m),
'meta' : {
'unit': '8-bit RGB',
'lattice': q['meta']['lattice'],
'description': 'Inverse Pole Figure (IPF) colors along sample direction ({} {} {})'.format(*m),
'creator': 'add_IPF_color'
}
}
def add_IPF_color(self, def add_IPF_color(self,
l: FloatSequence, l: FloatSequence,
q: str = 'O'): q: str = 'O'):
@ -916,20 +903,26 @@ class Result:
>>> r.add_IPF_color(np.array([0,1,1])) >>> r.add_IPF_color(np.array([0,1,1]))
""" """
self._add_generic_pointwise(self._add_IPF_color,{'q':q},{'l':l})
def IPF_color(l: FloatSequence, q: Dict[str, Any]) -> Dict[str, Any]:
m = util.scale_to_coprime(np.array(l))
lattice = q['meta']['lattice']
o = Orientation(rotation = q['data'],lattice=lattice)
return {
'data': np.uint8(o.IPF_color(l)*255),
'label': 'IPFcolor_({} {} {})'.format(*m),
'meta' : {
'unit': '8-bit RGB',
'lattice': q['meta']['lattice'],
'description': 'Inverse Pole Figure (IPF) colors along sample direction ({} {} {})'.format(*m),
'creator': 'add_IPF_color'
}
}
self._add_generic_pointwise(IPF_color,{'q':q},{'l':l})
@staticmethod
def _add_maximum_shear(T_sym: Dict[str, Any]) -> Dict[str, Any]:
return {
'data': mechanics.maximum_shear(T_sym['data']),
'label': f"max_shear({T_sym['label']})",
'meta': {
'unit': T_sym['meta']['unit'],
'description': f"maximum shear component of {T_sym['label']} ({T_sym['meta']['description']})",
'creator': 'add_maximum_shear'
}
}
def add_maximum_shear(self, T_sym: str): def add_maximum_shear(self, T_sym: str):
""" """
Add maximum shear components of symmetric tensor. Add maximum shear components of symmetric tensor.
@ -940,30 +933,20 @@ class Result:
Name of symmetric tensor dataset. Name of symmetric tensor dataset.
""" """
self._add_generic_pointwise(self._add_maximum_shear,{'T_sym':T_sym}) def maximum_shear(T_sym: Dict[str, Any]) -> Dict[str, Any]:
return {
'data': mechanics.maximum_shear(T_sym['data']),
'label': f"max_shear({T_sym['label']})",
'meta': {
'unit': T_sym['meta']['unit'],
'description': f"maximum shear component of {T_sym['label']} ({T_sym['meta']['description']})",
'creator': 'add_maximum_shear'
}
}
self._add_generic_pointwise(maximum_shear,{'T_sym':T_sym})
@staticmethod
def _add_equivalent_Mises(T_sym: Dict[str, Any], kind: str) -> Dict[str, Any]:
k = kind
if k is None:
if T_sym['meta']['unit'] == '1':
k = 'strain'
elif T_sym['meta']['unit'] == 'Pa':
k = 'stress'
if k not in ['stress', 'strain']:
raise ValueError(f'invalid von Mises kind "{kind}"')
return {
'data': (mechanics.equivalent_strain_Mises if k=='strain' else \
mechanics.equivalent_stress_Mises)(T_sym['data']),
'label': f"{T_sym['label']}_vM",
'meta': {
'unit': T_sym['meta']['unit'],
'description': f"Mises equivalent {k} of {T_sym['label']} ({T_sym['meta']['description']})",
'creator': 'add_Mises'
}
}
def add_equivalent_Mises(self, def add_equivalent_Mises(self,
T_sym: str, T_sym: str,
kind: Optional[str] = None): kind: Optional[str] = None):
@ -993,32 +976,30 @@ class Result:
>>> r.add_equivalent_Mises('epsilon_V^0.0(F)') >>> r.add_equivalent_Mises('epsilon_V^0.0(F)')
""" """
self._add_generic_pointwise(self._add_equivalent_Mises,{'T_sym':T_sym},{'kind':kind}) def equivalent_Mises(T_sym: Dict[str, Any], kind: str) -> Dict[str, Any]:
k = kind
if k is None:
if T_sym['meta']['unit'] == '1':
k = 'strain'
elif T_sym['meta']['unit'] == 'Pa':
k = 'stress'
if k not in ['stress', 'strain']:
raise ValueError(f'invalid von Mises kind "{kind}"')
return {
'data': (mechanics.equivalent_strain_Mises if k=='strain' else \
mechanics.equivalent_stress_Mises)(T_sym['data']),
'label': f"{T_sym['label']}_vM",
'meta': {
'unit': T_sym['meta']['unit'],
'description': f"Mises equivalent {k} of {T_sym['label']} ({T_sym['meta']['description']})",
'creator': 'add_Mises'
}
}
self._add_generic_pointwise(equivalent_Mises,{'T_sym':T_sym},{'kind':kind})
@staticmethod
def _add_norm(x: Dict[str, Any], ord: Union[int, float, Literal['fro', 'nuc']]) -> Dict[str, Any]:
o = ord
if len(x['data'].shape) == 2:
axis: Union[int, Tuple[int, int]] = 1
t = 'vector'
if o is None: o = 2
elif len(x['data'].shape) == 3:
axis = (1,2)
t = 'tensor'
if o is None: o = 'fro'
else:
raise ValueError(f'invalid shape of {x["label"]}')
return {
'data': np.linalg.norm(x['data'],ord=o,axis=axis,keepdims=True),
'label': f"|{x['label']}|_{o}",
'meta': {
'unit': x['meta']['unit'],
'description': f"{o}-norm of {t} {x['label']} ({x['meta']['description']})",
'creator': 'add_norm'
}
}
def add_norm(self, def add_norm(self,
x: str, x: str,
ord: Union[None, int, float, Literal['fro', 'nuc']] = None): ord: Union[None, int, float, Literal['fro', 'nuc']] = None):
@ -1033,22 +1014,32 @@ class Result:
Order of the norm. inf means NumPy's inf object. For details refer to numpy.linalg.norm. Order of the norm. inf means NumPy's inf object. For details refer to numpy.linalg.norm.
""" """
self._add_generic_pointwise(self._add_norm,{'x':x},{'ord':ord}) def norm(x: Dict[str, Any], ord: Union[int, float, Literal['fro', 'nuc']]) -> Dict[str, Any]:
o = ord
if len(x['data'].shape) == 2:
axis: Union[int, Tuple[int, int]] = 1
t = 'vector'
if o is None: o = 2
elif len(x['data'].shape) == 3:
axis = (1,2)
t = 'tensor'
if o is None: o = 'fro'
else:
raise ValueError(f'invalid shape of {x["label"]}')
return {
'data': np.linalg.norm(x['data'],ord=o,axis=axis,keepdims=True),
'label': f"|{x['label']}|_{o}",
'meta': {
'unit': x['meta']['unit'],
'description': f"{o}-norm of {t} {x['label']} ({x['meta']['description']})",
'creator': 'add_norm'
}
}
self._add_generic_pointwise(norm,{'x':x},{'ord':ord})
@staticmethod
def _add_stress_second_Piola_Kirchhoff(P: Dict[str, Any], F: Dict[str, Any]) -> Dict[str, Any]:
return {
'data': mechanics.stress_second_Piola_Kirchhoff(P['data'],F['data']),
'label': 'S',
'meta': {
'unit': P['meta']['unit'],
'description': "second Piola-Kirchhoff stress calculated "
f"from {P['label']} ({P['meta']['description']})"
f" and {F['label']} ({F['meta']['description']})",
'creator': 'add_stress_second_Piola_Kirchhoff'
}
}
def add_stress_second_Piola_Kirchhoff(self, def add_stress_second_Piola_Kirchhoff(self,
P: str = 'P', P: str = 'P',
F: str = 'F'): F: str = 'F'):
@ -1071,34 +1062,23 @@ class Result:
is taken into account. is taken into account.
""" """
self._add_generic_pointwise(self._add_stress_second_Piola_Kirchhoff,{'P':P,'F':F}) def stress_second_Piola_Kirchhoff(P: Dict[str, Any], F: Dict[str, Any]) -> Dict[str, Any]:
return {
'data': mechanics.stress_second_Piola_Kirchhoff(P['data'],F['data']),
'label': 'S',
'meta': {
'unit': P['meta']['unit'],
'description': "second Piola-Kirchhoff stress calculated "
f"from {P['label']} ({P['meta']['description']})"
f" and {F['label']} ({F['meta']['description']})",
'creator': 'add_stress_second_Piola_Kirchhoff'
}
}
self._add_generic_pointwise(stress_second_Piola_Kirchhoff,{'P':P,'F':F})
@staticmethod
def _add_pole(q: Dict[str, Any],
uvw: FloatSequence,
hkl: FloatSequence,
with_symmetry: bool,
normalize: bool) -> Dict[str, Any]:
c = q['meta']['c/a'] if 'c/a' in q['meta'] else 1
brackets = ['[]','()','⟨⟩','{}'][(uvw is None)*1+with_symmetry*2]
label = 'p^' + '{}{} {} {}{}'.format(brackets[0],
*(uvw if uvw else hkl),
brackets[-1],)
ori = Orientation(q['data'],lattice=q['meta']['lattice'],a=1,c=c)
return {
'data': ori.to_pole(uvw=uvw,hkl=hkl,with_symmetry=with_symmetry,normalize=normalize),
'label': label,
'meta' : {
'unit': '1',
'description': f'{"normalized " if normalize else ""}lab frame vector along lattice ' \
+ ('direction' if uvw is not None else 'plane') \
+ ('s' if with_symmetry else ''),
'creator': 'add_pole'
}
}
def add_pole(self, def add_pole(self,
q: str = 'O', q: str = 'O',
*, *,
@ -1124,22 +1104,33 @@ class Result:
Defaults to True. Defaults to True.
""" """
self._add_generic_pointwise(self._add_pole, def pole(q: Dict[str, Any],
{'q':q}, uvw: FloatSequence,
{'uvw':uvw,'hkl':hkl,'with_symmetry':with_symmetry,'normalize':normalize}) hkl: FloatSequence,
with_symmetry: bool,
normalize: bool) -> Dict[str, Any]:
c = q['meta']['c/a'] if 'c/a' in q['meta'] else 1
brackets = ['[]','()','⟨⟩','{}'][(uvw is None)*1+with_symmetry*2]
label = 'p^' + '{}{} {} {}{}'.format(brackets[0],
*(uvw if uvw else hkl),
brackets[-1],)
ori = Orientation(q['data'],lattice=q['meta']['lattice'],a=1,c=c)
return {
'data': ori.to_pole(uvw=uvw,hkl=hkl,with_symmetry=with_symmetry,normalize=normalize),
'label': label,
'meta' : {
'unit': '1',
'description': f'{"normalized " if normalize else ""}lab frame vector along lattice ' \
+ ('direction' if uvw is not None else 'plane') \
+ ('s' if with_symmetry else ''),
'creator': 'add_pole'
}
}
self._add_generic_pointwise(pole,{'q':q},{'uvw':uvw,'hkl':hkl,'with_symmetry':with_symmetry,'normalize':normalize})
@staticmethod
def _add_rotation(F: Dict[str, Any]) -> Dict[str, Any]:
return {
'data': mechanics.rotation(F['data']).as_matrix(),
'label': f"R({F['label']})",
'meta': {
'unit': F['meta']['unit'],
'description': f"rotational part of {F['label']} ({F['meta']['description']})",
'creator': 'add_rotation'
}
}
def add_rotation(self, F: str): def add_rotation(self, F: str):
""" """
Add rotational part of a deformation gradient. Add rotational part of a deformation gradient.
@ -1158,20 +1149,20 @@ class Result:
>>> r.add_rotation('F') >>> r.add_rotation('F')
""" """
self._add_generic_pointwise(self._add_rotation,{'F':F}) def rotation(F: Dict[str, Any]) -> Dict[str, Any]:
return {
'data': mechanics.rotation(F['data']).as_matrix(),
'label': f"R({F['label']})",
'meta': {
'unit': F['meta']['unit'],
'description': f"rotational part of {F['label']} ({F['meta']['description']})",
'creator': 'add_rotation'
}
}
self._add_generic_pointwise(rotation,{'F':F})
@staticmethod
def _add_spherical(T: Dict[str, Any]) -> Dict[str, Any]:
return {
'data': tensor.spherical(T['data'],False),
'label': f"p_{T['label']}",
'meta': {
'unit': T['meta']['unit'],
'description': f"spherical component of tensor {T['label']} ({T['meta']['description']})",
'creator': 'add_spherical'
}
}
def add_spherical(self, T: str): def add_spherical(self, T: str):
""" """
Add the spherical (hydrostatic) part of a tensor. Add the spherical (hydrostatic) part of a tensor.
@ -1190,22 +1181,20 @@ class Result:
>>> r.add_spherical('sigma') >>> r.add_spherical('sigma')
""" """
self._add_generic_pointwise(self._add_spherical,{'T':T}) def spherical(T: Dict[str, Any]) -> Dict[str, Any]:
return {
'data': tensor.spherical(T['data'],False),
'label': f"p_{T['label']}",
'meta': {
'unit': T['meta']['unit'],
'description': f"spherical component of tensor {T['label']} ({T['meta']['description']})",
'creator': 'add_spherical'
}
}
self._add_generic_pointwise(spherical,{'T':T})
@staticmethod
def _add_strain(F: Dict[str, Any], t: Literal['V', 'U'], m: float) -> Dict[str, Any]:
side = 'left' if t == 'V' else 'right'
return {
'data': mechanics.strain(F['data'],t,m),
'label': f"epsilon_{t}^{m}({F['label']})",
'meta': {
'unit': F['meta']['unit'],
'description': f'Seth-Hill strain tensor of order {m} based on {side} stretch tensor '+\
f"of {F['label']} ({F['meta']['description']})",
'creator': 'add_strain'
}
}
def add_strain(self, def add_strain(self,
F: str = 'F', F: str = 'F',
t: Literal['V', 'U'] = 'V', t: Literal['V', 'U'] = 'V',
@ -1266,21 +1255,22 @@ class Result:
| https://de.wikipedia.org/wiki/Verzerrungstensor | https://de.wikipedia.org/wiki/Verzerrungstensor
""" """
self._add_generic_pointwise(self._add_strain,{'F':F},{'t':t,'m':m}) def strain(F: Dict[str, Any], t: Literal['V', 'U'], m: float) -> Dict[str, Any]:
side = 'left' if t == 'V' else 'right'
return {
'data': mechanics.strain(F['data'],t,m),
'label': f"epsilon_{t}^{m}({F['label']})",
'meta': {
'unit': F['meta']['unit'],
'description': f'Seth-Hill strain tensor of order {m} based on {side} stretch tensor '+\
f"of {F['label']} ({F['meta']['description']})",
'creator': 'add_strain'
}
}
self._add_generic_pointwise(strain,{'F':F},{'t':t,'m':m})
@staticmethod
def _add_stretch_tensor(F: Dict[str, Any], t: str) -> Dict[str, Any]:
return {
'data': (mechanics.stretch_left if t.upper() == 'V' else mechanics.stretch_right)(F['data']),
'label': f"{t}({F['label']})",
'meta': {
'unit': F['meta']['unit'],
'description': f"{'left' if t.upper() == 'V' else 'right'} stretch tensor "\
+f"of {F['label']} ({F['meta']['description']})", # noqa
'creator': 'add_stretch_tensor'
}
}
def add_stretch_tensor(self, def add_stretch_tensor(self,
F: str = 'F', F: str = 'F',
t: Literal['V', 'U'] = 'V'): t: Literal['V', 'U'] = 'V'):
@ -1296,20 +1286,21 @@ class Result:
Defaults to 'V'. Defaults to 'V'.
""" """
self._add_generic_pointwise(self._add_stretch_tensor,{'F':F},{'t':t}) def stretch_tensor(F: Dict[str, Any], t: str) -> Dict[str, Any]:
return {
'data': (mechanics.stretch_left if t.upper() == 'V' else mechanics.stretch_right)(F['data']),
'label': f"{t}({F['label']})",
'meta': {
'unit': F['meta']['unit'],
'description': f"{'left' if t.upper() == 'V' else 'right'} stretch tensor "\
+f"of {F['label']} ({F['meta']['description']})", # noqa
'creator': 'add_stretch_tensor'
}
}
self._add_generic_pointwise(stretch_tensor,{'F':F},{'t':t})
@staticmethod
def _add_curl(f: Dict[str, Any], size: np.ndarray) -> Dict[str, Any]:
return {
'data': grid_filters.curl(size,f['data']),
'label': f"curl({f['label']})",
'meta': {
'unit': f['meta']['unit']+'/m',
'description': f"curl of {f['label']} ({f['meta']['description']})",
'creator': 'add_curl'
}
}
def add_curl(self, f: str): def add_curl(self, f: str):
""" """
Add curl of a field. Add curl of a field.
@ -1325,20 +1316,20 @@ class Result:
i.e. fields resulting from the grid solver. i.e. fields resulting from the grid solver.
""" """
self._add_generic_grid(self._add_curl,{'f':f},{'size':self.size}) def curl(f: Dict[str, Any], size: np.ndarray) -> Dict[str, Any]:
return {
'data': grid_filters.curl(size,f['data']),
'label': f"curl({f['label']})",
'meta': {
'unit': f['meta']['unit']+'/m',
'description': f"curl of {f['label']} ({f['meta']['description']})",
'creator': 'add_curl'
}
}
self._add_generic_grid(curl,{'f':f},{'size':self.size})
@staticmethod
def _add_divergence(f: Dict[str, Any], size: np.ndarray) -> Dict[str, Any]:
return {
'data': grid_filters.divergence(size,f['data']),
'label': f"divergence({f['label']})",
'meta': {
'unit': f['meta']['unit']+'/m',
'description': f"divergence of {f['label']} ({f['meta']['description']})",
'creator': 'add_divergence'
}
}
def add_divergence(self, f: str): def add_divergence(self, f: str):
""" """
Add divergence of a field. Add divergence of a field.
@ -1354,21 +1345,20 @@ class Result:
i.e. fields resulting from the grid solver. i.e. fields resulting from the grid solver.
""" """
self._add_generic_grid(self._add_divergence,{'f':f},{'size':self.size}) def divergence(f: Dict[str, Any], size: np.ndarray) -> Dict[str, Any]:
return {
'data': grid_filters.divergence(size,f['data']),
'label': f"divergence({f['label']})",
'meta': {
'unit': f['meta']['unit']+'/m',
'description': f"divergence of {f['label']} ({f['meta']['description']})",
'creator': 'add_divergence'
}
}
self._add_generic_grid(divergence,{'f':f},{'size':self.size})
@staticmethod
def _add_gradient(f: Dict[str, Any], size: np.ndarray) -> Dict[str, Any]:
return {
'data': grid_filters.gradient(size,f['data'] if len(f['data'].shape) == 4 else \
f['data'].reshape(f['data'].shape+(1,))),
'label': f"gradient({f['label']})",
'meta': {
'unit': f['meta']['unit']+'/m',
'description': f"gradient of {f['label']} ({f['meta']['description']})",
'creator': 'add_gradient'
}
}
def add_gradient(self, f: str): def add_gradient(self, f: str):
""" """
Add gradient of a field. Add gradient of a field.
@ -1384,7 +1374,19 @@ class Result:
i.e. fields resulting from the grid solver. i.e. fields resulting from the grid solver.
""" """
self._add_generic_grid(self._add_gradient,{'f':f},{'size':self.size}) def gradient(f: Dict[str, Any], size: np.ndarray) -> Dict[str, Any]:
return {
'data': grid_filters.gradient(size,f['data'] if len(f['data'].shape) == 4 else \
f['data'].reshape(f['data'].shape+(1,))),
'label': f"gradient({f['label']})",
'meta': {
'unit': f['meta']['unit']+'/m',
'description': f"gradient of {f['label']} ({f['meta']['description']})",
'creator': 'add_gradient'
}
}
self._add_generic_grid(gradient,{'f':f},{'size':self.size})
def _add_generic_grid(self, def _add_generic_grid(self,
@ -1446,29 +1448,6 @@ class Result:
f'damask.Result.{creator} v{damask.version}'.encode() f'damask.Result.{creator} v{damask.version}'.encode()
def _job_pointwise(self,
group: str,
callback: Callable,
datasets: Dict[str, str],
args: Dict[str, str],
lock: Lock) -> List[Union[None, Any]]:
"""Execute job for _add_generic_pointwise."""
try:
datasets_in = {}
lock.acquire()
with h5py.File(self.fname,'r') as f:
for arg,label in datasets.items():
loc = f[group+'/'+label]
datasets_in[arg]={'data' :loc[()],
'label':label,
'meta': {k:(v.decode() if not h5py3 and type(v) is bytes else v) \
for k,v in loc.attrs.items()}}
lock.release()
r = callback(**datasets_in,**args)
return [group,r]
except Exception as err:
print(f'Error during calculation: {err}.')
return [None,None]
def _add_generic_pointwise(self, def _add_generic_pointwise(self,
@ -1490,8 +1469,24 @@ class Result:
Arguments parsed to func. Arguments parsed to func.
""" """
pool = mp.Pool(int(os.environ.get('OMP_NUM_THREADS',4)))
lock = mp.Manager().Lock() def job_pointwise(group: str,
callback: Callable,
datasets: Dict[str, str],
args: Dict[str, str]) -> Union[None, Any]:
try:
datasets_in = {}
with h5py.File(self.fname,'r') as f:
for arg,label in datasets.items():
loc = f[group+'/'+label]
datasets_in[arg]={'data' :loc[()],
'label':label,
'meta': {k:(v.decode() if not h5py3 and type(v) is bytes else v) \
for k,v in loc.attrs.items()}}
return callback(**datasets_in,**args)
except Exception as err:
print(f'Error during calculation: {err}.')
return None
groups = [] groups = []
with h5py.File(self.fname,'r') as f: with h5py.File(self.fname,'r') as f:
@ -1506,12 +1501,10 @@ class Result:
print('No matching dataset found, no data was added.') print('No matching dataset found, no data was added.')
return return
default_arg = partial(self._job_pointwise,callback=func,datasets=datasets,args=args,lock=lock)
for group,result in util.show_progress(pool.imap_unordered(default_arg,groups),len(groups)):# type: ignore for group in util.show_progress(groups):
if not result: if not (result := job_pointwise(group, callback=func, datasets=datasets, args=args)): # type: ignore
continue continue
lock.acquire()
with h5py.File(self.fname, 'a') as f: with h5py.File(self.fname, 'a') as f:
try: try:
if not self._protected and '/'.join([group,result['label']]) in f: if not self._protected and '/'.join([group,result['label']]) in f:
@ -1543,10 +1536,6 @@ class Result:
except (OSError,RuntimeError) as err: except (OSError,RuntimeError) as err:
print(f'Could not add dataset: {err}.') print(f'Could not add dataset: {err}.')
lock.release()
pool.close()
pool.join()
def _mappings(self): def _mappings(self):
@ -2064,7 +2053,7 @@ class Result:
cfg_dir = (Path.cwd() if target_dir is None else Path(target_dir)) cfg_dir = (Path.cwd() if target_dir is None else Path(target_dir))
with h5py.File(self.fname,'r') as f_in: with h5py.File(self.fname,'r') as f_in:
f_in['setup'].visititems(partial(export, f_in['setup'].visititems(functools.partial(export,
output=output, output=output,
cfg_dir=cfg_dir, cfg_dir=cfg_dir,
overwrite=overwrite)) overwrite=overwrite))

89
python/damask/_rotation.py
View File

@ -375,6 +375,11 @@ class Rotation:
Return self@other. Return self@other.
Rotate vector, second-order tensor, or fourth-order tensor. Rotate vector, second-order tensor, or fourth-order tensor.
`other` is interpreted as an array of tensor quantities with the highest-possible order
considering the shape of `self`. Compatible innermost dimensions will blend.
For instance, shapes of (2,) and (3,3) for `self` and `other` prompt interpretation of
`other` as a second-rank tensor and result in (2,) rotated tensors, whereas
shapes of (2,1) and (3,3) for `self` and `other` result in (2,3) rotated vectors.
Parameters Parameters
---------- ----------
@ -386,29 +391,73 @@ class Rotation:
rotated : numpy.ndarray, shape (...,3), (...,3,3), or (...,3,3,3,3) rotated : numpy.ndarray, shape (...,3), (...,3,3), or (...,3,3,3,3)
Rotated vector or tensor, i.e. transformed to frame defined by rotation. Rotated vector or tensor, i.e. transformed to frame defined by rotation.
Examples
--------
All below examples rely on imported modules:
>>> import numpy as np
>>> import damask
Application of twelve (random) rotations to a set of five vectors.
>>> r = damask.Rotation.from_random(shape=(12))
>>> o = np.ones((5,3))
>>> (r@o).shape # (12) @ (5, 3)
(12,5, 3)
Application of a (random) rotation to all twelve second-rank tensors.
>>> r = damask.Rotation.from_random()
>>> o = np.ones((12,3,3))
>>> (r@o).shape # (1) @ (12, 3,3)
(12,3,3)
Application of twelve (random) rotations to the corresponding twelve second-rank tensors.
>>> r = damask.Rotation.from_random(shape=(12))
>>> o = np.ones((12,3,3))
>>> (r@o).shape # (12) @ (3,3)
(12,3,3)
Application of each of three (random) rotations to all three vectors.
>>> r = damask.Rotation.from_random(shape=(3))
>>> o = np.ones((3,3))
>>> (r[...,np.newaxis]@o[np.newaxis,...]).shape # (3,1) @ (1,3, 3)
(3,3,3)
Application of twelve (random) rotations to all twelve second-rank tensors.
>>> r = damask.Rotation.from_random(shape=(12))
>>> o = np.ones((12,3,3))
>>> (r@o[np.newaxis,...]).shape # (12) @ (1,12, 3,3)
(12,3,3,3)
""" """
if isinstance(other, np.ndarray): if isinstance(other, np.ndarray):
if self.shape + (3,) == other.shape: obs = util.shapeblender(self.shape,other.shape,keep_ones=False)[len(self.shape):]
q_m = self.quaternion[...,0] for l in [4,2,1]:
p_m = self.quaternion[...,1:] if obs[-l:] == l*(3,):
A = q_m**2 - np.einsum('...i,...i',p_m,p_m) bs = util.shapeblender(self.shape,other.shape[:-l],False)
B = 2. * np.einsum('...i,...i',p_m,other) self_ = self.broadcast_to(bs) if self.shape != bs else self
C = 2. * _P * q_m if l==1:
return np.block([(A * other[...,i]).reshape(self.shape+(1,)) + q_m = self_.quaternion[...,0]
(B * p_m[...,i]).reshape(self.shape+(1,)) + p_m = self_.quaternion[...,1:]
(C * ( p_m[...,(i+1)%3]*other[...,(i+2)%3]\ A = q_m**2 - np.einsum('...i,...i',p_m,p_m)
- p_m[...,(i+2)%3]*other[...,(i+1)%3])).reshape(self.shape+(1,)) B = 2. * np.einsum('...i,...i',p_m,other)
for i in [0,1,2]]) C = 2. * _P * q_m
if self.shape + (3,3) == other.shape: return np.block([(A * other[...,i]) +
R = self.as_matrix() (B * p_m[...,i]) +
return np.einsum('...im,...jn,...mn',R,R,other) (C * ( p_m[...,(i+1)%3]*other[...,(i+2)%3]
if self.shape + (3,3,3,3) == other.shape: - p_m[...,(i+2)%3]*other[...,(i+1)%3]))
R = self.as_matrix() for i in [0,1,2]]).reshape(bs+(3,),order='F')
return np.einsum('...im,...jn,...ko,...lp,...mnop',R,R,R,R,other) else:
else: return np.einsum({2: '...im,...jn,...mn',
raise ValueError('can only rotate vectors, second-order tensors, and fourth-order tensors') 4: '...im,...jn,...ko,...lp,...mnop'}[l],
*l*[self_.as_matrix()],
other)
raise ValueError('can only rotate vectors, second-order tensors, and fourth-order tensors')
elif isinstance(other, Rotation): elif isinstance(other, Rotation):
raise TypeError('use "R1*R2", i.e. multiplication, to compose rotations "R1" and "R2"') raise TypeError('use "R2*R1", i.e. multiplication, to compose rotations "R1" and "R2"')
else: else:
raise TypeError(f'cannot rotate "{type(other)}"') raise TypeError(f'cannot rotate "{type(other)}"')

63
python/damask/util.py
View File

@ -512,7 +512,8 @@ def shapeshifter(fro: _Tuple[int, ...],
return tuple(final_shape[::-1] if mode == 'left' else final_shape) return tuple(final_shape[::-1] if mode == 'left' else final_shape)
def shapeblender(a: _Tuple[int, ...], def shapeblender(a: _Tuple[int, ...],
b: _Tuple[int, ...]) -> _Tuple[int, ...]: b: _Tuple[int, ...],
keep_ones: bool = True) -> _Tuple[int, ...]:
""" """
Return a shape that overlaps the rightmost entries of 'a' with the leftmost of 'b'. Return a shape that overlaps the rightmost entries of 'a' with the leftmost of 'b'.
@ -522,6 +523,9 @@ def shapeblender(a: _Tuple[int, ...],
Shape of first array. Shape of first array.
b : tuple b : tuple
Shape of second array. Shape of second array.
keep_ones : bool, optional
Treat innermost '1's as literal value instead of dimensional placeholder.
Defaults to True.
Examples Examples
-------- --------
@ -531,13 +535,30 @@ def shapeblender(a: _Tuple[int, ...],
(1,2,3) (1,2,3)
>>> shapeblender((1,),(2,2,1)) >>> shapeblender((1,),(2,2,1))
(1,2,2,1) (1,2,2,1)
>>> shapeblender((1,),(2,2,1),False)
(2,2,1)
>>> shapeblender((3,2),(3,2)) >>> shapeblender((3,2),(3,2))
(3,2) (3,2)
""" """
i = min(len(a),len(b)) def is_broadcastable(a,b):
while i > 0 and a[-i:] != b[:i]: i -= 1 try:
return a + b[i:] _np.broadcast_shapes(a,b)
return True
except ValueError:
return False
a_,_b = a,b
if keep_ones:
i = min(len(a_),len(_b))
while i > 0 and a_[-i:] != _b[:i]: i -= 1
return a_ + _b[i:]
else:
a_ += max(0,len(_b)-len(a_))*(1,)
while not is_broadcastable(a_,_b):
a_ = a_ + ((1,) if len(a_)<=len(_b) else ())
_b = ((1,) if len(_b)<len(a_) else ()) + _b
return _np.broadcast_shapes(a_,_b)
def _docstringer(docstring: _Union[str, _Callable], def _docstringer(docstring: _Union[str, _Callable],
@ -698,7 +719,7 @@ def pass_on(keyword: str,
return wrapper return wrapper
return decorator return decorator
def DREAM3D_base_group(fname: _Union[str, _Path]) -> str: def DREAM3D_base_group(fname: _Union[str, _Path, _h5py.File]) -> str:
""" """
Determine the base group of a DREAM.3D file. Determine the base group of a DREAM.3D file.
@ -707,7 +728,7 @@ def DREAM3D_base_group(fname: _Union[str, _Path]) -> str:
Parameters Parameters
---------- ----------
fname : str or pathlib.Path fname : str, pathlib.Path, or _h5py.File
Filename of the DREAM.3D (HDF5) file. Filename of the DREAM.3D (HDF5) file.
Returns Returns
@ -716,15 +737,19 @@ def DREAM3D_base_group(fname: _Union[str, _Path]) -> str:
Path to the base group. Path to the base group.
""" """
with _h5py.File(_Path(fname).expanduser(),'r') as f: def get_base_group(f: _h5py.File) -> str:
base_group = f.visit(lambda path: path.rsplit('/',2)[0] if '_SIMPL_GEOMETRY/SPACING' in path else None) base_group = f.visit(lambda path: path.rsplit('/',2)[0] if '_SIMPL_GEOMETRY/SPACING' in path else None)
if base_group is None:
raise ValueError(f'could not determine base group in file "{fname}"')
return base_group
if base_group is None: if isinstance(fname,_h5py.File):
raise ValueError(f'could not determine base group in file "{fname}"') return get_base_group(fname)
return base_group with _h5py.File(_Path(fname).expanduser(),'r') as f:
return get_base_group(f)
def DREAM3D_cell_data_group(fname: _Union[str, _Path]) -> str: def DREAM3D_cell_data_group(fname: _Union[str, _Path, _h5py.File]) -> str:
""" """
Determine the cell data group of a DREAM.3D file. Determine the cell data group of a DREAM.3D file.
@ -734,7 +759,7 @@ def DREAM3D_cell_data_group(fname: _Union[str, _Path]) -> str:
Parameters Parameters
---------- ----------
fname : str or pathlib.Path fname : str, pathlib.Path, or h5py.File
Filename of the DREAM.3D (HDF5) file. Filename of the DREAM.3D (HDF5) file.
Returns Returns
@ -743,17 +768,21 @@ def DREAM3D_cell_data_group(fname: _Union[str, _Path]) -> str:
Path to the cell data group. Path to the cell data group.
""" """
base_group = DREAM3D_base_group(fname) def get_cell_data_group(f: _h5py.File) -> str:
with _h5py.File(_Path(fname).expanduser(),'r') as f: base_group = DREAM3D_base_group(f)
cells = tuple(f['/'.join([base_group,'_SIMPL_GEOMETRY','DIMENSIONS'])][()][::-1]) cells = tuple(f['/'.join([base_group,'_SIMPL_GEOMETRY','DIMENSIONS'])][()][::-1])
cell_data_group = f[base_group].visititems(lambda path,obj: path.split('/')[0] \ cell_data_group = f[base_group].visititems(lambda path,obj: path.split('/')[0] \
if isinstance(obj,_h5py._hl.dataset.Dataset) and _np.shape(obj)[:-1] == cells \ if isinstance(obj,_h5py._hl.dataset.Dataset) and _np.shape(obj)[:-1] == cells \
else None) else None)
if cell_data_group is None:
raise ValueError(f'could not determine cell-data group in file "{fname}/{base_group}"')
return cell_data_group
if cell_data_group is None: if isinstance(fname,_h5py.File):
raise ValueError(f'could not determine cell-data group in file "{fname}/{base_group}"') return get_cell_data_group(fname)
return cell_data_group with _h5py.File(_Path(fname).expanduser(),'r') as f:
return get_cell_data_group(f)
def Bravais_to_Miller(*, def Bravais_to_Miller(*,

25
python/tests/test_Orientation.py
View File

@ -162,7 +162,7 @@ class TestOrientation:
([np.arccos(3**(-.5)),np.pi/4,0],[0,0],[0,0,1],[0,0,1])]) ([np.arccos(3**(-.5)),np.pi/4,0],[0,0],[0,0,1],[0,0,1])])
def test_fiber_IPF(self,crystal,sample,direction,color): def test_fiber_IPF(self,crystal,sample,direction,color):
fiber = Orientation.from_fiber_component(crystal=crystal,sample=sample,family='cubic',shape=200) fiber = Orientation.from_fiber_component(crystal=crystal,sample=sample,family='cubic',shape=200)
print(np.allclose(fiber.IPF_color(direction),color)) assert np.allclose(fiber.IPF_color(direction),color)
@pytest.mark.parametrize('kwargs',[ @pytest.mark.parametrize('kwargs',[
@ -455,11 +455,9 @@ class TestOrientation:
p = Orientation.from_random(family=family,shape=right) p = Orientation.from_random(family=family,shape=right)
blend = util.shapeblender(o.shape,p.shape) blend = util.shapeblender(o.shape,p.shape)
for loc in np.random.randint(0,blend,(10,len(blend))): for loc in np.random.randint(0,blend,(10,len(blend))):
# print(f'{a}/{b} @ {loc}') l = () if left is None else tuple(np.minimum(np.array(left )-1,loc[:len(left)]))
# print(o[tuple(loc[:len(o.shape)])].disorientation(p[tuple(loc[-len(p.shape):])])) r = () if right is None else tuple(np.minimum(np.array(right)-1,loc[-len(right):]))
# print(o.disorientation(p)[tuple(loc)]) assert o[l].disorientation(p[r]).isclose(o.disorientation(p)[tuple(loc)])
assert o[tuple(loc[:len(o.shape)])].disorientation(p[tuple(loc[-len(p.shape):])]) \
.isclose(o.disorientation(p)[tuple(loc)])
@pytest.mark.parametrize('family',crystal_families) @pytest.mark.parametrize('family',crystal_families)
@pytest.mark.parametrize('left,right',[ @pytest.mark.parametrize('left,right',[
@ -467,13 +465,16 @@ class TestOrientation:
((2,2),(4,4)), ((2,2),(4,4)),
((3,1),(1,3)), ((3,1),(1,3)),
(None,(3,)), (None,(3,)),
(None,()),
]) ])
def test_IPF_color_blending(self,family,left,right): def test_IPF_color_blending(self,family,left,right):
o = Orientation.from_random(family=family,shape=left) o = Orientation.from_random(family=family,shape=left)
v = np.random.random(right+(3,)) v = np.random.random(right+(3,))
blend = util.shapeblender(o.shape,v.shape[:-1]) blend = util.shapeblender(o.shape,v.shape[:-1])
for loc in np.random.randint(0,blend,(10,len(blend))): for loc in np.random.randint(0,blend,(10,len(blend))):
assert np.allclose(o[tuple(loc[:len(o.shape)])].IPF_color(v[tuple(loc[-len(v.shape[:-1]):])]), l = () if left is None else tuple(np.minimum(np.array(left )-1,loc[:len(left)]))
r = () if right is None else tuple(np.minimum(np.array(right)-1,loc[-len(right):]))
assert np.allclose(o[l].IPF_color(v[r]),
o.IPF_color(v)[tuple(loc)]) o.IPF_color(v)[tuple(loc)])
@pytest.mark.parametrize('family',crystal_families) @pytest.mark.parametrize('family',crystal_families)
@ -488,7 +489,9 @@ class TestOrientation:
v = np.random.random(right+(3,)) v = np.random.random(right+(3,))
blend = util.shapeblender(o.shape,v.shape[:-1]) blend = util.shapeblender(o.shape,v.shape[:-1])
for loc in np.random.randint(0,blend,(10,len(blend))): for loc in np.random.randint(0,blend,(10,len(blend))):
assert np.allclose(o[tuple(loc[:len(o.shape)])].to_SST(v[tuple(loc[-len(v.shape[:-1]):])]), l = () if left is None else tuple(np.minimum(np.array(left )-1,loc[:len(left)]))
r = () if right is None else tuple(np.minimum(np.array(right)-1,loc[-len(right):]))
assert np.allclose(o[l].to_SST(v[r]),
o.to_SST(v)[tuple(loc)]) o.to_SST(v)[tuple(loc)])
@pytest.mark.parametrize('lattice,a,b,c,alpha,beta,gamma', @pytest.mark.parametrize('lattice,a,b,c,alpha,beta,gamma',
@ -514,8 +517,10 @@ class TestOrientation:
v = np.random.random(right+(3,)) v = np.random.random(right+(3,))
blend = util.shapeblender(o.shape,v.shape[:-1]) blend = util.shapeblender(o.shape,v.shape[:-1])
for loc in np.random.randint(0,blend,(10,len(blend))): for loc in np.random.randint(0,blend,(10,len(blend))):
assert np.allclose(o[tuple(loc[:len(o.shape)])].to_pole(uvw=v[tuple(loc[-len(v.shape[:-1]):])]), l = () if left is None else tuple(np.minimum(np.array(left )-1,loc[:len(left)]))
o.to_pole(uvw=v)[tuple(loc)]) r = () if right is None else tuple(np.minimum(np.array(right)-1,loc[-len(right):]))
assert np.allclose(o[l].to_pole(uvw=v[r]),
o.to_pole(uvw=v)[tuple(loc)])
def test_mul_invalid(self): def test_mul_invalid(self):
with pytest.raises(TypeError): with pytest.raises(TypeError):

2
python/tests/test_Result.py
View File

@ -326,7 +326,7 @@ class TestResult:
if shape == 'pseudo_scalar': default.add_calculation('#F#[:,0,0:1]','x','1','a pseudo scalar') if shape == 'pseudo_scalar': default.add_calculation('#F#[:,0,0:1]','x','1','a pseudo scalar')
if shape == 'scalar': default.add_calculation('#F#[:,0,0]','x','1','just a scalar') if shape == 'scalar': default.add_calculation('#F#[:,0,0]','x','1','just a scalar')
if shape == 'vector': default.add_calculation('#F#[:,:,1]','x','1','just a vector') if shape == 'vector': default.add_calculation('#F#[:,:,1]','x','1','just a vector')
x = default.place('x').reshape((np.product(default.cells),-1)) x = default.place('x').reshape((np.prod(default.cells),-1))
default.add_gradient('x') default.add_gradient('x')
in_file = default.place('gradient(x)') in_file = default.place('gradient(x)')
in_memory = grid_filters.gradient(default.size,x.reshape(tuple(default.cells)+x.shape[1:])).reshape(in_file.shape) in_memory = grid_filters.gradient(default.size,x.reshape(tuple(default.cells)+x.shape[1:])).reshape(in_file.shape)

2
python/tests/test_Rotation.py
View File

@ -1065,7 +1065,7 @@ class TestRotation:
@pytest.mark.parametrize('data',[np.random.rand(4), @pytest.mark.parametrize('data',[np.random.rand(4),
np.random.rand(3,2), np.random.rand(3,2),
np.random.rand(3,2,3,3)]) np.random.rand(3,3,3,1)])
def test_rotate_invalid_shape(self,data): def test_rotate_invalid_shape(self,data):
R = Rotation.from_random() R = Rotation.from_random()
with pytest.raises(ValueError): with pytest.raises(ValueError):

2
python/tests/test_grid_filters.py
View File

@ -398,7 +398,7 @@ class TestGridFilters:
np.arange(cells[1]), np.arange(cells[1]),
np.arange(cells[2]),indexing='ij')).reshape(tuple(cells)+(3,),order='F') np.arange(cells[2]),indexing='ij')).reshape(tuple(cells)+(3,),order='F')
x,y,z = map(np.random.randint,cells) x,y,z = map(np.random.randint,cells)
assert grid_filters.ravel_index(indices)[x,y,z] == np.arange(0,np.product(cells)).reshape(cells,order='F')[x,y,z] assert grid_filters.ravel_index(indices)[x,y,z] == np.arange(0,np.prod(cells)).reshape(cells,order='F')[x,y,z]
def test_unravel_index(self): def test_unravel_index(self):
cells = np.random.randint(8,32,(3)) cells = np.random.randint(8,32,(3))

42
python/tests/test_util.py
View File

@ -128,39 +128,47 @@ class TestUtil:
with pytest.raises(ValueError): with pytest.raises(ValueError):
util.shapeshifter(fro,to,mode) util.shapeshifter(fro,to,mode)
@pytest.mark.parametrize('a,b,answer', @pytest.mark.parametrize('a,b,ones,answer',
[ [
((),(1,),(1,)), ((),(1,),True,(1,)),
((1,),(),(1,)), ((1,),(),False,(1,)),
((1,),(7,),(1,7)), ((1,1),(7,),False,(1,7)),
((2,),(2,2),(2,2)), ((1,),(7,),False,(7,)),
((1,2),(2,2),(1,2,2)), ((1,),(7,),True,(1,7)),
((1,2,3),(2,3,4),(1,2,3,4)), ((2,),(2,2),False,(2,2)),
((1,2,3),(1,2,3),(1,2,3)), ((1,2),(2,2),False,(2,2)),
((1,1,2),(2,2),False,(1,2,2)),
((1,1,2),(2,2),True,(1,1,2,2)),
((1,2,3),(2,3,4),False,(1,2,3,4)),
((1,2,3),(1,2,3),False,(1,2,3)),
]) ])
def test_shapeblender(self,a,b,answer): def test_shapeblender(self,a,b,ones,answer):
assert util.shapeblender(a,b) == answer assert util.shapeblender(a,b,ones) == answer
@pytest.mark.parametrize('style',[util.emph,util.deemph,util.warn,util.strikeout]) @pytest.mark.parametrize('style',[util.emph,util.deemph,util.warn,util.strikeout])
def test_decorate(self,style): def test_decorate(self,style):
assert 'DAMASK' in style('DAMASK') assert 'DAMASK' in style('DAMASK')
@pytest.mark.parametrize('complete',[True,False]) @pytest.mark.parametrize('complete',[True,False])
def test_D3D_base_group(self,tmp_path,complete): @pytest.mark.parametrize('fhandle',[True,False])
def test_D3D_base_group(self,tmp_path,complete,fhandle):
base_group = ''.join(random.choices('DAMASK', k=10)) base_group = ''.join(random.choices('DAMASK', k=10))
with h5py.File(tmp_path/'base_group.dream3d','w') as f: with h5py.File(tmp_path/'base_group.dream3d','w') as f:
f.create_group('/'.join((base_group,'_SIMPL_GEOMETRY'))) f.create_group('/'.join((base_group,'_SIMPL_GEOMETRY')))
if complete: if complete:
f['/'.join((base_group,'_SIMPL_GEOMETRY'))].create_dataset('SPACING',data=np.ones(3)) f['/'.join((base_group,'_SIMPL_GEOMETRY'))].create_dataset('SPACING',data=np.ones(3))
fname = tmp_path/'base_group.dream3d'
if fhandle: fname = h5py.File(fname)
if complete: if complete:
assert base_group == util.DREAM3D_base_group(tmp_path/'base_group.dream3d') assert base_group == util.DREAM3D_base_group(fname)
else: else:
with pytest.raises(ValueError): with pytest.raises(ValueError):
util.DREAM3D_base_group(tmp_path/'base_group.dream3d') util.DREAM3D_base_group(fname)
@pytest.mark.parametrize('complete',[True,False]) @pytest.mark.parametrize('complete',[True,False])
def test_D3D_cell_data_group(self,tmp_path,complete): @pytest.mark.parametrize('fhandle',[True,False])
def test_D3D_cell_data_group(self,tmp_path,complete,fhandle):
base_group = ''.join(random.choices('DAMASK', k=10)) base_group = ''.join(random.choices('DAMASK', k=10))
cell_data_group = ''.join(random.choices('KULeuven', k=10)) cell_data_group = ''.join(random.choices('KULeuven', k=10))
cells = np.random.randint(1,50,3) cells = np.random.randint(1,50,3)
@ -172,11 +180,13 @@ class TestUtil:
if complete: if complete:
f['/'.join((base_group,cell_data_group))].create_dataset('data',shape=np.append(cells,1)) f['/'.join((base_group,cell_data_group))].create_dataset('data',shape=np.append(cells,1))
fname = tmp_path/'cell_data_group.dream3d'
if fhandle: fname = h5py.File(fname)
if complete: if complete:
assert cell_data_group == util.DREAM3D_cell_data_group(tmp_path/'cell_data_group.dream3d') assert cell_data_group == util.DREAM3D_cell_data_group(fname)
else: else:
with pytest.raises(ValueError): with pytest.raises(ValueError):
util.DREAM3D_cell_data_group(tmp_path/'cell_data_group.dream3d') util.DREAM3D_cell_data_group(fname)
@pytest.mark.parametrize('full,reduced',[({}, {}), @pytest.mark.parametrize('full,reduced',[({}, {}),