Merge branch 'development' into 'empty-table-init'

# Conflicts:
#   python/damask/_table.py
This commit is contained in:
Philip Eisenlohr 2022-05-25 13:29:23 +00:00
commit 4746ac890b
14 changed files with 295 additions and 207 deletions

@ -1 +1 @@
Subproject commit 1766c855e54668d6225c9b22db536be7052605bc
Subproject commit 0e82975b23a1bd4310c523388f1cadf1b8e03dd0

View File

@ -1 +1 @@
3.0.0-alpha6-371-g18d862cdb
3.0.0-alpha6-395-g9696c1967

View File

@ -482,7 +482,7 @@ class ConfigMaterial(Config):
"""
N,n,shaped = 1,1,{}
map_dim = {'O':-1,'F_i':-2}
map_dim = {'O':-1,'V_e':-2}
for k,v in kwargs.items():
shaped[k] = np.array(v)
s = shaped[k].shape[:map_dim.get(k,None)]
@ -494,12 +494,12 @@ class ConfigMaterial(Config):
if 'v' not in kwargs:
shaped['v'] = np.broadcast_to(1/n,(N,n))
map_shape = {'O':(N,n,4),'F_i':(N,n,3,3)}
map_shape = {'O':(N,n,4),'V_e':(N,n,3,3)}
for k,v in shaped.items():
target = map_shape.get(k,(N,n))
obj = np.broadcast_to(v.reshape(util.shapeshifter(v.shape, target, mode = 'right')), target)
for i in range(N):
if k in ['phase','O','v','F_i']:
if k in ['phase','O','v','V_e']:
for j in range(n):
mat[i]['constituents'][j][k] = obj[i,j].item() if isinstance(obj[i,j],np.generic) else obj[i,j]
else:

View File

@ -50,15 +50,16 @@ class Grid:
Coordinates of grid origin in meter. Defaults to [0.0,0.0,0.0].
initial_conditions : dictionary, optional
Labels and values of the inital conditions at each material point.
comments : (list of) str, optional
Comments, e.g. history of operations.
comments : str or iterable of str, optional
Additional, human-readable information, e.g. history of operations.
"""
self.material = material
self.size = size # type: ignore
self.origin = origin # type: ignore
self.initial_conditions = {} if initial_conditions is None else initial_conditions
self.comments = [] if comments is None else comments # type: ignore
comments_ = [comments] if isinstance(comments,str) else comments
self.comments = [] if comments_ is None else [str(c) for c in comments_]
def __repr__(self) -> str:
"""Give short human-readable summary."""

View File

