Merge branch 'python-documentation' into 'development'

improvements to python documentation

See merge request damask/DAMASK!371
This commit is contained in:
Philip Eisenlohr 2021-04-25 15:16:04 +00:00
commit b7777bc746
20 changed files with 989 additions and 278 deletions

@ -1 +1 @@
Subproject commit 23eac08c9f9638f8dae76710095222d00f948eec
Subproject commit 7f0594060779d9a8a4e774d558134309ab77b96e

View File

@ -5,6 +5,7 @@ references:
10.1016/j.ijplas.2020.102779
- K. Sedighiani et al.,
Mechanics of Materials, submitted
output: [rho_dip, rho_mob]
N_sl: [12, 12]
b_sl: [2.49e-10, 2.49e-10]
rho_mob_0: [2.81e12, 2.8e12]

View File

@ -1,11 +1,4 @@
"""
Tools for pre and post processing of DAMASK simulations.
Modules that contain only one class (of the same name),
are prefixed by a '_'. For example, '_colormap' contains
a class called 'Colormap' which is imported as 'damask.Colormap'.
"""
"""Tools for managing DAMASK simulations."""
from pathlib import Path as _Path
import re as _re
@ -15,7 +8,6 @@ with open(_Path(__file__).parent/_Path('VERSION')) as _f:
version = _re.sub(r'^v','',_f.readline().strip())
__version__ = version
# make classes directly accessible as damask.Class
from . import util # noqa
from . import seeds # noqa
from . import tensor # noqa
@ -23,6 +15,8 @@ from . import mechanics # noqa
from . import solver # noqa
from . import grid_filters # noqa
from . import lattice # noqa
#Modules that contain only one class (of the same name), are prefixed by a '_'.
#For example, '_colormap' containsa class called 'Colormap' which is imported as 'damask.Colormap'.
from ._rotation import Rotation # noqa
from ._orientation import Orientation # noqa
from ._table import Table # noqa

View File

@ -91,6 +91,16 @@ class Colormap(mpl.colors.ListedColormap):
- 'lab': CIE Lab.
- 'msh': Msh (for perceptual uniform interpolation).
Returns
-------
new : damask.Colormap
Colormap within given bounds.
Examples
--------
>>> import damask
>>> damask.Colormap.from_range((0,0,1),(0,0,0),'blue_to_black')
"""
low_high = np.vstack((low,high))
if model.lower() == 'rgb':
@ -150,6 +160,16 @@ class Colormap(mpl.colors.ListedColormap):
This parameter is not used for matplotlib colormaps
that are of type `ListedColormap`.
Returns
-------
new : damask.Colormap
Predefined colormap.
Examples
--------
>>> import damask
>>> damask.Colormap.from_predefined('strain')
"""
# matplotlib presets
try:
@ -220,6 +240,11 @@ class Colormap(mpl.colors.ListedColormap):
damask.Colormap
The reversed colormap.
Examples
--------
>>> import damask
>>> damask.Colormap.from_predefined('stress').reversed()
"""
rev = super(Colormap,self).reversed(name)
return Colormap(np.array(rev.colors),rev.name[:-4] if rev.name.endswith('_r_r') else rev.name)
@ -239,8 +264,8 @@ class Colormap(mpl.colors.ListedColormap):
Returns
-------
f
File handle
f : file object
File handle with write access.
"""
if fname is None:

View File

