DAMASK_EICMD/python/damask/_result.py

1356 lines
53 KiB
Python
Raw Normal View History

2020-06-26 15:15:06 +05:30
import multiprocessing as mp
import re
2021-04-04 23:14:06 +05:30
import fnmatch
import os
2021-04-05 11:23:19 +05:30
import copy
2020-05-25 22:20:31 +05:30
import datetime
import xml.etree.ElementTree as ET
import xml.dom.minidom
2020-06-03 14:13:07 +05:30
from pathlib import Path
from functools import partial
from collections import defaultdict
2021-03-31 18:49:23 +05:30
from collections.abc import Iterable
import h5py
2019-04-17 23:27:16 +05:30
import numpy as np
import numpy.ma as ma
from numpy.lib import recfunctions as rfn
2020-06-28 15:10:19 +05:30
import damask
2020-03-11 11:20:13 +05:30
from . import VTK
2020-02-16 00:39:24 +05:30
from . import Orientation
2020-02-21 22:17:47 +05:30
from . import grid_filters
2020-03-13 00:22:33 +05:30
from . import mechanics
from . import tensor
2020-03-13 00:22:33 +05:30
from . import util
h5py3 = h5py.__version__[0] == '3'
2021-04-03 14:38:22 +05:30
def _read(dataset):
2021-04-05 11:23:19 +05:30
"""Read a dataset and its metadata into a numpy.ndarray."""
2021-04-05 13:59:34 +05:30
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)
return np.array(dataset,dtype=dtype)
2021-03-31 15:35:51 +05:30
2021-04-03 14:38:22 +05:30
def _match(requested,existing):
"""Find matches among two sets of labels."""
2021-04-03 14:38:22 +05:30
def flatten_list(list_of_lists):
return [e for e_ in list_of_lists for e in e_]
if requested is True:
requested = '*'
elif requested is False or requested is None:
requested = []
requested_ = requested if hasattr(requested,'__iter__') and not isinstance(requested,str) else \
[requested]
2021-04-04 23:14:06 +05:30
return sorted(set(flatten_list([fnmatch.filter(existing,r) for r in requested_])),
2021-04-03 14:38:22 +05:30
key=util.natural_sort)
2021-04-06 21:09:44 +05:30
def _empty_like(dataset,N_materialpoints,fill_float,fill_int):
2021-04-05 11:23:19 +05:30
"""Create empty numpy.ma.MaskedArray."""
2021-04-04 16:55:16 +05:30
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,
mask = True)
2021-03-31 15:35:51 +05:30
class Result:
2019-09-20 01:02:15 +05:30
"""
Manipulate and read DADF5 files.
2019-09-16 08:49:14 +05:30
2020-03-19 12:15:31 +05:30
DADF5 (DAMASK HDF5) files contain DAMASK results.
2021-03-31 20:41:53 +05:30
The group/folder structure reflects the input data
in material.yaml.
2020-02-21 12:15:05 +05:30
"""
def __init__(self,fname):
"""
New result view bound to a HDF5 file.
2019-09-16 08:49:14 +05:30
2020-02-21 12:15:05 +05:30
Parameters
----------
2020-11-19 18:15:40 +05:30
fname : str or pathlib.Path
Name of the DADF5 file to be opened.
2019-09-20 01:02:15 +05:30
2020-02-21 12:15:05 +05:30
"""
with h5py.File(fname,'r') as f:
self.version_major = f.attrs['DADF5_version_major']
self.version_minor = f.attrs['DADF5_version_minor']
2021-03-25 23:52:59 +05:30
if self.version_major != 0 or not 7 <= self.version_minor <= 12:
2020-06-25 01:04:51 +05:30
raise TypeError(f'Unsupported DADF5 version {self.version_major}.{self.version_minor}')
2020-03-22 20:43:35 +05:30
2020-12-04 02:28:24 +05:30
self.structured = 'grid' in f['geometry'].attrs.keys() or \
'cells' in f['geometry'].attrs.keys()
2020-03-22 20:43:35 +05:30
if self.structured:
2020-12-04 02:28:24 +05:30
try:
self.cells = f['geometry'].attrs['cells']
except KeyError:
self.cells = f['geometry'].attrs['grid']
2020-03-22 20:43:35 +05:30
self.size = f['geometry'].attrs['size']
self.origin = f['geometry'].attrs['origin']
2020-03-22 20:43:35 +05:30
2021-03-25 23:52:59 +05:30
r=re.compile('inc[0-9]+' if self.version_minor < 12 else 'increment_[0-9]+')
2021-04-03 14:38:22 +05:30
self.increments = sorted([i for i in f.keys() if r.match(i)],key=util.natural_sort)
2021-04-06 21:09:44 +05:30
self.times = [round(f[i].attrs['time/s' if self.version_minor < 12 else
't/s'],12) for i in self.increments]
2020-03-22 20:43:35 +05:30
2021-03-25 23:52:59 +05:30
grp = 'mapping' if self.version_minor < 12 else 'cell_to'
2021-03-25 23:52:59 +05:30
self.N_materialpoints, self.N_constituents = np.shape(f[f'{grp}/phase'])
self.homogenizations = [m.decode() for m in np.unique(f[f'{grp}/homogenization']
['Name' if self.version_minor < 12 else 'label'])]
2021-04-03 14:38:22 +05:30
self.homogenizations = sorted(self.homogenizations,key=util.natural_sort)
2021-03-25 23:52:59 +05:30
self.phases = [c.decode() for c in np.unique(f[f'{grp}/phase']
['Name' if self.version_minor < 12 else 'label'])]
2021-04-03 14:38:22 +05:30
self.phases = sorted(self.phases,key=util.natural_sort)
2020-03-22 20:43:35 +05:30
self.fields = []
for c in self.phases:
self.fields += f['/'.join([self.increments[0],'phase',c])].keys()
for m in self.homogenizations:
self.fields += f['/'.join([self.increments[0],'homogenization',m])].keys()
2021-04-03 14:38:22 +05:30
self.fields = sorted(set(self.fields),key=util.natural_sort) # make unique
2021-01-13 19:27:58 +05:30
self.visible = {'increments': self.increments,
'phases': self.phases,
'homogenizations': self.homogenizations,
'fields': self.fields,
2021-01-13 19:27:58 +05:30
}
2020-06-03 14:13:07 +05:30
self.fname = Path(fname).absolute()
self._allow_modification = False
2021-04-05 11:23:19 +05:30
def __copy__(self):
"""Create deep copy."""
return copy.deepcopy(self)
copy = __copy__
def __repr__(self):
2020-11-05 11:45:59 +05:30
"""Show summary of file content."""
2021-01-13 19:27:58 +05:30
visible_increments = self.visible['increments']
2020-04-21 14:47:15 +05:30
2021-04-05 11:23:19 +05:30
first = self.view('increments',visible_increments[0:1]).list_data()
2020-04-21 14:47:15 +05:30
2021-04-05 11:23:19 +05:30
last = '' if len(visible_increments) < 2 else \
self.view('increments',visible_increments[-1:]).list_data()
2020-04-21 14:47:15 +05:30
2021-01-13 19:27:58 +05:30
in_between = '' if len(visible_increments) < 3 else \
''.join([f'\n{inc}\n ...\n' for inc in visible_increments[1:-1]])
2020-04-21 14:47:15 +05:30
return util.srepr(first + in_between + last)
2021-01-13 19:27:58 +05:30
def _manage_view(self,action,what,datasets):
2020-02-21 12:15:05 +05:30
"""
Manages the visibility of the groups.
2020-02-21 12:15:05 +05:30
Parameters
----------
2020-03-03 04:17:29 +05:30
action : str
2020-11-23 23:52:48 +05:30
Select from 'set', 'add', and 'del'.
2020-03-03 04:17:29 +05:30
what : str
Attribute to change (must be from self.visible).
2021-04-06 21:09:44 +05:30
datasets : (list of) int (for increments), (list of) float (for times), (list of) str, or bool
Name of datasets; supports '?' and '*' wildcards.
True is equivalent to '*', False is equivalent to [].
2020-03-03 04:17:29 +05:30
2020-02-21 12:15:05 +05:30
"""
# allow True/False and string arguments
2020-05-21 14:15:52 +05:30
if datasets is True:
2021-04-03 14:38:22 +05:30
datasets = '*'
elif datasets is False or datasets is None:
2020-03-03 04:17:29 +05:30
datasets = []
2021-03-30 23:27:49 +05:30
choice = list(datasets).copy() if hasattr(datasets,'__iter__') and not isinstance(datasets,str) else \
2020-03-03 18:37:02 +05:30
[datasets]
2021-03-25 23:52:59 +05:30
inc = 'inc' if self.version_minor < 12 else 'increment_' # compatibility hack
2020-03-03 18:37:02 +05:30
if what == 'increments':
2021-03-25 23:52:59 +05:30
choice = [c if isinstance(c,str) and c.startswith(inc) else
f'{inc}{c}' for c in choice]
2020-03-03 18:37:02 +05:30
elif what == 'times':
what = 'increments'
if choice == ['*']:
choice = self.increments
else:
iterator = map(float,choice)
choice = []
for c in iterator:
2020-03-22 21:33:28 +05:30
idx = np.searchsorted(self.times,c)
if idx >= len(self.times): continue
2020-03-03 18:37:02 +05:30
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])
2021-04-03 14:38:22 +05:30
valid = _match(choice,getattr(self,what))
2021-01-13 19:27:58 +05:30
existing = set(self.visible[what])
2020-02-21 12:15:05 +05:30
2021-04-05 11:23:19 +05:30
dup = self.copy()
2020-02-21 12:15:05 +05:30
if action == 'set':
2021-04-05 11:23:19 +05:30
dup.visible[what] = sorted(set(valid), key=util.natural_sort)
2020-02-21 12:15:05 +05:30
elif action == 'add':
2020-03-22 21:33:28 +05:30
add = existing.union(valid)
2021-04-05 11:23:19 +05:30
dup.visible[what] = sorted(add, key=util.natural_sort)
2020-02-21 12:15:05 +05:30
elif action == 'del':
2020-03-22 21:33:28 +05:30
diff = existing.difference(valid)
2021-04-05 11:23:19 +05:30
dup.visible[what] = sorted(diff, key=util.natural_sort)
return dup
2020-02-21 12:15:05 +05:30
def allow_modification(self):
2020-11-19 18:15:40 +05:30
"""Allow to overwrite existing data."""
print(util.warn('Warning: Modification of existing datasets allowed!'))
2021-04-05 11:23:19 +05:30
dup = self.copy()
dup._allow_modification = True
return dup
2020-06-01 15:03:22 +05:30
def disallow_modification(self):
"""Disallow to overwrite existing data (default case)."""
2021-04-05 11:23:19 +05:30
dup = self.copy()
dup._allow_modification = False
return dup
2021-03-25 23:52:59 +05:30
def increments_in_range(self,start,end):
2020-11-19 18:15:40 +05:30
"""
Select all increments within a given range.
Parameters
----------
start : int or str
Start increment.
end : int or str
End increment.
"""
2021-03-25 23:52:59 +05:30
# compatibility hack
ln = 3 if self.version_minor < 12 else 10
selected = []
2021-03-25 23:52:59 +05:30
for i,inc in enumerate([int(i[ln:]) for i in self.increments]):
s,e = map(lambda x: int(x[ln:] if isinstance(x,str) and x.startswith('inc') else x), (start,end))
if s <= inc <= e:
selected.append(self.increments[i])
return selected
def times_in_range(self,start,end):
2020-11-19 18:15:40 +05:30
"""
Select all increments within a given time range.
Parameters
----------
start : float
Time of start increment.
end : float
Time of end increment.
"""
selected = []
for i,time in enumerate(self.times):
if start <= time <= end:
2020-03-10 03:58:25 +05:30
selected.append(self.times[i])
return selected
2021-01-13 19:27:58 +05:30
def view(self,what,datasets):
2020-02-21 12:15:05 +05:30
"""
2021-01-13 19:27:58 +05:30
Set view.
2020-02-21 12:15:05 +05:30
Parameters
----------
2021-04-05 13:59:34 +05:30
what : {'increments', 'times', 'phases', 'homogenizations', 'fields'}
Attribute to change.
2021-04-06 21:09:44 +05:30
datasets : (list of) int (for increments), (list of) float (for times), (list of) str, or bool
Name of datasets; supports '?' and '*' wildcards.
True is equivalent to '*', False is equivalent to [].
2020-02-21 12:15:05 +05:30
"""
2021-04-05 11:23:19 +05:30
return self._manage_view('set',what,datasets)
2021-01-13 19:27:58 +05:30
def view_more(self,what,datasets):
2020-02-21 12:15:05 +05:30
"""
2021-01-13 19:27:58 +05:30
Add to view.
2020-02-21 12:15:05 +05:30
Parameters
----------
2021-04-05 13:59:34 +05:30
what : {'increments', 'times', 'phases', 'homogenizations', 'fields'}
Attribute to change.
2021-04-06 21:09:44 +05:30
datasets : (list of) int (for increments), (list of) float (for times), (list of) str, or bool
Name of datasets; supports '?' and '*' wildcards.
True is equivalent to '*', False is equivalent to [].
2020-02-21 12:15:05 +05:30
"""
2021-04-05 11:23:19 +05:30
return self._manage_view('add',what,datasets)
2021-01-13 19:27:58 +05:30
def view_less(self,what,datasets):
2020-02-21 12:15:05 +05:30
"""
2021-01-13 19:27:58 +05:30
Delete from view.
2020-02-21 12:15:05 +05:30
Parameters
----------
2021-04-05 13:59:34 +05:30
what : {'increments', 'times', 'phases', 'homogenizations', 'fields'}
Attribute to change.
2021-04-06 21:09:44 +05:30
datasets : (list of) int (for increments), (list of) float (for times), (list of) str, or bool
Name of datasets; supports '?' and '*' wildcards.
True is equivalent to '*', False is equivalent to [].
2019-10-19 16:40:46 +05:30
2020-02-21 12:15:05 +05:30
"""
2021-04-05 11:23:19 +05:30
return self._manage_view('del',what,datasets)
def rename(self,name_old,name_new):
"""
Rename dataset.
2020-06-25 01:04:51 +05:30
2020-06-02 01:43:01 +05:30
Parameters
----------
name_old : str
Name of the dataset to be renamed.
2020-06-02 01:43:01 +05:30
name_new : str
New name of the dataset.
"""
2021-04-04 23:14:06 +05:30
if not self._allow_modification:
2020-06-01 15:03:22 +05:30
raise PermissionError('Rename operation not permitted')
2021-04-04 23:14:06 +05:30
with h5py.File(self.fname,'a') as f:
for inc in self.visible['increments']:
2021-04-05 11:23:19 +05:30
for ty in ['phase','homogenization']:
for label in self.visible[ty+'s']:
2021-04-04 23:14:06 +05:30
for field in self.visible['fields']:
2021-04-05 11:23:19 +05:30
path_old = '/'.join([inc,ty,label,field,name_old])
path_new = '/'.join([inc,ty,label,field,name_new])
2021-04-04 23:14:06 +05:30
if path_old in f.keys():
f[path_new] = f[path_old]
f[path_new].attrs['renamed'] = f'original name: {name_old}' if h5py3 else \
f'original name: {name_old}'.encode()
del f[path_old]
2020-02-21 12:15:05 +05:30
def list_data(self):
"""Return information on all active datasets in the file."""
2021-03-25 23:52:59 +05:30
# compatibility hack
de = 'Description' if self.version_minor < 12 else 'description'
un = 'Unit' if self.version_minor < 12 else 'unit'
2021-04-05 11:23:19 +05:30
msg = ''
2020-02-21 12:15:05 +05:30
with h5py.File(self.fname,'r') as f:
2021-04-04 15:57:39 +05:30
for inc in self.visible['increments']:
2021-04-05 11:23:19 +05:30
msg = ''.join([msg,f'\n{inc} ({self.times[self.increments.index(inc)]}s)\n'])
for ty in ['phase','homogenization']:
msg = ' '.join([msg,f'{ty}\n'])
for label in self.visible[ty+'s']:
msg = ' '.join([msg,f'{label}\n'])
2021-04-04 15:57:39 +05:30
for field in self.visible['fields']:
2021-04-05 11:23:19 +05:30
msg = ' '.join([msg,f'{field}\n'])
for d in f['/'.join([inc,ty,label,field])].keys():
dataset = f['/'.join([inc,ty,label,field,d])]
2021-04-03 16:41:20 +05:30
unit = f' / {dataset.attrs[un]}' if h5py3 else \
f' / {dataset.attrs[un].decode()}'
description = dataset.attrs[de] if h5py3 else \
dataset.attrs[de].decode()
2021-04-05 11:23:19 +05:30
msg = ' '.join([msg,f'{d}{unit}: {description}\n'])
2020-02-21 12:15:05 +05:30
2021-04-05 11:23:19 +05:30
return msg
2020-02-21 12:15:05 +05:30
2020-07-03 10:59:31 +05:30
def enable_user_function(self,func):
globals()[func.__name__]=func
print(f'Function {func.__name__} enabled in add_calculation.')
2020-07-31 20:20:01 +05:30
@property
def coordinates0_point(self):
2020-02-21 12:15:05 +05:30
"""Return initial coordinates of the cell centers."""
if self.structured:
return grid_filters.coordinates0_point(self.cells,self.size,self.origin).reshape(-1,3,order='F')
2020-02-21 12:15:05 +05:30
else:
2020-02-21 22:17:47 +05:30
with h5py.File(self.fname,'r') as f:
return f['geometry/x_c'][()]
2020-02-21 12:15:05 +05:30
2020-07-31 20:20:01 +05:30
@property
def coordinates0_node(self):
2020-04-22 11:10:02 +05:30
"""Return initial coordinates of the cell centers."""
if self.structured:
return grid_filters.coordinates0_node(self.cells,self.size,self.origin).reshape(-1,3,order='F')
2020-04-22 11:10:02 +05:30
else:
with h5py.File(self.fname,'r') as f:
return f['geometry/x_n'][()]
@property
def geometry0(self):
if self.structured:
return VTK.from_rectilinear_grid(self.cells,self.size,self.origin)
else:
with h5py.File(self.fname,'r') as f:
return VTK.from_unstructured_grid(f['/geometry/x_n'][()],
f['/geometry/T_c'][()]-1,
f['/geometry/T_c'].attrs['VTK_TYPE'] if h5py3 else \
f['/geometry/T_c'].attrs['VTK_TYPE'].decode())
2020-02-21 12:15:05 +05:30
@staticmethod
def _add_absolute(x):
return {
'data': np.abs(x['data']),
2020-06-25 01:04:51 +05:30
'label': f'|{x["label"]}|',
'meta': {
2021-03-25 23:52:59 +05:30
'unit': x['meta']['unit'],
'description': f"absolute value of {x['label']} ({x['meta']['description']})",
'creator': 'add_absolute'
}
}
2020-02-21 12:15:05 +05:30
def add_absolute(self,x):
"""
Add absolute value.
Parameters
----------
x : str
2020-03-19 12:15:31 +05:30
Label of scalar, vector, or tensor dataset to take absolute value of.
2020-02-21 12:15:05 +05:30
"""
self._add_generic_pointwise(self._add_absolute,{'x':x})
2020-02-21 12:15:05 +05:30
@staticmethod
def _add_calculation(**kwargs):
formula = kwargs['formula']
for d in re.findall(r'#(.*?)#',formula):
2020-06-25 01:04:51 +05:30
formula = formula.replace(f'#{d}#',f"kwargs['{d}']['data']")
2020-02-21 12:15:05 +05:30
return {
'data': eval(formula),
'label': kwargs['label'],
'meta': {
2021-03-25 23:52:59 +05:30
'unit': kwargs['unit'],
'description': f"{kwargs['description']} (formula: {kwargs['formula']})",
'creator': 'add_calculation'
}
}
def add_calculation(self,label,formula,unit='n/a',description=None):
2020-02-21 12:15:05 +05:30
"""
Add result of a general formula.
Parameters
----------
label : str
Label of resulting dataset.
formula : str
2021-04-06 21:09:44 +05:30
Formula to calculate resulting dataset. Existing datasets are referenced by '#TheirLabel#'.
2020-02-21 12:15:05 +05:30
unit : str, optional
2020-03-19 12:15:31 +05:30
Physical unit of the result.
2020-02-21 12:15:05 +05:30
description : str, optional
2020-03-19 12:15:31 +05:30
Human-readable description of the result.
2020-02-21 12:15:05 +05:30
"""
dataset_mapping = {d:d for d in set(re.findall(r'#(.*?)#',formula))} # datasets used in the formula
2020-02-21 12:15:05 +05:30
args = {'formula':formula,'label':label,'unit':unit,'description':description}
self._add_generic_pointwise(self._add_calculation,dataset_mapping,args)
@staticmethod
2020-11-18 03:26:22 +05:30
def _add_stress_Cauchy(P,F):
return {
2020-11-18 03:26:22 +05:30
'data': mechanics.stress_Cauchy(P['data'],F['data']),
'label': 'sigma',
'meta': {
2021-03-25 23:52:59 +05:30
'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'
}
}
2020-11-18 03:26:22 +05:30
def add_stress_Cauchy(self,P='P',F='F'):
2020-02-21 12:15:05 +05:30
"""
Add Cauchy stress calculated from first Piola-Kirchhoff stress and deformation gradient.
2020-02-21 12:15:05 +05:30
Parameters
----------
P : str, optional
2021-04-06 21:09:44 +05:30
Label of the dataset containing the first Piola-Kirchhoff stress. Defaults to 'P'.
2020-02-21 12:15:05 +05:30
F : str, optional
2021-04-06 21:09:44 +05:30
Label of the dataset containing the deformation gradient. Defaults to 'F'.
2020-02-21 12:15:05 +05:30
"""
2020-11-18 03:26:22 +05:30
self._add_generic_pointwise(self._add_stress_Cauchy,{'P':P,'F':F})
2020-02-21 12:15:05 +05:30
@staticmethod
def _add_determinant(T):
return {
'data': np.linalg.det(T['data']),
2020-06-25 01:04:51 +05:30
'label': f"det({T['label']})",
'meta': {
2021-03-25 23:52:59 +05:30
'unit': T['meta']['unit'],
'description': f"determinant of tensor {T['label']} ({T['meta']['description']})",
'creator': 'add_determinant'
}
}
def add_determinant(self,T):
2020-02-21 12:15:05 +05:30
"""
Add the determinant of a tensor.
Parameters
----------
T : str
2020-03-19 12:15:31 +05:30
Label of tensor dataset.
2020-02-21 12:15:05 +05:30
"""
self._add_generic_pointwise(self._add_determinant,{'T':T})
2020-02-21 12:15:05 +05:30
@staticmethod
def _add_deviator(T):
return {
2020-11-19 19:08:54 +05:30
'data': tensor.deviatoric(T['data']),
2020-06-25 01:04:51 +05:30
'label': f"s_{T['label']}",
'meta': {
2021-03-25 23:52:59 +05:30
'unit': T['meta']['unit'],
'description': f"deviator of tensor {T['label']} ({T['meta']['description']})",
'creator': 'add_deviator'
}
}
def add_deviator(self,T):
2020-02-21 12:15:05 +05:30
"""
Add the deviatoric part of a tensor.
Parameters
----------
T : str
2020-03-19 12:15:31 +05:30
Label of tensor dataset.
2020-02-21 12:15:05 +05:30
"""
self._add_generic_pointwise(self._add_deviator,{'T':T})
2020-02-21 12:15:05 +05:30
@staticmethod
def _add_eigenvalue(T_sym,eigenvalue):
if eigenvalue == 'max':
2021-03-25 23:52:59 +05:30
label,p = 'maximum',2
elif eigenvalue == 'mid':
2021-03-25 23:52:59 +05:30
label,p = 'intermediate',1
elif eigenvalue == 'min':
2021-03-25 23:52:59 +05:30
label,p = 'minimum',0
return {
'data': tensor.eigenvalues(T_sym['data'])[:,p],
2020-06-25 01:04:51 +05:30
'label': f"lambda_{eigenvalue}({T_sym['label']})",
'meta' : {
2021-03-25 23:52:59 +05:30
'unit': T_sym['meta']['unit'],
'description': f"{label} eigenvalue of {T_sym['label']} ({T_sym['meta']['description']})",
'creator': 'add_eigenvalue'
}
}
def add_eigenvalue(self,T_sym,eigenvalue='max'):
2020-02-21 12:15:05 +05:30
"""
Add eigenvalues of symmetric tensor.
Parameters
----------
2020-03-03 03:44:59 +05:30
T_sym : str
2020-03-19 12:15:31 +05:30
Label of symmetric tensor dataset.
eigenvalue : str, optional
2021-04-06 21:09:44 +05:30
Eigenvalue. Select from 'max', 'mid', 'min'. Defaults to 'max'.
2020-02-21 12:15:05 +05:30
"""
self._add_generic_pointwise(self._add_eigenvalue,{'T_sym':T_sym},{'eigenvalue':eigenvalue})
2020-02-21 12:15:05 +05:30
@staticmethod
def _add_eigenvector(T_sym,eigenvalue):
if eigenvalue == 'max':
label,p = 'maximum',2
elif eigenvalue == 'mid':
label,p = 'intermediate',1
elif eigenvalue == 'min':
label,p = 'minimum',0
return {
'data': tensor.eigenvectors(T_sym['data'])[:,p],
2020-06-25 01:04:51 +05:30
'label': f"v_{eigenvalue}({T_sym['label']})",
'meta' : {
2021-03-25 23:52:59 +05:30
'unit': '1',
'description': f"eigenvector corresponding to {label} eigenvalue"
f" of {T_sym['label']} ({T_sym['meta']['description']})",
'creator': 'add_eigenvector'
}
2020-06-25 01:04:51 +05:30
}
def add_eigenvector(self,T_sym,eigenvalue='max'):
2020-02-21 12:15:05 +05:30
"""
Add eigenvector of symmetric tensor.
2020-02-21 12:15:05 +05:30
Parameters
----------
2020-03-03 03:44:59 +05:30
T_sym : str
2020-03-19 12:15:31 +05:30
Label of symmetric tensor dataset.
eigenvalue : str, optional
2021-04-06 21:09:44 +05:30
Eigenvalue to which the eigenvector corresponds.
Select from 'max', 'mid', 'min'. Defaults to 'max'.
2020-02-21 12:15:05 +05:30
"""
self._add_generic_pointwise(self._add_eigenvector,{'T_sym':T_sym},{'eigenvalue':eigenvalue})
2020-02-21 12:15:05 +05:30
@staticmethod
2020-12-02 19:15:47 +05:30
def _add_IPF_color(l,q):
m = util.scale_to_coprime(np.array(l))
try:
2021-03-25 23:52:59 +05:30
lattice = {'fcc':'cF','bcc':'cI','hex':'hP'}[q['meta']['lattice']]
except KeyError:
2021-03-25 23:52:59 +05:30
lattice = q['meta']['lattice']
try:
o = Orientation(rotation = (rfn.structured_to_unstructured(q['data'])),lattice=lattice)
except ValueError:
o = Orientation(rotation = q['data'],lattice=lattice)
return {
'data': np.uint8(o.IPF_color(l)*255),
'label': 'IPFcolor_({} {} {})'.format(*m),
'meta' : {
2021-03-25 23:52:59 +05:30
'unit': '8-bit RGB',
'lattice': q['meta']['lattice'],
'description': 'Inverse Pole Figure (IPF) colors along sample direction ({} {} {})'.format(*m),
2021-03-25 23:52:59 +05:30
'creator': 'add_IPF_color'
}
}
2020-12-02 19:15:47 +05:30
def add_IPF_color(self,l,q='O'):
2020-02-21 12:15:05 +05:30
"""
Add RGB color tuple of inverse pole figure (IPF) color.
Parameters
----------
2020-02-21 22:12:01 +05:30
l : numpy.array of shape (3)
2020-03-19 12:15:31 +05:30
Lab frame direction for inverse pole figure.
2020-12-02 19:15:47 +05:30
q : str
Label of the dataset containing the crystallographic orientation as quaternions.
Defaults to 'O'.
2020-02-21 12:15:05 +05:30
"""
2020-07-01 01:13:57 +05:30
self._add_generic_pointwise(self._add_IPF_color,{'q':q},{'l':l})
2020-02-21 12:15:05 +05:30
@staticmethod
2020-03-03 03:44:59 +05:30
def _add_maximum_shear(T_sym):
return {
2020-03-03 03:44:59 +05:30
'data': mechanics.maximum_shear(T_sym['data']),
2020-06-25 01:04:51 +05:30
'label': f"max_shear({T_sym['label']})",
'meta': {
2021-03-25 23:52:59 +05:30
'unit': T_sym['meta']['unit'],
'description': f"maximum shear component of {T_sym['label']} ({T_sym['meta']['description']})",
'creator': 'add_maximum_shear'
}
2020-02-21 12:15:05 +05:30
}
2020-03-03 03:44:59 +05:30
def add_maximum_shear(self,T_sym):
2020-02-21 12:15:05 +05:30
"""
Add maximum shear components of symmetric tensor.
Parameters
----------
2020-03-03 03:44:59 +05:30
T_sym : str
2020-03-19 12:15:31 +05:30
Label of symmetric tensor dataset.
2020-02-21 12:15:05 +05:30
"""
2020-03-03 03:44:59 +05:30
self._add_generic_pointwise(self._add_maximum_shear,{'T_sym':T_sym})
2020-02-21 12:15:05 +05:30
@staticmethod
2020-11-18 03:26:22 +05:30
def _add_equivalent_Mises(T_sym,kind):
2020-11-06 04:17:37 +05:30
k = kind
if k is None:
2021-03-25 23:52:59 +05:30
if T_sym['meta']['unit'] == '1':
k = 'strain'
2021-03-25 23:52:59 +05:30
elif T_sym['meta']['unit'] == 'Pa':
k = 'stress'
if k not in ['stress', 'strain']:
2021-04-05 19:11:28 +05:30
raise ValueError(f'invalid von Mises kind {kind}')
2020-02-21 12:15:05 +05:30
return {
2020-11-18 03:26:22 +05:30
'data': (mechanics.equivalent_strain_Mises if k=='strain' else \
mechanics.equivalent_stress_Mises)(T_sym['data']),
2020-06-25 01:04:51 +05:30
'label': f"{T_sym['label']}_vM",
'meta': {
2021-03-25 23:52:59 +05:30
'unit': T_sym['meta']['unit'],
'description': f"Mises equivalent {k} of {T_sym['label']} ({T_sym['meta']['description']})",
'creator': 'add_Mises'
}
}
2020-11-18 03:26:22 +05:30
def add_equivalent_Mises(self,T_sym,kind=None):
2020-02-21 12:15:05 +05:30
"""
Add the equivalent Mises stress or strain of a symmetric tensor.
Parameters
----------
2020-03-03 03:44:59 +05:30
T_sym : str
2020-03-19 12:15:31 +05:30
Label of symmetric tensorial stress or strain dataset.
kind : {'stress', 'strain', None}, optional
Kind of the von Mises equivalent. Defaults to None, in which case
2021-04-06 21:09:44 +05:30
it is selected based on the unit of the dataset ('1' -> strain, 'Pa' -> stress).
2020-02-21 12:15:05 +05:30
"""
2020-11-18 03:26:22 +05:30
self._add_generic_pointwise(self._add_equivalent_Mises,{'T_sym':T_sym},{'kind':kind})
2020-02-21 12:15:05 +05:30
@staticmethod
def _add_norm(x,ord):
o = ord
if len(x['data'].shape) == 2:
axis = 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
2020-02-21 12:15:05 +05:30
return {
'data': np.linalg.norm(x['data'],ord=o,axis=axis,keepdims=True),
2020-06-25 01:04:51 +05:30
'label': f"|{x['label']}|_{o}",
'meta': {
2021-03-25 23:52:59 +05:30
'unit': x['meta']['unit'],
'description': f"{o}-norm of {t} {x['label']} ({x['meta']['description']})",
'creator': 'add_norm'
}
}
2020-02-21 12:15:05 +05:30
def add_norm(self,x,ord=None):
"""
Add the norm of vector or tensor.
Parameters
----------
x : str
2020-03-19 12:15:31 +05:30
Label of vector or tensor dataset.
2021-04-06 21:09:44 +05:30
ord : {non-zero int, inf, -inf, 'fro', 'nuc'}, optional
2020-03-19 12:15:31 +05:30
Order of the norm. inf means NumPys inf object. For details refer to numpy.linalg.norm.
2020-02-21 12:15:05 +05:30
"""
self._add_generic_pointwise(self._add_norm,{'x':x},{'ord':ord})
2020-02-21 12:15:05 +05:30
@staticmethod
2020-11-18 03:26:22 +05:30
def _add_stress_second_Piola_Kirchhoff(P,F):
return {
2020-11-18 03:26:22 +05:30
'data': mechanics.stress_second_Piola_Kirchhoff(P['data'],F['data']),
'label': 'S',
'meta': {
2021-03-25 23:52:59 +05:30
'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'
}
}
2020-11-18 03:26:22 +05:30
def add_stress_second_Piola_Kirchhoff(self,P='P',F='F'):
2020-02-21 12:15:05 +05:30
"""
2020-06-25 01:04:51 +05:30
Add second Piola-Kirchhoff stress calculated from first Piola-Kirchhoff stress and deformation gradient.
2020-02-21 12:15:05 +05:30
Parameters
----------
P : str, optional
2021-04-06 21:09:44 +05:30
Label of first Piola-Kirchhoff stress dataset. Defaults to 'P'.
2020-02-21 12:15:05 +05:30
F : str, optional
2021-04-06 21:09:44 +05:30
Label of deformation gradient dataset. Defaults to 'F'.
2020-02-21 12:15:05 +05:30
"""
2020-11-18 03:26:22 +05:30
self._add_generic_pointwise(self._add_stress_second_Piola_Kirchhoff,{'P':P,'F':F})
2020-02-21 12:15:05 +05:30
# The add_pole functionality needs discussion.
# The new Crystal object can perform such a calculation but the outcome depends on the lattice parameters
# as well as on whether a direction or plane is concerned (see the DAMASK_examples/pole_figure notebook).
# Below code appears to be too simplistic.
# @staticmethod
# def _add_pole(q,p,polar):
# pole = np.array(p)
# unit_pole = pole/np.linalg.norm(pole)
# m = util.scale_to_coprime(pole)
# rot = Rotation(q['data'].view(np.double).reshape(-1,4))
#
# rotatedPole = rot @ np.broadcast_to(unit_pole,rot.shape+(3,)) # rotate pole according to crystal orientation
# xy = rotatedPole[:,0:2]/(1.+abs(unit_pole[2])) # stereographic projection
# coords = xy if not polar else \
# np.block([np.sqrt(xy[:,0:1]*xy[:,0:1]+xy[:,1:2]*xy[:,1:2]),np.arctan2(xy[:,1:2],xy[:,0:1])])
# return {
# 'data': coords,
# 'label': 'p^{}_[{} {} {})'.format(u'rφ' if polar else 'xy',*m),
# 'meta' : {
2021-03-25 23:52:59 +05:30
# 'unit': '1',
# 'description': '{} coordinates of stereographic projection of pole (direction/plane) in crystal frame'\
# .format('Polar' if polar else 'Cartesian'),
2021-03-25 23:52:59 +05:30
# 'creator': 'add_pole'
# }
# }
# def add_pole(self,q,p,polar=False):
# """
# Add coordinates of stereographic projection of given pole in crystal frame.
#
# Parameters
# ----------
# q : str
# Label of the dataset containing the crystallographic orientation as quaternions.
# p : numpy.array of shape (3)
# Crystallographic direction or plane.
# polar : bool, optional
# Give pole in polar coordinates. Defaults to False.
#
# """
# self._add_generic_pointwise(self._add_pole,{'q':q},{'p':p,'polar':polar})
2020-02-21 12:15:05 +05:30
@staticmethod
2020-11-20 03:16:52 +05:30
def _add_rotation(F):
return {
2020-11-20 03:16:52 +05:30
'data': mechanics.rotation(F['data']).as_matrix(),
2020-06-25 01:04:51 +05:30
'label': f"R({F['label']})",
'meta': {
2021-03-25 23:52:59 +05:30
'unit': F['meta']['unit'],
'description': f"rotational part of {F['label']} ({F['meta']['description']})",
'creator': 'add_rotation'
}
2020-02-21 12:15:05 +05:30
}
2020-11-20 03:16:52 +05:30
def add_rotation(self,F):
2020-02-21 12:15:05 +05:30
"""
Add rotational part of a deformation gradient.
Parameters
----------
F : str, optional
2020-03-19 12:15:31 +05:30
Label of deformation gradient dataset.
2020-02-21 12:15:05 +05:30
"""
2020-11-20 03:16:52 +05:30
self._add_generic_pointwise(self._add_rotation,{'F':F})
2020-02-21 12:15:05 +05:30
@staticmethod
def _add_spherical(T):
return {
2020-11-19 19:08:54 +05:30
'data': tensor.spherical(T['data'],False),
2020-06-25 01:04:51 +05:30
'label': f"p_{T['label']}",
'meta': {
2021-03-25 23:52:59 +05:30
'unit': T['meta']['unit'],
'description': f"spherical component of tensor {T['label']} ({T['meta']['description']})",
'creator': 'add_spherical'
}
}
def add_spherical(self,T):
2020-02-21 12:15:05 +05:30
"""
Add the spherical (hydrostatic) part of a tensor.
Parameters
----------
T : str
2020-03-19 12:15:31 +05:30
Label of tensor dataset.
2020-02-21 12:15:05 +05:30
"""
self._add_generic_pointwise(self._add_spherical,{'T':T})
2020-02-21 12:15:05 +05:30
@staticmethod
2020-11-16 05:42:23 +05:30
def _add_strain(F,t,m):
return {
2020-11-16 05:42:23 +05:30
'data': mechanics.strain(F['data'],t,m),
2020-06-25 01:04:51 +05:30
'label': f"epsilon_{t}^{m}({F['label']})",
'meta': {
2021-03-25 23:52:59 +05:30
'unit': F['meta']['unit'],
'description': f"strain tensor of {F['label']} ({F['meta']['description']})",
'creator': 'add_strain'
}
}
2020-11-16 05:42:23 +05:30
def add_strain(self,F='F',t='V',m=0.0):
2020-02-21 12:15:05 +05:30
"""
Add strain tensor of a deformation gradient.
2020-02-21 12:15:05 +05:30
2021-04-06 21:09:44 +05:30
For details, see damask.mechanics.strain.
2020-02-21 12:15:05 +05:30
Parameters
----------
F : str, optional
2021-04-06 21:09:44 +05:30
Label of deformation gradient dataset. Defaults to 'F'.
t : {'V', 'U'}, optional
Type of the polar decomposition, 'V' for left stretch tensor and 'U' for right stretch tensor.
Defaults to 'V'.
2020-02-21 12:15:05 +05:30
m : float, optional
2021-04-06 21:09:44 +05:30
Order of the strain calculation. Defaults to 0.0.
2020-02-21 12:15:05 +05:30
"""
2020-11-16 05:42:23 +05:30
self._add_generic_pointwise(self._add_strain,{'F':F},{'t':t,'m':m})
2020-02-21 12:15:05 +05:30
@staticmethod
def _add_stretch_tensor(F,t):
return {
'data': (mechanics.stretch_left if t.upper() == 'V' else mechanics.stretch_right)(F['data']),
2020-06-25 01:04:51 +05:30
'label': f"{t}({F['label']})",
'meta': {
2021-03-25 23:52:59 +05:30
'unit': F['meta']['unit'],
'description': '{} stretch tensor of {} ({})'.format('left' if t.upper() == 'V' else 'right',
F['label'],F['meta']['description']),
'creator': 'add_stretch_tensor'
}
}
def add_stretch_tensor(self,F='F',t='V'):
2020-02-21 12:15:05 +05:30
"""
Add stretch tensor of a deformation gradient.
2020-02-21 12:15:05 +05:30
Parameters
----------
F : str, optional
2021-04-06 21:09:44 +05:30
Label of deformation gradient dataset. Defaults to 'F'.
t : {'V', 'U'}, optional
Type of the polar decomposition, 'V' for left stretch tensor and 'U' for right stretch tensor.
Defaults to 'V'.
2020-02-21 12:15:05 +05:30
"""
self._add_generic_pointwise(self._add_stretch_tensor,{'F':F},{'t':t})
def _job(self,group,func,datasets,args,lock):
2020-02-22 04:29:33 +05:30
"""Execute job for _add_generic_pointwise."""
try:
2020-02-22 03:46:25 +05:30
datasets_in = {}
lock.acquire()
with h5py.File(self.fname,'r') as f:
2020-03-12 03:05:58 +05:30
for arg,label in datasets.items():
loc = f[group+'/'+label]
datasets_in[arg]={'data' :loc[()],
'label':label,
'meta': {k:(v if h5py3 else v.decode()) for k,v in loc.attrs.items()}}
2020-02-22 03:46:25 +05:30
lock.release()
r = func(**datasets_in,**args)
return [group,r]
except Exception as err:
2020-06-25 01:04:51 +05:30
print(f'Error during calculation: {err}.')
return None
def _add_generic_pointwise(self,func,datasets,args={}):
2020-02-22 03:46:25 +05:30
"""
General function to add pointwise data.
2020-02-22 03:46:25 +05:30
Parameters
----------
func : function
Callback function that calculates a new dataset from one or
more datasets per HDF5 group.
2020-02-22 03:46:25 +05:30
datasets : dictionary
2021-04-06 21:09:44 +05:30
Details of the datasets to be used:
{arg (name to which the data is passed in func): label (in HDF5 file)}.
2020-02-22 03:46:25 +05:30
args : dictionary, optional
2020-03-19 12:15:31 +05:30
Arguments parsed to func.
2020-02-22 04:29:33 +05:30
2020-02-22 03:46:25 +05:30
"""
chunk_size = 1024**2//8
2021-01-13 17:18:20 +05:30
pool = mp.Pool(int(os.environ.get('OMP_NUM_THREADS',1)))
2020-06-26 15:15:06 +05:30
lock = mp.Manager().Lock()
2020-02-22 03:46:25 +05:30
2021-04-04 22:02:17 +05:30
groups = []
with h5py.File(self.fname,'r') as f:
for inc in self.visible['increments']:
2021-04-05 11:23:19 +05:30
for ty in ['phase','homogenization']:
for label in self.visible[ty+'s']:
2021-04-04 22:02:17 +05:30
for field in self.visible['fields']:
2021-04-05 11:23:19 +05:30
group = '/'.join([inc,ty,label,field])
2021-04-04 22:02:17 +05:30
if set(datasets.values()).issubset(f[group].keys()): groups.append(group)
if len(groups) == 0:
print('No matching dataset found, no data was added.')
return
2020-02-22 03:46:25 +05:30
default_arg = partial(self._job,func=func,datasets=datasets,args=args,lock=lock)
for result in util.show_progress(pool.imap_unordered(default_arg,groups),len(groups)):
if not result:
continue
2020-02-22 03:46:25 +05:30
lock.acquire()
with h5py.File(self.fname, 'a') as f:
try:
if self._allow_modification and result[0]+'/'+result[1]['label'] in f:
dataset = f[result[0]+'/'+result[1]['label']]
dataset[...] = result[1]['data']
2021-03-25 23:52:59 +05:30
dataset.attrs['overwritten'] = True
else:
if result[1]['data'].size >= chunk_size*2:
shape = result[1]['data'].shape
chunks = (chunk_size//np.prod(shape[1:]),)+shape[1:]
dataset = f[result[0]].create_dataset(result[1]['label'],data=result[1]['data'],
maxshape=shape, chunks=chunks,
compression='gzip', compression_opts=6,
shuffle=True,fletcher32=True)
else:
dataset = f[result[0]].create_dataset(result[1]['label'],data=result[1]['data'])
2020-05-25 22:20:31 +05:30
now = datetime.datetime.now().astimezone()
2021-03-25 23:52:59 +05:30
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()
2020-02-22 03:46:25 +05:30
for l,v in result[1]['meta'].items():
2021-03-25 23:52:59 +05:30
dataset.attrs[l.lower()]=v if h5py3 else v.encode()
creator = dataset.attrs['creator'] if h5py3 else \
dataset.attrs['creator'].decode()
2021-04-03 11:38:48 +05:30
dataset.attrs['creator'] = f'damask.Result.{creator} v{damask.version}' if h5py3 else \
f'damask.Result.{creator} v{damask.version}'.encode()
except (OSError,RuntimeError) as err:
2020-06-25 01:04:51 +05:30
print(f'Could not add dataset: {err}.')
2020-02-22 03:46:25 +05:30
lock.release()
pool.close()
pool.join()
2020-02-21 12:15:05 +05:30
def save_XDMF(self,output='*'):
"""
Write XDMF file to directly visualize data in DADF5 file.
2021-04-05 13:59:34 +05:30
Parameters
----------
2021-04-06 21:09:44 +05:30
output : (list of) str
Labels of the datasets to read.
Defaults to '*', in which case all datasets are considered.
2021-04-05 13:59:34 +05:30
"""
2021-04-05 13:59:34 +05:30
u = 'Unit' if self.version_minor < 12 else 'unit' # compatibility hack
if self.N_constituents != 1 or len(self.phases) != 1 or not self.structured:
2021-04-06 21:09:44 +05:30
raise TypeError('XDMF output requires structured grid with single phase and single constituent.')
attribute_type_map = defaultdict(lambda:'Matrix', ( ((),'Scalar'), ((3,),'Vector'), ((3,3),'Tensor')) )
def number_type_map(dtype):
if dtype in np.sctypes['int']: return 'Int'
if dtype in np.sctypes['uint']: return 'UInt'
if dtype in np.sctypes['float']: return 'Float'
2021-04-06 21:09:44 +05:30
xdmf = ET.Element('Xdmf')
xdmf.attrib={'Version': '2.0',
'xmlns:xi': 'http://www.w3.org/2001/XInclude'}
2021-04-06 21:09:44 +05:30
domain = ET.SubElement(xdmf, 'Domain')
collection = ET.SubElement(domain, 'Grid')
collection.attrib={'GridType': 'Collection',
'CollectionType': 'Temporal'}
time = ET.SubElement(collection, 'Time')
time.attrib={'TimeType': 'List'}
time_data = ET.SubElement(time, 'DataItem')
2021-04-05 13:59:34 +05:30
times = [self.times[self.increments.index(i)] for i in self.visible['increments']]
time_data.attrib={'Format': 'XML',
'NumberType': 'Float',
2021-04-05 13:59:34 +05:30
'Dimensions': f'{len(times)}'}
time_data.text = ' '.join(map(str,times))
attributes = []
data_items = []
2021-04-05 19:28:10 +05:30
with h5py.File(self.fname,'r') as f:
for inc in self.visible['increments']:
2021-04-06 21:09:44 +05:30
grid = ET.SubElement(collection,'Grid')
2021-04-05 19:28:10 +05:30
grid.attrib = {'GridType': 'Uniform',
'Name': inc}
2021-04-06 21:09:44 +05:30
topology = ET.SubElement(grid, 'Topology')
topology.attrib = {'TopologyType': '3DCoRectMesh',
'Dimensions': '{} {} {}'.format(*(self.cells+1))}
2021-04-06 21:09:44 +05:30
geometry = ET.SubElement(grid, 'Geometry')
geometry.attrib = {'GeometryType':'Origin_DxDyDz'}
2021-04-06 21:09:44 +05:30
origin = ET.SubElement(geometry, 'DataItem')
origin.attrib = {'Format': 'XML',
'NumberType': 'Float',
'Dimensions': '3'}
origin.text = "{} {} {}".format(*self.origin)
2021-04-06 21:09:44 +05:30
delta = ET.SubElement(geometry, 'DataItem')
delta.attrib = {'Format': 'XML',
'NumberType': 'Float',
'Dimensions': '3'}
2021-04-05 19:28:10 +05:30
delta.text="{} {} {}".format(*(self.size/self.cells))
attributes.append(ET.SubElement(grid, 'Attribute'))
2021-04-06 21:09:44 +05:30
attributes[-1].attrib = {'Name': 'u / m',
'Center': 'Node',
'AttributeType': 'Vector'}
data_items.append(ET.SubElement(attributes[-1], 'DataItem'))
2021-04-06 21:09:44 +05:30
data_items[-1].attrib = {'Format': 'HDF',
'Precision': '8',
'Dimensions': '{} {} {} 3'.format(*(self.cells+1))}
data_items[-1].text = f'{os.path.split(self.fname)[1]}:/{inc}/geometry/u_n'
2021-04-05 11:23:19 +05:30
for ty in ['phase','homogenization']:
for label in self.visible[ty+'s']:
for field in self.visible['fields']:
for out in _match(output,f['/'.join([inc,ty,label,field])].keys()):
2021-04-05 11:23:19 +05:30
name = '/'.join([inc,ty,label,field,out])
shape = f[name].shape[1:]
dtype = f[name].dtype
2021-03-25 23:52:59 +05:30
unit = f[name].attrs[u] if h5py3 else f[name].attrs[u].decode()
attributes.append(ET.SubElement(grid, 'Attribute'))
2021-04-06 21:09:44 +05:30
attributes[-1].attrib = {'Name': name.split('/',2)[2]+f' / {unit}',
'Center': 'Cell',
'AttributeType': attribute_type_map[shape]}
data_items.append(ET.SubElement(attributes[-1], 'DataItem'))
2021-04-06 21:09:44 +05:30
data_items[-1].attrib = {'Format': 'HDF',
'NumberType': number_type_map(dtype),
'Precision': f'{dtype.itemsize}',
'Dimensions': '{} {} {} {}'.format(*self.cells,1 if shape == () else
np.prod(shape))}
data_items[-1].text = f'{os.path.split(self.fname)[1]}:{name}'
with open(self.fname.with_suffix('.xdmf').name,'w',newline='\n') as f:
f.write(xml.dom.minidom.parseString(ET.tostring(xdmf).decode()).toprettyxml())
2021-04-05 19:11:28 +05:30
def _mappings(self):
2021-04-06 21:09:44 +05:30
grp = 'mapping' if self.version_minor < 12 else 'cell_to' # compatibility hack
name = 'Name' if self.version_minor < 12 else 'label' # compatibility hack
member = 'member' if self.version_minor < 12 else 'entry' # compatibility hack
2021-04-05 19:11:28 +05:30
with h5py.File(self.fname,'r') as f:
at_cell_ph = []
in_data_ph = []
for c in range(self.N_constituents):
at_cell_ph.append({label: np.where(f['/'.join([grp,'phase'])][:,c][name] == label.encode())[0] \
for label in self.visible['phases']})
in_data_ph.append({label: f['/'.join([grp,'phase'])][member][at_cell_ph[c][label]][:,c] \
for label in self.visible['phases']})
at_cell_ho = {label: np.where(f['/'.join([grp,'homogenization'])][:][name] == label.encode())[0] \
for label in self.visible['homogenizations']}
in_data_ho = {label: f['/'.join([grp,'homogenization'])][member][at_cell_ho[label]] \
for label in self.visible['homogenizations']}
return at_cell_ph,in_data_ph,at_cell_ho,in_data_ho
2021-04-05 13:59:34 +05:30
def save_VTK(self,output='*',mode='cell',constituents=None,fill_float=np.nan,fill_int=0,parallel=True):
2020-02-21 12:15:05 +05:30
"""
2021-04-06 21:09:44 +05:30
Export to VTK cell/point data.
2020-02-21 12:15:05 +05:30
Parameters
----------
2021-04-06 21:09:44 +05:30
output : (list of) str, optional
Labels of the datasets to place.
Defaults to '*', in which case all datasets are exported.
mode : {'cell', 'point'}
2020-03-19 12:15:31 +05:30
Export in cell format or point format.
Defaults to 'cell'.
2021-04-06 21:09:44 +05:30
constituents : (list of) int, optional
Constituents to consider.
Defaults to None, in which case all constituents are considered.
2021-04-02 15:51:27 +05:30
fill_float : float
Fill value for non-existent entries of floating point type.
2021-04-06 21:09:44 +05:30
Defaults to NaN.
2021-04-02 15:51:27 +05:30
fill_int : int
Fill value for non-existent entries of integer type.
Defaults to 0.
2021-04-05 13:59:34 +05:30
parallel : bool
Write out VTK files in parallel in a separate background process.
Defaults to True.
2020-02-21 12:15:05 +05:30
"""
if mode.lower()=='cell':
v = self.geometry0
2020-02-21 23:22:58 +05:30
elif mode.lower()=='point':
v = VTK.from_poly_data(self.coordinates0_point)
2020-03-12 03:05:58 +05:30
2021-04-05 13:59:34 +05:30
ln = 3 if self.version_minor < 12 else 10 # compatibility hack
2021-03-25 23:52:59 +05:30
N_digits = int(np.floor(np.log10(max(1,int(self.increments[-1][ln:])))))+1
2020-02-21 12:15:05 +05:30
2021-04-02 15:51:27 +05:30
constituents_ = constituents if isinstance(constituents,Iterable) else \
(range(self.N_constituents) if constituents is None else [constituents])
suffixes = [''] if self.N_constituents == 1 or isinstance(constituents,int) else \
[f'#{c}' for c in constituents_]
2021-04-05 19:11:28 +05:30
at_cell_ph,in_data_ph,at_cell_ho,in_data_ho = self._mappings()
2021-04-02 15:51:27 +05:30
with h5py.File(self.fname,'r') as f:
for inc in util.show_progress(self.visible['increments']):
u = _read(f['/'.join([inc,'geometry','u_n' if mode.lower() == 'cell' else 'u_p'])])
2021-04-02 15:51:27 +05:30
v.add(u,'u')
for ty in ['phase','homogenization']:
for field in self.visible['fields']:
outs = {}
2021-04-02 15:51:27 +05:30
for label in self.visible[ty+'s']:
if field not in f['/'.join([inc,ty,label])].keys(): continue
2021-04-02 15:51:27 +05:30
for out in _match(output,f['/'.join([inc,ty,label,field])].keys()):
data = ma.array(_read(f['/'.join([inc,ty,label,field,out])]))
2021-04-02 15:51:27 +05:30
if ty == 'phase':
if out+suffixes[0] not in outs.keys():
for c,suffix in zip(constituents_,suffixes):
outs[out+suffix] = \
2021-04-06 21:09:44 +05:30
_empty_like(data,self.N_materialpoints,fill_float,fill_int)
2021-04-02 15:51:27 +05:30
for c,suffix in zip(constituents_,suffixes):
outs[out+suffix][at_cell_ph[c][label]] = data[in_data_ph[c][label]]
if ty == 'homogenization':
if out not in outs.keys():
2021-04-06 21:09:44 +05:30
outs[out] = _empty_like(data,self.N_materialpoints,fill_float,fill_int)
2020-03-22 20:43:35 +05:30
2021-04-02 15:51:27 +05:30
outs[out][at_cell_ho[label]] = data[in_data_ho[label]]
2020-03-22 20:43:35 +05:30
2021-04-02 15:51:27 +05:30
for label,dataset in outs.items():
v.add(dataset,' / '.join(['/'.join([ty,field,label]),dataset.dtype.metadata['unit']]))
2020-03-22 20:43:35 +05:30
2021-04-05 13:59:34 +05:30
v.save(f'{self.fname.stem}_inc{inc[ln:].zfill(N_digits)}',parallel=parallel)
2021-04-06 21:09:44 +05:30
def get(self,output='*',flatten=True,prune=True):
2021-03-31 20:41:53 +05:30
"""
2021-04-06 21:09:44 +05:30
Collect data per phase/homogenization reflecting the group/folder structure in the DADF5 file.
2021-03-31 20:41:53 +05:30
Parameters
----------
2021-04-06 21:09:44 +05:30
output : (list of) str
Labels of the datasets to read.
Defaults to '*', in which case all datasets are read.
flatten : bool
Remove singular levels of the folder hierarchy.
This might be beneficial in case of single increment,
phase/homogenization, or field. Defaults to True.
prune : bool
Remove branches with no data. Defaults to True.
2021-03-31 20:41:53 +05:30
2021-04-06 21:09:44 +05:30
Returns
-------
data : dict of numpy.ndarray
Datasets structured by phase/homogenization and according to selected view.
2021-03-31 20:41:53 +05:30
"""
r = {}
2021-03-31 18:49:23 +05:30
with h5py.File(self.fname,'r') as f:
2021-03-30 23:11:36 +05:30
for inc in util.show_progress(self.visible['increments']):
2021-03-31 15:35:51 +05:30
r[inc] = {'phase':{},'homogenization':{},'geometry':{}}
2021-03-31 18:49:23 +05:30
for out in _match(output,f['/'.join([inc,'geometry'])].keys()):
r[inc]['geometry'][out] = _read(f['/'.join([inc,'geometry',out])])
2021-04-01 18:26:17 +05:30
for ty in ['phase','homogenization']:
2021-04-01 19:22:43 +05:30
for label in self.visible[ty+'s']:
r[inc][ty][label] = {}
for field in _match(self.visible['fields'],f['/'.join([inc,ty,label])].keys()):
2021-04-01 19:22:43 +05:30
r[inc][ty][label][field] = {}
for out in _match(output,f['/'.join([inc,ty,label,field])].keys()):
r[inc][ty][label][field][out] = _read(f['/'.join([inc,ty,label,field,out])])
2021-03-30 23:11:36 +05:30
if prune: r = util.dict_prune(r)
if flatten: r = util.dict_flatten(r)
2021-04-01 18:26:17 +05:30
return r
2021-04-06 21:09:44 +05:30
def place(self,output='*',flatten=True,prune=True,constituents=None,fill_float=np.nan,fill_int=0):
2021-03-31 20:41:53 +05:30
"""
2021-04-06 21:09:44 +05:30
Merge data into spatial order that is compatible with the damask.VTK geometry representation.
2021-03-31 20:41:53 +05:30
The returned data structure reflects the group/folder structure
in the DADF5 file.
2021-03-31 20:41:53 +05:30
Multi-phase data is fused into a single output.
2021-04-02 15:51:27 +05:30
`place` is equivalent to `read` if only one phase/homogenization
and one constituent is present.
2021-03-31 20:41:53 +05:30
Parameters
----------
2021-04-06 21:09:44 +05:30
output : (list of) str, optional
Labels of the datasets to place.
Defaults to '*', in which case all datasets are placed.
flatten : bool
Remove singular levels of the folder hierarchy.
2021-04-06 21:09:44 +05:30
This might be beneficial in case of single increment or field.
Defaults to True.
prune : bool
Remove branches with no data. Defaults to True.
2021-04-06 21:09:44 +05:30
constituents : (list of) int, optional
Constituents to consider.
Defaults to 'None', in which case all constituents are considered.
2021-04-01 03:32:51 +05:30
fill_float : float
Fill value for non-existent entries of floating point type.
2021-04-06 21:09:44 +05:30
Defaults to NaN.
2021-04-01 03:32:51 +05:30
fill_int : int
Fill value for non-existent entries of integer type.
2021-04-01 03:32:51 +05:30
Defaults to 0.
2021-04-01 19:22:43 +05:30
2021-04-08 18:28:22 +05:30
Returns
-------
data : dict of numpy.ma.MaskedArray
Datasets structured by spatial position and according to selected view.
2021-03-31 20:41:53 +05:30
"""
r = {}
2021-04-02 11:17:03 +05:30
constituents_ = constituents if isinstance(constituents,Iterable) else \
(range(self.N_constituents) if constituents is None else [constituents])
2021-04-01 19:22:43 +05:30
suffixes = [''] if self.N_constituents == 1 or isinstance(constituents,int) else \
2021-03-31 20:41:53 +05:30
[f'#{c}' for c in constituents_]
2021-04-05 19:11:28 +05:30
at_cell_ph,in_data_ph,at_cell_ho,in_data_ho = self._mappings()
with h5py.File(self.fname,'r') as f:
2021-03-31 11:24:23 +05:30
for inc in util.show_progress(self.visible['increments']):
2021-03-31 15:35:51 +05:30
r[inc] = {'phase':{},'homogenization':{},'geometry':{}}
2021-03-31 18:49:23 +05:30
for out in _match(output,f['/'.join([inc,'geometry'])].keys()):
2021-04-06 01:48:18 +05:30
r[inc]['geometry'][out] = ma.array(_read(f['/'.join([inc,'geometry',out])]),fill_value = fill_float)
2021-03-31 11:24:23 +05:30
2021-04-01 19:22:43 +05:30
for ty in ['phase','homogenization']:
for label in self.visible[ty+'s']:
for field in _match(self.visible['fields'],f['/'.join([inc,ty,label])].keys()):
2021-04-01 19:22:43 +05:30
if field not in r[inc][ty].keys():
r[inc][ty][field] = {}
for out in _match(output,f['/'.join([inc,ty,label,field])].keys()):
data = ma.array(_read(f['/'.join([inc,ty,label,field,out])]))
2021-04-01 19:22:43 +05:30
if ty == 'phase':
if out+suffixes[0] not in r[inc][ty][field].keys():
2021-04-02 11:17:03 +05:30
for c,suffix in zip(constituents_,suffixes):
2021-04-01 19:22:43 +05:30
r[inc][ty][field][out+suffix] = \
2021-04-06 21:09:44 +05:30
_empty_like(data,self.N_materialpoints,fill_float,fill_int)
2021-04-01 19:22:43 +05:30
for c,suffix in zip(constituents_,suffixes):
2021-04-01 19:22:43 +05:30
r[inc][ty][field][out+suffix][at_cell_ph[c][label]] = data[in_data_ph[c][label]]
if ty == 'homogenization':
if out not in r[inc][ty][field].keys():
r[inc][ty][field][out] = \
2021-04-06 21:09:44 +05:30
_empty_like(data,self.N_materialpoints,fill_float,fill_int)
2021-04-01 19:22:43 +05:30
r[inc][ty][field][out][at_cell_ho[label]] = data[in_data_ho[label]]
2021-03-31 18:49:23 +05:30
if prune: r = util.dict_prune(r)
if flatten: r = util.dict_flatten(r)
2021-04-01 19:22:43 +05:30
return r