@ -1,4 +1,5 @@
import multiprocessing as mp
from multiprocessing.synchronize import Lock
import re
import fnmatch
import os
@ -10,6 +11,7 @@ from pathlib import Path
from functools import partial
from collections import defaultdict
from collections.abc import Iterable
from typing import Union, Callable, Any, Sequence, Literal, Dict, List, Tuple
import h5py
import numpy as np
@ -22,19 +24,21 @@ from . import grid_filters
from . import mechanics
from . import tensor
from . import util
from ._typehints import FloatSequence, IntSequence
h5py3 = h5py.__version__[0] == '3'
chunk_size = 1024**2//8 # for compression in HDF5
prefix_inc = 'increment_'
def _read(dataset):
def _read(dataset: h5py._hl.dataset.Dataset) -> np.ndarray:
"""Read a dataset and its metadata into a numpy.ndarray."""
metadata = {k:(v.decode() if not h5py3 and type(v) is bytes else v) for k,v in dataset.attrs.items()}
dtype = np.dtype(dataset.dtype,metadata=metadata)
dtype = np.dtype(dataset.dtype,metadata=metadata) # type: ignore
return np.array(dataset,dtype=dtype)
def _match(requested,existing):
def _match(requested,
existing: h5py._hl.base.KeysViewHDF5) -> List[Any]:
"""Find matches among two sets of labels."""
def flatten_list(list_of_lists):
return [e for e_ in list_of_lists for e in e_]
@ -50,7 +54,10 @@ def _match(requested,existing):
return sorted(set(flatten_list([fnmatch.filter(existing,r) for r in requested_])),
key=util.natural_sort)
def _empty_like(dataset,N_materialpoints,fill_float,fill_int):
def _empty_like(dataset: np.ma.core.MaskedArray,
N_materialpoints: int,
fill_float: float,
fill_int: int) -> np.ma.core.MaskedArray:
"""Create empty numpy.ma.MaskedArray."""
return ma.array(np.empty((N_materialpoints,)+dataset.shape[1:],dataset.dtype),
fill_value = fill_float if dataset.dtype in np.sctypes['float'] else fill_int,
@ -84,7 +91,7 @@ class Result:
"""
def __init__(self,fname):
def __init__(self, fname: Union[str, Path]):
"""
New result view bound to a HDF5 file.
@ -102,7 +109,7 @@ class Result:
if self.version_major != 0 or not 12 <= self.version_minor <= 14:
raise TypeError(f'unsupported DADF5 version "{self.version_major}.{self.version_minor}"')
if self.version_major == 0 and self.version_minor < 14:
self.export_setup = None
self.export_setup = None # type: ignore
self.structured = 'cells' in f['geometry'].attrs.keys()
@ -111,7 +118,7 @@ class Result:
self.size = f['geometry'].attrs['size']
self.origin = f['geometry'].attrs['origin']
else:
self.add_curl = self.add_divergence = self.add_gradient = None
self.add_curl = self.add_divergence = self.add_gradient = None # type: ignore
r = re.compile(rf'{prefix_inc}([0-9]+)')
self.increments = sorted([i for i in f.keys() if r.match(i)],key=util.natural_sort)
@ -126,7 +133,7 @@ class Result:
self.phase = f['cell_to/phase']['label'].astype('str')
self.phases = sorted(np.unique(self.phase),key=util.natural_sort)
self.fields = []
self.fields: List[str] = []
for c in self.phases:
self.fields += f['/'.join([self.increments[0],'phase',c])].keys()
for m in self.homogenizations:
@ -144,14 +151,14 @@ class Result:
self._protected = True
def __copy__(self):
def __copy__(self) -> "Result":
"""Create deep copy."""
return copy.deepcopy(self)
copy = __copy__
def __repr__(self):
def __repr__(self) -> str:
"""Give short human-readable summary."""
with h5py.File(self.fname,'r') as f:
header = [f'Created by {f.attrs["creator"]}',
@ -171,12 +178,12 @@ class Result:
def _manage_view(self,
action,
increments=None,
times=None,
phases=None,
homogenizations=None,
fields=None):
action: Literal['set', 'add', 'del'],
increments: Union[int, Sequence[int], str, Sequence[str], bool] = None,
times: Union[float, Sequence[float], str, Sequence[str], bool] = None,
phases: Union[str, Sequence[str], bool] = None,
homogenizations: Union[str, Sequence[str], bool] = None,
fields: Union[str, Sequence[str], bool] = None) -> "Result":
"""
Manages the visibility of the groups.
@ -204,8 +211,7 @@ class Result:
datasets = '*'
elif datasets is False:
datasets = []
choice = list(datasets).copy() if hasattr(datasets,'__iter__') and not isinstance(datasets,str) else \
[datasets]
choice = [datasets] if not hasattr(datasets,'__iter__') or isinstance(datasets,str) else list(datasets) # type: ignore
if what == 'increments':
choice = [c if isinstance(c,str) and c.startswith(prefix_inc) else
@ -216,7 +222,7 @@ class Result:
if choice == ['*']:
choice = self.increments
else:
iterator = map(float,choice)
iterator = map(float,choice) # type: ignore
choice = []
for c in iterator:
idx = np.searchsorted(self.times,c)
@ -224,7 +230,7 @@ class Result:
if np.isclose(c,self.times[idx]):
choice.append(self.increments[idx])
elif np.isclose(c,self.times[idx+1]):
choice.append(self.increments[idx+1])
choice.append(self.increments[idx+1]) # type: ignore
valid = _match(choice,getattr(self,what))
existing = set(self.visible[what])
@ -241,7 +247,9 @@ class Result:
return dup
def increments_in_range(self,start=None,end=None):
def increments_in_range(self,
start: Union[str, int] = None,
end: Union[str, int] = None) -> Sequence[int]:
"""
Get all increments within a given range.
@ -263,7 +271,9 @@ class Result:
self.incs[-1] if end is None else end))
return [i for i in self.incs if s <= i <= e]
def times_in_range(self,start=None,end=None):
def times_in_range(self,
start: float = None,
end: float = None) -> Sequence[int]:
"""
Get all increments within a given time range.
@ -286,12 +296,12 @@ class Result:
def view(self,*,
increments=None,
times=None,
phases=None,
homogenizations=None,
fields=None,
protected=None):
increments: Union[int, Sequence[int], str, Sequence[str], bool] = None,
times: Union[float, Sequence[float], str, Sequence[str], bool] = None,
phases: Union[str, Sequence[str], bool] = None,
homogenizations: Union[str, Sequence[str], bool] = None,
fields: Union[str, Sequence[str], bool] = None,
protected: bool = None) -> "Result":
"""
Set view.
@ -343,11 +353,11 @@ class Result:
def view_more(self,*,
increments=None,
times=None,
phases=None,
homogenizations=None,
fields=None):
increments: Union[int, Sequence[int], str, Sequence[str], bool] = None,
times: Union[float, Sequence[float], str, Sequence[str], bool] = None,
phases: Union[str, Sequence[str], bool] = None,
homogenizations: Union[str, Sequence[str], bool] = None,
fields: Union[str, Sequence[str], bool] = None) -> "Result":
"""
Add to view.
@ -386,11 +396,11 @@ class Result:
def view_less(self,*,
increments=None,
times=None,
phases=None,
homogenizations=None,
fields=None):
increments: Union[int, Sequence[int], str, Sequence[str], bool] = None,
times: Union[float, Sequence[float], str, Sequence[str], bool] = None,
phases: Union[str, Sequence[str], bool] = None,
homogenizations: Union[str, Sequence[str], bool] = None,
fields: Union[str, Sequence[str], bool] = None) -> "Result":
"""
Remove from view.
@ -427,7 +437,9 @@ class Result:
return self._manage_view('del',increments,times,phases,homogenizations,fields)
def rename(self,name_src,name_dst):
def rename(self,
name_src: str,
name_dst: str):
"""
Rename/move datasets (within the same group/folder).
@ -468,7 +480,7 @@ class Result:
del f[path_src]
def remove(self,name):
def remove(self, name: str):
"""
Remove/delete datasets.
@ -502,7 +514,7 @@ class Result:
if path in f.keys(): del f[path]
def list_data(self):
def list_data(self) -> List[str]:
"""
Collect information on all active datasets in the file.
@ -533,7 +545,8 @@ class Result:
return msg
def enable_user_function(self,func):
def enable_user_function(self,
func: Callable):
globals()[func.__name__]=func
print(f'Function {func.__name__} enabled in add_calculation.')
@ -544,7 +557,7 @@ class Result:
@property
def coordinates0_point(self):
def coordinates0_point(self) -> np.ndarray:
"""Initial/undeformed cell center coordinates."""
if self.structured:
return grid_filters.coordinates0_point(self.cells,self.size,self.origin).reshape(-1,3,order='F')
@ -553,7 +566,7 @@ class Result:
return f['geometry/x_p'][()]
@property
def coordinates0_node(self):
def coordinates0_node(self) -> np.ndarray:
"""Initial/undeformed nodal coordinates."""
if self.structured:
return grid_filters.coordinates0_node(self.cells,self.size,self.origin).reshape(-1,3,order='F')
@ -562,7 +575,7 @@ class Result:
return f['geometry/x_n'][()]
@property
def geometry0(self):
def geometry0(self) -> VTK:
"""Initial/undeformed geometry."""
if self.structured:
return VTK.from_image_data(self.cells,self.size,self.origin)
@ -575,7 +588,7 @@ class Result:
@staticmethod
def _add_absolute(x):
def _add_absolute(x: Dict[str, Any]) -> Dict[str, Any]:
return {
'data': np.abs(x['data']),
'label': f'|{x["label"]}|',
@ -585,7 +598,7 @@ class Result:
'creator': 'add_absolute'
}
}
def add_absolute(self,x):
def add_absolute(self, x: str):
"""
Add absolute value.
@ -599,7 +612,7 @@ class Result:
@staticmethod
def _add_calculation(**kwargs):
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']")
@ -617,7 +630,11 @@ class Result:
'creator': 'add_calculation'
}
}
def add_calculation(self,formula,name,unit='n/a',description=None):
def add_calculation(self,
formula: str,
name: str,
unit: str = 'n/a',
description: str = None):
"""
Add result of a general formula.
@ -667,7 +684,7 @@ class Result:
@staticmethod
def _add_stress_Cauchy(P,F):
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',
@ -679,7 +696,9 @@ class Result:
'creator': 'add_stress_Cauchy'
}
}
def add_stress_Cauchy(self,P='P',F='F'):
def add_stress_Cauchy(self,
P: str = 'P',
F: str = 'F'):
"""
Add Cauchy stress calculated from first Piola-Kirchhoff stress and deformation gradient.
@ -695,7 +714,7 @@ class Result:
@staticmethod
def _add_determinant(T):
def _add_determinant(T: Dict[str, Any]) -> Dict[str, Any]:
return {
'data': np.linalg.det(T['data']),
'label': f"det({T['label']})",
@ -705,7 +724,7 @@ class Result:
'creator': 'add_determinant'
}
}
def add_determinant(self,T):
def add_determinant(self, T: str):
"""
Add the determinant of a tensor.
@ -727,7 +746,7 @@ class Result:
@staticmethod
def _add_deviator(T):
def _add_deviator(T: Dict[str, Any]) -> Dict[str, Any]:
return {
'data': tensor.deviatoric(T['data']),
'label': f"s_{T['label']}",
@ -737,7 +756,7 @@ class Result:
'creator': 'add_deviator'
}
}
def add_deviator(self,T):
def add_deviator(self, T: str):
"""
Add the deviatoric part of a tensor.
@ -759,13 +778,15 @@ class Result:
@staticmethod
def _add_eigenvalue(T_sym,eigenvalue):
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],
@ -776,7 +797,9 @@ class Result:
'creator': 'add_eigenvalue'
}
}
def add_eigenvalue(self,T_sym,eigenvalue='max'):
def add_eigenvalue(self,
T_sym: str,
eigenvalue: Literal['max', 'mid', 'min'] = 'max'):
"""
Add eigenvalues of symmetric tensor.
@ -800,13 +823,16 @@ class Result:
@staticmethod
def _add_eigenvector(T_sym,eigenvalue):
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']})",
@ -817,7 +843,9 @@ class Result:
'creator': 'add_eigenvector'
}
}
def add_eigenvector(self,T_sym,eigenvalue='max'):
def add_eigenvector(self,
T_sym: str,
eigenvalue: Literal['max', 'mid', 'min'] = 'max'):
"""
Add eigenvector of symmetric tensor.
@ -834,7 +862,7 @@ class Result:
@staticmethod
def _add_IPF_color(l,q):
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)
@ -849,7 +877,9 @@ class Result:
'creator': 'add_IPF_color'
}
}
def add_IPF_color(self,l,q='O'):
def add_IPF_color(self,
l: FloatSequence,
q: str = 'O'):
"""
Add RGB color tuple of inverse pole figure (IPF) color.
@ -874,7 +904,7 @@ class Result:
@staticmethod
def _add_maximum_shear(T_sym):
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']})",
@ -884,7 +914,7 @@ class Result:
'creator': 'add_maximum_shear'
}
}
def add_maximum_shear(self,T_sym):
def add_maximum_shear(self, T_sym: str):
"""
Add maximum shear components of symmetric tensor.
@ -898,7 +928,7 @@ class Result:
@staticmethod
def _add_equivalent_Mises(T_sym,kind):
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':
@ -918,7 +948,9 @@ class Result:
'creator': 'add_Mises'
}
}
def add_equivalent_Mises(self,T_sym,kind=None):
def add_equivalent_Mises(self,
T_sym: str,
kind: str = None):
"""
Add the equivalent Mises stress or strain of a symmetric tensor.
@ -949,10 +981,10 @@ class Result:
@staticmethod
def _add_norm(x,ord):
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 = 1
axis: Union[int, Tuple[int, int]] = 1
t = 'vector'
if o is None: o = 2
elif len(x['data'].shape) == 3:
@ -971,7 +1003,9 @@ class Result:
'creator': 'add_norm'
}
}
def add_norm(self,x,ord=None):
def add_norm(self,
x: str,
ord: Union[int, float, Literal['fro', 'nuc']] = None):
"""
Add the norm of vector or tensor.
@ -987,7 +1021,7 @@ class Result:
@staticmethod
def _add_stress_second_Piola_Kirchhoff(P,F):
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',
@ -999,7 +1033,9 @@ class Result:
'creator': 'add_stress_second_Piola_Kirchhoff'
}
}
def add_stress_second_Piola_Kirchhoff(self,P='P',F='F'):
def add_stress_second_Piola_Kirchhoff(self,
P: str = 'P',
F: str = 'F'):
"""
Add second Piola-Kirchhoff stress calculated from first Piola-Kirchhoff stress and deformation gradient.
@ -1023,7 +1059,11 @@ class Result:
@staticmethod
def _add_pole(q,uvw,hkl,with_symmetry,normalize):
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],
@ -1042,7 +1082,13 @@ class Result:
'creator': 'add_pole'
}
}
def add_pole(self,q='O',*,uvw=None,hkl=None,with_symmetry=False,normalize=True):
def add_pole(self,
q: str = 'O',
*,
uvw: FloatSequence = None,
hkl: FloatSequence = None,
with_symmetry: bool = False,
normalize: bool = True):
"""
Add lab frame vector along lattice direction [uvw] or plane normal (hkl).
@ -1067,7 +1113,7 @@ class Result:
@staticmethod
def _add_rotation(F):
def _add_rotation(F: Dict[str, Any]) -> Dict[str, Any]:
return {
'data': mechanics.rotation(F['data']).as_matrix(),
'label': f"R({F['label']})",
@ -1077,7 +1123,7 @@ class Result:
'creator': 'add_rotation'
}
}
def add_rotation(self,F):
def add_rotation(self, F: str):
"""
Add rotational part of a deformation gradient.
@ -1099,7 +1145,7 @@ class Result:
@staticmethod
def _add_spherical(T):
def _add_spherical(T: Dict[str, Any]) -> Dict[str, Any]:
return {
'data': tensor.spherical(T['data'],False),
'label': f"p_{T['label']}",
@ -1109,7 +1155,7 @@ class Result:
'creator': 'add_spherical'
}
}
def add_spherical(self,T):
def add_spherical(self, T: str):
"""
Add the spherical (hydrostatic) part of a tensor.
@ -1131,7 +1177,7 @@ class Result:
@staticmethod
def _add_strain(F,t,m):
def _add_strain(F: Dict[str, Any], t: Literal['V', 'U'], m: float) -> Dict[str, Any]:
return {
'data': mechanics.strain(F['data'],t,m),
'label': f"epsilon_{t}^{m}({F['label']})",
@ -1141,7 +1187,10 @@ class Result:
'creator': 'add_strain'
}
}
def add_strain(self,F='F',t='V',m=0.0):
def add_strain(self,
F: str = 'F',
t: Literal['V', 'U'] = 'V',
m: float = 0.0):
"""
Add strain tensor of a deformation gradient.
@ -1177,7 +1226,7 @@ class Result:
@staticmethod
def _add_stretch_tensor(F,t):
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']})",
@ -1188,7 +1237,9 @@ class Result:
'creator': 'add_stretch_tensor'
}
}
def add_stretch_tensor(self,F='F',t='V'):
def add_stretch_tensor(self,
F: str = 'F',
t: Literal['V', 'U'] = 'V'):
"""
Add stretch tensor of a deformation gradient.
@ -1205,7 +1256,7 @@ class Result:
@staticmethod
def _add_curl(f,size):
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']})",
@ -1215,7 +1266,7 @@ class Result:
'creator': 'add_curl'
}
}
def add_curl(self,f):
def add_curl(self, f: str):
"""
Add curl of a field.
@ -1234,7 +1285,7 @@ class Result:
@staticmethod
def _add_divergence(f,size):
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']})",
@ -1244,7 +1295,7 @@ class Result:
'creator': 'add_divergence'
}
}
def add_divergence(self,f):
def add_divergence(self, f: str):
"""
Add divergence of a field.
@ -1263,7 +1314,7 @@ class Result:
@staticmethod
def _add_gradient(f,size):
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,))),
@ -1274,7 +1325,7 @@ class Result:
'creator': 'add_gradient'
}
}
def add_gradient(self,f):
def add_gradient(self, f: str):
"""
Add gradient of a field.
@ -1292,7 +1343,11 @@ class Result:
self._add_generic_grid(self._add_gradient,{'f':f},{'size':self.size})
def _add_generic_grid(self,func,datasets,args={},constituents=None):
def _add_generic_grid(self,
func: Callable,
datasets: Dict[str, str],
args: Dict[str, str] = {},
constituents = None):
"""
General function to add data on a regular grid.
@ -1313,11 +1368,13 @@ class Result:
at_cell_ph,in_data_ph,at_cell_ho,in_data_ho = self._mappings()
increments = self.place(list(datasets.values()),False)
if not increments: raise RuntimeError("received invalid dataset")
with h5py.File(self.fname, 'a') as f:
for increment in self.place(datasets.values(),False).items():
for increment in increments.items():
for ty in increment[1].items():
for field in ty[1].items():
d = list(field[1].values())[0]
d: np.ma.MaskedArray = list(field[1].values())[0]
if np.any(d.mask): continue
dataset = {'f':{'data':np.reshape(d.data,tuple(self.cells)+d.data.shape[1:]),
'label':list(datasets.values())[0],
@ -1331,21 +1388,26 @@ class Result:
result1 = result[at_cell_ho[x]]
path = '/'.join(['/',increment[0],ty[0],x,field[0]])
dataset = f[path].create_dataset(r['label'],data=result1)
h5_dataset = f[path].create_dataset(r['label'],data=result1)
now = datetime.datetime.now().astimezone()
dataset.attrs['created'] = now.strftime('%Y-%m-%d %H:%M:%S%z') if h5py3 else \
h5_dataset.attrs['created'] = now.strftime('%Y-%m-%d %H:%M:%S%z') if h5py3 else \
now.strftime('%Y-%m-%d %H:%M:%S%z').encode()
for l,v in r['meta'].items():
dataset.attrs[l.lower()]=v if h5py3 else v.encode()
creator = dataset.attrs['creator'] if h5py3 else \
dataset.attrs['creator'].decode()
dataset.attrs['creator'] = f'damask.Result.{creator} v{damask.version}' if h5py3 else \
h5_dataset.attrs[l.lower()]=v if h5py3 else v.encode()
creator = h5_dataset.attrs['creator'] if h5py3 else \
h5_dataset.attrs['creator'].decode()
h5_dataset.attrs['creator'] = f'damask.Result.{creator} v{damask.version}' if h5py3 else \
f'damask.Result.{creator} v{damask.version}'.encode()
def _job_pointwise(self,group,func,datasets,args,lock):
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 = {}
@ -1358,19 +1420,23 @@ class Result:
'meta': {k:(v.decode() if not h5py3 and type(v) is bytes else v) \
for k,v in loc.attrs.items()}}
lock.release()
r = func(**datasets_in,**args)
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,func,datasets,args={}):
def _add_generic_pointwise(self,
func: Callable,
datasets: Dict[str, Any],
args: Dict[str, Any] = {}):
"""
General function to add pointwise data.
Parameters
----------
func : function
callback : function
Callback function that calculates a new dataset from one or
more datasets per HDF5 group.
datasets : dictionary
@ -1396,9 +1462,9 @@ class Result:
print('No matching dataset found, no data was added.')
return
default_arg = partial(self._job_pointwise,func=func,datasets=datasets,args=args,lock=lock)
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)):
for group,result in util.show_progress(pool.imap_unordered(default_arg,groups),len(groups)):# type: ignore
if not result:
continue
lock.acquire()
@ -1410,15 +1476,14 @@ class Result:
dataset.attrs['overwritten'] = True
else:
shape = result['data'].shape
if result['data'].size >= chunk_size*2:
if compress := (result['data'].size >= chunk_size*2):
chunks = (chunk_size//np.prod(shape[1:]),)+shape[1:]
compression = ('gzip',6)
else:
chunks = shape
compression = (None,None)
dataset = f[group].create_dataset(result['label'],data=result['data'],
maxshape=shape, chunks=chunks,
compression=compression[0], compression_opts=compression[1],
compression = 'gzip' if compress else None,
compression_opts = 6 if compress else None,
shuffle=True,fletcher32=True)
now = datetime.datetime.now().astimezone()
@ -1440,7 +1505,8 @@ class Result:
pool.join()
def export_XDMF(self,output='*'):
def export_XDMF(self,
output: Union[str, List[str]] = '*'):
"""
Write XDMF file to directly visualize data from DADF5 file.
@ -1574,7 +1640,13 @@ class Result:
return at_cell_ph,in_data_ph,at_cell_ho,in_data_ho
def export_VTK(self,output='*',mode='cell',constituents=None,fill_float=np.nan,fill_int=0,parallel=True):
def export_VTK(self,
output: Union[str,list] = '*',
mode: str = 'cell',
constituents: IntSequence = None,
fill_float: float = np.nan,
fill_int: int = 0,
parallel: bool = True):
"""
Export to VTK cell/point data.
@ -1612,12 +1684,12 @@ class Result:
else:
raise ValueError(f'invalid mode "{mode}"')
v.comments = util.execution_stamp('Result','export_VTK')
v.comments = [util.execution_stamp('Result','export_VTK')]
N_digits = int(np.floor(np.log10(max(1,self.incs[-1]))))+1
constituents_ = constituents if isinstance(constituents,Iterable) else \
(range(self.N_constituents) if constituents is None else [constituents])
(range(self.N_constituents) if constituents is None else [constituents]) # type: ignore
suffixes = [''] if self.N_constituents == 1 or isinstance(constituents,int) else \
[f'#{c}' for c in constituents_]
@ -1637,7 +1709,7 @@ class Result:
for ty in ['phase','homogenization']:
for field in self.visible['fields']:
outs = {}
outs: Dict[str, np.ma.core.MaskedArray] = {}
for label in self.visible[ty+'s']:
if field not in f['/'.join([inc,ty,label])].keys(): continue
@ -1665,7 +1737,10 @@ class Result:
v.save(f'{self.fname.stem}_inc{inc[10:].zfill(N_digits)}',parallel=parallel)
def get(self,output='*',flatten=True,prune=True):
def get(self,
output: Union[str, List[str]] = '*',
flatten: bool = True,
prune: bool = True):
"""
Collect data per phase/homogenization reflecting the group/folder structure in the DADF5 file.
@ -1687,7 +1762,7 @@ class Result:
Datasets structured by phase/homogenization and according to selected view.
"""
r = {}
r = {} # type: ignore
with h5py.File(self.fname,'r') as f:
for inc in util.show_progress(self.visible['increments']):
@ -1710,7 +1785,13 @@ class Result:
return None if (type(r) == dict and r == {}) else r
def place(self,output='*',flatten=True,prune=True,constituents=None,fill_float=np.nan,fill_int=0):
def place(self,
output: Union[str, List[str]] = '*',
flatten: bool = True,
prune: bool = True,
constituents: IntSequence = None,
fill_float: float = np.nan,
fill_int: int = 0):
"""
Merge data into spatial order that is compatible with the damask.VTK geometry representation.
@ -1748,10 +1829,10 @@ class Result:
Datasets structured by spatial position and according to selected view.
"""
r = {}
r = {} # type: ignore
constituents_ = constituents if isinstance(constituents,Iterable) else \
(range(self.N_constituents) if constituents is None else [constituents])
constituents_ = list(map(int,constituents)) if isinstance(constituents,Iterable) else \
(range(self.N_constituents) if constituents is None else [constituents]) # type: ignore
suffixes = [''] if self.N_constituents == 1 or isinstance(constituents,int) else \
[f'#{c}' for c in constituents_]
@ -1797,7 +1878,9 @@ class Result:
return None if (type(r) == dict and r == {}) else r
def export_setup(self,output='*',overwrite=False):
def export_setup(self,
output: Union[str, List[str]] = '*',
overwrite: bool = False):
"""
Export configuration files.
@ -1811,7 +1894,7 @@ class Result:
Defaults to False.
"""
def export(name,obj,output,overwrite):
def export(name: str, obj: Union[h5py.Dataset,h5py.Group], output: Union[str,List[str]], overwrite: bool):
if type(obj) == h5py.Dataset and _match(output,[name]):
d = obj.attrs['description'] if h5py3 else obj.attrs['description'].decode()
if not Path(name).exists() or overwrite:

View File

@ -1,6 +1,6 @@
import re
import copy
from typing import Union, Tuple, List
from typing import Union, Tuple, List, Iterable
import pandas as pd
import numpy as np
@ -12,9 +12,9 @@ class Table:
"""Manipulate multi-dimensional spreadsheet-like data."""
def __init__(self,
shapes: dict = {},
data: np.ndarray = None,
comments: Union[str, list] = None):
shapes: dict,
data: np.ndarray,
comments: Union[str, Iterable[str]] = None):
"""
New spreadsheet.
@ -30,7 +30,7 @@ class Table:
"""
comments_ = [comments] if isinstance(comments,str) else comments
self.comments = [] if comments_ is None else [c for c in comments_]
self.comments = [] if comments_ is None else [str(c) for c in comments_]
self.shapes = { k:(v,) if isinstance(v,(np.int64,np.int32,int)) else v for k,v in shapes.items() }
self.data = pd.DataFrame(data=data)
self._relabel('uniform')
@ -434,8 +434,8 @@ class Table:
def rename(self,
old: Union[str, List[str]],
new: Union[str, List[str]],
old: Union[str, Iterable[str]],
new: Union[str, Iterable[str]],
info: str = None) -> 'Table':
"""
Rename column data.