@ -61,6 +61,11 @@ class Config(dict):
other : damask.Config or dict
Key-value pairs that update self.
Returns
-------
updated : damask.Config
Updated configuration.
"""
duplicate = self.copy()
duplicate.update(other)
@ -81,6 +86,11 @@ class Config(dict):
keys : iterable or scalar
Label of the key(s) to remove.
Returns
-------
updated : damask.Config
Updated configuration.
"""
duplicate = self.copy()
for k in keys if isinstance(keys, Iterable) and not isinstance(keys, str) else [keys]:
@ -98,6 +108,11 @@ class Config(dict):
fname : file, str, or pathlib.Path
Filename or file for writing.
Returns
-------
loaded : damask.Config
Configuration from file.
"""
try:
fhandle = open(fname)

View File

@ -54,60 +54,17 @@ class ConfigMaterial(Config):
Parameters
----------
fname : file, str, or pathlib.Path, optional
Filename or file for writing. Defaults to 'material.yaml'.
Filename or file to read from. Defaults to 'material.yaml'.
Returns
-------
loaded : damask.ConfigMaterial
Material configuration from file.
"""
return super(ConfigMaterial,cls).load(fname)
@staticmethod
def from_table(table,**kwargs):
"""
Generate from an ASCII table.
Parameters
----------
table : damask.Table
Table that contains material information.
**kwargs
Keyword arguments where the key is the name and the value specifies
the label of the data column in the table.
Examples
--------
>>> import damask
>>> import damask.ConfigMaterial as cm
>>> t = damask.Table.load('small.txt')
>>> t
pos pos pos qu qu qu qu phase homog
0 0 0 0 0.19 0.8 0.24 -0.51 Aluminum SX
1 1 0 0 0.8 0.19 0.24 -0.51 Steel SX
1 1 1 0 0.8 0.19 0.24 -0.51 Steel SX
>>> cm.from_table(t,O='qu',phase='phase',homogenization='homog')
material:
- constituents:
- O: [0.19, 0.8, 0.24, -0.51]
v: 1.0
phase: Aluminum
homogenization: SX
- constituents:
- O: [0.8, 0.19, 0.24, -0.51]
v: 1.0
phase: Steel
homogenization: SX
homogenization: {}
phase: {}
"""
kwargs_ = {k:table.get(v) for k,v in kwargs.items()}
_,idx = np.unique(np.hstack(list(kwargs_.values())),return_index=True,axis=0)
idx = np.sort(idx)
kwargs_ = {k:np.atleast_1d(v[idx].squeeze()) for k,v in kwargs_.items()}
return ConfigMaterial().material_add(**kwargs_)
@staticmethod
def load_DREAM3D(fname,
grain_data=None,cell_data=None,cell_ensemble_data='CellEnsembleData',
@ -151,6 +108,11 @@ class ConfigMaterial(Config):
and grain- or cell-wise data. Defaults to None, in which case
it is set as the path that contains _SIMPL_GEOMETRY/SPACING.
Returns
-------
loaded : damask.ConfigMaterial
Material configuration from file.
"""
b = util.DREAM3D_base_group(fname) if base_group is None else base_group
c = util.DREAM3D_cell_data_group(fname) if cell_data is None else cell_data
@ -181,9 +143,74 @@ class ConfigMaterial(Config):
return base_config.material_add(**constituent,homogenization='direct')
@staticmethod
def from_table(table,**kwargs):
"""
Generate from an ASCII table.
Parameters
----------
table : damask.Table
Table that contains material information.
**kwargs
Keyword arguments where the key is the name and the value specifies
the label of the data column in the table.
Returns
-------
new : damask.ConfigMaterial
Material configuration from values in table.
Examples
--------
>>> import damask
>>> import damask.ConfigMaterial as cm
>>> t = damask.Table.load('small.txt')
>>> t
pos pos pos qu qu qu qu phase homog
0 0 0 0 0.19 0.8 0.24 -0.51 Aluminum SX
1 1 0 0 0.8 0.19 0.24 -0.51 Steel SX
1 1 1 0 0.8 0.19 0.24 -0.51 Steel SX
>>> cm.from_table(t,O='qu',phase='phase',homogenization='homog')
material:
- constituents:
- O: [0.19, 0.8, 0.24, -0.51]
v: 1.0
phase: Aluminum
homogenization: SX
- constituents:
- O: [0.8, 0.19, 0.24, -0.51]
v: 1.0
phase: Steel
homogenization: SX
homogenization: {}
phase: {}
"""
kwargs_ = {k:table.get(v) for k,v in kwargs.items()}
_,idx = np.unique(np.hstack(list(kwargs_.values())),return_index=True,axis=0)
idx = np.sort(idx)
kwargs_ = {k:np.atleast_1d(v[idx].squeeze()) for k,v in kwargs_.items()}
return ConfigMaterial().material_add(**kwargs_)
@property
def is_complete(self):
"""Check for completeness."""
"""
Check for completeness.
Only the general file layout is considered.
This check does not consider whether parameters for
a particular phase/homogenization model are missing.
Returns
-------
complete : bool
Whether the material.yaml definition is complete.
"""
ok = True
for top_level in ['homogenization','phase','material']:
ok &= top_level in self
@ -236,7 +263,19 @@ class ConfigMaterial(Config):
@property
def is_valid(self):
"""Check for valid content."""
"""
Check for valid content.
Only the generic file content is considered.
This check does not consider whether parameters for a
particular phase/homogenization mode are out of bounds.
Returns
-------
valid : bool
Whether the material.yaml definition is valid.
"""
ok = True
if 'phase' in self:
@ -282,7 +321,7 @@ class ConfigMaterial(Config):
Returns
-------
cfg : damask.ConfigMaterial
updated : damask.ConfigMaterial
Updated material configuration.
"""
@ -311,7 +350,7 @@ class ConfigMaterial(Config):
Returns
-------
cfg : damask.ConfigMaterial
updated : damask.ConfigMaterial
Updated material configuration.
"""
@ -336,53 +375,61 @@ class ConfigMaterial(Config):
Returns
-------
cfg : damask.ConfigMaterial
updated : damask.ConfigMaterial
Updated material configuration.
Examples
--------
Create a dual-phase steel microstructure for micromechanical simulations:
>>> import numpy as np
>>> import damask
>>> m = damask.ConfigMaterial().material_add(phase = ['Aluminum','Steel'],
... O = damask.Rotation.from_random(2),
... homogenization = 'SX')
>>> m = damask.ConfigMaterial()
>>> m = m.material_add(phase = ['Ferrite','Martensite'],
... O = damask.Rotation.from_random(2),
... homogenization = 'SX')
>>> m
material:
- constituents:
- O: [0.577764, -0.146299, -0.617669, 0.513010]
v: 1.0
phase: Aluminum
phase: Ferrite
homogenization: SX
- constituents:
- O: [0.184176, 0.340305, 0.737247, 0.553840]
v: 1.0
phase: Steel
phase: Martensite
homogenization: SX
homogenization: {}
phase: {}
>>> m = damask.ConfigMaterial().material_add(phase = np.array(['Austenite','Martensite']).reshape(1,2),
... O = damask.Rotation.from_random((2,2)),
... v = np.array([0.2,0.8]).reshape(1,2),
... homogenization = ['A','B'])
Create a duplex stainless steel microstructure for forming simulations:
>>> import numpy as np
>>> import damask
>>> m = damask.ConfigMaterial()
>>> m = m.material_add(phase = np.array(['Austenite','Ferrite']).reshape(1,2),
... O = damask.Rotation.from_random((2,2)),
... v = np.array([0.2,0.8]).reshape(1,2),
... homogenization = 'Taylor')
>>> m
material:
- constituents:
- phase: Austenite
O: [0.659802978293224, 0.6953785848195171, 0.22426295326327111, -0.17554139512785227]
v: 0.2
- phase: Martensite
- phase: Ferrite
O: [0.49356745891301596, 0.2841806579193434, -0.7487679215072818, -0.339085707289975]
v: 0.8
homogenization: A
homogenization: Taylor
- constituents:
- phase: Austenite
O: [0.26542221365204055, 0.7268854930702071, 0.4474726435701472, -0.44828201137283735]
v: 0.2
- phase: Martensite
- phase: Ferrite
O: [0.6545817158479885, -0.08004812803625233, -0.6226561293931374, 0.4212059104577611]
v: 0.8
homogenization: B
homogenization: Taylor
homogenization: {}
phase: {}

View File

@ -161,6 +161,11 @@ class Grid:
Grid file to read. Valid extension is .vtr, which will be appended
if not given.
Returns
-------
loaded : damask.Grid
Grid-based geometry from file.
"""
v = VTK.load(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr')
comments = v.get_comments()
@ -190,6 +195,11 @@ class Grid:
fname : str, pathlib.Path, or file handle
Geometry file to read.
Returns
-------
loaded : damask.Grid
Grid-based geometry from file.
"""
warnings.warn('Support for ASCII-based geom format will be removed in DAMASK 3.1.0', DeprecationWarning,2)
try:
@ -281,6 +291,10 @@ class Grid:
and grain- or cell-wise data. Defaults to None, in which case
it is set as the path that contains _SIMPL_GEOMETRY/SPACING.
Returns
-------
loaded : damask.Grid
Grid-based geometry from file.
"""
b = util.DREAM3D_base_group(fname) if base_group is None else base_group
@ -319,6 +333,11 @@ class Grid:
Label(s) of the columns containing the material definition.
Each unique combination of values results in one material ID.
Returns
-------
new : damask.Grid
Grid-based geometry from values in table.
"""
cells,size,origin = grid_filters.cellsSizeOrigin_coordinates0_point(table.get(coordinates))
@ -356,6 +375,11 @@ class Grid:
periodic : Boolean, optional
Assume grid to be periodic. Defaults to True.
Returns
-------
new : damask.Grid
Grid-based geometry from tessellation.
"""
if periodic:
weights_p = np.tile(weights,27) # Laguerre weights (1,2,3,1,2,3,...,1,2,3)
@ -405,6 +429,11 @@ class Grid:
periodic : Boolean, optional
Assume grid to be periodic. Defaults to True.
Returns
-------
new : damask.Grid
Grid-based geometry from tessellation.
"""
coords = grid_filters.coordinates0_point(cells,size).reshape(-1,3)
KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds)
@ -478,6 +507,11 @@ class Grid:
materials : (int, int), optional
Material IDs. Defaults to (1,2).
Returns
-------
new : damask.Grid
Grid-based geometry from definition of minimal surface.
Notes
-----
The following triply-periodic minimal surfaces are implemented:
@ -598,6 +632,11 @@ class Grid:
periodic : Boolean, optional
Assume grid to be periodic. Defaults to True.
Returns
-------
updated : damask.Grid
Updated grid-based geometry.
"""
# radius and center
r = np.array(dimension)/2.0*self.size/self.cells if np.array(dimension).dtype in np.sctypes['int'] else \
@ -638,6 +677,11 @@ class Grid:
reflect : bool, optional
Reflect (include) outermost layers. Defaults to False.
Returns
-------
updated : damask.Grid
Updated grid-based geometry.
"""
valid = ['x','y','z']
if not set(directions).issubset(valid):
@ -670,6 +714,11 @@ class Grid:
Direction(s) along which the grid is flipped.
Valid entries are 'x', 'y', 'z'.
Returns
-------
updated : damask.Grid
Updated grid-based geometry.
"""
valid = ['x','y','z']
if not set(directions).issubset(valid):
@ -695,6 +744,11 @@ class Grid:
periodic : Boolean, optional
Assume grid to be periodic. Defaults to True.
Returns
-------
updated : damask.Grid
Updated grid-based geometry.
"""
return Grid(material = ndimage.interpolation.zoom(
self.material,
@ -723,6 +777,11 @@ class Grid:
periodic : Boolean, optional
Assume grid to be periodic. Defaults to True.
Returns
-------
updated : damask.Grid
Updated grid-based geometry.
"""
def mostFrequent(arr,selection=None):
me = arr[arr.size//2]
@ -746,7 +805,15 @@ class Grid:
def renumber(self):
"""Renumber sorted material indices as 0,...,N-1."""
"""
Renumber sorted material indices as 0,...,N-1.
Returns
-------
updated : damask.Grid
Updated grid-based geometry.
"""
_,renumbered = np.unique(self.material,return_inverse=True)
return Grid(material = renumbered.reshape(self.cells),
@ -767,6 +834,11 @@ class Grid:
fill : int or float, optional
Material index to fill the corners. Defaults to material.max() + 1.
Returns
-------
updated : damask.Grid
Updated grid-based geometry.
"""
if fill is None: fill = np.nanmax(self.material) + 1
dtype = float if isinstance(fill,float) or self.material.dtype in np.sctypes['float'] else int
@ -802,6 +874,11 @@ class Grid:
fill : int or float, optional
Material index to fill the background. Defaults to material.max() + 1.
Returns
-------
updated : damask.Grid
Updated grid-based geometry.
"""
if offset is None: offset = 0
if fill is None: fill = np.nanmax(self.material) + 1
@ -834,6 +911,11 @@ class Grid:
to_material : iterable of ints
New material indices.
Returns
-------
updated : damask.Grid
Updated grid-based geometry.
"""
def mp(entry,mapper):
return mapper[entry] if entry in mapper else entry
@ -849,7 +931,15 @@ class Grid:
def sort(self):
"""Sort material indices such that min(material) is located at (0,0,0)."""
"""
Sort material indices such that min(material) is located at (0,0,0).
Returns
-------
updated : damask.Grid
Updated grid-based geometry.
"""
a = self.material.flatten(order='F')
from_ma = pd.unique(a)
sort_idx = np.argsort(from_ma)
@ -884,6 +974,11 @@ class Grid:
periodic : Boolean, optional
Assume grid to be periodic. Defaults to True.
Returns
-------
updated : damask.Grid
Updated grid-based geometry.
"""
def tainted_neighborhood(stencil,trigger):
@ -917,6 +1012,11 @@ class Grid:
Direction(s) along which the boundaries are determined.
Valid entries are 'x', 'y', 'z'. Defaults to 'xyz'.
Returns
-------
grain_boundaries : damask.VTK
VTK-based geometry of grain boundary network.
"""
valid = ['x','y','z']
if not set(directions).issubset(valid):

View File

@ -8,15 +8,15 @@ from . import tensor
_parameter_doc = \
"""lattice : str
Either a crystal family out of [triclinic, monoclinic, orthorhombic, tetragonal, hexagonal, cubic]
or a Bravais lattice out of [aP, mP, mS, oP, oS, oI, oF, tP, tI, hP, cP, cI, cF].
When specifying a Bravais lattice, additional lattice parameters might be required:
Either a crystal family out of {triclinic, monoclinic, orthorhombic, tetragonal, hexagonal, cubic}
or a Bravais lattice out of {aP, mP, mS, oP, oS, oI, oF, tP, tI, hP, cP, cI, cF}.
When specifying a Bravais lattice, additional lattice parameters might be required.
a : float, optional
Length of lattice parameter "a".
Length of lattice parameter 'a'.
b : float, optional
Length of lattice parameter "b".
Length of lattice parameter 'b'.
c : float, optional
Length of lattice parameter "c".
Length of lattice parameter 'c'.
alpha : float, optional
Angle between b and c lattice basis.
beta : float, optional

View File

@ -34,7 +34,7 @@ def _read(dataset):
return np.array(dataset,dtype=dtype)
def _match(requested,existing):
"""Find matches among two sets of labels."""
"""Find matches among two sets of names."""
def flatten_list(list_of_lists):
return [e for e_ in list_of_lists for e in e_]
@ -57,11 +57,30 @@ def _empty_like(dataset,N_materialpoints,fill_float,fill_int):
class Result:
"""
Manipulate and read DADF5 files.
Add data to and export data from a DADF5 file.
A DADF5 (DAMASK HDF5) file contains DAMASK results.
Its group/folder structure reflects the layout in material.yaml.
This class provides a customizable view on the DADF5 file.
Upon initialization, all attributes are visible.
Derived quantities are added to the file and existing data is
exported based on the current view.
Examples
--------
Open 'my_file.hdf5', which is assumed to contain deformation gradient 'F'
and first Piola-Kirchhoff stress 'P', add the Mises equivalent of the
Cauchy stress, and export it to VTK (file) and numpy.ndarray (memory).
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r.add_Cauchy()
>>> r.add_equivalent_Mises('sigma')
>>> r.save_VTK()
>>> r_last = r.view('increments',-1)
>>> sigma_vM_last = r_last.get('sigma_vM')
DADF5 (DAMASK HDF5) files contain DAMASK results.
The group/folder structure reflects the input data
in material.yaml.
"""
def __init__(self,fname):
@ -163,6 +182,11 @@ class Result:
Name of datasets; supports '?' and '*' wildcards.
True is equivalent to '*', False is equivalent to [].
Returns
-------
view : damask.Result
Modified or new view on the DADF5 file.
"""
# allow True/False and string arguments
if datasets is True:
@ -176,6 +200,7 @@ class Result:
if what == 'increments':
choice = [c if isinstance(c,str) and c.startswith(inc) else
f'{inc}{c}' for c in choice]
if datasets == -1: choice = [self.increments[-1]]
elif what == 'times':
what = 'increments'
if choice == ['*']:
@ -207,15 +232,31 @@ class Result:
return dup
def allow_modification(self):
"""Allow to overwrite existing data."""
def modification_enable(self):
"""
Allow modification of existing data.
Returns
-------
modified_view : damask.Result
View without write-protection of existing data.
"""
print(util.warn('Warning: Modification of existing datasets allowed!'))
dup = self.copy()
dup._allow_modification = True
return dup
def disallow_modification(self):
"""Disallow to overwrite existing data (default case)."""
def modification_disable(self):
"""
Prevent modification of existing data (default case).
Returns
-------
modified_view : damask.Result
View with write-protection of existing data.
"""
dup = self.copy()
dup._allow_modification = False
return dup
@ -223,7 +264,7 @@ class Result:
def increments_in_range(self,start,end):
"""
Select all increments within a given range.
Get all increments within a given range.
Parameters
----------
@ -232,6 +273,11 @@ class Result:
end : int or str
End increment.
Returns
-------
increments : list of ints
Increment number of all increments within the given bounds.
"""
# compatibility hack
ln = 3 if self.version_minor < 12 else 10
@ -239,13 +285,13 @@ class Result:
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])
selected.append(int(self.increments[i].split('_')[1]))
return selected
def times_in_range(self,start,end):
"""
Select all increments within a given time range.
Get all increments within a given time range.
Parameters
----------
@ -254,6 +300,11 @@ class Result:
end : float
Time of end increment.
Returns
-------
times : list of float
Simulation time of all increments within the given bounds.
"""
selected = []
for i,time in enumerate(self.times):
@ -274,6 +325,25 @@ class Result:
Name of datasets; supports '?' and '*' wildcards.
True is equivalent to '*', False is equivalent to [].
Returns
-------
view : damask.Result
View with only the selected attributes being visible.
Examples
--------
Get a view that shows only results from the initial configuration:
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r_first = r.view('increment',0)
Get a view that shows all results between simulation times of 10 to 40:
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r_t10to40 = r.view('times',r.times_in_range(10.0,40.0))
"""
return self._manage_view('set',what,datasets)
@ -290,13 +360,27 @@ class Result:
Name of datasets; supports '?' and '*' wildcards.
True is equivalent to '*', False is equivalent to [].
Returns
-------
modified_view : damask.Result
View with additional visible attributes.
Examples
--------
Get a view that shows only results from first and last increment:
>>> import damask
>>> r_empty = damask.Result('my_file.hdf5').view('increments',False)
>>> r_first = r_empty.view_more('increments',0)
>>> r_first_and_last = r.first.view_more('increments',-1)
"""
return self._manage_view('add',what,datasets)
def view_less(self,what,datasets):
"""
Delete from view.
Remove from view.
Parameters
----------
@ -306,37 +390,98 @@ class Result:
Name of datasets; supports '?' and '*' wildcards.
True is equivalent to '*', False is equivalent to [].
Returns
-------
modified_view : damask.Result
View with fewer visible attributes.
Examples
--------
Get a view that omits the undeformed configuration:
>>> import damask
>>> r_all = damask.Result('my_file.hdf5')
>>> r_deformed = r_all.view_less('increments',0)
"""
return self._manage_view('del',what,datasets)
def rename(self,name_old,name_new):
def rename(self,name_src,name_dst):
"""
Rename dataset.
Rename/move datasets (within the same group/folder).
This operation is discouraged because the history of the
data becomes untracable and scientific integrity cannot be
ensured.
Parameters
----------
name_old : str
Name of the dataset to be renamed.
name_new : str
New name of the dataset.
name_src : str
Name of the datasets to be renamed.
name_dst : str
New name of the datasets.
Examples
--------
Rename datasets containing the deformation gradient from 'F' to 'def_grad':
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r_unprotected = r.modification_enable()
>>> r_unprotected.rename('F','def_grad')
"""
if not self._allow_modification:
raise PermissionError('Rename operation not permitted')
raise PermissionError('Renaming datasets not permitted')
with h5py.File(self.fname,'a') as f:
for inc in self.visible['increments']:
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()):
path_old = '/'.join([inc,ty,label,field,name_old])
path_new = '/'.join([inc,ty,label,field,name_new])
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]
path_src = '/'.join([inc,ty,label,field,name_src])
path_dst = '/'.join([inc,ty,label,field,name_dst])
if path_src in f.keys():
f[path_dst] = f[path_src]
f[path_dst].attrs['renamed'] = f'original name: {name_src}' if h5py3 else \
f'original name: {name_src}'.encode()
del f[path_src]
def remove(self,name):
"""
Remove/delete datasets.
This operation is discouraged because the history of the
data becomes untracable and scientific integrity cannot be
ensured.
Parameters
----------
name : str
Name of the datasets to be deleted.
Examples
--------
Delete the deformation gradient 'F':
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r_unprotected = r.modification_enable()
>>> r_unprotected.remove('F')
"""
if not self._allow_modification:
raise PermissionError('Removing datasets not permitted')
with h5py.File(self.fname,'a') as f:
for inc in self.visible['increments']:
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()):
path = '/'.join([inc,ty,label,field,name])
if path in f.keys(): del f[path]
def list_data(self):
@ -372,7 +517,7 @@ class Result:
@property
def coordinates0_point(self):
"""Return initial coordinates of the cell centers."""
"""Initial/undeformed cell center coordinates."""
if self.structured:
return grid_filters.coordinates0_point(self.cells,self.size,self.origin).reshape(-1,3,order='F')
else:
@ -381,7 +526,7 @@ class Result:
@property
def coordinates0_node(self):
"""Return initial coordinates of the cell centers."""
"""Initial/undeformed nodal coordinates."""
if self.structured:
return grid_filters.coordinates0_node(self.cells,self.size,self.origin).reshape(-1,3,order='F')
else:
@ -390,6 +535,7 @@ class Result:
@property
def geometry0(self):
"""Initial/undeformed geometry."""
if self.structured:
return VTK.from_rectilinear_grid(self.cells,self.size,self.origin)
else:
@ -418,7 +564,7 @@ class Result:
Parameters
----------
x : str
Label of scalar, vector, or tensor dataset to take absolute value of.
Name of scalar, vector, or tensor dataset to take absolute value of.
"""
self._add_generic_pointwise(self._add_absolute,{'x':x})
@ -439,24 +585,52 @@ class Result:
'creator': 'add_calculation'
}
}
def add_calculation(self,label,formula,unit='n/a',description=None):
def add_calculation(self,formula,name,unit='n/a',description=None):
"""
Add result of a general formula.
Parameters
----------
label : str
Label of resulting dataset.
formula : str
Formula to calculate resulting dataset. Existing datasets are referenced by '#TheirLabel#'.
Formula to calculate resulting dataset.
Existing datasets are referenced by '#TheirName#'.
name : str
Name of resulting dataset.
unit : str, optional
Physical unit of the result.
description : str, optional
Human-readable description of the result.
Examples
--------
Add total dislocation density, i.e. the sum of mobile dislocation
density 'rho_mob' and dislocation dipole density 'rho_dip' over
all slip systems:
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r.add_calculation('np.sum(#rho_mob#,axis=1)','rho_mob_total',
... '1/m²','total mobile dislocation density')
>>> r.add_calculation(''np.sum(#rho_dip#,axis=1)',rho_dip_total',
... '1/m²','total dislocation dipole density')
>>> r.add_calculation('#rho_dip_total#+#rho_mob_total','rho_total',
... '1/m²','total dislocation density')
Add Mises equivalent of the Cauchy stress without storage of
intermediate results. Define a user function for better readability:
>>> import damask
>>> def equivalent_stress(F,P):
... sigma = damask.mechanics.stress_Cauchy(F=F,P=P)
... return damask.mechanics.equivalent_stress_Mises(sigma)
>>> r = damask.Result('my_file.hdf5')
>>> r.enable_user_function(equivalent_stress)
>>> r.add_calculation('equivalent_stress(#F#,#P#)','sigma_vM','Pa',
... 'Mises equivalent of the Cauchy stress')
"""
dataset_mapping = {d:d for d in set(re.findall(r'#(.*?)#',formula))} # datasets used in the formula
args = {'formula':formula,'label':label,'unit':unit,'description':description}
args = {'formula':formula,'label':name,'unit':unit,'description':description}
self._add_generic_pointwise(self._add_calculation,dataset_mapping,args)
@ -480,9 +654,9 @@ class Result:
Parameters
----------
P : str, optional
Label of the dataset containing the first Piola-Kirchhoff stress. Defaults to 'P'.
Name of the dataset containing the first Piola-Kirchhoff stress. Defaults to 'P'.
F : str, optional
Label of the dataset containing the deformation gradient. Defaults to 'F'.
Name of the dataset containing the deformation gradient. Defaults to 'F'.
"""
self._add_generic_pointwise(self._add_stress_Cauchy,{'P':P,'F':F})
@ -506,7 +680,15 @@ class Result:
Parameters
----------
T : str
Label of tensor dataset.
Name of tensor dataset.
Examples
--------
Add the determinant of plastic deformation gradient 'F_p':
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r.add_determinant('F_p')
"""
self._add_generic_pointwise(self._add_determinant,{'T':T})
@ -530,7 +712,7 @@ class Result:
Parameters
----------
T : str
Label of tensor dataset.
Name of tensor dataset.
"""
self._add_generic_pointwise(self._add_deviator,{'T':T})
@ -561,7 +743,7 @@ class Result:
Parameters
----------
T_sym : str
Label of symmetric tensor dataset.
Name of symmetric tensor dataset.
eigenvalue : str, optional
Eigenvalue. Select from 'max', 'mid', 'min'. Defaults to 'max'.
@ -594,7 +776,7 @@ class Result:
Parameters
----------
T_sym : str
Label of symmetric tensor dataset.
Name of symmetric tensor dataset.
eigenvalue : str, optional
Eigenvalue to which the eigenvector corresponds.
Select from 'max', 'mid', 'min'. Defaults to 'max'.
@ -634,9 +816,17 @@ class Result:
l : numpy.array of shape (3)
Lab frame direction for inverse pole figure.
q : str
Label of the dataset containing the crystallographic orientation as quaternions.
Name of the dataset containing the crystallographic orientation as quaternions.
Defaults to 'O'.
Examples
--------
Add the IPF color along [0,1,1] for orientation 'O':
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r.add_IPF_color(np.array([0,1,1]))
"""
self._add_generic_pointwise(self._add_IPF_color,{'q':q},{'l':l})
@ -659,7 +849,7 @@ class Result:
Parameters
----------
T_sym : str
Label of symmetric tensor dataset.
Name of symmetric tensor dataset.
"""
self._add_generic_pointwise(self._add_maximum_shear,{'T_sym':T_sym})
@ -693,11 +883,25 @@ class Result:
Parameters
----------
T_sym : str
Label of symmetric tensorial stress or strain dataset.
Name of symmetric tensorial stress or strain dataset.
kind : {'stress', 'strain', None}, optional
Kind of the von Mises equivalent. Defaults to None, in which case
it is selected based on the unit of the dataset ('1' -> strain, 'Pa' -> stress).
Examples
--------
Add the Mises equivalent of the Cauchy stress 'sigma':
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r.add_equivalent_Mises('sigma')
Add the Mises equivalent of the spatial logarithmic strain 'epsilon_V^0.0(F)':
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r.add_equivalent_Mises('epsilon_V^0.0(F)')
"""
self._add_generic_pointwise(self._add_equivalent_Mises,{'T_sym':T_sym},{'kind':kind})
@ -732,7 +936,7 @@ class Result:
Parameters
----------
x : str
Label of vector or tensor dataset.
Name of vector or tensor dataset.
ord : {non-zero int, inf, -inf, 'fro', 'nuc'}, optional
Order of the norm. inf means NumPys inf object. For details refer to numpy.linalg.norm.
@ -760,9 +964,9 @@ class Result:
Parameters
----------
P : str, optional
Label of first Piola-Kirchhoff stress dataset. Defaults to 'P'.
Name of first Piola-Kirchhoff stress dataset. Defaults to 'P'.
F : str, optional
Label of deformation gradient dataset. Defaults to 'F'.
Name of deformation gradient dataset. Defaults to 'F'.
"""
self._add_generic_pointwise(self._add_stress_second_Piola_Kirchhoff,{'P':P,'F':F})
@ -801,7 +1005,7 @@ class Result:
# Parameters
# ----------
# q : str
# Label of the dataset containing the crystallographic orientation as quaternions.
# Name of the dataset containing the crystallographic orientation as quaternions.
# p : numpy.array of shape (3)
# Crystallographic direction or plane.
# polar : bool, optional
@ -828,8 +1032,16 @@ class Result:
Parameters
----------
F : str, optional
Label of deformation gradient dataset.
F : str
Name of deformation gradient dataset.
Examples
--------
Add the rotational part of deformation gradient 'F':
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r.add_rotation('F')
"""
self._add_generic_pointwise(self._add_rotation,{'F':F})
@ -853,7 +1065,15 @@ class Result:
Parameters
----------
T : str
Label of tensor dataset.
Name of tensor dataset.
Examples
--------
Add the hydrostatic part of the Cauchy stress 'sigma':
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r.add_spherical('sigma')
"""
self._add_generic_pointwise(self._add_spherical,{'T':T})
@ -879,13 +1099,28 @@ class Result:
Parameters
----------
F : str, optional
Label of deformation gradient dataset. Defaults to 'F'.
Name 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'.
m : float, optional
Order of the strain calculation. Defaults to 0.0.
Examples
--------
Add the Biot strain based on the deformation gradient 'F':
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r.strain(t='U',m=0.5)
Add the plastic Euler-Almansi strain based on the
plastic deformation gradient 'F_p':
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r.strain('F_p','V',-1)
"""
self._add_generic_pointwise(self._add_strain,{'F':F},{'t':t,'m':m})
@ -909,7 +1144,7 @@ class Result:
Parameters
----------
F : str, optional
Label of deformation gradient dataset. Defaults to 'F'.
Name 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'.
@ -1016,10 +1251,14 @@ class Result:
"""
Write XDMF file to directly visualize data in DADF5 file.
The XDMF format is only supported for structured grids
with single phase and single constituent.
For other cases use `save_VTK`.
Parameters
----------
output : (list of) str
Labels of the datasets to read.
Names of the datasets included in the XDMF file.
Defaults to '*', in which case all datasets are considered.
"""
@ -1149,10 +1388,16 @@ class Result:
"""
Export to VTK cell/point data.
One VTK file per visible increment is created.
For cell data, the VTK format is a rectilinear grid (.vtr) for
grid-based simulations and an unstructured grid (.vtu) for
mesh-baed simulations. For point data, the VTK format is poly
data (.vtp).
Parameters
----------
output : (list of) str, optional
Labels of the datasets to place.
Names of the datasets to export to the VTK file.
Defaults to '*', in which case all datasets are exported.
mode : {'cell', 'point'}
Export in cell format or point format.
@ -1231,7 +1476,7 @@ class Result:
Parameters
----------
output : (list of) str
Labels of the datasets to read.
Names of the datasets to read.
Defaults to '*', in which case all datasets are read.
flatten : bool
Remove singular levels of the folder hierarchy.
@ -1277,13 +1522,13 @@ class Result:
in the DADF5 file.
Multi-phase data is fused into a single output.
`place` is equivalent to `read` if only one phase/homogenization
`place` is equivalent to `get` if only one phase/homogenization
and one constituent is present.
Parameters
----------
output : (list of) str, optional
Labels of the datasets to place.
Names of the datasets to read.
Defaults to '*', in which case all datasets are placed.
flatten : bool
Remove singular levels of the folder hierarchy.
@ -1293,7 +1538,7 @@ class Result:
Remove branches with no data. Defaults to True.
constituents : (list of) int, optional
Constituents to consider.
Defaults to 'None', in which case all constituents are considered.
Defaults to None, in which case all constituents are considered.
fill_float : float
Fill value for non-existent entries of floating point type.
Defaults to NaN.

View File

@ -24,7 +24,7 @@ class Rotation:
- Euler angle triplets are implemented using the Bunge convention,
with angular ranges of [0,2π], [0,π], [0,2π].
- The rotation angle ω is limited to the interval [0,π].
- The real part of a quaternion is positive, Re(q) > 0
- The real part of a quaternion is positive, Re(q) 0
- P = -1 (as default).
Examples
@ -357,7 +357,7 @@ class Rotation:
Parameters
----------
other : Rotation or list of Rotations.
other : damask.Rotation
"""
return self.copy(rotation=np.vstack(tuple(map(lambda x:x.quaternion,
@ -365,12 +365,28 @@ class Rotation:
def flatten(self,order = 'C'):
"""Flatten array."""
"""
Flatten array.
Returns
-------
flattened : damask.Rotation
Rotation flattened to single dimension.
"""
return self.copy(rotation=self.quaternion.reshape((-1,4),order=order))
def reshape(self,shape,order = 'C'):
"""Reshape array."""
"""
Reshape array.
Returns
-------
reshaped : damask.Rotation
Rotation of given shape.
"""
if isinstance(shape,(int,np.integer)): shape = (shape,)
return self.copy(rotation=self.quaternion.reshape(tuple(shape)+(4,),order=order))
@ -387,6 +403,11 @@ class Rotation:
Where to preferentially locate missing dimensions.
Either 'left' or 'right' (default).
Returns
-------
broadcasted : damask.Rotation
Rotation broadcasted to given shape.
"""
if isinstance(shape,(int,np.integer)): shape = (shape,)
return self.copy(rotation=np.broadcast_to(self.quaternion.reshape(util.shapeshifter(self.shape,shape,mode)+(4,)),
@ -404,7 +425,7 @@ class Rotation:
Returns
-------
average : Rotation
average : damask.Rotation
Weighted average of original Rotation field.
References
@ -438,9 +459,14 @@ class Rotation:
Parameters
----------
other : Rotation
other : damask.Rotation
Rotation to which the misorientation is computed.
Returns
-------
g : damask.Rotation
Misorientation.
"""
return other*~self

View File

@ -49,7 +49,7 @@ class Table:
Returns
-------
slice : Table
slice : damask.Table
Sliced part of the Table.
Examples
@ -157,7 +157,7 @@ class Table:
Parameters
----------
other : Table
other : damask.Table
Table to compare against.
rtol : float, optional
Relative tolerance of equality.
@ -185,7 +185,7 @@ class Table:
Parameters
----------
other : Table
other : damask.Table
Table to compare against.
rtol : float, optional
Relative tolerance of equality.
@ -223,6 +223,11 @@ class Table:
fname : file, str, or pathlib.Path
Filename or file for reading.
Returns
-------
loaded : damask.Table
Table data from file.
"""
try:
f = open(fname)
@ -275,6 +280,11 @@ class Table:
fname : file, str, or pathlib.Path
Filename or file for reading.
Returns
-------
loaded : damask.Table
Table data from file.
"""
try:
f = open(fname)
@ -334,14 +344,14 @@ class Table:
----------
label : str
Column label.
data : np.ndarray
data : numpy.ndarray
New data.
info : str, optional
Human-readable information about the new data.
Returns
-------
table : Table
updated : damask.Table
Updated table.
"""
@ -367,14 +377,14 @@ class Table:
----------
label : str
Column label.
data : np.ndarray
data : numpy.ndarray
Modified data.
info : str, optional
Human-readable information about the modified data.
Returns
-------
table : Table
udated : damask.Table
Updated table.
"""
@ -402,7 +412,7 @@ class Table:
Returns
-------
table : Table
udated : damask.Table
Updated table.
"""
@ -425,7 +435,7 @@ class Table:
Returns
-------
table : Table
udated : damask.Table
Updated table.
"""
@ -451,7 +461,7 @@ class Table:
Returns
-------
table : Table
udated : damask.Table
Updated table.
"""
@ -479,13 +489,13 @@ class Table:
Parameters
----------
other : Table
other : damask.Table
Table to append.
Returns
-------
table : Table
Concatenated table.
udated : damask.Table
Updated table.
"""
if self.shapes != other.shapes or not self.data.columns.equals(other.data.columns):
@ -504,13 +514,13 @@ class Table:
Parameters
----------
other : Table
other : damask.Table
Table to join.
Returns
-------
table : Table
Joined table.
udated : damask.Table
Updated table.
"""
if set(self.shapes) & set(other.shapes) or self.data.shape[0] != other.data.shape[0]:

View File

@ -3,7 +3,6 @@ import multiprocessing as mp
from pathlib import Path
import numpy as np
import numpy.ma as ma
import vtk
from vtk.util.numpy_support import numpy_to_vtk as np_to_vtk
from vtk.util.numpy_support import numpy_to_vtkIdTypeArray as np_to_vtkIdTypeArray
@ -47,9 +46,14 @@ class VTK:
grid : iterable of int, len (3)
Number of cells along each dimension.
size : iterable of float, len (3)
Physical lengths along each dimension.
Physical length along each dimension.
origin : iterable of float, len (3), optional
Spatial origin coordinates.
Coordinates of grid origin.
Returns
-------
new : damask.VTK
VTK-based geometry without nodal or cell data.
"""
vtk_data = vtk.vtkRectilinearGrid()
@ -68,7 +72,7 @@ class VTK:
"""
Create VTK of type vtk.vtkUnstructuredGrid.
This is the common type for FEM solver results.
This is the common type for mesh solver results.
Parameters
----------
@ -80,6 +84,11 @@ class VTK:
cell_type : str
Name of the vtk.vtkCell subclass. Tested for TRIANGLE, QUAD, TETRA, and HEXAHEDRON.
Returns
-------
new : damask.VTK
VTK-based geometry without nodal or cell data.
"""
vtk_nodes = vtk.vtkPoints()
vtk_nodes.SetData(np_to_vtk(nodes))
@ -110,6 +119,11 @@ class VTK:
points : numpy.ndarray of shape (:,3)
Spatial position of the points.
Returns
-------
new : damask.VTK
VTK-based geometry without nodal or cell data.
"""
N = points.shape[0]
vtk_points = vtk.vtkPoints()
@ -137,25 +151,30 @@ class VTK:
fname : str or pathlib.Path
Filename for reading. Valid extensions are .vtr, .vtu, .vtp, and .vtk.
dataset_type : str, optional
Name of the vtk.vtkDataSet subclass when opening a .vtk file. Valid types are vtkRectilinearGrid,
vtkUnstructuredGrid, and vtkPolyData.
Name of the vtk.vtkDataSet subclass when opening a .vtk file.
Valid types are vtkRectilinearGrid, vtkUnstructuredGrid, and vtkPolyData.
Returns
-------
loaded : damask.VTK
VTK-based geometry from file.
"""
if not os.path.isfile(fname): # vtk has a strange error handling
raise FileNotFoundError(f'no such file: {fname}')
raise FileNotFoundError(f'No such file: {fname}')
ext = Path(fname).suffix
if ext == '.vtk' or dataset_type is not None:
reader = vtk.vtkGenericDataObjectReader()
reader.SetFileName(str(fname))
if dataset_type is None:
raise TypeError('Dataset type for *.vtk file not given.')
elif dataset_type.lower().endswith('rectilineargrid'):
elif dataset_type.lower().endswith(('rectilineargrid','rectilinear_grid')):
reader.Update()
vtk_data = reader.GetRectilinearGridOutput()
elif dataset_type.lower().endswith('unstructuredgrid'):
elif dataset_type.lower().endswith(('unstructuredgrid','unstructured_grid')):
reader.Update()
vtk_data = reader.GetUnstructuredGridOutput()
elif dataset_type.lower().endswith('polydata'):
elif dataset_type.lower().endswith(('polydata','poly_data')):
reader.Update()
vtk_data = reader.GetPolyDataOutput()
else:
@ -246,7 +265,7 @@ class VTK:
raise ValueError('No label defined for numpy.ndarray')
N_data = data.shape[0]
data_ = np.where(data.mask,data.fill_value,data) if isinstance(data,ma.MaskedArray) else\
data_ = np.where(data.mask,data.fill_value,data) if isinstance(data,np.ma.MaskedArray) else\
data
d = np_to_vtk((data_.astype(np.single) if data_.dtype in [np.double, np.longdouble] else
data_).reshape(N_data,-1),deep=True) # avoid large files
@ -277,6 +296,11 @@ class VTK:
label : str
Data label.
Returns
-------
data : numpy.ndarray
Data stored under the given label.
"""
cell_data = self.vtk_data.GetCellData()
for a in range(cell_data.GetNumberOfArrays()):

View File

@ -1,15 +1,14 @@
"""
Filters for operations on regular grids.
Notes
-----
The grids are defined as (x,y,z,...) where x is fastest and z is slowest.
This convention is consistent with the layout in grid vtr files.
When converting to/from a plain list (e.g. storage in ASCII table),
the following operations are required for tensorial data:
D3 = D1.reshape(cells+(-1,),order='F').reshape(cells+(3,3))
D1 = D3.reshape(cells+(-1,)).reshape(-1,9,order='F')
- D3 = D1.reshape(cells+(-1,),order='F').reshape(cells+(3,3))
- D1 = D3.reshape(cells+(-1,)).reshape(-1,9,order='F')
"""
from scipy import spatial as _spatial
@ -29,10 +28,12 @@ def _ks(size,cells,first_order=False):
Correction for first order derivatives, defaults to False.
"""
k_sk = _np.where(_np.arange(cells[0])>cells[0]//2,_np.arange(cells[0])-cells[0],_np.arange(cells[0]))/size[0]
k_sk = _np.where(_np.arange(cells[0])>cells[0]//2,
_np.arange(cells[0])-cells[0],_np.arange(cells[0]))/size[0]
if cells[0]%2 == 0 and first_order: k_sk[cells[0]//2] = 0 # Nyquist freq=0 for even cells (Johnson, MIT, 2011)
k_sj = _np.where(_np.arange(cells[1])>cells[1]//2,_np.arange(cells[1])-cells[1],_np.arange(cells[1]))/size[1]
k_sj = _np.where(_np.arange(cells[1])>cells[1]//2,
_np.arange(cells[1])-cells[1],_np.arange(cells[1]))/size[1]
if cells[1]%2 == 0 and first_order: k_sj[cells[1]//2] = 0 # Nyquist freq=0 for even cells (Johnson, MIT, 2011)
k_si = _np.arange(cells[2]//2+1)/size[2]
@ -40,74 +41,89 @@ def _ks(size,cells,first_order=False):
return _np.stack(_np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij'), axis=-1)
def curl(size,field):
"""
def curl(size,f):
u"""
Calculate curl of a vector or tensor field in Fourier space.
Parameters
----------
size : numpy.ndarray of shape (3)
Physical size of the periodic field.
field : numpy.ndarray of shape (:,:,:,3) or (:,:,:,3,3)
f : numpy.ndarray of shape (:,:,:,3) or (:,:,:,3,3)
Periodic field of which the curl is calculated.
Returns
-------
× f : numpy.ndarray
Curl of f.
"""
n = _np.prod(field.shape[3:])
k_s = _ks(size,field.shape[:3],True)
n = _np.prod(f.shape[3:])
k_s = _ks(size,f.shape[:3],True)
e = _np.zeros((3, 3, 3))
e[0, 1, 2] = e[1, 2, 0] = e[2, 0, 1] = +1.0 # Levi-Civita symbol
e[0, 2, 1] = e[2, 1, 0] = e[1, 0, 2] = -1.0
field_fourier = _np.fft.rfftn(field,axes=(0,1,2))
curl_ = (_np.einsum('slm,ijkl,ijkm ->ijks', e,k_s,field_fourier)*2.0j*_np.pi if n == 3 else # vector, 3 -> 3
_np.einsum('slm,ijkl,ijknm->ijksn',e,k_s,field_fourier)*2.0j*_np.pi) # tensor, 3x3 -> 3x3
f_fourier = _np.fft.rfftn(f,axes=(0,1,2))
curl_ = (_np.einsum('slm,ijkl,ijkm ->ijks', e,k_s,f_fourier)*2.0j*_np.pi if n == 3 else # vector, 3 -> 3
_np.einsum('slm,ijkl,ijknm->ijksn',e,k_s,f_fourier)*2.0j*_np.pi) # tensor, 3x3 -> 3x3
return _np.fft.irfftn(curl_,axes=(0,1,2),s=field.shape[:3])
return _np.fft.irfftn(curl_,axes=(0,1,2),s=f.shape[:3])
def divergence(size,field):
"""
def divergence(size,f):
u"""
Calculate divergence of a vector or tensor field in Fourier space.
Parameters
----------
size : numpy.ndarray of shape (3)
Physical size of the periodic field.
field : numpy.ndarray of shape (:,:,:,3) or (:,:,:,3,3)
f : numpy.ndarray of shape (:,:,:,3) or (:,:,:,3,3)
Periodic field of which the divergence is calculated.
Returns
-------
· f : numpy.ndarray
Divergence of f.
"""
n = _np.prod(field.shape[3:])
k_s = _ks(size,field.shape[:3],True)
n = _np.prod(f.shape[3:])
k_s = _ks(size,f.shape[:3],True)
field_fourier = _np.fft.rfftn(field,axes=(0,1,2))
div_ = (_np.einsum('ijkl,ijkl ->ijk', k_s,field_fourier)*2.0j*_np.pi if n == 3 else # vector, 3 -> 1
_np.einsum('ijkm,ijklm->ijkl',k_s,field_fourier)*2.0j*_np.pi) # tensor, 3x3 -> 3
f_fourier = _np.fft.rfftn(f,axes=(0,1,2))
div_ = (_np.einsum('ijkl,ijkl ->ijk', k_s,f_fourier)*2.0j*_np.pi if n == 3 else # vector, 3 -> 1
_np.einsum('ijkm,ijklm->ijkl',k_s,f_fourier)*2.0j*_np.pi) # tensor, 3x3 -> 3
return _np.fft.irfftn(div_,axes=(0,1,2),s=field.shape[:3])
return _np.fft.irfftn(div_,axes=(0,1,2),s=f.shape[:3])
def gradient(size,field):
"""
Calculate gradient of a scalar or vector field in Fourier space.
def gradient(size,f):
u"""
Calculate gradient of a scalar or vector fieldin Fourier space.
Parameters
----------
size : numpy.ndarray of shape (3)
Physical size of the periodic field.
field : numpy.ndarray of shape (:,:,:,1) or (:,:,:,3)
f : numpy.ndarray of shape (:,:,:,1) or (:,:,:,3)
Periodic field of which the gradient is calculated.
Returns
-------
f : numpy.ndarray
Divergence of f.
"""
n = _np.prod(field.shape[3:])
k_s = _ks(size,field.shape[:3],True)
n = _np.prod(f.shape[3:])
k_s = _ks(size,f.shape[:3],True)
field_fourier = _np.fft.rfftn(field,axes=(0,1,2))
grad_ = (_np.einsum('ijkl,ijkm->ijkm', field_fourier,k_s)*2.0j*_np.pi if n == 1 else # scalar, 1 -> 3
_np.einsum('ijkl,ijkm->ijklm',field_fourier,k_s)*2.0j*_np.pi) # vector, 3 -> 3x3
f_fourier = _np.fft.rfftn(f,axes=(0,1,2))
grad_ = (_np.einsum('ijkl,ijkm->ijkm', f_fourier,k_s)*2.0j*_np.pi if n == 1 else # scalar, 1 -> 3
_np.einsum('ijkl,ijkm->ijklm',f_fourier,k_s)*2.0j*_np.pi) # vector, 3 -> 3x3
return _np.fft.irfftn(grad_,axes=(0,1,2),s=field.shape[:3])
return _np.fft.irfftn(grad_,axes=(0,1,2),s=f.shape[:3])
def coordinates0_point(cells,size,origin=_np.zeros(3)):
@ -123,6 +139,11 @@ def coordinates0_point(cells,size,origin=_np.zeros(3)):
origin : numpy.ndarray, optional
Physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
Returns
-------
x_p_0 : numpy.ndarray
Undeformed cell center coordinates.
"""
start = origin + size/cells*.5
end = origin + size - size/cells*.5
@ -144,6 +165,11 @@ def displacement_fluct_point(size,F):
F : numpy.ndarray
Deformation gradient field.
Returns
-------
u_p_fluct : numpy.ndarray
Fluctuating part of the cell center displacements.
"""
integrator = 0.5j*size/_np.pi
@ -171,6 +197,11 @@ def displacement_avg_point(size,F):
F : numpy.ndarray
Deformation gradient field.
Returns
-------
u_p_avg : numpy.ndarray
Average part of the cell center displacements.
"""
F_avg = _np.average(F,axis=(0,1,2))
return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),coordinates0_point(F.shape[:3],size))
@ -187,6 +218,11 @@ def displacement_point(size,F):
F : numpy.ndarray
Deformation gradient field.
Returns
-------
u_p : numpy.ndarray
Cell center displacements.
"""
return displacement_avg_point(size,F) + displacement_fluct_point(size,F)
@ -204,6 +240,11 @@ def coordinates_point(size,F,origin=_np.zeros(3)):
origin : numpy.ndarray of shape (3), optional
Physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
Returns
-------
x_p : numpy.ndarray
Cell center coordinates.
"""
return coordinates0_point(F.shape[:3],size,origin) + displacement_point(size,F)
@ -220,6 +261,11 @@ def cellsSizeOrigin_coordinates0_point(coordinates0,ordered=True):
Expect coordinates0 data to be ordered (x fast, z slow).
Defaults to True.
Returns
-------
cells, size, origin : Three numpy.ndarray, each of shape (3)
Information to reconstruct grid.
"""
coords = [_np.unique(coordinates0[:,i]) for i in range(3)]
mincorner = _np.array(list(map(min,coords)))
@ -252,19 +298,6 @@ def cellsSizeOrigin_coordinates0_point(coordinates0,ordered=True):
return (cells,size,origin)
def coordinates0_check(coordinates0):
"""
Check whether coordinates lie on a regular grid.
Parameters
----------
coordinates0 : numpy.ndarray
Array of undeformed cell coordinates.
"""
cellsSizeOrigin_coordinates0_point(coordinates0,ordered=True)
def coordinates0_node(cells,size,origin=_np.zeros(3)):
"""
Nodal positions (undeformed).
@ -278,6 +311,11 @@ def coordinates0_node(cells,size,origin=_np.zeros(3)):
origin : numpy.ndarray of shape (3), optional
Physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
Returns
-------
x_n_0 : numpy.ndarray
Undeformed nodal coordinates.
"""
return _np.stack(_np.meshgrid(_np.linspace(origin[0],size[0]+origin[0],cells[0]+1),
_np.linspace(origin[1],size[1]+origin[1],cells[1]+1),
@ -296,6 +334,11 @@ def displacement_fluct_node(size,F):
F : numpy.ndarray
Deformation gradient field.
Returns
-------
u_n_fluct : numpy.ndarray
Fluctuating part of the nodal displacements.
"""
return point_to_node(displacement_fluct_point(size,F))
@ -311,6 +354,11 @@ def displacement_avg_node(size,F):
F : numpy.ndarray
Deformation gradient field.
Returns
-------
u_n_avg : numpy.ndarray
Average part of the nodal displacements.
"""
F_avg = _np.average(F,axis=(0,1,2))
return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),coordinates0_node(F.shape[:3],size))
@ -327,6 +375,11 @@ def displacement_node(size,F):
F : numpy.ndarray
Deformation gradient field.
Returns
-------
u_p : numpy.ndarray
Nodal displacements.
"""
return displacement_avg_node(size,F) + displacement_fluct_node(size,F)
@ -344,28 +397,15 @@ def coordinates_node(size,F,origin=_np.zeros(3)):
origin : numpy.ndarray of shape (3), optional
Physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
Returns
-------
x_n : numpy.ndarray
Nodal coordinates.
"""
return coordinates0_node(F.shape[:3],size,origin) + displacement_node(size,F)
def point_to_node(cell_data):
"""Interpolate periodic point data to nodal data."""
n = ( cell_data + _np.roll(cell_data,1,(0,1,2))
+ _np.roll(cell_data,1,(0,)) + _np.roll(cell_data,1,(1,)) + _np.roll(cell_data,1,(2,))
+ _np.roll(cell_data,1,(0,1)) + _np.roll(cell_data,1,(1,2)) + _np.roll(cell_data,1,(2,0)))*0.125
return _np.pad(n,((0,1),(0,1),(0,1))+((0,0),)*len(cell_data.shape[3:]),mode='wrap')
def node_2_point(node_data):
"""Interpolate periodic nodal data to point data."""
c = ( node_data + _np.roll(node_data,1,(0,1,2))
+ _np.roll(node_data,1,(0,)) + _np.roll(node_data,1,(1,)) + _np.roll(node_data,1,(2,))
+ _np.roll(node_data,1,(0,1)) + _np.roll(node_data,1,(1,2)) + _np.roll(node_data,1,(2,0)))*0.125
return c[1:,1:,1:]
def cellsSizeOrigin_coordinates0_node(coordinates0,ordered=True):
"""
Return grid 'DNA', i.e. cells, size, and origin from 1D array of nodal positions.
@ -378,6 +418,11 @@ def cellsSizeOrigin_coordinates0_node(coordinates0,ordered=True):
Expect coordinates0 data to be ordered (x fast, z slow).
Defaults to True.
Returns
-------
cells, size, origin : Three numpy.ndarray, each of shape (3)
Information to reconstruct grid.
"""
coords = [_np.unique(coordinates0[:,i]) for i in range(3)]
mincorner = _np.array(list(map(min,coords)))
@ -402,6 +447,72 @@ def cellsSizeOrigin_coordinates0_node(coordinates0,ordered=True):
return (cells,size,origin)
def point_to_node(cell_data):
"""
Interpolate periodic point data to nodal data.
Parameters
----------
cell_data : numpy.ndarray of shape (:,:,:,...)
Data defined on the cell centers of a periodic grid.
Returns
-------
node_data : numpy.ndarray of shape (:,:,:,...)
Data defined on the nodes of a periodic grid.
"""
n = ( cell_data + _np.roll(cell_data,1,(0,1,2))
+ _np.roll(cell_data,1,(0,)) + _np.roll(cell_data,1,(1,)) + _np.roll(cell_data,1,(2,))
+ _np.roll(cell_data,1,(0,1)) + _np.roll(cell_data,1,(1,2)) + _np.roll(cell_data,1,(2,0)))*0.125
return _np.pad(n,((0,1),(0,1),(0,1))+((0,0),)*len(cell_data.shape[3:]),mode='wrap')
def node_to_point(node_data):
"""
Interpolate periodic nodal data to point data.
Parameters
----------
node_data : numpy.ndarray of shape (:,:,:,...)
Data defined on the nodes of a periodic grid.
Returns
-------
cell_data : numpy.ndarray of shape (:,:,:,...)
Data defined on the cell centers of a periodic grid.
"""
c = ( node_data + _np.roll(node_data,1,(0,1,2))
+ _np.roll(node_data,1,(0,)) + _np.roll(node_data,1,(1,)) + _np.roll(node_data,1,(2,))
+ _np.roll(node_data,1,(0,1)) + _np.roll(node_data,1,(1,2)) + _np.roll(node_data,1,(2,0)))*0.125
return c[1:,1:,1:]
def coordinates0_valid(coordinates0):
"""
Check whether coordinates form a regular grid.
Parameters
----------
coordinates0 : numpy.ndarray
Array of undeformed cell coordinates.
Returns
-------
valid : bool
Whether the coordinates form a regular grid.
"""
try:
cellsSizeOrigin_coordinates0_point(coordinates0,ordered=True)
return True
except ValueError:
return False
def regrid(size,F,cells):
"""
Return mapping from coordinates in deformed configuration to a regular grid.

View File

@ -1,4 +1,9 @@
"""Finite-strain continuum mechanics."""
"""
Finite-strain continuum mechanics.
All routines operate on numpy.ndarrays of shape (...,3,3).
"""
from . import tensor as _tensor
from . import _rotation
@ -154,7 +159,6 @@ def strain(F,t,m):
return eps
def stress_Cauchy(P,F):
"""
Calculate the Cauchy stress (true stress).

View File

@ -24,6 +24,11 @@ def from_random(size,N_seeds,cells=None,rng_seed=None):
A seed to initialize the BitGenerator. Defaults to None.
If None, then fresh, unpredictable entropy will be pulled from the OS.
Returns
-------
coords : numpy.ndarray of shape (N_seeds,3)
Seed coordinates in 3D space.
"""
rng = _np.random.default_rng(rng_seed)
if cells is None:
@ -56,6 +61,11 @@ def from_Poisson_disc(size,N_seeds,N_candidates,distance,periodic=True,rng_seed=
A seed to initialize the BitGenerator. Defaults to None.
If None, then fresh, unpredictable entropy will be pulled from the OS.
Returns
-------
coords : numpy.ndarray of shape (N_seeds,3)
Seed coordinates in 3D space.
"""
rng = _np.random.default_rng(rng_seed)
coords = _np.empty((N_seeds,3))
@ -94,6 +104,11 @@ def from_grid(grid,selection=None,invert=False,average=False,periodic=True):
periodic : boolean, optional
Center of gravity with periodic boundaries.
Returns
-------
coords, materials : numpy.ndarray of shape (:,3), numpy.ndarray of shape (:)
Seed coordinates in 3D space, material IDs.
"""
material = grid.material.reshape((-1,1),order='F')
mask = _np.full(grid.cells.prod(),True,dtype=bool) if selection is None else \

View File

@ -1,3 +1,3 @@
"""Tools to control the various solvers."""
"""Run simulations directly from Python."""
from ._marc import Marc # noqa

View File

@ -1,10 +1,7 @@
"""
Tensor operations.
Tensor mathematics.
Notes
-----
This is not a tensor class, but a collection of routines
to operate on numpy.ndarrays of shape (...,3,3).
All routines operate on numpy.ndarrays of shape (...,3,3).
"""

View File

@ -49,7 +49,7 @@ _colors = {
####################################################################################################
def srepr(arg,glue = '\n'):
r"""
Join arguments with glue string.
Join items with glue string.
Parameters
----------
@ -58,6 +58,11 @@ def srepr(arg,glue = '\n'):
glue : str, optional
Glue used for joining operation. Defaults to \n.
Returns
-------
joined : str
String representation of the joined items.
"""
if (not hasattr(arg, 'strip') and
(hasattr(arg, '__getitem__') or
@ -68,19 +73,71 @@ def srepr(arg,glue = '\n'):
def emph(what):
"""Formats string with emphasis."""
"""
Format with emphasis.
Parameters
----------
what : object with __repr__ or iterable of objects with __repr__.
Message to format.
Returns
-------
formatted : str
Formatted string representation of the joined items.
"""
return _colors['bold']+srepr(what)+_colors['end_color']
def deemph(what):
"""Formats string with deemphasis."""
"""
Format with deemphasis.
Parameters
----------
what : object with __repr__ or iterable of objects with __repr__.
Message to format.
Returns
-------
formatted : str
Formatted string representation of the joined items.
"""
return _colors['dim']+srepr(what)+_colors['end_color']
def warn(what):
"""Formats string for warning."""
"""
Format for warning.
Parameters
----------
what : object with __repr__ or iterable of objects with __repr__.
Message to format.
Returns
-------
formatted : str
Formatted string representation of the joined items.
"""
return _colors['warning']+emph(what)+_colors['end_color']
def strikeout(what):
"""Formats string as strikeout."""
"""
Format as strikeout.
Parameters
----------
what : object with __repr__ or iterable of objects with __repr__.
Message to format.
Returns
-------
formatted : str
Formatted string representation of the joined items.
"""
return _colors['crossout']+srepr(what)+_colors['end_color']
@ -97,6 +154,11 @@ def execute(cmd,wd='./',env=None):
env : dict, optional
Environment for execution.
Returns
-------
stdout, stderr : str
Output of the executed command.
"""
print(f"executing '{cmd}' in '{wd}'")
process = subprocess.run(shlex.split(cmd),
@ -157,6 +219,11 @@ def scale_to_coprime(v):
v : numpy.ndarray of shape (:)
Vector to scale.
Returns
-------
m : numpy.ndarray of shape (:)
Vector scaled to co-prime numbers.
"""
MAX_DENOMINATOR = 1000000
@ -371,9 +438,14 @@ def DREAM3D_base_group(fname):
Parameters
----------
fname : str
fname : str or pathlib.Path
Filename of the DREAM.3D (HDF5) file.
Returns
-------
path : str
Path to the base group.
"""
with h5py.File(fname,'r') as f:
base_group = f.visit(lambda path: path.rsplit('/',2)[0] if '_SIMPL_GEOMETRY/SPACING' in path else None)
@ -393,9 +465,14 @@ def DREAM3D_cell_data_group(fname):
Parameters
----------
fname : str
fname : str or pathlib.Path
Filename of the DREAM.3D (HDF5) file.
Returns
-------
path : str
Path to the cell data group.
"""
base_group = DREAM3D_base_group(fname)
with h5py.File(fname,'r') as f:

View File

@ -110,14 +110,14 @@ class TestResult:
def test_add_calculation(self,default,tmp_path,mode):
if mode == 'direct':
default.add_calculation('x','2.0*np.abs(#F#)-1.0','-','my notes')
default.add_calculation('2.0*np.abs(#F#)-1.0','x','-','my notes')
else:
with open(tmp_path/'f.py','w') as f:
f.write("import numpy as np\ndef my_func(field):\n return 2.0*np.abs(field)-1.0\n")
sys.path.insert(0,str(tmp_path))
import f
default.enable_user_function(f.my_func)
default.add_calculation('x','my_func(#F#)','-','my notes')
default.add_calculation('my_func(#F#)','x','-','my notes')
in_memory = 2.0*np.abs(default.place('F'))-1.0
in_file = default.place('x')
@ -193,14 +193,14 @@ class TestResult:
def test_add_Mises_invalid(self,default):
default.add_stress_Cauchy('P','F')
default.add_calculation('sigma_y','#sigma#',unit='y')
default.add_calculation('#sigma#','sigma_y',unit='y')
default.add_equivalent_Mises('sigma_y')
assert default.get('sigma_y_vM') is None
def test_add_Mises_stress_strain(self,default):
default.add_stress_Cauchy('P','F')
default.add_calculation('sigma_y','#sigma#',unit='y')
default.add_calculation('sigma_x','#sigma#',unit='x')
default.add_calculation('#sigma#','sigma_y',unit='y')
default.add_calculation('#sigma#','sigma_x',unit='x')
default.add_equivalent_Mises('sigma_y',kind='strain')
default.add_equivalent_Mises('sigma_x',kind='stress')
assert not np.allclose(default.place('sigma_y_vM'),default.place('sigma_x_vM'))
@ -271,7 +271,7 @@ class TestResult:
@pytest.mark.parametrize('overwrite',['off','on'])
def test_add_overwrite(self,default,overwrite):
last = default.view('times',default.times_in_range(0,np.inf)[-1])
last = default.view('increments',-1)
last.add_stress_Cauchy()
@ -279,13 +279,13 @@ class TestResult:
created_first = datetime.strptime(created_first,'%Y-%m-%d %H:%M:%S%z')
if overwrite == 'on':
last = last.allow_modification()
last = last.modification_enable()
else:
last = last.disallow_modification()
last = last.modification_disable()
time.sleep(2.)
try:
last.add_calculation('sigma','#sigma#*0.0+311.','not the Cauchy stress')
last.add_calculation('#sigma#*0.0+311.','sigma','not the Cauchy stress')
except ValueError:
pass
@ -301,14 +301,24 @@ class TestResult:
def test_rename(self,default,allowed):
if allowed == 'on':
F = default.place('F')
default = default.allow_modification()
default = default.modification_enable()
default.rename('F','new_name')
assert np.all(F == default.place('new_name'))
default = default.disallow_modification()
default = default.modification_disable()
with pytest.raises(PermissionError):
default.rename('P','another_new_name')
@pytest.mark.parametrize('allowed',['off','on'])
def test_remove(self,default,allowed):
if allowed == 'on':
unsafe = default.modification_enable()
unsafe.remove('F')
assert unsafe.get('F') is None
else:
with pytest.raises(PermissionError):
default.remove('F')
@pytest.mark.parametrize('mode',['cell','node'])
def test_coordinates(self,default,mode):
if mode == 'cell':
@ -352,7 +362,7 @@ class TestResult:
def test_XDMF(self,tmp_path,single_phase,update,ref_path):
for shape in [('scalar',()),('vector',(3,)),('tensor',(3,3)),('matrix',(12,))]:
for dtype in ['f4','f8','i1','i2','i4','i8','u1','u2','u4','u8']:
single_phase.add_calculation(f'{shape[0]}_{dtype}',f"np.ones(np.shape(#F#)[0:1]+{shape[1]},'{dtype}')")
single_phase.add_calculation(f"np.ones(np.shape(#F#)[0:1]+{shape[1]},'{dtype}')",f'{shape[0]}_{dtype}')
fname = os.path.splitext(os.path.basename(single_phase.fname))[0]+'.xdmf'
os.chdir(tmp_path)
single_phase.save_XDMF()

View File

@ -60,7 +60,7 @@ class TestGridFilters:
cell_field_x = np.interp(coordinates0_point_x,coordinates_node_x,node_field_x,period=np.pi*2.)
cell_field = np.broadcast_to(cell_field_x.reshape(-1,1,1),cells)
assert np.allclose(cell_field,grid_filters.node_2_point(node_field))
assert np.allclose(cell_field,grid_filters.node_to_point(node_field))
@pytest.mark.parametrize('mode',['point','node'])
def test_coordinates0_origin(self,mode):
@ -93,14 +93,24 @@ class TestGridFilters:
F = np.broadcast_to(np.random.random((3,3)), tuple(cells)+(3,3))
assert np.allclose(function(size,F),0.0)
@pytest.mark.parametrize('function',[grid_filters.coordinates0_check,
grid_filters.cellsSizeOrigin_coordinates0_node,
grid_filters.cellsSizeOrigin_coordinates0_point])
@pytest.mark.parametrize('function',[grid_filters.cellsSizeOrigin_coordinates0_point,
grid_filters.cellsSizeOrigin_coordinates0_node])
def test_invalid_coordinates(self,function):
invalid_coordinates = np.random.random((np.random.randint(12,52),3))
with pytest.raises(ValueError):
function(invalid_coordinates)
@pytest.mark.parametrize('function',[grid_filters.coordinates0_point,
grid_filters.coordinates0_node])
def test_valid_coordinates_check(self,function):
valid_coordinates = function(np.random.randint(4,10,(3)),np.random.rand(3))
assert grid_filters.coordinates0_valid(valid_coordinates.reshape(-1,3,order='F'))
def test_invalid_coordinates_check(self):
invalid_coordinates = np.random.random((np.random.randint(12,52),3))
assert not grid_filters.coordinates0_valid(invalid_coordinates)
@pytest.mark.parametrize('function',[grid_filters.cellsSizeOrigin_coordinates0_node,
grid_filters.cellsSizeOrigin_coordinates0_point])
def test_uneven_spaced_coordinates(self,function):