DAMASK_EICMD/python/damask/_result.py

1267 lines
49 KiB
Python
Raw Normal View History

2020-06-26 15:15:06 +05:30
import multiprocessing as mp
import re
import inspect
import glob
import os
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
import h5py
2019-04-17 23:27:16 +05:30
import numpy as np
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-03-13 00:22:33 +05:30
from . import Table
2020-02-21 11:51:45 +05:30
from . import Rotation
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 util
class Result:
2019-09-20 01:02:15 +05:30
"""
2020-02-21 12:15:05 +05:30
Read and write to 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.
2020-02-21 12:15:05 +05:30
"""
def __init__(self,fname):
"""
Open an existing DADF5 file.
2019-09-16 08:49:14 +05:30
2020-02-21 12:15:05 +05:30
Parameters
----------
fname : str
2020-06-01 15:03:22 +05:30
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:
2020-03-22 20:43:35 +05:30
try:
self.version_major = f.attrs['DADF5_version_major']
self.version_minor = f.attrs['DADF5_version_minor']
except KeyError:
self.version_major = f.attrs['DADF5-major']
self.version_minor = f.attrs['DADF5-minor']
if self.version_major != 0 or not 2 <= self.version_minor <= 6:
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
self.structured = 'grid' in f['geometry'].attrs.keys()
if self.structured:
self.grid = f['geometry'].attrs['grid']
self.size = f['geometry'].attrs['size']
self.origin = f['geometry'].attrs['origin'] if self.version_major == 0 and self.version_minor >= 5 else \
np.zeros(3)
r=re.compile('inc[0-9]+')
increments_unsorted = {int(i[3:]):i for i in f.keys() if r.match(i)}
self.increments = [increments_unsorted[i] for i in sorted(increments_unsorted)]
self.times = [round(f[i].attrs['time/s'],12) for i in self.increments]
self.Nmaterialpoints, self.Nconstituents = np.shape(f['mapping/cellResults/constituent'])
self.materialpoints = [m.decode() for m in np.unique(f['mapping/cellResults/materialpoint']['Name'])]
self.constituents = [c.decode() for c in np.unique(f['mapping/cellResults/constituent'] ['Name'])]
# faster, but does not work with (deprecated) DADF5_postResults
#self.materialpoints = [m for m in f['inc0/materialpoint']]
#self.constituents = [c for c in f['inc0/constituent']]
2020-03-22 20:43:35 +05:30
2020-03-22 21:33:28 +05:30
self.con_physics = []
2020-03-22 20:43:35 +05:30
for c in self.constituents:
self.con_physics += f['/'.join([self.increments[0],'constituent',c])].keys()
2020-04-21 14:47:15 +05:30
self.con_physics = list(set(self.con_physics)) # make unique
2020-03-22 20:43:35 +05:30
self.mat_physics = []
for m in self.materialpoints:
self.mat_physics += f['/'.join([self.increments[0],'materialpoint',m])].keys()
2020-04-21 14:47:15 +05:30
self.mat_physics = list(set(self.mat_physics)) # make unique
2020-03-22 21:33:28 +05:30
self.selection = {'increments': self.increments,
'constituents': self.constituents,'materialpoints': self.materialpoints,
'con_physics': self.con_physics, 'mat_physics': self.mat_physics
}
2020-06-03 14:13:07 +05:30
self.fname = Path(fname).absolute()
self._allow_modification = False
def __repr__(self):
"""Show selected data."""
2020-03-19 16:00:36 +05:30
all_selected_increments = self.selection['increments']
2020-04-21 14:47:15 +05:30
2020-03-19 16:00:36 +05:30
self.pick('increments',all_selected_increments[0:1])
first = self.list_data()
2020-04-21 14:47:15 +05:30
2020-03-19 16:00:36 +05:30
self.pick('increments',all_selected_increments[-1:])
2020-04-21 14:47:15 +05:30
last = '' if len(all_selected_increments) < 2 else self.list_data()
2020-03-19 16:00:36 +05:30
self.pick('increments',all_selected_increments)
2020-04-21 14:47:15 +05:30
in_between = '' if len(all_selected_increments) < 3 else \
2020-06-25 01:04:51 +05:30
''.join([f'\n{inc}\n ...\n' for inc in all_selected_increments[1:-2]])
2020-04-21 14:47:15 +05:30
return util.srepr(first + in_between + last)
2020-03-03 04:17:29 +05:30
def _manage_selection(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-03-19 12:15:31 +05:30
select from 'set', 'add', and 'del'
2020-03-03 04:17:29 +05:30
what : str
2020-03-19 12:15:31 +05:30
attribute to change (must be from self.selection)
datasets : list of str or bool
2020-03-19 12:15:31 +05:30
name of datasets as list, 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:
2020-03-03 04:17:29 +05:30
datasets = ['*']
2020-02-21 12:15:05 +05:30
elif datasets is False:
2020-03-03 04:17:29 +05:30
datasets = []
2020-03-03 18:37:02 +05:30
choice = datasets if hasattr(datasets,'__iter__') and not isinstance(datasets,str) else \
[datasets]
if what == 'increments':
2020-03-03 18:54:27 +05:30
choice = [c if isinstance(c,str) and c.startswith('inc') else
2020-06-25 01:04:51 +05:30
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])
2020-02-21 12:15:05 +05:30
valid = [e for e_ in [glob.fnmatch.filter(getattr(self,what),s) for s in choice] for e in e_]
existing = set(self.selection[what])
2020-02-21 12:15:05 +05:30
if action == 'set':
2020-03-03 04:17:29 +05:30
self.selection[what] = valid
2020-02-21 12:15:05 +05:30
elif action == 'add':
2020-03-22 21:33:28 +05:30
add = existing.union(valid)
add_sorted = sorted(add, key=lambda x: int("".join([i for i in x if i.isdigit()])))
2020-03-03 04:17:29 +05:30
self.selection[what] = add_sorted
2020-02-21 12:15:05 +05:30
elif action == 'del':
2020-03-22 21:33:28 +05:30
diff = existing.difference(valid)
diff_sorted = sorted(diff, key=lambda x: int("".join([i for i in x if i.isdigit()])))
2020-03-03 04:17:29 +05:30
self.selection[what] = diff_sorted
2020-02-21 12:15:05 +05:30
def allow_modification(self):
print(util.bcolors().WARNING+util.bcolors().BOLD+
'Warning: Modification of existing datasets allowed!'+
util.bcolors().ENDC)
self._allow_modification = True
2020-06-01 15:03:22 +05:30
def disallow_modification(self):
self._allow_modification = False
def incs_in_range(self,start,end):
selected = []
for i,inc in enumerate([int(i[3:]) for i in self.increments]):
s,e = map(lambda x: int(x[3:] 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):
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
2020-03-10 03:58:25 +05:30
def iterate(self,what):
2020-02-21 12:15:05 +05:30
"""
2020-03-03 04:17:29 +05:30
Iterate over selection items by setting each one selected.
2020-02-21 12:15:05 +05:30
Parameters
----------
what : str
2020-03-19 12:15:31 +05:30
attribute to change (must be from self.selection)
2020-02-21 12:15:05 +05:30
"""
datasets = self.selection[what]
2020-03-19 16:00:36 +05:30
last_selection = datasets.copy()
2020-02-21 12:15:05 +05:30
for dataset in datasets:
2020-03-19 16:00:36 +05:30
if last_selection != self.selection[what]:
2020-03-03 11:19:46 +05:30
self._manage_selection('set',what,datasets)
raise Exception
self._manage_selection('set',what,dataset)
2020-03-19 16:00:36 +05:30
last_selection = self.selection[what]
2020-03-03 11:19:46 +05:30
yield dataset
2020-03-03 04:17:29 +05:30
self._manage_selection('set',what,datasets)
2020-03-03 04:17:29 +05:30
def pick(self,what,datasets):
2020-02-21 12:15:05 +05:30
"""
2020-03-03 04:17:29 +05:30
Set selection.
2020-02-21 12:15:05 +05:30
Parameters
----------
2020-03-03 04:17:29 +05:30
what : str
2020-03-19 12:15:31 +05:30
attribute to change (must be from self.selection)
datasets : list of str or bool
2020-03-19 12:15:31 +05:30
name of datasets as list, supports ? and * wildcards.
True is equivalent to [*], False is equivalent to []
2020-02-21 12:15:05 +05:30
"""
2020-03-03 04:17:29 +05:30
self._manage_selection('set',what,datasets)
2020-03-03 04:17:29 +05:30
def pick_more(self,what,datasets):
2020-02-21 12:15:05 +05:30
"""
2020-03-03 04:17:29 +05:30
Add to selection.
2020-02-21 12:15:05 +05:30
Parameters
----------
2020-03-03 04:17:29 +05:30
what : str
2020-03-19 12:15:31 +05:30
attribute to change (must be from self.selection)
datasets : list of str or bool
2020-03-19 12:15:31 +05:30
name of datasets as list, supports ? and * wildcards.
True is equivalent to [*], False is equivalent to []
2020-02-21 12:15:05 +05:30
"""
2020-03-03 04:17:29 +05:30
self._manage_selection('add',what,datasets)
2020-03-03 04:17:29 +05:30
def pick_less(self,what,datasets):
2020-02-21 12:15:05 +05:30
"""
2020-03-03 04:17:29 +05:30
Delete from selection.
2020-02-21 12:15:05 +05:30
Parameters
----------
2020-03-03 04:17:29 +05:30
what : str
2020-03-19 12:15:31 +05:30
attribute to change (must be from self.selection)
datasets : list of str or bool
2020-03-19 12:15:31 +05:30
name of datasets as list, 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
"""
2020-03-03 04:17:29 +05:30
self._manage_selection('del',what,datasets)
def rename(self,name_old,name_new):
"""
Rename datasets.
2020-06-25 01:04:51 +05:30
2020-06-02 01:43:01 +05:30
Parameters
----------
name_old : str
name of the datasets to be renamed
name_new : str
new name of the datasets
"""
if self._allow_modification:
with h5py.File(self.fname,'a') as f:
for path_old in self.get_dataset_location(name_old):
path_new = os.path.join(os.path.dirname(path_old),name_new)
f[path_new] = f[path_old]
f[path_new].attrs['Renamed'] = 'Original name: {}'.encode()
del f[path_old]
else:
2020-06-01 15:03:22 +05:30
raise PermissionError('Rename operation not permitted')
2020-03-10 03:58:25 +05:30
# def datamerger(regular expression to filter groups into one copy)
def place(self,datasets,component=0,tagged=False,split=True):
"""
Distribute datasets onto geometry and return Table or (split) dictionary of Tables.
2020-03-11 11:20:13 +05:30
2020-03-10 03:58:25 +05:30
Must not mix nodal end cell data.
2020-03-11 11:20:13 +05:30
2020-03-10 03:58:25 +05:30
Only data within
- inc?????/constituent/*_*/*
- inc?????/materialpoint/*_*/*
- inc?????/geometry/*
are considered.
2020-03-11 11:20:13 +05:30
2020-03-10 03:58:25 +05:30
Parameters
----------
datasets : iterable or str
component : int
2020-03-19 12:15:31 +05:30
homogenization component to consider for constituent data
tagged : bool
2020-03-19 12:15:31 +05:30
tag Table.column name with '#component'
defaults to False
split : bool
2020-03-19 12:15:31 +05:30
split Table by increment and return dictionary of Tables
defaults to True
2020-03-10 03:58:25 +05:30
"""
sets = datasets if hasattr(datasets,'__iter__') and not isinstance(datasets,str) \
else [datasets]
tag = f'#{component}' if tagged else ''
tbl = {} if split else None
inGeom = {}
inData = {}
with h5py.File(self.fname,'r') as f:
for dataset in sets:
for group in self.groups_with_datasets(dataset):
path = os.path.join(group,dataset)
inc,prop,name,cat,item = (path.split('/') + ['']*5)[:5]
key = '/'.join([prop,name+tag])
if key not in inGeom:
if prop == 'geometry':
inGeom[key] = inData[key] = np.arange(self.Nmaterialpoints)
elif prop == 'constituent':
inGeom[key] = np.where(f['mapping/cellResults/constituent'][:,component]['Name'] == str.encode(name))[0]
inData[key] = f['mapping/cellResults/constituent'][inGeom[key],component]['Position']
else:
inGeom[key] = np.where(f['mapping/cellResults/materialpoint']['Name'] == str.encode(name))[0]
inData[key] = f['mapping/cellResults/materialpoint'][inGeom[key].tolist()]['Position']
shape = np.shape(f[path])
data = np.full((self.Nmaterialpoints,) + (shape[1:] if len(shape)>1 else (1,)),
2020-03-22 21:33:28 +05:30
np.nan,
dtype=np.dtype(f[path]))
2020-03-10 03:58:25 +05:30
data[inGeom[key]] = (f[path] if len(shape)>1 else np.expand_dims(f[path],1))[inData[key]]
path = (os.path.join(*([prop,name]+([cat] if cat else [])+([item] if item else []))) if split else path)+tag
if split:
try:
tbl[inc].add(path,data)
except KeyError:
2020-03-13 00:22:33 +05:30
tbl[inc] = Table(data.reshape(self.Nmaterialpoints,-1),{path:data.shape[1:]})
2020-03-10 03:58:25 +05:30
else:
try:
tbl.add(path,data)
except AttributeError:
2020-03-13 00:22:33 +05:30
tbl = Table(data.reshape(self.Nmaterialpoints,-1),{path:data.shape[1:]})
2020-03-10 03:58:25 +05:30
return tbl
2020-02-21 12:15:05 +05:30
def groups_with_datasets(self,datasets):
"""
Return groups that contain all requested datasets.
Only groups within
2020-05-13 14:39:37 +05:30
- inc*/constituent/*/*
- inc*/materialpoint/*/*
- inc*/geometry/*
are considered as they contain user-relevant data.
2020-02-21 12:15:05 +05:30
Single strings will be treated as list with one entry.
2020-02-21 12:15:05 +05:30
Wild card matching is allowed, but the number of arguments need to fit.
2020-02-21 12:15:05 +05:30
Parameters
----------
datasets : iterable or str or bool
2020-02-21 12:15:05 +05:30
Examples
--------
2020-03-19 12:15:31 +05:30
datasets = False matches no group
datasets = True matches all groups
datasets = ['F','P'] matches a group with ['F','P','sigma']
datasets = ['*','P'] matches a group with ['F','P']
datasets = ['*'] does not match a group with ['F','P','sigma']
datasets = ['*','*'] does not match a group with ['F','P','sigma']
datasets = ['*','*','*'] matches a group with ['F','P','sigma']
2020-02-21 12:15:05 +05:30
"""
if datasets is False: return []
2020-03-03 18:37:02 +05:30
sets = datasets if isinstance(datasets,bool) or (hasattr(datasets,'__iter__') and not isinstance(datasets,str)) else \
[datasets]
2019-09-12 06:27:24 +05:30
2020-02-21 12:15:05 +05:30
groups = []
2020-02-21 12:15:05 +05:30
with h5py.File(self.fname,'r') as f:
2020-03-10 03:58:25 +05:30
for i in self.iterate('increments'):
2020-03-03 18:54:27 +05:30
for o,p in zip(['constituents','materialpoints'],['con_physics','mat_physics']):
2020-03-22 20:43:35 +05:30
for oo in self.iterate(o):
for pp in self.iterate(p):
group = '/'.join([i,o[:-1],oo,pp]) # o[:-1]: plural/singular issue
if sets is True:
groups.append(group)
else:
match = [e for e_ in [glob.fnmatch.filter(f[group].keys(),s) for s in sets] for e in e_]
2020-03-22 21:33:28 +05:30
if len(set(match)) == len(sets): groups.append(group)
2020-02-21 12:15:05 +05:30
return groups
def list_data(self):
"""Return information on all active datasets in the file."""
message = ''
with h5py.File(self.fname,'r') as f:
2020-03-10 03:58:25 +05:30
for i in self.iterate('increments'):
2020-06-25 01:04:51 +05:30
message += f'\n{i} ({self.times[self.increments.index(i)]}s)\n'
2020-03-03 18:54:27 +05:30
for o,p in zip(['constituents','materialpoints'],['con_physics','mat_physics']):
2020-06-26 15:15:06 +05:30
message += f' {o[:-1]}\n'
2020-03-22 20:43:35 +05:30
for oo in self.iterate(o):
2020-06-26 15:15:06 +05:30
message += f' {oo}\n'
2020-03-22 20:43:35 +05:30
for pp in self.iterate(p):
2020-06-26 15:15:06 +05:30
message += f' {pp}\n'
2020-03-22 20:43:35 +05:30
group = '/'.join([i,o[:-1],oo,pp]) # o[:-1]: plural/singular issue
for d in f[group].keys():
try:
dataset = f['/'.join([group,d])]
2020-06-25 01:04:51 +05:30
unit = f" / {dataset.attrs['Unit'].decode()}" if 'Unit' in dataset.attrs else ''
description = dataset.attrs['Description'].decode()
2020-06-26 15:15:06 +05:30
message += f' {d}{unit}: {description}\n'
2020-03-22 20:43:35 +05:30
except KeyError:
pass
2020-02-21 12:15:05 +05:30
return message
def get_dataset_location(self,label):
"""Return the location of all active datasets with given label."""
path = []
with h5py.File(self.fname,'r') as f:
2020-03-10 03:58:25 +05:30
for i in self.iterate('increments'):
2020-03-03 18:54:27 +05:30
k = '/'.join([i,'geometry',label])
try:
2020-02-21 12:15:05 +05:30
f[k]
path.append(k)
2020-03-22 20:43:35 +05:30
except KeyError:
2020-02-21 12:15:05 +05:30
pass
2020-03-03 18:54:27 +05:30
for o,p in zip(['constituents','materialpoints'],['con_physics','mat_physics']):
2020-03-10 03:58:25 +05:30
for oo in self.iterate(o):
for pp in self.iterate(p):
2020-03-03 18:54:27 +05:30
k = '/'.join([i,o[:-1],oo,pp,label])
try:
f[k]
path.append(k)
2020-03-22 20:43:35 +05:30
except KeyError:
2020-03-03 18:54:27 +05:30
pass
2020-02-21 12:15:05 +05:30
return path
def get_constituent_ID(self,c=0):
"""Pointwise constituent ID."""
with h5py.File(self.fname,'r') as f:
2020-03-03 18:54:27 +05:30
names = f['/mapping/cellResults/constituent']['Name'][:,c].astype('str')
2020-02-21 12:15:05 +05:30
return np.array([int(n.split('_')[0]) for n in names.tolist()],dtype=np.int32)
2020-02-21 12:15:05 +05:30
def get_crystal_structure(self): # ToDo: extension to multi constituents/phase
"""Info about the crystal structure."""
with h5py.File(self.fname,'r') as f:
2020-03-03 18:54:27 +05:30
return f[self.get_dataset_location('orientation')[0]].attrs['Lattice'].astype('str') # np.bytes_ to string
2020-02-21 12:15:05 +05:30
def read_dataset(self,path,c=0,plain=False):
"""
Dataset for all points/cells.
2020-02-21 12:15:05 +05:30
If more than one path is given, the dataset is composed of the individual contributions.
"""
with h5py.File(self.fname,'r') as f:
2020-03-22 20:43:35 +05:30
shape = (self.Nmaterialpoints,) + np.shape(f[path[0]])[1:]
if len(shape) == 1: shape = shape +(1,)
dataset = np.full(shape,np.nan,dtype=np.dtype(f[path[0]]))
for pa in path:
label = pa.split('/')[2]
2020-03-22 21:33:28 +05:30
if pa.split('/')[1] == 'geometry':
2020-03-22 20:43:35 +05:30
dataset = np.array(f[pa])
continue
p = np.where(f['mapping/cellResults/constituent'][:,c]['Name'] == str.encode(label))[0]
if len(p)>0:
u = (f['mapping/cellResults/constituent']['Position'][p,c])
a = np.array(f[pa])
if len(a.shape) == 1:
a=a.reshape([a.shape[0],1])
dataset[p,:] = a[u,:]
p = np.where(f['mapping/cellResults/materialpoint']['Name'] == str.encode(label))[0]
if len(p)>0:
u = (f['mapping/cellResults/materialpoint']['Position'][p.tolist()])
a = np.array(f[pa])
if len(a.shape) == 1:
a=a.reshape([a.shape[0],1])
dataset[p,:] = a[u,:]
2020-02-21 12:15:05 +05:30
if plain and dataset.dtype.names is not None:
2020-03-03 18:54:27 +05:30
return dataset.view(('float64',len(dataset.dtype.names)))
2020-02-21 12:15:05 +05:30
else:
2020-03-03 18:54:27 +05:30
return dataset
2020-02-21 12:15:05 +05:30
def cell_coordinates(self):
"""Return initial coordinates of the cell centers."""
if self.structured:
return grid_filters.cell_coord0(self.grid,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-04-22 11:10:02 +05:30
def node_coordinates(self):
"""Return initial coordinates of the cell centers."""
if self.structured:
return grid_filters.node_coord0(self.grid,self.size,self.origin).reshape(-1,3,order='F')
else:
with h5py.File(self.fname,'r') as f:
return f['geometry/x_n'][()]
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': {
'Unit': x['meta']['Unit'],
2020-06-25 01:04:51 +05:30
'Description': f"Absolute value of {x['label']} ({x['meta']['Description']})",
2020-05-25 23:29:29 +05:30
'Creator': inspect.stack()[0][3][1:]
}
}
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': {
'Unit': kwargs['unit'],
2020-06-25 01:04:51 +05:30
'Description': f"{kwargs['description']} (formula: {kwargs['formula']})",
2020-05-25 23:29:29 +05:30
'Creator': inspect.stack()[0][3][1:]
}
}
def add_calculation(self,label,formula,unit='n/a',description=None,vectorized=True):
2020-02-21 12:15:05 +05:30
"""
Add result of a general formula.
Parameters
----------
label : str
Label of resulting dataset.
formula : str
2020-03-19 12:15:31 +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
vectorized : bool, optional
2020-03-19 12:15:31 +05:30
Indicate whether the formula can be used in vectorized form. Defaults to True.
2020-02-21 12:15:05 +05:30
"""
if not vectorized:
2020-03-22 20:43:35 +05:30
raise NotImplementedError
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
def _add_Cauchy(P,F):
return {
'data': mechanics.Cauchy(P['data'],F['data']),
'label': 'sigma',
'meta': {
'Unit': P['meta']['Unit'],
2020-06-25 01:04:51 +05:30
'Description': "Cauchy stress calculated "
f"from {P['label']} ({P['meta']['Description']})"
f" and {F['label']} ({F['meta']['Description']})",
2020-05-25 23:29:29 +05:30
'Creator': inspect.stack()[0][3][1:]
}
}
def add_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
2020-03-19 12:15:31 +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
2020-03-19 12:15:31 +05:30
Label of the dataset containing the deformation gradient. Defaults to F.
2020-02-21 12:15:05 +05:30
"""
self._add_generic_pointwise(self._add_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': {
'Unit': T['meta']['Unit'],
2020-06-25 01:04:51 +05:30
'Description': f"Determinant of tensor {T['label']} ({T['meta']['Description']})",
2020-05-25 23:29:29 +05:30
'Creator': inspect.stack()[0][3][1:]
}
}
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 {
'data': mechanics.deviatoric_part(T['data']),
2020-06-25 01:04:51 +05:30
'label': f"s_{T['label']}",
'meta': {
'Unit': T['meta']['Unit'],
2020-06-25 01:04:51 +05:30
'Description': f"Deviator of tensor {T['label']} ({T['meta']['Description']})",
2020-05-25 23:29:29 +05:30
'Creator': inspect.stack()[0][3][1:]
}
}
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':
label,p = 'Maximum',2
elif eigenvalue == 'mid':
label,p = 'Intermediate',1
elif eigenvalue == 'min':
label,p = 'Minimum',0
return {
'data': mechanics.eigenvalues(T_sym['data'])[:,p],
2020-06-25 01:04:51 +05:30
'label': f"lambda_{eigenvalue}({T_sym['label']})",
'meta' : {
2020-03-03 03:44:59 +05:30
'Unit': T_sym['meta']['Unit'],
2020-06-25 01:04:51 +05:30
'Description': f"{label} eigenvalue of {T_sym['label']} ({T_sym['meta']['Description']})",
2020-05-25 23:29:29 +05:30
'Creator': inspect.stack()[0][3][1:]
}
}
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
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': mechanics.eigenvectors(T_sym['data'])[:,p],
2020-06-25 01:04:51 +05:30
'label': f"v_{eigenvalue}({T_sym['label']})",
'meta' : {
'Unit': '1',
2020-06-25 01:04:51 +05:30
'Description': f"Eigenvector corresponding to {label} eigenvalue"
f" of {T_sym['label']} ({T_sym['meta']['Description']})",
2020-05-25 23:29:29 +05:30
'Creator': inspect.stack()[0][3][1:]
}
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
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-07-01 01:13:57 +05:30
def _add_IPF_color(q,l):
m = util.scale_to_coprime(np.array(l))
2020-02-21 12:15:05 +05:30
o = Orientation(Rotation(rfn.structured_to_unstructured(q['data'])),
lattice = q['meta']['Lattice'])
return {
'data': np.uint8(o.IPF_color(l)*255),
'label': 'IPFcolor_[{} {} {}]'.format(*m),
'meta' : {
'Unit': '8-bit RGB',
'Lattice': q['meta']['Lattice'],
2020-03-19 16:00:36 +05:30
'Description': 'Inverse Pole Figure (IPF) colors along sample direction [{} {} {}]'.format(*m),
2020-05-25 23:29:29 +05:30
'Creator': inspect.stack()[0][3][1:]
}
}
2020-07-01 01:13:57 +05:30
def add_IPF_color(self,q,l):
2020-02-21 12:15:05 +05:30
"""
Add RGB color tuple of inverse pole figure (IPF) color.
Parameters
----------
q : str
2020-03-19 12:15:31 +05:30
Label of the dataset containing the crystallographic orientation as quaternions.
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-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': {
2020-03-03 03:44:59 +05:30
'Unit': T_sym['meta']['Unit'],
2020-06-25 01:04:51 +05:30
'Description': f"Maximum shear component of {T_sym['label']} ({T_sym['meta']['Description']})",
2020-05-25 23:29:29 +05:30
'Creator': inspect.stack()[0][3][1:]
}
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-03-03 03:44:59 +05:30
def _add_Mises(T_sym):
t = 'strain' if T_sym['meta']['Unit'] == '1' else \
'stress'
2020-02-21 12:15:05 +05:30
return {
2020-06-25 01:04:51 +05:30
'data': (mechanics.Mises_strain if t=='strain' else mechanics.Mises_stress)(T_sym['data']),
'label': f"{T_sym['label']}_vM",
'meta': {
2020-03-03 03:44:59 +05:30
'Unit': T_sym['meta']['Unit'],
2020-06-25 01:04:51 +05:30
'Description': f"Mises equivalent {t} of {T_sym['label']} ({T_sym['meta']['Description']})",
2020-05-25 23:29:29 +05:30
'Creator': inspect.stack()[0][3][1:]
}
}
2020-03-03 03:44:59 +05:30
def add_Mises(self,T_sym):
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.
2020-02-21 12:15:05 +05:30
"""
2020-03-03 03:44:59 +05:30
self._add_generic_pointwise(self._add_Mises,{'T_sym':T_sym})
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': {
'Unit': x['meta']['Unit'],
2020-06-25 01:04:51 +05:30
'Description': f"{o}-norm of {t} {x['label']} ({x['meta']['Description']})",
2020-05-25 23:29:29 +05:30
'Creator': inspect.stack()[0][3][1:]
}
}
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.
2020-02-21 12:15:05 +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
def _add_PK2(P,F):
return {
'data': mechanics.PK2(P['data'],F['data']),
'label': 'S',
'meta': {
'Unit': P['meta']['Unit'],
2020-06-25 01:04:51 +05:30
'Description': "2. Piola-Kirchhoff stress calculated "
f"from {P['label']} ({P['meta']['Description']})"
f" and {F['label']} ({F['meta']['Description']})",
2020-05-25 23:29:29 +05:30
'Creator': inspect.stack()[0][3][1:]
}
}
def add_PK2(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
2020-06-25 01:04:51 +05:30
Label of first Piola-Kirchhoff stress dataset. Defaults to P.
2020-02-21 12:15:05 +05:30
F : str, optional
2020-03-19 12:15:31 +05:30
Label of deformation gradient dataset. Defaults to F.
2020-02-21 12:15:05 +05:30
"""
self._add_generic_pointwise(self._add_PK2,{'P':P,'F':F})
2020-02-21 12:15:05 +05:30
@staticmethod
def _add_pole(q,p,polar):
pole = np.array(p)
unit_pole = pole/np.linalg.norm(pole)
m = util.scale_to_coprime(pole)
2020-05-23 12:08:57 +05:30
rot = Rotation(q['data'].view(np.double).reshape(-1,4))
2020-02-21 12:15:05 +05:30
2020-05-23 12:08:57 +05:30
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'' if polar else 'xy',*m),
'meta' : {
2020-05-25 23:29:29 +05:30
'Unit': '1',
'Description': '{} coordinates of stereographic projection of pole (direction/plane) in crystal frame'\
.format('Polar' if polar else 'Cartesian'),
2020-05-25 23:29:29 +05:30
'Creator': inspect.stack()[0][3][1:]
}
}
2020-02-21 12:15:05 +05:30
def add_pole(self,q,p,polar=False):
"""
Add coordinates of stereographic projection of given pole in crystal frame.
Parameters
----------
q : str
2020-03-19 12:15:31 +05:30
Label of the dataset containing the crystallographic orientation as quaternions.
2020-02-21 12:15:05 +05:30
p : numpy.array of shape (3)
2020-03-19 12:15:31 +05:30
Crystallographic direction or plane.
2020-02-21 12:15:05 +05:30
polar : bool, optional
2020-03-19 12:15:31 +05:30
Give pole in polar coordinates. Defaults to False.
2020-02-21 12:15:05 +05:30
"""
self._add_generic_pointwise(self._add_pole,{'q':q},{'p':p,'polar':polar})
2020-02-21 12:15:05 +05:30
@staticmethod
def _add_rotational_part(F):
return {
'data': mechanics.rotational_part(F['data']),
2020-06-25 01:04:51 +05:30
'label': f"R({F['label']})",
'meta': {
'Unit': F['meta']['Unit'],
2020-06-25 01:04:51 +05:30
'Description': f"Rotational part of {F['label']} ({F['meta']['Description']})",
2020-05-25 23:29:29 +05:30
'Creator': inspect.stack()[0][3][1:]
}
2020-02-21 12:15:05 +05:30
}
def add_rotational_part(self,F):
"""
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
"""
self._add_generic_pointwise(self._add_rotational_part,{'F':F})
2020-02-21 12:15:05 +05:30
@staticmethod
def _add_spherical(T):
return {
'data': mechanics.spherical_part(T['data']),
2020-06-25 01:04:51 +05:30
'label': f"p_{T['label']}",
'meta': {
'Unit': T['meta']['Unit'],
2020-06-25 01:04:51 +05:30
'Description': f"Spherical component of tensor {T['label']} ({T['meta']['Description']})",
2020-05-25 23:29:29 +05:30
'Creator': inspect.stack()[0][3][1:]
}
}
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
def _add_strain_tensor(F,t,m):
return {
'data': mechanics.strain_tensor(F['data'],t,m),
2020-06-25 01:04:51 +05:30
'label': f"epsilon_{t}^{m}({F['label']})",
'meta': {
'Unit': F['meta']['Unit'],
2020-06-25 01:04:51 +05:30
'Description': f"Strain tensor of {F['label']} ({F['meta']['Description']})",
2020-05-25 23:29:29 +05:30
'Creator': inspect.stack()[0][3][1:]
}
}
def add_strain_tensor(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
For details refer to damask.mechanics.strain_tensor
Parameters
----------
F : str, optional
2020-03-19 12:15:31 +05:30
Label of deformation gradient dataset. Defaults to F.
2020-02-21 12:15:05 +05:30
t : {V, U}, optional
2020-03-19 12:15:31 +05:30
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
2020-03-19 12:15:31 +05:30
Order of the strain calculation. Defaults to 0.0.
2020-02-21 12:15:05 +05:30
"""
self._add_generic_pointwise(self._add_strain_tensor,{'F':F},{'t':t,'m':m})
2020-02-21 12:15:05 +05:30
@staticmethod
def _add_stretch_tensor(F,t):
return {
2020-06-25 01:04:51 +05:30
'data': (mechanics.left_stretch if t.upper() == 'V' else mechanics.right_stretch)(F['data']),
'label': f"{t}({F['label']})",
'meta': {
'Unit': F['meta']['Unit'],
2020-06-25 01:04:51 +05:30
'Description': '{} stretch tensor of {} ({})'.format('Left' if t.upper() == 'V' else 'Right',
F['label'],F['meta']['Description']),
2020-05-25 23:29:29 +05:30
'Creator': inspect.stack()[0][3][1:]
}
}
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
2020-03-19 12:15:31 +05:30
Label of deformation gradient dataset. Defaults to F.
2020-02-21 12:15:05 +05:30
t : {V, U}, optional
2020-03-19 12:15:31 +05:30
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.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
Details of the datasets to be used: label (in HDF5 file) and
arg (argument to which the data is parsed in func).
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
"""
num_threads = damask.environment.options['DAMASK_NUM_THREADS']
2020-06-26 15:15:06 +05:30
pool = mp.Pool(int(num_threads) if num_threads is not None else None)
lock = mp.Manager().Lock()
2020-02-22 03:46:25 +05:30
groups = self.groups_with_datasets(datasets.values())
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']
dataset.attrs['Overwritten'] = 'Yes'.encode()
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()
dataset.attrs['Created'] = 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():
dataset.attrs[l]=v.encode()
2020-06-28 15:10:19 +05:30
creator = f"damask.Result.{dataset.attrs['Creator'].decode()} v{damask.version}"
dataset.attrs['Creator'] = creator.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
2020-05-22 22:34:02 +05:30
def write_XDMF(self):
"""
Write XDMF file to directly visualize data in DADF5 file.
This works only for scalar, 3-vector and 3x3-tensor data.
Selection is not taken into account.
"""
if len(self.constituents) != 1 or not self.structured:
2020-05-23 19:24:28 +05:30
raise NotImplementedError('XDMF only available for grid results with 1 constituent.')
xdmf=ET.Element('Xdmf')
xdmf.attrib={'Version': '2.0',
'xmlns:xi': 'http://www.w3.org/2001/XInclude'}
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')
time_data.attrib={'Format': 'XML',
'NumberType': 'Float',
2020-06-25 01:04:51 +05:30
'Dimensions': f'{len(self.times)}'}
time_data.text = ' '.join(map(str,self.times))
attributes = []
data_items = []
for inc in self.increments:
grid=ET.SubElement(collection,'Grid')
grid.attrib = {'GridType': 'Uniform',
'Name': inc}
topology=ET.SubElement(grid, 'Topology')
topology.attrib={'TopologyType': '3DCoRectMesh',
'Dimensions': '{} {} {}'.format(*self.grid+1)}
geometry=ET.SubElement(grid, 'Geometry')
geometry.attrib={'GeometryType':'Origin_DxDyDz'}
origin=ET.SubElement(geometry, 'DataItem')
origin.attrib={'Format': 'XML',
'NumberType': 'Float',
'Dimensions': '3'}
origin.text="{} {} {}".format(*self.origin)
delta=ET.SubElement(geometry, 'DataItem')
delta.attrib={'Format': 'XML',
'NumberType': 'Float',
'Dimensions': '3'}
delta.text="{} {} {}".format(*(self.size/self.grid))
with h5py.File(self.fname,'r') as f:
attributes.append(ET.SubElement(grid, 'Attribute'))
attributes[-1].attrib={'Name': 'u',
'Center': 'Node',
'AttributeType': 'Vector'}
data_items.append(ET.SubElement(attributes[-1], 'DataItem'))
data_items[-1].attrib={'Format': 'HDF',
'Precision': '8',
'Dimensions': '{} {} {} 3'.format(*(self.grid+1))}
2020-06-25 01:04:51 +05:30
data_items[-1].text=f'{os.path.split(self.fname)[1]}:/{inc}/geometry/u_n'
for o,p in zip(['constituents','materialpoints'],['con_physics','mat_physics']):
for oo in getattr(self,o):
for pp in getattr(self,p):
g = '/'.join([inc,o[:-1],oo,pp])
for l in f[g]:
name = '/'.join([g,l])
shape = f[name].shape[1:]
dtype = f[name].dtype
prec = f[name].dtype.itemsize
if (shape not in [(1,), (3,), (3,3)]) or dtype != np.float64: continue
attributes.append(ET.SubElement(grid, 'Attribute'))
2020-06-25 01:04:51 +05:30
attributes[-1].attrib={'Name': name.split('/',2)[2],
'Center': 'Cell',
'AttributeType': 'Tensor'}
data_items.append(ET.SubElement(attributes[-1], 'DataItem'))
data_items[-1].attrib={'Format': 'HDF',
'NumberType': 'Float',
2020-06-25 01:04:51 +05:30
'Precision': f'{prec}',
'Dimensions': '{} {} {} {}'.format(*self.grid,np.prod(shape))}
2020-06-25 01:04:51 +05:30
data_items[-1].text=f'{os.path.split(self.fname)[1]}:{name}'
2020-06-03 14:13:07 +05:30
with open(self.fname.with_suffix('.xdmf').name,'w') as f:
f.write(xml.dom.minidom.parseString(ET.tostring(xdmf).decode()).toprettyxml())
2020-03-19 16:00:36 +05:30
def to_vtk(self,labels=[],mode='cell'):
2020-02-21 12:15:05 +05:30
"""
Export to vtk cell/point data.
Parameters
----------
2020-03-19 16:00:36 +05:30
labels : str or list of, optional
2020-03-19 12:15:31 +05:30
Labels of the datasets to be exported.
mode : str, either 'cell' or 'point'
2020-03-19 12:15:31 +05:30
Export in cell format or point format.
Defaults to 'cell'.
2020-02-21 12:15:05 +05:30
"""
if mode.lower()=='cell':
2020-02-21 12:15:05 +05:30
2020-03-22 20:43:35 +05:30
if self.structured:
v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin)
else:
with h5py.File(self.fname,'r') as f:
v = VTK.from_unstructuredGrid(f['/geometry/x_n'][()],
f['/geometry/T_c'][()]-1,
f['/geometry/T_c'].attrs['VTK_TYPE'].decode())
2020-02-21 23:22:58 +05:30
elif mode.lower()=='point':
2020-03-12 03:05:58 +05:30
v = VTK.from_polyData(self.cell_coordinates())
N_digits = int(np.floor(np.log10(max(1,int(self.increments[-1][3:])))))+1
2020-02-21 12:15:05 +05:30
2020-03-22 20:43:35 +05:30
for inc in util.show_progress(self.iterate('increments'),len(self.selection['increments'])):
materialpoints_backup = self.selection['materialpoints'].copy()
self.pick('materialpoints',False)
for label in (labels if isinstance(labels,list) else [labels]):
for p in self.iterate('con_physics'):
if p != 'generic':
for c in self.iterate('constituents'):
x = self.get_dataset_location(label)
if len(x) == 0:
continue
array = self.read_dataset(x,0)
v.add(array,'1_'+x[0].split('/',1)[1]) #ToDo: hard coded 1!
else:
x = self.get_dataset_location(label)
if len(x) == 0:
continue
array = self.read_dataset(x,0)
ph_name = re.compile(r'(?<=(constituent\/))(.*?)(?=(generic))') # identify phase name
dset_name = '1_' + re.sub(ph_name,r'',x[0].split('/',1)[1]) # removing phase name
v.add(array,dset_name)
self.pick('materialpoints',materialpoints_backup)
constituents_backup = self.selection['constituents'].copy()
self.pick('constituents',False)
for label in (labels if isinstance(labels,list) else [labels]):
for p in self.iterate('mat_physics'):
if p != 'generic':
for m in self.iterate('materialpoints'):
x = self.get_dataset_location(label)
if len(x) == 0:
continue
array = self.read_dataset(x,0)
v.add(array,'1_'+x[0].split('/',1)[1]) #ToDo: why 1_?
else:
x = self.get_dataset_location(label)
if len(x) == 0:
continue
array = self.read_dataset(x,0)
v.add(array,'1_'+x[0].split('/',1)[1])
self.pick('constituents',constituents_backup)
u = self.read_dataset(self.get_dataset_location('u_n' if mode.lower() == 'cell' else 'u_p'))
v.add(u,'u')
2020-06-25 01:04:51 +05:30
v.write(f'{self.fname.stem}_inc{inc[3:].zfill(N_digits)}')