View File

@ -103,7 +103,7 @@ class VTK:
Parameters
----------
comments : str or list of str
comments : str or sequence of str
Comments.
"""
@ -153,11 +153,11 @@ class VTK:
Parameters
----------
cells : iterable of int, len (3)
cells : sequence of int, len (3)
Number of cells along each dimension.
size : iterable of float, len (3)
size : sequence of float, len (3)
Physical length along each dimension.
origin : iterable of float, len (3), optional
origin : sequence of float, len (3), optional
Coordinates of grid origin.
Returns
@ -257,7 +257,7 @@ class VTK:
Parameters
----------
grid : iterables of floats, len (3)
grid : sequence of sequences of floats, len (3)
Grid coordinates along x, y, and z directions.
Returns

View File

@ -113,15 +113,15 @@ class TestConfigMaterial:
@pytest.mark.parametrize('N,n,kw',[
(1,1,{'phase':'Gold',
'O':[1,0,0,0],
'F_i':np.eye(3),
'V_e':np.eye(3),
'homogenization':'SX'}),
(3,1,{'phase':'Gold',
'O':Rotation.from_random(3),
'F_i':np.broadcast_to(np.eye(3),(3,3,3)),
'V_e':np.broadcast_to(np.eye(3),(3,3,3)),
'homogenization':'SX'}),
(2,3,{'phase':np.broadcast_to(['a','b','c'],(2,3)),
'O':Rotation.from_random((2,3)),
'F_i':np.broadcast_to(np.eye(3),(2,3,3,3)),
'V_e':np.broadcast_to(np.eye(3),(2,3,3,3)),
'homogenization':['SX','PX']}),
])
def test_material_add(self,kw,N,n):

View File

@ -48,7 +48,7 @@ module HDF5_utilities
!> @details for parallel IO, all dimension except for the last need to match
!--------------------------------------------------------------------------------------------------
interface HDF5_write
#if defined(__GFORTRAN__) && __GNUC__<11
#if defined(__GFORTRAN__)
module procedure HDF5_write_real1
module procedure HDF5_write_real2
module procedure HDF5_write_real3
@ -1214,7 +1214,7 @@ subroutine HDF5_read_int7(dataset,loc_id,datasetName,parallel)
end subroutine HDF5_read_int7
#if defined(__GFORTRAN__) && __GNUC__<11
#if defined(__GFORTRAN__)
!--------------------------------------------------------------------------------------------------
!> @brief write dataset of type real with 1 dimension
@ -1631,7 +1631,7 @@ subroutine HDF5_write_str(dataset,loc_id,datasetName)
end subroutine HDF5_write_str
#if defined(__GFORTRAN__) && __GNUC__<11
#if defined(__GFORTRAN__)
!--------------------------------------------------------------------------------------------------
!> @brief write dataset of type integer with 1 dimension

View File

@ -428,6 +428,8 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
msg = 'too many systems requested'
case (146)
msg = 'number of values does not match'
case (147)
msg = 'V_e needs to be symmetric'
case (148)
msg = 'Nconstituents mismatch between homogenization and material'

View File

@ -27,7 +27,7 @@ module material
type(tRotationContainer), dimension(:), allocatable, public, protected :: material_O_0
type(tTensorContainer), dimension(:), allocatable, public, protected :: material_F_i_0
type(tTensorContainer), dimension(:), allocatable, public, protected :: material_V_e_0
integer, dimension(:), allocatable, public, protected :: &
homogenization_Nconstituents !< number of grains in each homogenization
@ -133,7 +133,7 @@ subroutine parse()
allocate(material_v(homogenization_maxNconstituents,discretization_Ncells),source=0.0_pReal)
allocate(material_O_0(materials%length))
allocate(material_F_i_0(materials%length))
allocate(material_V_e_0(materials%length))
allocate(ho_of(materials%length))
allocate(ph_of(materials%length,homogenization_maxNconstituents),source=-1)
@ -154,7 +154,7 @@ subroutine parse()
if (constituents%length /= homogenization%get_asInt('N_constituents')) call IO_error(148)
allocate(material_O_0(ma)%data(constituents%length))
allocate(material_F_i_0(ma)%data(1:3,1:3,constituents%length))
allocate(material_V_e_0(ma)%data(1:3,1:3,constituents%length))
do co = 1, constituents%length
constituent => constituents%get(co)
@ -162,7 +162,9 @@ subroutine parse()
ph_of(ma,co) = phases%getIndex(constituent%get_asString('phase'))
call material_O_0(ma)%data(co)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4))
material_F_i_0(ma)%data(1:3,1:3,co) = constituent%get_as2dFloat('F_i',defaultVal=math_I3,requiredShape=[3,3])
material_V_e_0(ma)%data(1:3,1:3,co) = constituent%get_as2dFloat('V_e',defaultVal=math_I3,requiredShape=[3,3])
if (any(dNeq(material_V_e_0(ma)%data(1:3,1:3,co),transpose(material_V_e_0(ma)%data(1:3,1:3,co))))) &
call IO_error(147)
end do
if (dNeq(sum(v_of(ma,:)),1.0_pReal,1.e-9_pReal)) call IO_error(153,ext_msg='constituent')

View File

@ -210,7 +210,6 @@ module subroutine mechanical_init(phases)
Nmembers
class(tNode), pointer :: &
num_crystallite, &
phase, &
mech
@ -239,11 +238,8 @@ module subroutine mechanical_init(phases)
allocate(phase_mechanical_Fe(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Fi(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Fi0(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Fp(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Fp0(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_F(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_F0(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Li(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_Li0(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_Lp(ph)%data(3,3,Nmembers),source=0.0_pReal)
@ -265,23 +261,19 @@ module subroutine mechanical_init(phases)
ma = discretization_materialAt((ce-1)/discretization_nIPs+1)
do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
ph = material_phaseID(co,ce)
phase_mechanical_Fi0(ph)%data(1:3,1:3,material_phaseEntry(co,ce)) = material_F_i_0(ma)%data(1:3,1:3,co)
en = material_phaseEntry(co,ce)
phase_mechanical_F(ph)%data(1:3,1:3,en) = math_I3
phase_mechanical_Fp(ph)%data(1:3,1:3,en) = material_O_0(ma)%data(co)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
phase_mechanical_Fe(ph)%data(1:3,1:3,en) = matmul(material_V_e_0(ma)%data(1:3,1:3,co), &
transpose(phase_mechanical_Fp(ph)%data(1:3,1:3,en)))
phase_mechanical_Fi(ph)%data(1:3,1:3,en) = material_O_0(ma)%data(co)%rotate(math_inv33(material_V_e_0(ma)%data(1:3,1:3,co)))
enddo
enddo
do ph = 1, phases%length
do en = 1, count(material_phaseID == ph)
phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = phase_O_0(ph)%data(en)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) &
/ math_det33(phase_mechanical_Fp0(ph)%data(1:3,1:3,en))**(1.0_pReal/3.0_pReal)
phase_mechanical_F0(ph)%data(1:3,1:3,en) = math_I3
phase_mechanical_Fe(ph)%data(1:3,1:3,en) = math_inv33(matmul(phase_mechanical_Fi0(ph)%data(1:3,1:3,en), &
phase_mechanical_Fp0(ph)%data(1:3,1:3,en))) ! assuming that euler angles are given in internal strain free configuration
enddo
phase_mechanical_Fp(ph)%data = phase_mechanical_Fp0(ph)%data
phase_mechanical_Fi(ph)%data = phase_mechanical_Fi0(ph)%data
phase_mechanical_F(ph)%data = phase_mechanical_F0(ph)%data
phase_mechanical_F0(ph)%data = phase_mechanical_F(ph)%data
phase_mechanical_Fp0(ph)%data = phase_mechanical_Fp(ph)%data
phase_mechanical_Fi0(ph)%data = phase_mechanical_Fi(ph)%data
enddo
@ -318,7 +310,6 @@ module subroutine mechanical_init(phases)
end select
call eigen_init(phases)
@ -1244,6 +1235,9 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
end function phase_mechanical_dPdF
!--------------------------------------------------------------------------------------------------
!< @brief Write restart information to file.
!--------------------------------------------------------------------------------------------------
module subroutine mechanical_restartWrite(groupHandle,ph)
integer(HID_T), intent(in) :: groupHandle
@ -1251,36 +1245,41 @@ module subroutine mechanical_restartWrite(groupHandle,ph)
call HDF5_write(plasticState(ph)%state,groupHandle,'omega_plastic')
call HDF5_write(phase_mechanical_Fi(ph)%data,groupHandle,'F_i')
call HDF5_write(phase_mechanical_Li(ph)%data,groupHandle,'L_i')
call HDF5_write(phase_mechanical_Lp(ph)%data,groupHandle,'L_p')
call HDF5_write(phase_mechanical_Fp(ph)%data,groupHandle,'F_p')
call HDF5_write(phase_mechanical_S(ph)%data,groupHandle,'S')
call HDF5_write(phase_mechanical_F(ph)%data,groupHandle,'F')
call HDF5_write(phase_mechanical_Fp(ph)%data,groupHandle,'F_p')
call HDF5_write(phase_mechanical_Fi(ph)%data,groupHandle,'F_i')
call HDF5_write(phase_mechanical_Lp(ph)%data,groupHandle,'L_p')
call HDF5_write(phase_mechanical_Li(ph)%data,groupHandle,'L_i')
end subroutine mechanical_restartWrite
!--------------------------------------------------------------------------------------------------
!< @brief Read restart information from file.
!--------------------------------------------------------------------------------------------------
module subroutine mechanical_restartRead(groupHandle,ph)
integer(HID_T), intent(in) :: groupHandle
integer, intent(in) :: ph
integer :: en
call HDF5_read(plasticState(ph)%state0,groupHandle,'omega_plastic')
call HDF5_read(phase_mechanical_Fi0(ph)%data,groupHandle,'F_i')
call HDF5_read(phase_mechanical_Li0(ph)%data,groupHandle,'L_i')
call HDF5_read(phase_mechanical_Lp0(ph)%data,groupHandle,'L_p')
call HDF5_read(phase_mechanical_Fp0(ph)%data,groupHandle,'F_p')
call HDF5_read(phase_mechanical_S0(ph)%data,groupHandle,'S')
call HDF5_read(phase_mechanical_F0(ph)%data,groupHandle,'F')
call HDF5_read(phase_mechanical_Fp0(ph)%data,groupHandle,'F_p')
call HDF5_read(phase_mechanical_Fi0(ph)%data,groupHandle,'F_i')
call HDF5_read(phase_mechanical_Lp0(ph)%data,groupHandle,'L_p')
call HDF5_read(phase_mechanical_Li0(ph)%data,groupHandle,'L_i')
end subroutine mechanical_restartRead
!----------------------------------------------------------------------------------------------
!< @brief Get first Piola-Kichhoff stress (for use by non-mech physics)
!----------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
!< @brief Get first Piola-Kichhoff stress (for use by non-mech physics).
!--------------------------------------------------------------------------------------------------
module function mechanical_S(ph,en) result(S)
integer, intent(in) :: ph,en
@ -1292,9 +1291,9 @@ module function mechanical_S(ph,en) result(S)
end function mechanical_S
!----------------------------------------------------------------------------------------------
!< @brief Get plastic velocity gradient (for use by non-mech physics)
!----------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
!< @brief Get plastic velocity gradient (for use by non-mech physics).
!--------------------------------------------------------------------------------------------------
module function mechanical_L_p(ph,en) result(L_p)
integer, intent(in) :: ph,en
@ -1306,9 +1305,9 @@ module function mechanical_L_p(ph,en) result(L_p)
end function mechanical_L_p
!----------------------------------------------------------------------------------------------
!< @brief Get elastic deformation gradient (for use by non-mech physics)
!----------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
!< @brief Get elastic deformation gradient (for use by non-mech physics).
!--------------------------------------------------------------------------------------------------
module function mechanical_F_e(ph,en) result(F_e)
integer, intent(in) :: ph,en
@ -1320,9 +1319,9 @@ module function mechanical_F_e(ph,en) result(F_e)
end function mechanical_F_e
!----------------------------------------------------------------------------------------------
!< @brief Get second Piola-Kichhoff stress (for use by homogenization)
!----------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
!< @brief Get second Piola-Kichhoff stress (for use by homogenization).
!--------------------------------------------------------------------------------------------------
module function phase_P(co,ce) result(P)
integer, intent(in) :: co, ce
@ -1334,9 +1333,9 @@ module function phase_P(co,ce) result(P)
end function phase_P
!----------------------------------------------------------------------------------------------
!< @brief Get deformation gradient (for use by homogenization)
!----------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
!< @brief Get deformation gradient (for use by homogenization).
!--------------------------------------------------------------------------------------------------
module function phase_F(co,ce) result(F)
integer, intent(in) :: co, ce
@ -1348,9 +1347,9 @@ module function phase_F(co,ce) result(F)
end function phase_F
!----------------------------------------------------------------------------------------------
!< @brief Set deformation gradient (for use by homogenization)
!----------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
!< @brief Set deformation gradient (for use by homogenization).
!--------------------------------------------------------------------------------------------------
module subroutine phase_set_F(F,co,ce)
real(pReal), dimension(3,3), intent(in) :: F

View File

@ -323,17 +323,16 @@ pure function rotTensor2(self,T,active) result(tRot)
logical :: passive
if (present(active)) then
passive = .not. active
else
passive = .true.
endif
if (passive) then
tRot = matmul(matmul(self%asMatrix(),T),transpose(self%asMatrix()))
else
tRot = matmul(matmul(transpose(self%asMatrix()),T),self%asMatrix())
endif
tRot = merge(matmul(matmul(self%asMatrix(),T),transpose(self%asMatrix())), &
matmul(matmul(transpose(self%asMatrix()),T),self%asMatrix()), &
passive)
end function rotTensor2
@ -450,6 +449,7 @@ pure function qu2om(qu) result(om)
om(1,3) = 2.0_pReal*(qu(2)*qu(4)+qu(1)*qu(3))
if (sign(1.0_pReal,P) < 0.0_pReal) om = transpose(om)
om = om/math_det33(om)**(1.0_pReal/3.0_pReal)
end function qu2om
@ -575,7 +575,7 @@ end function qu2cu
!---------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief convert rotation matrix to cubochoric
!> @brief convert rotation matrix to unit quaternion
!> @details the original formulation (direct conversion) had (numerical?) issues
!---------------------------------------------------------------------------------------------------
pure function om2qu(om) result(qu)
@ -602,14 +602,15 @@ pure function om2qu(om) result(qu)
endif
endif
if(sign(1.0_pReal,qu(1))<0.0_pReal) qu =-1.0_pReal * qu
qu = qu*[1.0_pReal,P,P,P]
qu(2:4) = merge(qu(2:4),qu(2:4)*P,dEq0(qu(2:4)))
qu = qu/norm2(qu)
end function om2qu
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief orientation matrix to Euler angles
!> @brief convert orientation matrix to Euler angles
!> @details Two step check for special cases to avoid invalid operations (not needed for python)
!---------------------------------------------------------------------------------------------------
pure function om2eu(om) result(eu)
@ -1334,8 +1335,8 @@ pure function cu2ho(cu) result(ho)
! transform to sphere grid (inverse Lambert)
! [note that there is no need to worry about dividing by zero, since XYZ(3) can not become zero]
c = sum(T**2)
s = PI * c/(24.0*XYZ(3)**2)
c = sqrt(PI) * c / sqrt(24.0_pReal) / XYZ(3)
s = c * PI/(24.0*XYZ(3)**2)
c = c * sqrt(PI/24.0_pReal) / XYZ(3)
q = sqrt( 1.0 - s )
LamXYZ = [ T(order(2)) * q, T(order(1)) * q, PREF * XYZ(3) - c ]
end if special