diff --git a/.gitattributes b/.gitattributes index c2785cf15..1437f7412 100644 --- a/.gitattributes +++ b/.gitattributes @@ -9,10 +9,11 @@ *.hdf5 binary *.pdf binary *.dream3d binary +*.pbz2 binary # ignore files from MSC.Marc in language statistics installation/mods_MarcMentat/20*/* linguist-vendored -src/marc/include/* linguist-vendored +src/Marc/include/* linguist-vendored # ignore reference files for tests in language statistics python/tests/reference/* linguist-vendored diff --git a/PRIVATE b/PRIVATE index f1ac733d5..129812414 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit f1ac733d5eb90de1cdcaf79a261157c9034f8136 +Subproject commit 1298124143e7e2901d0b9c2e79ab6388cb78a1e3 diff --git a/VERSION b/VERSION index e5c6789ca..6855a7dca 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v3.0.0-alpha2-753-g565dab120 +v3.0.0-alpha2-829-g73b07eda4 diff --git a/python/damask/_configmaterial.py b/python/damask/_configmaterial.py index dedf4373b..557907594 100644 --- a/python/damask/_configmaterial.py +++ b/python/damask/_configmaterial.py @@ -1,5 +1,3 @@ -import os - import numpy as np import h5py @@ -159,18 +157,18 @@ class ConfigMaterial(Config): f = h5py.File(fname,'r') if grain_data is None: - phase = f[os.path.join(b,c,phases)][()].flatten() - O = Rotation.from_Euler_angles(f[os.path.join(b,c,Euler_angles)]).as_quaternion().reshape(-1,4) # noqa + phase = f['/'.join([b,c,phases])][()].flatten() + O = Rotation.from_Euler_angles(f['/'.join([b,c,Euler_angles])]).as_quaternion().reshape(-1,4) # noqa _,idx = np.unique(np.hstack([O,phase.reshape(-1,1)]),return_index=True,axis=0) idx = np.sort(idx) else: - phase = f[os.path.join(b,grain_data,phases)][()] - O = Rotation.from_Euler_angles(f[os.path.join(b,grain_data,Euler_angles)]).as_quaternion() # noqa + phase = f['/'.join([b,grain_data,phases])][()] + O = Rotation.from_Euler_angles(f['/'.join([b,grain_data,Euler_angles])]).as_quaternion() # noqa idx = np.arange(phase.size) if cell_ensemble_data is not None and phase_names is not None: try: - names = np.array([s.decode() for s in f[os.path.join(b,cell_ensemble_data,phase_names)]]) + names = np.array([s.decode() for s in f['/'.join([b,cell_ensemble_data,phase_names])]]) phase = names[phase] except KeyError: pass diff --git a/python/damask/_grid.py b/python/damask/_grid.py index ed44989d7..308be5c80 100644 --- a/python/damask/_grid.py +++ b/python/damask/_grid.py @@ -69,9 +69,9 @@ class Grid: copy = __copy__ - def diff(self,other): + def __eq__(self,other): """ - Report property differences of self relative to other. + Test equality of other. Parameters ---------- @@ -79,28 +79,10 @@ class Grid: Grid to compare self against. """ - message = [] - if np.any(other.cells != self.cells): - message.append(util.deemph(f'cells a b c: {util.srepr(other.cells," x ")}')) - message.append(util.emph( f'cells a b c: {util.srepr( self.cells," x ")}')) - - if not np.allclose(other.size,self.size): - message.append(util.deemph(f'size x y z: {util.srepr(other.size," x ")}')) - message.append(util.emph( f'size x y z: {util.srepr( self.size," x ")}')) - - if not np.allclose(other.origin,self.origin): - message.append(util.deemph(f'origin x y z: {util.srepr(other.origin," ")}')) - message.append(util.emph( f'origin x y z: {util.srepr( self.origin," ")}')) - - if other.N_materials != self.N_materials: - message.append(util.deemph(f'# materials: {other.N_materials}')) - message.append(util.emph( f'# materials: { self.N_materials}')) - - if np.nanmax(other.material) != np.nanmax(self.material): - message.append(util.deemph(f'max material: {np.nanmax(other.material)}')) - message.append(util.emph( f'max material: {np.nanmax( self.material)}')) - - return util.return_message(message) + return (np.allclose(other.size,self.size) + and np.allclose(other.origin,self.origin) + and np.all(other.cells == self.cells) + and np.all(other.material == self.material)) @property @@ -270,7 +252,7 @@ class Grid: """ Load DREAM.3D (HDF5) file. - Data in DREAM.3D files can be stored per cell ('CellData') and/or + Data in DREAM.3D files can be stored per cell ('CellData') and/or per grain ('Grain Data'). Per default, cell-wise data is assumed. damask.ConfigMaterial.load_DREAM3D gives the corresponding material definition. @@ -305,18 +287,18 @@ class Grid: c = util.DREAM3D_cell_data_group(fname) if cell_data is None else cell_data f = h5py.File(fname, 'r') - cells = f[os.path.join(b,'_SIMPL_GEOMETRY','DIMENSIONS')][()] - size = f[os.path.join(b,'_SIMPL_GEOMETRY','SPACING')] * cells - origin = f[os.path.join(b,'_SIMPL_GEOMETRY','ORIGIN')][()] + cells = f['/'.join([b,'_SIMPL_GEOMETRY','DIMENSIONS'])][()] + size = f['/'.join([b,'_SIMPL_GEOMETRY','SPACING'])] * cells + origin = f['/'.join([b,'_SIMPL_GEOMETRY','ORIGIN'])][()] if feature_IDs is None: - phase = f[os.path.join(b,c,phases)][()].reshape(-1,1) - O = Rotation.from_Euler_angles(f[os.path.join(b,c,Euler_angles)]).as_quaternion().reshape(-1,4) # noqa + phase = f['/'.join([b,c,phases])][()].reshape(-1,1) + O = Rotation.from_Euler_angles(f['/'.join([b,c,Euler_angles])]).as_quaternion().reshape(-1,4) # noqa unique,unique_inverse = np.unique(np.hstack([O,phase]),return_inverse=True,axis=0) ma = np.arange(cells.prod()) if len(unique) == cells.prod() else \ np.arange(unique.size)[np.argsort(pd.unique(unique_inverse))][unique_inverse] else: - ma = f[os.path.join(b,c,feature_IDs)][()].flatten() + ma = f['/'.join([b,c,feature_IDs])][()].flatten() return Grid(ma.reshape(cells,order='F'),size,origin,util.execution_stamp('Grid','load_DREAM3D')) diff --git a/python/damask/_result.py b/python/damask/_result.py index b2a0dd4a1..89c704b94 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -1,21 +1,23 @@ import multiprocessing as mp import re -import glob +import fnmatch import os +import copy import datetime import xml.etree.ElementTree as ET import xml.dom.minidom from pathlib import Path from functools import partial from collections import defaultdict +from collections.abc import Iterable import h5py import numpy as np +import numpy.ma as ma from numpy.lib import recfunctions as rfn import damask from . import VTK -from . import Table from . import Orientation from . import grid_filters from . import mechanics @@ -24,11 +26,42 @@ from . import util h5py3 = h5py.__version__[0] == '3' + +def _read(dataset): + """Read a dataset and its metadata into a numpy.ndarray.""" + metadata = {k:(v.decode() if not h5py3 and type(v) is bytes else v) for k,v in dataset.attrs.items()} + dtype = np.dtype(dataset.dtype,metadata=metadata) + return np.array(dataset,dtype=dtype) + +def _match(requested,existing): + """Find matches among two sets of labels.""" + def flatten_list(list_of_lists): + return [e for e_ in list_of_lists for e in e_] + + if requested is True: + requested = '*' + elif requested is False or requested is None: + requested = [] + + requested_ = requested if hasattr(requested,'__iter__') and not isinstance(requested,str) else \ + [requested] + + return sorted(set(flatten_list([fnmatch.filter(existing,r) for r in requested_])), + key=util.natural_sort) + +def _empty_like(dataset,N_materialpoints,fill_float,fill_int): + """Create empty numpy.ma.MaskedArray.""" + return ma.array(np.empty((N_materialpoints,)+dataset.shape[1:],dataset.dtype), + fill_value = fill_float if dataset.dtype in np.sctypes['float'] else fill_int, + mask = True) + class Result: """ Manipulate and read DADF5 files. DADF5 (DAMASK HDF5) files contain DAMASK results. + The group/folder structure reflects the input data + in material.yaml. """ def __init__(self,fname): @@ -61,10 +94,9 @@ class Result: self.origin = f['geometry'].attrs['origin'] r=re.compile('inc[0-9]+' if self.version_minor < 12 else 'increment_[0-9]+') - increments_unsorted = {int(i[10:]):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] if self.version_minor < 12 else \ - [round(f[i].attrs['t/s'],12) for i in self.increments] + self.increments = sorted([i for i in f.keys() if r.match(i)],key=util.natural_sort) + self.times = [round(f[i].attrs['time/s' if self.version_minor < 12 else + 't/s'],12) for i in self.increments] grp = 'mapping' if self.version_minor < 12 else 'cell_to' @@ -72,24 +104,22 @@ class Result: self.homogenizations = [m.decode() for m in np.unique(f[f'{grp}/homogenization'] ['Name' if self.version_minor < 12 else 'label'])] + self.homogenizations = sorted(self.homogenizations,key=util.natural_sort) self.phases = [c.decode() for c in np.unique(f[f'{grp}/phase'] ['Name' if self.version_minor < 12 else 'label'])] + self.phases = sorted(self.phases,key=util.natural_sort) - self.out_type_ph = [] + self.fields = [] for c in self.phases: - self.out_type_ph += f['/'.join([self.increments[0],'phase',c])].keys() - self.out_type_ph = list(set(self.out_type_ph)) # make unique - - self.out_type_ho = [] + self.fields += f['/'.join([self.increments[0],'phase',c])].keys() for m in self.homogenizations: - self.out_type_ho += f['/'.join([self.increments[0],'homogenization',m])].keys() - self.out_type_ho = list(set(self.out_type_ho)) # make unique + self.fields += f['/'.join([self.increments[0],'homogenization',m])].keys() + self.fields = sorted(set(self.fields),key=util.natural_sort) # make unique self.visible = {'increments': self.increments, 'phases': self.phases, 'homogenizations': self.homogenizations, - 'out_type_ph': self.out_type_ph, - 'out_type_ho': self.out_type_ho + 'fields': self.fields, } self.fname = Path(fname).absolute() @@ -97,17 +127,21 @@ class Result: self._allow_modification = False + def __copy__(self): + """Create deep copy.""" + return copy.deepcopy(self) + + copy = __copy__ + + def __repr__(self): """Show summary of file content.""" visible_increments = self.visible['increments'] - self.view('increments',visible_increments[0:1]) - first = self.list_data() + first = self.view('increments',visible_increments[0:1]).list_data() - self.view('increments',visible_increments[-1:]) - last = '' if len(visible_increments) < 2 else self.list_data() - - self.view('increments',visible_increments) + last = '' if len(visible_increments) < 2 else \ + self.view('increments',visible_increments[-1:]).list_data() in_between = '' if len(visible_increments) < 3 else \ ''.join([f'\n{inc}\n ...\n' for inc in visible_increments[1:-1]]) @@ -125,19 +159,15 @@ class Result: Select from 'set', 'add', and 'del'. what : str Attribute to change (must be from self.visible). - datasets : str, int, list of str, list of int, or bool - Name of datasets as list; supports ? and * wildcards. - True is equivalent to [*], False is equivalent to []. + datasets : (list of) int (for increments), (list of) float (for times), (list of) str, or bool + Name of datasets; supports '?' and '*' wildcards. + True is equivalent to '*', False is equivalent to []. """ - def natural_sort(key): - convert = lambda text: int(text) if text.isdigit() else text - return [ convert(c) for c in re.split('([0-9]+)', key) ] - # allow True/False and string arguments if datasets is True: - datasets = ['*'] - elif datasets is False: + datasets = '*' + elif datasets is False or datasets is None: datasets = [] choice = list(datasets).copy() if hasattr(datasets,'__iter__') and not isinstance(datasets,str) else \ [datasets] @@ -161,53 +191,34 @@ class Result: elif np.isclose(c,self.times[idx+1]): choice.append(self.increments[idx+1]) - valid = [e for e_ in [glob.fnmatch.filter(getattr(self,what),s) for s in choice] for e in e_] + valid = _match(choice,getattr(self,what)) existing = set(self.visible[what]) + dup = self.copy() if action == 'set': - self.visible[what] = valid + dup.visible[what] = sorted(set(valid), key=util.natural_sort) elif action == 'add': add = existing.union(valid) - add_sorted = sorted(add, key=natural_sort) - self.visible[what] = add_sorted + dup.visible[what] = sorted(add, key=util.natural_sort) elif action == 'del': diff = existing.difference(valid) - diff_sorted = sorted(diff, key=natural_sort) - self.visible[what] = diff_sorted + dup.visible[what] = sorted(diff, key=util.natural_sort) - - def _get_attribute(self,path,attr): - """ - Get the attribute of a dataset. - - Parameters - ---------- - Path : str - Path to the dataset. - attr : str - Name of the attribute to get. - - Returns - ------- - attr at path, str or None. - The requested attribute, None if not found. - - """ - with h5py.File(self.fname,'r') as f: - try: - return f[path].attrs[attr] if h5py3 else f[path].attrs[attr].decode() - except KeyError: - return None + return dup def allow_modification(self): """Allow to overwrite existing data.""" print(util.warn('Warning: Modification of existing datasets allowed!')) - self._allow_modification = True + dup = self.copy() + dup._allow_modification = True + return dup def disallow_modification(self): """Disallow to overwrite existing data (default case).""" - self._allow_modification = False + dup = self.copy() + dup._allow_modification = False + return dup def increments_in_range(self,start,end): @@ -251,42 +262,20 @@ class Result: return selected - def iterate(self,what): - """ - Iterate over visible items and view them independently. - - Parameters - ---------- - what : str - Attribute to change (must be from self.visible). - - """ - datasets = self.visible[what] - last_view = datasets.copy() - for dataset in datasets: - if last_view != self.visible[what]: - self._manage_view('set',what,datasets) - raise Exception - self._manage_view('set',what,dataset) - last_view = self.visible[what] - yield dataset - self._manage_view('set',what,datasets) - - def view(self,what,datasets): """ Set view. Parameters ---------- - what : str - Attribute to change (must be from self.visible). - datasets : list of str or bool - Name of datasets as list; supports ? and * wildcards. - True is equivalent to [*], False is equivalent to []. + what : {'increments', 'times', 'phases', 'homogenizations', 'fields'} + Attribute to change. + datasets : (list of) int (for increments), (list of) float (for times), (list of) str, or bool + Name of datasets; supports '?' and '*' wildcards. + True is equivalent to '*', False is equivalent to []. """ - self._manage_view('set',what,datasets) + return self._manage_view('set',what,datasets) def view_more(self,what,datasets): @@ -295,14 +284,14 @@ class Result: Parameters ---------- - what : str - Attribute to change (must be from self.visible). - datasets : list of str or bool - Name of datasets as list; supports ? and * wildcards. - True is equivalent to [*], False is equivalent to []. + what : {'increments', 'times', 'phases', 'homogenizations', 'fields'} + Attribute to change. + datasets : (list of) int (for increments), (list of) float (for times), (list of) str, or bool + Name of datasets; supports '?' and '*' wildcards. + True is equivalent to '*', False is equivalent to []. """ - self._manage_view('add',what,datasets) + return self._manage_view('add',what,datasets) def view_less(self,what,datasets): @@ -311,14 +300,14 @@ class Result: Parameters ---------- - what : str - Attribute to change (must be from self.visible). - datasets : list of str or bool - Name of datasets as list; supports ? and * wildcards. - True is equivalent to [*], False is equivalent to []. + what : {'increments', 'times', 'phases', 'homogenizations', 'fields'} + Attribute to change. + datasets : (list of) int (for increments), (list of) float (for times), (list of) str, or bool + Name of datasets; supports '?' and '*' wildcards. + True is equivalent to '*', False is equivalent to []. """ - self._manage_view('del',what,datasets) + return self._manage_view('del',what,datasets) def rename(self,name_old,name_new): @@ -333,137 +322,21 @@ class Result: New name of the dataset. """ - 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'] = f'Original name: {name_old}' if h5py3 else \ - f'Original name: {name_old}'.encode() - del f[path_old] - else: + if not self._allow_modification: raise PermissionError('Rename operation not permitted') - - def place(self,datasets,constituent=0,tagged=False,split=True): - """ - Distribute datasets onto geometry and return Table or (split) dictionary of Tables. - - Must not mix nodal end cell data. - - Only data within - - inc*/phase/*/* - - inc*/homogenization/*/* - - inc*/geometry/* - are considered. - - Parameters - ---------- - datasets : iterable or str - constituent : int - Constituent to consider for phase data. - tagged : bool - Tag Table.column name with '#constituent'. - Defaults to False. - split : bool - Split Table by increment and return dictionary of Tables. - Defaults to True. - - """ - sets = datasets if hasattr(datasets,'__iter__') and not isinstance(datasets,str) else \ - [datasets] - tag = f'#{constituent}' if tagged else '' - tbl = {} if split else None - inGeom = {} - inData = {} - # compatibility hack - name = 'Name' if self.version_minor < 12 else 'label' - member = 'Position' if self.version_minor < 12 else 'entry' - grp = 'mapping' if self.version_minor < 12 else 'cell_to' - with h5py.File(self.fname,'r') as f: - for dataset in sets: - for group in self.groups_with_datasets(dataset): - 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.N_materialpoints) - elif prop == 'phase': - inGeom[key] = np.where(f[f'{grp}/phase'][:,constituent][name] == str.encode(name))[0] - inData[key] = f[f'{grp}/phase'][inGeom[key],constituent][member] - elif prop == 'homogenization': - inGeom[key] = np.where(f[f'{grp}/homogenization'][name] == str.encode(name))[0] - inData[key] = f[f'{grp}/homogenization'][inGeom[key].tolist()][member] - shape = np.shape(f[path]) - data = np.full((self.N_materialpoints,) + (shape[1:] if len(shape)>1 else (1,)), - np.nan, - dtype=np.dtype(f[path])) - data[inGeom[key]] = (f[path] if len(shape)>1 else np.expand_dims(f[path],1))[inData[key]] - path = ('/'.join([prop,name]+([cat] if cat else [])+([item] if item else [])) if split else path)+tag - if split: - try: - tbl[inc] = tbl[inc].add(path,data) - except KeyError: - tbl[inc] = Table(data.reshape(self.N_materialpoints,-1),{path:data.shape[1:]}) - else: - try: - tbl = tbl.add(path,data) - except AttributeError: - tbl = Table(data.reshape(self.N_materialpoints,-1),{path:data.shape[1:]}) - - return tbl - - - def groups_with_datasets(self,datasets): - """ - Return groups that contain all requested datasets. - - Only groups within - - inc*/phase/*/* - - inc*/homogenization/*/* - - inc*/geometry/* - - are considered as they contain user-relevant data. - Single strings will be treated as list with one entry. - - Wild card matching is allowed, but the number of arguments needs to fit. - - Parameters - ---------- - datasets : iterable or str or bool - - Examples - -------- - 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'] - - """ - if datasets is False: return [] - - sets = datasets if isinstance(datasets,bool) or (hasattr(datasets,'__iter__') and not isinstance(datasets,str)) else \ - [datasets] - - groups = [] - - with h5py.File(self.fname,'r') as f: - for i in self.iterate('increments'): - for o,p in zip(['phases','homogenizations'],['out_type_ph','out_type_ho']): - 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: - if group in f.keys(): - match = [e for e_ in [glob.fnmatch.filter(f[group].keys(),s) for s in sets] for e in e_] - if len(set(match)) == len(sets): groups.append(group) - return groups + 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 self.visible['fields']: + 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] def list_data(self): @@ -471,54 +344,25 @@ class Result: # compatibility hack de = 'Description' if self.version_minor < 12 else 'description' un = 'Unit' if self.version_minor < 12 else 'unit' - message = '' + msg = '' with h5py.File(self.fname,'r') as f: - for i in self.iterate('increments'): - message += f'\n{i} ({self.times[self.increments.index(i)]}s)\n' - for o,p in zip(['phases','homogenizations'],['out_type_ph','out_type_ho']): - message += f' {o[:-1]}\n' - for oo in self.iterate(o): - message += f' {oo}\n' - for pp in self.iterate(p): - message += f' {pp}\n' - group = '/'.join([i,o[:-1],oo,pp]) # o[:-1]: plural/singular issue - for d in f[group].keys(): - try: - dataset = f['/'.join([group,d])] - if un in dataset.attrs: - unit = f" / {dataset.attrs[un]}" if h5py3 else \ - f" / {dataset.attrs[un].decode()}" - else: - unit = '' - description = dataset.attrs[de] if h5py3 else \ - dataset.attrs[de].decode() - message += f' {d}{unit}: {description}\n' - except KeyError: - pass - return message + for inc in self.visible['increments']: + msg = ''.join([msg,f'\n{inc} ({self.times[self.increments.index(inc)]}s)\n']) + for ty in ['phase','homogenization']: + msg = ' '.join([msg,f'{ty}\n']) + for label in self.visible[ty+'s']: + msg = ' '.join([msg,f'{label}\n']) + for field in self.visible['fields']: + msg = ' '.join([msg,f'{field}\n']) + for d in f['/'.join([inc,ty,label,field])].keys(): + dataset = f['/'.join([inc,ty,label,field,d])] + unit = f' / {dataset.attrs[un]}' if h5py3 else \ + f' / {dataset.attrs[un].decode()}' + description = dataset.attrs[de] if h5py3 else \ + dataset.attrs[de].decode() + msg = ' '.join([msg,f'{d}{unit}: {description}\n']) - - 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: - for i in self.iterate('increments'): - k = '/'.join([i,'geometry',label]) - try: - f[k] - path.append(k) - except KeyError: - pass - for o,p in zip(['phases','homogenizations'],['out_type_ph','out_type_ho']): - for oo in self.iterate(o): - for pp in self.iterate(p): - k = '/'.join([i,o[:-1],oo,pp,label]) - try: - f[k] - path.append(k) - except KeyError: - pass - return path + return msg def enable_user_function(self,func): @@ -526,60 +370,6 @@ class Result: print(f'Function {func.__name__} enabled in add_calculation.') - def read_dataset(self,path,c=0,plain=False): - """ - Dataset for all points/cells. - - If more than one path is given, the dataset is composed of the individual contributions. - - Parameters - ---------- - path : list of strings - The name of the datasets to consider. - c : int, optional - The constituent to consider. Defaults to 0. - plain: boolean, optional - Convert into plain numpy datatype. - Only relevant for compound datatype, e.g. the orientation. - Defaults to False. - - """ - # compatibility hack - name = 'Name' if self.version_minor < 12 else 'label' - member = 'Position' if self.version_minor < 12 else 'entry' - grp = 'mapping' if self.version_minor < 12 else 'cell_to' - with h5py.File(self.fname,'r') as f: - shape = (self.N_materialpoints,) + 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] - - if pa.split('/')[1] == 'geometry': - dataset = np.array(f[pa]) - continue - - p = np.where(f[f'{grp}/phase'][:,c][name] == str.encode(label))[0] - if len(p)>0: - u = (f[f'{grp}/phase'][member][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[f'{grp}/homogenization'][name] == str.encode(label))[0] - if len(p)>0: - u = (f[f'{grp}/homogenization'][member][p.tolist()]) - a = np.array(f[pa]) - if len(a.shape) == 1: - a=a.reshape([a.shape[0],1]) - dataset[p,:] = a[u,:] - - if plain and dataset.dtype.names is not None: - return dataset.view(('float64',len(dataset.dtype.names))) - else: - return dataset - @property def coordinates0_point(self): """Return initial coordinates of the cell centers.""" @@ -598,6 +388,17 @@ class Result: with h5py.File(self.fname,'r') as f: return f['geometry/x_n'][()] + @property + def geometry0(self): + if self.structured: + return VTK.from_rectilinear_grid(self.cells,self.size,self.origin) + else: + with h5py.File(self.fname,'r') as f: + return VTK.from_unstructured_grid(f['/geometry/x_n'][()], + f['/geometry/T_c'][()]-1, + f['/geometry/T_c'].attrs['VTK_TYPE'] if h5py3 else \ + f['/geometry/T_c'].attrs['VTK_TYPE'].decode()) + @staticmethod def _add_absolute(x): @@ -647,7 +448,7 @@ class Result: 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 '#TheirLabel#'. unit : str, optional Physical unit of the result. description : str, optional @@ -679,9 +480,9 @@ class Result: Parameters ---------- P : str, optional - Label of the dataset containing the first Piola-Kirchhoff stress. Defaults to ‘P’. + Label 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’. + Label of the dataset containing the deformation gradient. Defaults to 'F'. """ self._add_generic_pointwise(self._add_stress_Cauchy,{'P':P,'F':F}) @@ -762,7 +563,7 @@ class Result: T_sym : str Label of symmetric tensor dataset. eigenvalue : str, optional - Eigenvalue. Select from ‘max’, ‘mid’, ‘min’. Defaults to ‘max’. + Eigenvalue. Select from 'max', 'mid', 'min'. Defaults to 'max'. """ self._add_generic_pointwise(self._add_eigenvalue,{'T_sym':T_sym},{'eigenvalue':eigenvalue}) @@ -795,8 +596,8 @@ class Result: T_sym : str Label of symmetric tensor dataset. eigenvalue : str, optional - Eigenvalue to which the eigenvector corresponds. Select from - ‘max’, ‘mid’, ‘min’. Defaults to ‘max’. + Eigenvalue to which the eigenvector corresponds. + Select from 'max', 'mid', 'min'. Defaults to 'max'. """ self._add_generic_pointwise(self._add_eigenvector,{'T_sym':T_sym},{'eigenvalue':eigenvalue}) @@ -816,11 +617,11 @@ class Result: return { 'data': np.uint8(o.IPF_color(l)*255), - 'label': 'IPFcolor_[{} {} {}]'.format(*m), + 'label': 'IPFcolor_({} {} {})'.format(*m), 'meta' : { 'unit': '8-bit RGB', 'lattice': q['meta']['lattice'], - 'description': 'Inverse Pole Figure (IPF) colors along sample direction [{} {} {}]'.format(*m), + 'description': 'Inverse Pole Figure (IPF) colors along sample direction ({} {} {})'.format(*m), 'creator': 'add_IPF_color' } } @@ -873,7 +674,7 @@ class Result: elif T_sym['meta']['unit'] == 'Pa': k = 'stress' if k not in ['stress', 'strain']: - raise ValueError('invalid von Mises kind {kind}') + raise ValueError(f'invalid von Mises kind {kind}') return { 'data': (mechanics.equivalent_strain_Mises if k=='strain' else \ @@ -895,7 +696,7 @@ class Result: Label of symmetric tensorial stress or strain dataset. kind : {'stress', 'strain', None}, optional Kind of the von Mises equivalent. Defaults to None, in which case - it is selected based on the unit of the dataset ('1' -> strain, 'Pa' -> stress'). + it is selected based on the unit of the dataset ('1' -> strain, 'Pa' -> stress). """ self._add_generic_pointwise(self._add_equivalent_Mises,{'T_sym':T_sym},{'kind':kind}) @@ -932,7 +733,7 @@ class Result: ---------- x : str Label of vector or tensor dataset. - ord : {non-zero int, inf, -inf, ‘fro’, ‘nuc’}, optional + ord : {non-zero int, inf, -inf, 'fro', 'nuc'}, optional Order of the norm. inf means NumPy’s inf object. For details refer to numpy.linalg.norm. """ @@ -959,9 +760,9 @@ class Result: Parameters ---------- P : str, optional - Label of first Piola-Kirchhoff stress dataset. Defaults to ‘P’. + Label of first Piola-Kirchhoff stress dataset. Defaults to 'P'. F : str, optional - Label of deformation gradient dataset. Defaults to ‘F’. + Label of deformation gradient dataset. Defaults to 'F'. """ self._add_generic_pointwise(self._add_stress_second_Piola_Kirchhoff,{'P':P,'F':F}) @@ -1073,17 +874,17 @@ class Result: """ Add strain tensor of a deformation gradient. - For details refer to damask.mechanics.strain + For details, see damask.mechanics.strain. Parameters ---------- F : str, optional - Label of deformation gradient dataset. Defaults to ‘F’. - t : {‘V’, ‘U’}, optional - Type of the polar decomposition, ‘V’ for left stretch tensor and ‘U’ for right stretch tensor. - Defaults to ‘V’. + Label of deformation gradient dataset. Defaults to 'F'. + t : {'V', 'U'}, optional + Type of the polar decomposition, 'V' for left stretch tensor and 'U' for right stretch tensor. + Defaults to 'V'. m : float, optional - Order of the strain calculation. Defaults to ‘0.0’. + Order of the strain calculation. Defaults to 0.0. """ self._add_generic_pointwise(self._add_strain,{'F':F},{'t':t,'m':m}) @@ -1108,10 +909,10 @@ class Result: Parameters ---------- F : str, optional - Label of deformation gradient dataset. Defaults to ‘F’. - t : {‘V’, ‘U’}, optional - Type of the polar decomposition, ‘V’ for left stretch tensor and ‘U’ for right stretch tensor. - Defaults to ‘V’. + Label of deformation gradient dataset. Defaults to 'F'. + t : {'V', 'U'}, optional + Type of the polar decomposition, 'V' for left stretch tensor and 'U' for right stretch tensor. + Defaults to 'V'. """ self._add_generic_pointwise(self._add_stretch_tensor,{'F':F},{'t':t}) @@ -1146,8 +947,8 @@ class Result: Callback function that calculates a new dataset from one or more datasets per HDF5 group. datasets : dictionary - Details of the datasets to be used: label (in HDF5 file) and - arg (argument to which the data is parsed in func). + Details of the datasets to be used: + {arg (name to which the data is passed in func): label (in HDF5 file)}. args : dictionary, optional Arguments parsed to func. @@ -1156,7 +957,15 @@ class Result: pool = mp.Pool(int(os.environ.get('OMP_NUM_THREADS',1))) lock = mp.Manager().Lock() - groups = self.groups_with_datasets(datasets.values()) + groups = [] + with h5py.File(self.fname,'r') as f: + for inc in self.visible['increments']: + for ty in ['phase','homogenization']: + for label in self.visible[ty+'s']: + for field in self.visible['fields']: + group = '/'.join([inc,ty,label,field]) + if set(datasets.values()).issubset(f[group].keys()): groups.append(group) + if len(groups) == 0: print('No matching dataset found, no data was added.') return @@ -1192,8 +1001,8 @@ class Result: dataset.attrs[l.lower()]=v if h5py3 else v.encode() creator = dataset.attrs['creator'] if h5py3 else \ dataset.attrs['creator'].decode() - dataset.attrs['creator'] = f"damask.Result.{creator} v{damask.version}" if h5py3 else \ - f"damask.Result.{creator} v{damask.version}".encode() + dataset.attrs['creator'] = f'damask.Result.{creator} v{damask.version}' if h5py3 else \ + f'damask.Result.{creator} v{damask.version}'.encode() except (OSError,RuntimeError) as err: print(f'Could not add dataset: {err}.') @@ -1203,17 +1012,20 @@ class Result: pool.join() - def save_XDMF(self): + def save_XDMF(self,output='*'): """ Write XDMF file to directly visualize data in DADF5 file. - The view is not taken into account, i.e. the content of the - whole file will be included. + Parameters + ---------- + output : (list of) str + Labels of the datasets to read. + Defaults to '*', in which case all datasets are considered. + """ - # compatibility hack - u = 'Unit' if self.version_minor < 12 else 'unit' + u = 'Unit' if self.version_minor < 12 else 'unit' # compatibility hack if self.N_constituents != 1 or len(self.phases) != 1 or not self.structured: - raise TypeError('XDMF output requires homogeneous grid') + raise TypeError('XDMF output requires structured grid with single phase and single constituent.') attribute_type_map = defaultdict(lambda:'Matrix', ( ((),'Scalar'), ((3,),'Vector'), ((3,3),'Tensor')) ) @@ -1224,11 +1036,11 @@ class Result: if dtype in np.sctypes['float']: return 'Float' - xdmf=ET.Element('Xdmf') + xdmf = ET.Element('Xdmf') xdmf.attrib={'Version': '2.0', 'xmlns:xi': 'http://www.w3.org/2001/XInclude'} - domain=ET.SubElement(xdmf, 'Domain') + domain = ET.SubElement(xdmf, 'Domain') collection = ET.SubElement(domain, 'Grid') collection.attrib={'GridType': 'Collection', @@ -1238,150 +1050,306 @@ class Result: time.attrib={'TimeType': 'List'} time_data = ET.SubElement(time, 'DataItem') + times = [self.times[self.increments.index(i)] for i in self.visible['increments']] time_data.attrib={'Format': 'XML', 'NumberType': 'Float', - 'Dimensions': f'{len(self.times)}'} - time_data.text = ' '.join(map(str,self.times)) + 'Dimensions': f'{len(times)}'} + time_data.text = ' '.join(map(str,times)) attributes = [] data_items = [] - for inc in self.increments: + with h5py.File(self.fname,'r') as f: + for inc in self.visible['increments']: - grid=ET.SubElement(collection,'Grid') - grid.attrib = {'GridType': 'Uniform', - 'Name': inc} + grid = ET.SubElement(collection,'Grid') + grid.attrib = {'GridType': 'Uniform', + 'Name': inc} - topology=ET.SubElement(grid, 'Topology') - topology.attrib={'TopologyType': '3DCoRectMesh', - 'Dimensions': '{} {} {}'.format(*self.cells+1)} + topology = ET.SubElement(grid, 'Topology') + topology.attrib = {'TopologyType': '3DCoRectMesh', + 'Dimensions': '{} {} {}'.format(*(self.cells+1))} - geometry=ET.SubElement(grid, 'Geometry') - geometry.attrib={'GeometryType':'Origin_DxDyDz'} + 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) + 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.cells)) + delta = ET.SubElement(geometry, 'DataItem') + delta.attrib = {'Format': 'XML', + 'NumberType': 'Float', + 'Dimensions': '3'} + delta.text="{} {} {}".format(*(self.size/self.cells)) - - with h5py.File(self.fname,'r') as f: attributes.append(ET.SubElement(grid, 'Attribute')) - attributes[-1].attrib={'Name': 'u / m', - 'Center': 'Node', - 'AttributeType': 'Vector'} + attributes[-1].attrib = {'Name': 'u / m', + 'Center': 'Node', + 'AttributeType': 'Vector'} data_items.append(ET.SubElement(attributes[-1], 'DataItem')) - data_items[-1].attrib={'Format': 'HDF', - 'Precision': '8', - 'Dimensions': '{} {} {} 3'.format(*(self.cells+1))} - data_items[-1].text=f'{os.path.split(self.fname)[1]}:/{inc}/geometry/u_n' + data_items[-1].attrib = {'Format': 'HDF', + 'Precision': '8', + 'Dimensions': '{} {} {} 3'.format(*(self.cells+1))} + data_items[-1].text = f'{os.path.split(self.fname)[1]}:/{inc}/geometry/u_n' - for o,p in zip(['phases','homogenizations'],['out_type_ph','out_type_ho']): - 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]) + for ty in ['phase','homogenization']: + for label in self.visible[ty+'s']: + for field in self.visible['fields']: + for out in _match(output,f['/'.join([inc,ty,label,field])].keys()): + name = '/'.join([inc,ty,label,field,out]) shape = f[name].shape[1:] dtype = f[name].dtype - if dtype not in np.sctypes['int']+np.sctypes['uint']+np.sctypes['float']: continue unit = f[name].attrs[u] if h5py3 else f[name].attrs[u].decode() attributes.append(ET.SubElement(grid, 'Attribute')) - attributes[-1].attrib={'Name': name.split('/',2)[2]+f' / {unit}', - 'Center': 'Cell', - 'AttributeType': attribute_type_map[shape]} + attributes[-1].attrib = {'Name': name.split('/',2)[2]+f' / {unit}', + 'Center': 'Cell', + 'AttributeType': attribute_type_map[shape]} data_items.append(ET.SubElement(attributes[-1], 'DataItem')) - data_items[-1].attrib={'Format': 'HDF', - 'NumberType': number_type_map(dtype), - 'Precision': f'{dtype.itemsize}', - 'Dimensions': '{} {} {} {}'.format(*self.cells,1 if shape == () else - np.prod(shape))} - data_items[-1].text=f'{os.path.split(self.fname)[1]}:{name}' + data_items[-1].attrib = {'Format': 'HDF', + 'NumberType': number_type_map(dtype), + 'Precision': f'{dtype.itemsize}', + 'Dimensions': '{} {} {} {}'.format(*self.cells,1 if shape == () else + np.prod(shape))} + data_items[-1].text = f'{os.path.split(self.fname)[1]}:{name}' with open(self.fname.with_suffix('.xdmf').name,'w',newline='\n') as f: f.write(xml.dom.minidom.parseString(ET.tostring(xdmf).decode()).toprettyxml()) - def save_VTK(self,labels=[],mode='cell'): + def _mappings(self): + grp = 'mapping' if self.version_minor < 12 else 'cell_to' # compatibility hack + name = 'Name' if self.version_minor < 12 else 'label' # compatibility hack + member = 'member' if self.version_minor < 12 else 'entry' # compatibility hack + + with h5py.File(self.fname,'r') as f: + + at_cell_ph = [] + in_data_ph = [] + for c in range(self.N_constituents): + at_cell_ph.append({label: np.where(f['/'.join([grp,'phase'])][:,c][name] == label.encode())[0] \ + for label in self.visible['phases']}) + in_data_ph.append({label: f['/'.join([grp,'phase'])][member][at_cell_ph[c][label]][:,c] \ + for label in self.visible['phases']}) + + at_cell_ho = {label: np.where(f['/'.join([grp,'homogenization'])][:][name] == label.encode())[0] \ + for label in self.visible['homogenizations']} + in_data_ho = {label: f['/'.join([grp,'homogenization'])][member][at_cell_ho[label]] \ + for label in self.visible['homogenizations']} + + return at_cell_ph,in_data_ph,at_cell_ho,in_data_ho + + + def save_VTK(self,output='*',mode='cell',constituents=None,fill_float=np.nan,fill_int=0,parallel=True): """ - Export to vtk cell/point data. + Export to VTK cell/point data. Parameters ---------- - labels : str or list of, optional - Labels of the datasets to be exported. - mode : str, either 'cell' or 'point' + output : (list of) str, optional + Labels of the datasets to place. + Defaults to '*', in which case all datasets are exported. + mode : {'cell', 'point'} Export in cell format or point format. Defaults to 'cell'. + constituents : (list of) int, optional + Constituents to consider. + 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. + fill_int : int + Fill value for non-existent entries of integer type. + Defaults to 0. + parallel : bool + Write out VTK files in parallel in a separate background process. + Defaults to True. """ if mode.lower()=='cell': - - if self.structured: - v = VTK.from_rectilinear_grid(self.cells,self.size,self.origin) - else: - with h5py.File(self.fname,'r') as f: - v = VTK.from_unstructured_grid(f['/geometry/x_n'][()], - f['/geometry/T_c'][()]-1, - f['/geometry/T_c'].attrs['VTK_TYPE'] if h5py3 else \ - f['/geometry/T_c'].attrs['VTK_TYPE'].decode()) + v = self.geometry0 elif mode.lower()=='point': v = VTK.from_poly_data(self.coordinates0_point) - # compatibility hack - ln = 3 if self.version_minor < 12 else 10 - + ln = 3 if self.version_minor < 12 else 10 # compatibility hack N_digits = int(np.floor(np.log10(max(1,int(self.increments[-1][ln:])))))+1 - for inc in util.show_progress(self.iterate('increments'),len(self.visible['increments'])): + constituents_ = constituents if isinstance(constituents,Iterable) else \ + (range(self.N_constituents) if constituents is None else [constituents]) - viewed_backup_ho = self.visible['homogenizations'].copy() - self.view('homogenizations',False) - for label in (labels if isinstance(labels,list) else [labels]): - for o in self.iterate('out_type_ph'): - for c in range(self.N_constituents): - prefix = '' if self.N_constituents == 1 else f'constituent{c}/' - if o not in ['mechanics', 'mechanical']: # compatibility hack - for _ in self.iterate('phases'): - path = self.get_dataset_location(label) - if len(path) == 0: - continue - array = self.read_dataset(path,c) - v.add(array,prefix+path[0].split('/',1)[1]+f' / {self._get_attribute(path[0],"unit")}') - else: - paths = self.get_dataset_location(label) - if len(paths) == 0: - continue - array = self.read_dataset(paths,c) - if self.version_minor < 12: - ph_name = re.compile(r'(?<=(phase\/))(.*?)(?=(mechanics))') # identify phase name - else: - ph_name = re.compile(r'(?<=(phase\/))(.*?)(?=(mechanical))') # identify phase name - dset_name = prefix+re.sub(ph_name,r'',paths[0].split('/',1)[1]) # remove phase name - v.add(array,dset_name+f' / {self._get_attribute(paths[0],"unit")}') - self.view('homogenizations',viewed_backup_ho) + suffixes = [''] if self.N_constituents == 1 or isinstance(constituents,int) else \ + [f'#{c}' for c in constituents_] - viewed_backup_ph = self.visible['phases'].copy() - self.view('phases',False) - for label in (labels if isinstance(labels,list) else [labels]): - for _ in self.iterate('out_type_ho'): - paths = self.get_dataset_location(label) - if len(paths) == 0: - continue - array = self.read_dataset(paths) - v.add(array,paths[0].split('/',1)[1]+f' / {self._get_attribute(paths[0],"unit")}') - self.view('phases',viewed_backup_ph) + at_cell_ph,in_data_ph,at_cell_ho,in_data_ho = self._mappings() - u = self.read_dataset(self.get_dataset_location('u_n' if mode.lower() == 'cell' else 'u_p')) - v.add(u,'u') + with h5py.File(self.fname,'r') as f: - v.save(f'{self.fname.stem}_inc{inc[ln:].zfill(N_digits)}') + for inc in util.show_progress(self.visible['increments']): + + u = _read(f['/'.join([inc,'geometry','u_n' if mode.lower() == 'cell' else 'u_p'])]) + v.add(u,'u') + + for ty in ['phase','homogenization']: + for field in self.visible['fields']: + outs = {} + for label in self.visible[ty+'s']: + if field not in f['/'.join([inc,ty,label])].keys(): continue + + for out in _match(output,f['/'.join([inc,ty,label,field])].keys()): + data = ma.array(_read(f['/'.join([inc,ty,label,field,out])])) + + if ty == 'phase': + if out+suffixes[0] not in outs.keys(): + for c,suffix in zip(constituents_,suffixes): + outs[out+suffix] = \ + _empty_like(data,self.N_materialpoints,fill_float,fill_int) + + for c,suffix in zip(constituents_,suffixes): + outs[out+suffix][at_cell_ph[c][label]] = data[in_data_ph[c][label]] + + if ty == 'homogenization': + if out not in outs.keys(): + outs[out] = _empty_like(data,self.N_materialpoints,fill_float,fill_int) + + outs[out][at_cell_ho[label]] = data[in_data_ho[label]] + + for label,dataset in outs.items(): + v.add(dataset,' / '.join(['/'.join([ty,field,label]),dataset.dtype.metadata['unit']])) + + v.save(f'{self.fname.stem}_inc{inc[ln:].zfill(N_digits)}',parallel=parallel) + + + def get(self,output='*',flatten=True,prune=True): + """ + Collect data per phase/homogenization reflecting the group/folder structure in the DADF5 file. + + Parameters + ---------- + output : (list of) str + Labels of the datasets to read. + Defaults to '*', in which case all datasets are read. + flatten : bool + Remove singular levels of the folder hierarchy. + This might be beneficial in case of single increment, + phase/homogenization, or field. Defaults to True. + prune : bool + Remove branches with no data. Defaults to True. + + Returns + ------- + data : dict of numpy.ndarray + Datasets structured by phase/homogenization and according to selected view. + + """ + r = {} + + with h5py.File(self.fname,'r') as f: + for inc in util.show_progress(self.visible['increments']): + r[inc] = {'phase':{},'homogenization':{},'geometry':{}} + + for out in _match(output,f['/'.join([inc,'geometry'])].keys()): + r[inc]['geometry'][out] = _read(f['/'.join([inc,'geometry',out])]) + + for ty in ['phase','homogenization']: + for label in self.visible[ty+'s']: + r[inc][ty][label] = {} + for field in _match(self.visible['fields'],f['/'.join([inc,ty,label])].keys()): + r[inc][ty][label][field] = {} + for out in _match(output,f['/'.join([inc,ty,label,field])].keys()): + r[inc][ty][label][field][out] = _read(f['/'.join([inc,ty,label,field,out])]) + + if prune: r = util.dict_prune(r) + if flatten: r = util.dict_flatten(r) + + return r + + + def place(self,output='*',flatten=True,prune=True,constituents=None,fill_float=np.nan,fill_int=0): + """ + Merge data into spatial order that is compatible with the damask.VTK geometry representation. + + The returned data structure reflects the group/folder structure + in the DADF5 file. + + Multi-phase data is fused into a single output. + `place` is equivalent to `read` if only one phase/homogenization + and one constituent is present. + + Parameters + ---------- + output : (list of) str, optional + Labels of the datasets to place. + Defaults to '*', in which case all datasets are placed. + flatten : bool + Remove singular levels of the folder hierarchy. + This might be beneficial in case of single increment or field. + Defaults to True. + prune : bool + 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. + fill_float : float + Fill value for non-existent entries of floating point type. + Defaults to NaN. + fill_int : int + Fill value for non-existent entries of integer type. + Defaults to 0. + + Returns + ------- + data : dict of numpy.ma.MaskedArray + Datasets structured by spatial position and according to selected view. + + """ + r = {} + + constituents_ = constituents if isinstance(constituents,Iterable) else \ + (range(self.N_constituents) if constituents is None else [constituents]) + + suffixes = [''] if self.N_constituents == 1 or isinstance(constituents,int) else \ + [f'#{c}' for c in constituents_] + + at_cell_ph,in_data_ph,at_cell_ho,in_data_ho = self._mappings() + + with h5py.File(self.fname,'r') as f: + + for inc in util.show_progress(self.visible['increments']): + r[inc] = {'phase':{},'homogenization':{},'geometry':{}} + + for out in _match(output,f['/'.join([inc,'geometry'])].keys()): + r[inc]['geometry'][out] = ma.array(_read(f['/'.join([inc,'geometry',out])]),fill_value = fill_float) + + 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()): + if field not in r[inc][ty].keys(): + r[inc][ty][field] = {} + + for out in _match(output,f['/'.join([inc,ty,label,field])].keys()): + data = ma.array(_read(f['/'.join([inc,ty,label,field,out])])) + + if ty == 'phase': + if out+suffixes[0] not in r[inc][ty][field].keys(): + for c,suffix in zip(constituents_,suffixes): + r[inc][ty][field][out+suffix] = \ + _empty_like(data,self.N_materialpoints,fill_float,fill_int) + + for c,suffix in zip(constituents_,suffixes): + r[inc][ty][field][out+suffix][at_cell_ph[c][label]] = data[in_data_ph[c][label]] + + if ty == 'homogenization': + if out not in r[inc][ty][field].keys(): + r[inc][ty][field][out] = \ + _empty_like(data,self.N_materialpoints,fill_float,fill_int) + + r[inc][ty][field][out][at_cell_ho[label]] = data[in_data_ho[label]] + + if prune: r = util.dict_prune(r) + if flatten: r = util.dict_flatten(r) + + return r diff --git a/python/damask/_table.py b/python/damask/_table.py index e5f436f75..01d8c9863 100644 --- a/python/damask/_table.py +++ b/python/damask/_table.py @@ -230,7 +230,6 @@ class Table: f = fname f.seek(0) - f.seek(0) comments = [] line = f.readline().strip() while line.startswith('#'): @@ -515,7 +514,7 @@ class Table: """ if set(self.shapes) & set(other.shapes) or self.data.shape[0] != other.data.shape[0]: - raise KeyError('Dublicated keys or row count mismatch') + raise KeyError('Duplicated keys or row count mismatch') else: dup = self.copy() dup.data = dup.data.join(other.data) diff --git a/python/damask/_vtk.py b/python/damask/_vtk.py index ebd2f39a6..340bd1ea3 100644 --- a/python/damask/_vtk.py +++ b/python/damask/_vtk.py @@ -2,8 +2,8 @@ import os import multiprocessing as mp from pathlib import Path -import pandas as pd 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 @@ -224,14 +224,14 @@ class VTK: # Check https://blog.kitware.com/ghost-and-blanking-visibility-changes/ for missing data - # Needs support for pd.DataFrame and/or table + # Needs support for damask.Table def add(self,data,label=None): """ Add data to either cells or points. Parameters ---------- - data : numpy.ndarray + data : numpy.ndarray or numpy.ma.MaskedArray Data to add. First dimension needs to match either number of cells or number of points. label : str @@ -246,8 +246,10 @@ class VTK: raise ValueError('No label defined for numpy.ndarray') N_data = data.shape[0] - 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 + data_ = np.where(data.mask,data.fill_value,data) if isinstance(data,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 d.SetName(label) if N_data == N_points: @@ -256,8 +258,6 @@ class VTK: self.vtk_data.GetCellData().AddArray(d) else: raise ValueError(f'Cell / point count ({N_cells} / {N_points}) differs from data ({N_data}).') - elif isinstance(data,pd.DataFrame): - raise NotImplementedError('pd.DataFrame') elif isinstance(data,Table): raise NotImplementedError('damask.Table') else: diff --git a/python/damask/util.py b/python/damask/util.py index 902e0409c..3fe30320a 100644 --- a/python/damask/util.py +++ b/python/damask/util.py @@ -17,15 +17,16 @@ __all__=[ 'srepr', 'emph','deemph','warn','strikeout', 'execute', + 'natural_sort', 'show_progress', 'scale_to_coprime', 'project_stereographic', 'hybrid_IA', - 'return_message', 'execution_stamp', 'shapeshifter', 'shapeblender', 'extend_docstring', 'extended_docstring', - 'DREAM3D_base_group', 'DREAM3D_cell_data_group' + 'DREAM3D_base_group', 'DREAM3D_cell_data_group', + 'dict_prune', 'dict_flatten' ] # https://svn.blender.org/svnroot/bf-blender/trunk/blender/build_files/scons/tools/bcolors.py @@ -58,11 +59,12 @@ def srepr(arg,glue = '\n'): Glue used for joining operation. Defaults to \n. """ - if (not hasattr(arg, "strip") and - (hasattr(arg, "__getitem__") or - hasattr(arg, "__iter__"))): + if (not hasattr(arg, 'strip') and + (hasattr(arg, '__getitem__') or + hasattr(arg, '__iter__'))): return glue.join(str(x) for x in arg) - return arg if isinstance(arg,str) else repr(arg) + else: + return arg if isinstance(arg,str) else repr(arg) def emph(what): @@ -112,6 +114,11 @@ def execute(cmd,wd='./',env=None): return process.stdout, process.stderr +def natural_sort(key): + convert = lambda text: int(text) if text.isdigit() else text + return [ convert(c) for c in re.split('([0-9]+)', key) ] + + def show_progress(iterable,N_iter=None,prefix='',bar_length=50): """ Decorate a loop with a status bar. @@ -130,11 +137,15 @@ def show_progress(iterable,N_iter=None,prefix='',bar_length=50): Character length of bar. Defaults to 50. """ - status = _ProgressBar(N_iter if N_iter else len(iterable),prefix,bar_length) + if N_iter == 1 or (hasattr(iterable,'__len__') and len(iterable) == 1): + for item in iterable: + yield item + else: + status = _ProgressBar(N_iter if N_iter is not None else len(iterable),prefix,bar_length) - for i,item in enumerate(iterable): - yield item - status.update(i) + for i,item in enumerate(iterable): + yield item + status.update(i) def scale_to_coprime(v): @@ -388,7 +399,7 @@ def DREAM3D_cell_data_group(fname): """ base_group = DREAM3D_base_group(fname) with h5py.File(fname,'r') as f: - cells = tuple(f[os.path.join(base_group,'_SIMPL_GEOMETRY','DIMENSIONS')][()][::-1]) + cells = tuple(f['/'.join([base_group,'_SIMPL_GEOMETRY','DIMENSIONS'])][()][::-1]) cell_data_group = f[base_group].visititems(lambda path,obj: path.split('/')[0] \ if isinstance(obj,h5py._hl.dataset.Dataset) and np.shape(obj)[:-1] == cells \ else None) @@ -399,29 +410,59 @@ def DREAM3D_cell_data_group(fname): return cell_data_group +def dict_prune(d): + """ + Recursively remove empty dictionaries. + + Parameters + ---------- + d : dict + Dictionary to prune. + + Returns + ------- + pruned : dict + Pruned dictionary. + + """ + # https://stackoverflow.com/questions/48151953 + new = {} + for k,v in d.items(): + if isinstance(v, dict): + v = dict_prune(v) + if not isinstance(v,dict) or v != {}: + new[k] = v + return new + + +def dict_flatten(d): + """ + Recursively remove keys of single-entry dictionaries. + + Parameters + ---------- + d : dict + Dictionary to flatten. + + Returns + ------- + flattened : dict + Flattened dictionary. + + """ + if isinstance(d,dict) and len(d) == 1: + entry = d[list(d.keys())[0]] + new = dict_flatten(entry.copy()) if isinstance(entry,dict) else entry + else: + new = {k: (dict_flatten(v) if isinstance(v, dict) else v) for k,v in d.items()} + + return new + + + #################################################################################################### # Classes #################################################################################################### -class return_message: - """Object with formatted return message.""" - - def __init__(self,message): - """ - Set return message. - - Parameters - ---------- - message : str or list of str - message for output to screen - - """ - self.message = message - - def __repr__(self): - """Return message suitable for interactive shells.""" - return srepr(self.message) - - class _ProgressBar: """ Report progress of an interation as a status bar. diff --git a/python/tests/reference/ConfigMaterial/measured.material_yaml b/python/tests/reference/ConfigMaterial/measured.material.yaml similarity index 100% rename from python/tests/reference/ConfigMaterial/measured.material_yaml rename to python/tests/reference/ConfigMaterial/measured.material.yaml diff --git a/python/tests/reference/Result/12grains6x7x8_tensionY.yaml b/python/tests/reference/Result/12grains6x7x8.material.yaml similarity index 100% rename from python/tests/reference/Result/12grains6x7x8_tensionY.yaml rename to python/tests/reference/Result/12grains6x7x8.material.yaml diff --git a/python/tests/reference/Result/4grains2x4x3.material.yaml b/python/tests/reference/Result/4grains2x4x3.material.yaml new file mode 100644 index 000000000..4cff681c0 --- /dev/null +++ b/python/tests/reference/Result/4grains2x4x3.material.yaml @@ -0,0 +1,678 @@ +material: + - constituents: + - phase: A + O: [0.44698096036898677, -0.035696623617731814, 0.050331402131209624, -0.8924127532086364] + v: 0.125 + - phase: A + O: [0.6226205515354057, 0.4257792272120092, -0.6365250917893261, 0.16090837766592916] + v: 0.125 + - phase: A + O: [0.4319357791504236, 0.5804552577310466, 0.4870786336250474, 0.48913963356904305] + v: 0.125 + - phase: A + O: [0.6712122425649047, -0.5315470950319967, 0.11982471915778993, -0.5025672570639611] + v: 0.125 + - phase: A + O: [0.03979536866120591, -0.01919843067074171, -0.5618628311727827, 0.8260495795286167] + v: 0.125 + - phase: A + O: [0.6951890945163803, -0.46325436790720337, 0.03393768848173515, -0.5485943371753935] + v: 0.125 + - phase: A + O: [0.3489469292795834, -0.5959323325634578, -0.6067706771532161, 0.3936115355256417] + v: 0.125 + - phase: A + O: [0.39500485674202723, -0.07573446369604235, -0.015761384830863343, -0.9154163167144754] + v: 0.125 + homogenization: RGC + - constituents: + - phase: B + O: [0.9221278604171202, 0.32701451341562027, 0.15952215192749652, -0.1315081750406033] + v: 0.125 + - phase: B + O: [0.8117843306311869, 0.4102779561929707, -0.18584663447455768, 0.37167085930736465] + v: 0.125 + - phase: B + O: [0.6548302235705483, 0.6510678661349768, -0.3832730015357094, 0.02024396894884073] + v: 0.125 + - phase: B + O: [0.04255709858343937, -0.3005730931619659, -0.3060125731417352, -0.9023308783957146] + v: 0.125 + - phase: B + O: [0.6960740097259569, 0.16329320418261428, 0.6991228405747408, 0.006599715032396515] + v: 0.125 + - phase: B + O: [0.48818278044780217, -0.664294296073913, 0.2679908412891581, 0.4985695238008905] + v: 0.125 + - phase: B + O: [0.545426202615872, 0.179298582176746, -0.7847831868864313, -0.23340442478628137] + v: 0.125 + - phase: B + O: [0.8566603715426528, 0.26967937521434965, -0.03159131558814077, 0.43864339866435076] + v: 0.125 + homogenization: RGC + - constituents: + - phase: C + O: [0.140304599816322, -0.017971728003341014, -0.9895572782263894, 0.027713342853867628] + v: 0.125 + - phase: C + O: [0.8096406385129331, 0.2395956967578664, -0.25531574170054405, -0.47105181307727034] + v: 0.125 + - phase: C + O: [0.15570676199790476, -0.7501738841096782, 0.2328253832435299, -0.5989882209070811] + v: 0.125 + - phase: C + O: [0.4595649941522109, 0.3331746723324176, 0.23957383093695328, 0.7876541331042811] + v: 0.125 + - phase: C + O: [0.2969510261483044, 0.0027640644222379847, -0.8133130954773664, 0.5003341450894221] + v: 0.125 + - phase: C + O: [0.114170281047612, -0.9068850613039258, 0.3701597984948465, -0.16585040273553361] + v: 0.125 + - phase: C + O: [0.29658840046433166, -0.5331902828876833, 0.7365409179409275, 0.2919776004129383] + v: 0.125 + - phase: C + O: [0.931315147629592, -0.2017743466108175, 0.30303719188324046, 0.010376376100021569] + v: 0.125 + homogenization: RGC + - constituents: + - phase: A + O: [0.20901169650798976, 0.5778057521557359, -0.3628020593274987, 0.7005920990464581] + v: 0.125 + - phase: B + O: [0.6025573623737024, 0.7345291803280383, 0.3086529661391562, 0.04609614722911203] + v: 0.125 + - phase: C + O: [0.6402291970447783, -0.5402560454335561, -0.5285008279381848, -0.13753856002062462] + v: 0.125 + - phase: A + O: [0.12506307996119548, 0.5336909245447153, 0.8344001590799459, 0.05752910234470594] + v: 0.125 + - phase: B + O: [0.5504712846116686, -0.6241210958640075, 0.42061504390637533, 0.3612993320712448] + v: 0.125 + - phase: C + O: [0.3135636745691277, 0.6574027353927014, 0.6044800436379738, -0.32265049563317466] + v: 0.125 + - phase: A + O: [0.1619452502554878, -0.4269267792831907, -0.8080415555188633, 0.3722581169097929] + v: 0.125 + - phase: A + O: [0.17044595054343245, 0.8646695166660998, -0.30004868613437863, -0.365055599656809] + v: 0.125 + homogenization: RGC + - constituents: + - phase: B + O: [0.43255174377606787, 0.2026672781261794, 0.35000490906160386, 0.8058048938583007] + v: 0.125 + - phase: B + O: [0.19591912697365416, -0.8181375036268687, 0.3260880890153691, -0.4311998133665892] + v: 0.125 + - phase: A + O: [0.1373248436058589, -0.43524771427953446, -0.8885320171891298, -0.047033700395389226] + v: 0.125 + - phase: A + O: [0.250519065956112, -0.7636612069861932, 0.5672425862778043, -0.17971534951065127] + v: 0.125 + - phase: B + O: [0.600521695167639, -0.2273230963128902, 0.7563911883237204, -0.12478090295368073] + v: 0.125 + - phase: A + O: [0.4281846857945264, 0.5284504248063085, 0.21428541450090705, 0.7010561921167584] + v: 0.125 + - phase: A + O: [0.34007599923699156, 0.5500172687924186, 0.6138467460818403, -0.45279298923219496] + v: 0.125 + - phase: A + O: [0.8581382626679844, -0.20034169756111123, 0.3134907520728513, -0.35381559424127107] + v: 0.125 + homogenization: RGC + - constituents: + - phase: C + O: [0.18637512817696594, 0.9354389080016575, 0.2967480269088709, -0.04646471262558266] + v: 0.125 + - phase: C + O: [0.33425800495572805, 0.24948308477875414, 0.7297193930640264, 0.5417927499686225] + v: 0.125 + - phase: A + O: [0.11918214701666907, -0.4185869380028542, 0.618905484805619, -0.6538628235673086] + v: 0.125 + - phase: A + O: [0.40732205687168244, -0.6166051531766635, -0.6737097202702335, -0.0014282419993970551] + v: 0.125 + - phase: A + O: [0.68397144415648, 0.6848389352552103, 0.24111983008651763, 0.07099242125789083] + v: 0.125 + - phase: C + O: [0.05719438496521803, 0.2575447319897859, 0.5020591650211859, -0.8236116245968056] + v: 0.125 + - phase: A + O: [0.10565206937493385, -0.1109750475948143, 0.8781814355546643, 0.45312199824690885] + v: 0.125 + - phase: A + O: [0.22821192086563635, -0.9595048060506683, -0.07745760541549483, -0.1458429487626446] + v: 0.125 + homogenization: RGC + - constituents: + - phase: A + O: [0.5745354940858776, 0.6366614824775295, 0.3307586622977252, -0.39391601907009377] + v: 0.125 + - phase: B + O: [0.3505113084912295, -0.6046500894621007, 0.6274019946374426, -0.34337563841687757] + v: 0.125 + - phase: A + O: [0.030949493058962222, 0.4416700940450431, 0.8819578143863066, 0.16161704906526728] + v: 0.125 + - phase: A + O: [0.13563849174261922, 0.2687388407244776, 0.809782239640578, -0.5036212459840637] + v: 0.125 + - phase: A + O: [0.2305505587522992, -0.6837828098041563, 0.3813348695443346, 0.577815910255975] + v: 0.125 + - phase: A + O: [0.10699977862890839, 0.8478562366641894, 0.26664808851893185, -0.4456339823355066] + v: 0.125 + - phase: A + O: [0.9203073703838264, 0.35592932457392046, 0.15147737103605818, -0.058337517855697775] + v: 0.125 + - phase: B + O: [0.8544895693446002, 0.43485422381677186, -0.2696115267291114, 0.08977195867747481] + v: 0.125 + homogenization: RGC + - constituents: + - phase: A + O: [0.22196756947117874, 0.940293527559, -0.2573416285961528, -0.018808676859083388] + v: 0.125 + - phase: A + O: [0.28669138183611736, -0.5847463393505065, 0.6896415010701208, 0.3166612862331461] + v: 0.125 + - phase: A + O: [0.18214681844078473, 0.07037956730772883, 0.6195004937983107, -0.7603212421214639] + v: 0.125 + - phase: A + O: [0.5744613003376511, 0.17468957401795682, -0.7944671360487727, 0.09110289173380101] + v: 0.125 + - phase: C + O: [0.4329446256398051, -0.7351014072423584, 0.5116672951703971, -0.10188940697110307] + v: 0.125 + - phase: C + O: [0.2946245183318559, -0.8241105780568372, 0.40376354741611825, -0.2664829189845001] + v: 0.125 + - phase: A + O: [0.22212591383189398, -0.06959072743195248, -0.8166434230649855, -0.5281199945320578] + v: 0.125 + - phase: A + O: [0.038576201076459135, 0.9876698316278196, 0.06851090552108202, -0.13537516843004982] + v: 0.125 + homogenization: RGC + - constituents: + - phase: C + O: [0.004892880719601024, -0.02167422880443418, 0.3339942463386166, 0.9423131809205983] + v: 0.125 + - phase: A + O: [0.3719521264362446, -0.4055315244649798, -0.39177426332273346, 0.7373660725193388] + v: 0.125 + - phase: A + O: [0.010947854228877327, 0.16142752376310676, -0.19859031041582603, -0.9666349816080735] + v: 0.125 + - phase: B + O: [0.9101129345691923, -0.17592965401443023, -0.30944622338368544, 0.21209959453471436] + v: 0.125 + - phase: A + O: [0.5280075626141939, -0.4828698476824812, -0.1580687189774305, 0.6804843893155447] + v: 0.125 + - phase: A + O: [0.8213575048478461, 0.18229514326534527, -0.35486977770450967, -0.4076858727549186] + v: 0.125 + - phase: A + O: [0.16635046542653203, -0.8410358537737095, -0.47641083611371204, 0.19498443669415633] + v: 0.125 + - phase: B + O: [0.5406080589325533, 0.5405688205487984, -0.4980750162732086, -0.4092060056158767] + v: 0.125 + homogenization: RGC + - constituents: + - phase: A + O: [0.41778322805525603, 0.23584193505422144, 0.8768801833720732, -0.03028035724639885] + v: 0.125 + - phase: A + O: [0.8727428843707065, 0.0764211249039828, 0.3927224252414984, -0.2797298092108618] + v: 0.125 + - phase: A + O: [0.8766166377340643, -0.2644877070285106, 0.022928146762306964, -0.40132757613285325] + v: 0.125 + - phase: A + O: [0.6643649540361558, -0.04353445577200276, -0.2630693201131682, -0.6982252443333509] + v: 0.125 + - phase: A + O: [0.43011162904952843, -0.7133156350005553, 0.22368449997894274, -0.5061126711408104] + v: 0.125 + - phase: A + O: [0.42119993016818996, -0.2666762728079305, 0.7798784815692744, 0.37850223028772895] + v: 0.125 + - phase: A + O: [0.17194138099812353, 0.48707029434265314, 0.28918455812138427, -0.8059596647559721] + v: 0.125 + - phase: A + O: [0.7411029340880614, -0.6424443807902435, 0.06269338980849715, 0.18466509565000905] + v: 0.125 + homogenization: RGC + - constituents: + - phase: A + O: [0.2238963963710994, -0.0650586111271435, -0.4249089470195672, -0.8746943280672199] + v: 0.125 + - phase: A + O: [0.11834025356863993, -0.8856449810417544, 0.417012782917471, 0.16651994122112387] + v: 0.125 + - phase: A + O: [0.11716309944373668, 0.0038756181092136927, -0.7000966032560261, -0.7043596622623864] + v: 0.125 + - phase: A + O: [0.45702665319142677, -0.16783958203387148, -0.15365667123574342, 0.8598523945190182] + v: 0.125 + - phase: A + O: [0.5197700015250588, 0.41758011415389684, -0.06820554688133636, 0.7421684425738383] + v: 0.125 + - phase: A + O: [0.675852597832269, -0.5552996378928718, 0.4213174293100185, 0.23949363648960756] + v: 0.125 + - phase: A + O: [0.5936583232008587, 0.6945289973655595, -0.40553899121228437, -0.027154994370431035] + v: 0.125 + - phase: A + O: [0.1865986177378815, 0.3352356341235895, 0.35399611880081633, -0.8529271793922534] + v: 0.125 + homogenization: RGC + - constituents: + - phase: A + O: [0.04264635665959766, -0.25662858077651657, 0.9118921459098731, -0.3174520027030545] + v: 0.125 + - phase: A + O: [0.8821155978581948, 0.3761807268996467, 0.1631361963837017, 0.23183337584133834] + v: 0.125 + - phase: A + O: [0.162154410316598, -0.8783087060961443, 0.0990284993753847, -0.4387175860642615] + v: 0.125 + - phase: A + O: [0.5420568096555354, 0.5088379086577298, 0.34515629818391275, -0.572822422433749] + v: 0.125 + - phase: A + O: [0.4535105167294993, -0.351824950997862, 0.5869707747562448, -0.5709752399650516] + v: 0.125 + - phase: A + O: [0.8098666527705685, 0.36698878904202453, -0.4428476393801139, 0.11541751055678014] + v: 0.125 + - phase: A + O: [0.2576165318782538, 0.537377864604212, -0.6660199792713344, 0.44863809506978913] + v: 0.125 + - phase: A + O: [0.665871020774723, -0.24779933206030424, -0.6964987066828571, 0.10050286718299539] + v: 0.125 + homogenization: RGC + - constituents: + - phase: A + O: [0.5998961636312351, -0.5496518322385892, 0.10155928244475286, 0.5724449041843198] + v: 0.125 + - phase: A + O: [0.7039639404137236, 0.03679879173607231, 0.58207433274894, 0.4053024681380869] + v: 0.125 + - phase: A + O: [0.40054268633071116, 0.37951399049616763, 0.5867236594596857, -0.5926972539795397] + v: 0.125 + - phase: A + O: [0.9355573844663049, 0.24834752099784146, 0.19019564545030732, 0.16395580391231754] + v: 0.125 + - phase: A + O: [0.7776851244299001, 0.5132289287723679, -0.3500697782737797, 0.0961928492714775] + v: 0.125 + - phase: A + O: [0.14184146561701433, -0.12106060201428649, 0.06736696083733706, 0.9801464287845448] + v: 0.125 + - phase: A + O: [0.35378845612452303, 0.6454299208817863, -0.6758747727133072, 0.03804257027716481] + v: 0.125 + - phase: A + O: [0.5450137677982273, -0.6067354536251314, -0.5773779505037945, 0.03829862264786783] + v: 0.125 + homogenization: RGC + - constituents: + - phase: A + O: [0.8960165930555501, 0.24799888375444512, -0.3529216677489853, 0.10534284531447227] + v: 0.125 + - phase: A + O: [0.018356029815615817, -0.2889337455527527, 0.9520448827103584, 0.09894891689798985] + v: 0.125 + - phase: A + O: [0.853080700130541, -0.3850382673796441, 0.3507108308253732, 0.03163486778610022] + v: 0.125 + - phase: A + O: [0.3115000788768896, -0.36279737244767013, 0.6484519555216468, -0.5923308440263011] + v: 0.125 + - phase: A + O: [0.44453200197416054, 0.5366087739174971, -0.4155131463908387, -0.5846290688564767] + v: 0.125 + - phase: A + O: [0.6294708415113138, -0.5730458394479108, 0.3971797612768184, -0.34297691294104227] + v: 0.125 + - phase: A + O: [0.7012827551365688, -0.6926674077585473, -0.1570646347205316, -0.06119689614043776] + v: 0.125 + - phase: A + O: [0.344127722728084, -0.02877253391263471, -0.12761292039202274, 0.9297651285627186] + v: 0.125 + homogenization: RGC + - constituents: + - phase: A + O: [0.6978484686818835, 0.22305447962726752, 0.026182556001199425, -0.6801240237175888] + v: 0.125 + - phase: A + O: [0.9628977633574072, 0.16069766981498335, 0.2167078702347067, -0.006469560701840579] + v: 0.125 + - phase: A + O: [0.5032342094162214, -0.32097872787787546, 0.8022984104181043, -0.006726616067118365] + v: 0.125 + - phase: A + O: [0.8004358150893949, 0.0875805156227694, -0.3568768378095974, -0.47357267851983204] + v: 0.125 + - phase: A + O: [0.5880748138588325, 0.45495016865260346, 0.5556078663040557, -0.37214010298397254] + v: 0.125 + - phase: A + O: [0.0242946564098417, 0.6005127022483121, 0.7620365834257689, -0.24102802664656872] + v: 0.125 + - phase: A + O: [0.34550850696510604, 0.7299878005482168, -0.29949045445085315, -0.5079834154363129] + v: 0.125 + - phase: A + O: [0.10265979489126274, 0.34339772969920185, 0.0799881993576351, 0.9301294822302113] + v: 0.125 + homogenization: RGC + - constituents: + - phase: A + O: [0.5849010547063512, -0.46890126598541493, -0.141130615667434, -0.6466100125129551] + v: 0.125 + - phase: A + O: [0.3232591731614976, -0.2858051445366468, 0.8855713416105193, -0.17199513144701628] + v: 0.125 + - phase: A + O: [0.7530648506011821, -0.3316323026169325, -0.30797159948874325, -0.477563441396382] + v: 0.125 + - phase: A + O: [0.5841938587340145, -0.5195581188838846, 0.45610755437383976, 0.4251385601923404] + v: 0.125 + - phase: A + O: [0.5699287377636442, -0.7161652717835199, -0.40080048103879345, -0.04058955236816678] + v: 0.125 + - phase: A + O: [0.05536551119448681, -0.40478746018638956, -0.6031583160532334, -0.6850414717532457] + v: 0.125 + - phase: A + O: [0.2227977302479141, -0.48882932798587353, -0.15904170725500447, -0.8283192590122909] + v: 0.125 + - phase: A + O: [0.05275211816268626, -0.223660120175789, 0.6943402501868806, 0.6819713935662708] + v: 0.125 + homogenization: RGC + - constituents: + - phase: A + O: [0.7102673449281939, 0.43290229875199215, 0.4268837400347354, 0.35480441225815024] + v: 0.125 + - phase: A + O: [0.2541648349387051, 0.3725842249682419, -0.13402399303768983, -0.8823937903655196] + v: 0.125 + - phase: A + O: [0.4091596391203177, 0.10748924891263324, 0.2807193957630742, -0.8615283349522196] + v: 0.125 + - phase: A + O: [0.16364908456142396, -0.9217193970202214, 0.11860582627366922, 0.33103623404821947] + v: 0.125 + - phase: A + O: [0.80868984786943, -0.5530964468638372, 0.18400280484186798, -0.07904440669549186] + v: 0.125 + - phase: A + O: [0.785539568992473, -0.42139749663876175, 0.03451711851851934, 0.45184101618034067] + v: 0.125 + - phase: A + O: [0.07920711616878613, -0.6177377682037178, -0.7622215143727179, 0.1764784562213612] + v: 0.125 + - phase: A + O: [0.6079704911212058, -0.06435941038671157, 0.3854845277687917, 0.6911088388028229] + v: 0.125 + homogenization: RGC + - constituents: + - phase: A + O: [0.6650964387751211, -0.34961073112696767, -0.4316460762174112, -0.4990999185490129] + v: 0.125 + - phase: A + O: [0.1643442326674122, 0.21205555067771598, 0.18237328729946675, -0.9459193415378059] + v: 0.125 + - phase: A + O: [0.5079103592861302, -0.4522879915621313, 0.7008603696203147, 0.21507529359320504] + v: 0.125 + - phase: A + O: [0.4588629674996609, 0.7103482146539822, -0.14094577428597305, -0.5147664321867083] + v: 0.125 + - phase: A + O: [0.044305288641285946, -0.3875768281048735, -0.7143901476515312, -0.5809199261972351] + v: 0.125 + - phase: A + O: [0.8455643302581954, -0.4342346154649184, -0.29235576797838947, -0.10483018199359342] + v: 0.125 + - phase: A + O: [0.5790795198800359, 0.5645942182636775, -0.26425181035873874, 0.5254248367567554] + v: 0.125 + - phase: A + O: [0.6475654656469717, 0.04725837905512413, 0.4718561325845725, -0.5964707901086468] + v: 0.125 + homogenization: RGC + - constituents: + - phase: A + O: [0.9334013387182871, -0.3485154208652402, -0.04006434439320742, 0.07545721043330728] + v: 0.125 + - phase: A + O: [0.38741862953083567, -0.9054265066714982, 0.12281627778824236, 0.12257980428821813] + v: 0.125 + - phase: A + O: [0.8578234268453049, 0.25286192172706, 0.21735917871049656, -0.3910943675459601] + v: 0.125 + - phase: A + O: [0.4841852989554279, -0.5324884938481191, 0.6128392956763806, -0.32626461326610684] + v: 0.125 + - phase: A + O: [0.24479559088002964, -0.39790158622306016, -0.6885075670842895, 0.5547132380199182] + v: 0.125 + - phase: A + O: [0.006342119419357721, 0.8016474931043842, -0.09875135732320732, 0.589550034982232] + v: 0.125 + - phase: A + O: [0.5815233878895042, -0.41292157227416415, -0.6523344380499, -0.2564880219859532] + v: 0.125 + - phase: A + O: [0.08441325811399321, -0.9287379499833768, -0.04047714959468877, 0.3587224867163255] + v: 0.125 + homogenization: RGC + - constituents: + - phase: A + O: [0.9272909114072972, 0.22505658735105372, -0.12087397257258231, 0.27362489080097285] + v: 0.125 + - phase: A + O: [0.22477798467053944, -0.16730490431874273, 0.7619888558500476, 0.5838295214860949] + v: 0.125 + - phase: A + O: [0.6944164222400087, -0.3089622138004769, 0.47430811205294454, -0.44425217816873547] + v: 0.125 + - phase: A + O: [0.9319464859005648, -0.33314612927635295, 0.09454266632036588, 0.10747598899664906] + v: 0.125 + - phase: A + O: [0.7631364561921935, 0.6348137033139695, -0.07779160762534851, 0.09264327875398014] + v: 0.125 + - phase: A + O: [0.30213352773158436, 0.7607414448583625, -0.5489414690920735, -0.1692662075144206] + v: 0.125 + - phase: A + O: [0.44393342791682283, -0.20509632166957376, 0.6929260547202925, 0.529822699688679] + v: 0.125 + - phase: A + O: [0.5553300723266665, 0.8100017529476545, -0.18277171291698383, 0.045827633026130514] + v: 0.125 + homogenization: RGC + - constituents: + - phase: A + O: [0.3369089353527777, 0.11156804371056815, -0.8851691518465128, -0.30086627182417736] + v: 0.125 + - phase: A + O: [0.1904884889827469, -0.2955994647758964, -0.5609071888421532, 0.7494786304455029] + v: 0.125 + - phase: A + O: [0.14357497570029057, -0.8835693969815234, -0.07717841484904021, -0.4390157620766682] + v: 0.125 + - phase: A + O: [0.2844319189561915, -0.26616939939342976, -0.09666838434612286, 0.9159189689996324] + v: 0.125 + - phase: A + O: [0.3986740568181236, -0.009203742364828175, 0.3059441820684439, 0.8645070531841439] + v: 0.125 + - phase: A + O: [0.6689177202684424, -0.4911686233022882, 0.011604901743798118, -0.5578241597938561] + v: 0.125 + - phase: A + O: [0.24260666144446064, -0.7362708394584294, 0.555535585659367, 0.3007116091075576] + v: 0.125 + - phase: A + O: [0.8264968442974631, 0.08346800006449286, -0.5223744566842243, 0.19251230177687406] + v: 0.125 + homogenization: RGC + - constituents: + - phase: A + O: [0.5964950823065459, -0.03914167887903209, -0.764986315653944, 0.23970290490697482] + v: 0.125 + - phase: A + O: [0.8080255887350464, 0.369717513486498, 0.4541656797460721, 0.06432063052809132] + v: 0.125 + - phase: A + O: [0.9384315445539018, 0.21987525747326847, 0.26439717033746396, 0.033094465621676714] + v: 0.125 + - phase: A + O: [0.3398146205195601, 0.19417091506154938, 0.7316894460339152, 0.5580808489707294] + v: 0.125 + - phase: A + O: [0.3246633791904426, 0.47411988229633434, 0.8166638950637336, 0.05351737963768399] + v: 0.125 + - phase: A + O: [0.3795740054462351, -0.8471170612518093, 0.01606588144238899, 0.37156176657331014] + v: 0.125 + - phase: A + O: [0.3324564774831089, 0.5840719896629913, 0.33174639891221797, 0.6620248698345199] + v: 0.125 + - phase: A + O: [0.5142455012854315, 0.2115524033587584, -0.25751707686577474, -0.7902417985422785] + v: 0.125 + homogenization: RGC + - constituents: + - phase: A + O: [0.5917389153130538, 0.014487837662305996, 0.1546249763526937, -0.7910286185416618] + v: 0.125 + - phase: A + O: [0.13783983489759696, 0.19317471156196087, 0.9651184311351051, -0.11058989380440512] + v: 0.125 + - phase: A + O: [0.8586698401906224, 0.4868991116453938, -0.04218722674591906, 0.15438781857849304] + v: 0.125 + - phase: A + O: [0.9890970857322728, 0.05214494425397627, -0.1362936500688951, -0.019796482909146578] + v: 0.125 + - phase: A + O: [0.9020138640071081, 0.25707037984798153, 0.30534577613336306, -0.16446813047303388] + v: 0.125 + - phase: A + O: [0.3175523590958934, -0.563302190026314, -0.6464135983869942, -0.40496987760149405] + v: 0.125 + - phase: A + O: [0.752248294923397, -0.30191235437425856, 0.44098110520789696, -0.3853661867764942] + v: 0.125 + - phase: A + O: [0.07541991372724928, 0.21981465995106916, -0.8076581447316371, 0.5419240473835979] + v: 0.125 + homogenization: RGC + - constituents: + - phase: A + O: [0.3735354245257942, -0.19362809638145217, 0.7860424876438061, 0.45289806196843757] + v: 0.125 + - phase: A + O: [0.0320727653734856, -0.8866094001869422, 0.46133859224884993, 0.007862094078374177] + v: 0.125 + - phase: A + O: [0.4437387774582536, 0.2752317472571834, 0.5589728539449029, -0.6441216742466462] + v: 0.125 + - phase: A + O: [0.4189270336917096, -0.7997419507282093, 0.08375652034255085, 0.4217793238031132] + v: 0.125 + - phase: A + O: [0.07774372790271536, 0.8494512032072281, -0.23549059135564554, 0.4657603971191082] + v: 0.125 + - phase: A + O: [0.7403304241587872, 0.5263542470691401, -0.1337642666383515, -0.39619337529526266] + v: 0.125 + - phase: A + O: [0.3862324493674717, -0.7790044035666689, 0.2926767318730596, 0.39789064439798927] + v: 0.125 + - phase: A + O: [0.12717246708576896, 0.2790094483039878, 0.8667779585068207, -0.39328979394229346] + v: 0.125 + homogenization: RGC + +homogenization: + RGC: + N_constituents: 8 + mechanical: + type: RGC + D_alpha: [4.0e-06, 4.0e-06, 2.0e-06] + a_g: [0.0, 0.0, 0.0] + c_alpha: 2.0 + cluster_size: [2, 2, 2] + output: [M, Delta_V, avg_a_dot, max_a_dot] + xi_alpha: 10.0 + +phase: + A: + lattice: cF + mechanical: + output: [F, F_e, F_p, L_p] + elastic: {C_11: 106.75e9, C_12: 60.41e9, C_44: 28.34e9, type: hooke} + plastic: + N_sl: [12] + a_sl: 2.25 + atol_xi: 1.0 + dot_gamma_0_sl: 0.001 + h_0_sl_sl: 75e6 + h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4] + n_sl: 20 + output: [xi_sl] + type: phenopowerlaw + xi_0_sl: [31e6] + xi_inf_sl: [63e6] + B: + lattice: cI + mechanical: + output: [F, P, F_e, F_p, L_p, O] + elastic: {C_11: 106.75e9, C_12: 60.41e9, C_44: 28.34e9, type: hooke} + plastic: + N_sl: [12] + a_sl: 2.25 + atol_xi: 1.0 + dot_gamma_0_sl: 0.001 + h_0_sl_sl: 75e6 + h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4] + n_sl: 20 + output: [xi_sl] + type: phenopowerlaw + xi_0_sl: [31e6] + xi_inf_sl: [63e6] + C: + lattice: cI + mechanical: + output: [F] + elastic: {C_11: 106.75e9, C_12: 60.41e9, C_44: 28.34e9, type: hooke} diff --git a/python/tests/reference/Result/4grains2x4x3.vtr b/python/tests/reference/Result/4grains2x4x3.vtr new file mode 100644 index 000000000..a769a8e55 --- /dev/null +++ b/python/tests/reference/Result/4grains2x4x3.vtr @@ -0,0 +1,30 @@ + + + + + + AAAAAACAAAAAAAAA + + + + + + + + AQAAAACAAADAAAAAIgAAAA==eF5jZIAARjSamQAfXRwdoMszoYnD+DBAyBx0GgYADegAHw== + + + + + AQAAAACAAAAYAAAAEQAAAA==eF5jYEAGD+wh9Ad7AA0uAk8= + + + AQAAAACAAAAoAAAAFwAAAA==eF5jYEAGF+wh9AMo/QJKf7AHADzEBIU= + + + AQAAAACAAAAgAAAAGAAAAA==eF5jYICAUDC4ag+hn9pDRD/YAwBmSwdk + + + + + diff --git a/python/tests/reference/Result/4grains2x4x3_compressionY.hdf5 b/python/tests/reference/Result/4grains2x4x3_compressionY.hdf5 new file mode 100644 index 000000000..c3520c27e Binary files /dev/null and b/python/tests/reference/Result/4grains2x4x3_compressionY.hdf5 differ diff --git a/python/tests/reference/Result/6grains6x7x8_single_phase_tensionY.yaml b/python/tests/reference/Result/6grains6x7x8_single_phase.material.yaml similarity index 100% rename from python/tests/reference/Result/6grains6x7x8_single_phase_tensionY.yaml rename to python/tests/reference/Result/6grains6x7x8_single_phase.material.yaml diff --git a/python/tests/reference/Result/compressionY.yaml b/python/tests/reference/Result/compressionY.yaml new file mode 100644 index 000000000..cfe67e692 --- /dev/null +++ b/python/tests/reference/Result/compressionY.yaml @@ -0,0 +1,17 @@ +--- +solver: + mechanical: spectral_basic + +loadstep: + - boundary_conditions: + mechanical: + dot_F: [x, 0, 0, + 0, -1.0e-3, 0, + 0, 0, x] + P: [0, x, x, + x, x, x, + x, x, 0] + discretization: + t: 5 + N: 10 + f_out: 2 diff --git a/python/tests/reference/Result/get/test_get[0].pbz2 b/python/tests/reference/Result/get/test_get[0].pbz2 new file mode 100644 index 000000000..e4ef0bde9 Binary files /dev/null and b/python/tests/reference/Result/get/test_get[0].pbz2 differ diff --git a/python/tests/reference/Result/get/test_get[1].pbz2 b/python/tests/reference/Result/get/test_get[1].pbz2 new file mode 100644 index 000000000..c1963d4ea Binary files /dev/null and b/python/tests/reference/Result/get/test_get[1].pbz2 differ diff --git a/python/tests/reference/Result/get/test_get[2].pbz2 b/python/tests/reference/Result/get/test_get[2].pbz2 new file mode 100644 index 000000000..c8398ab2f Binary files /dev/null and b/python/tests/reference/Result/get/test_get[2].pbz2 differ diff --git a/python/tests/reference/Result/get/test_get[3].pbz2 b/python/tests/reference/Result/get/test_get[3].pbz2 new file mode 100644 index 000000000..d2eb08545 Binary files /dev/null and b/python/tests/reference/Result/get/test_get[3].pbz2 differ diff --git a/python/tests/reference/Result/get/test_get[4].pbz2 b/python/tests/reference/Result/get/test_get[4].pbz2 new file mode 100644 index 000000000..c10fa916b Binary files /dev/null and b/python/tests/reference/Result/get/test_get[4].pbz2 differ diff --git a/python/tests/reference/Result/get/test_get[5].pbz2 b/python/tests/reference/Result/get/test_get[5].pbz2 new file mode 100644 index 000000000..c1963d4ea Binary files /dev/null and b/python/tests/reference/Result/get/test_get[5].pbz2 differ diff --git a/python/tests/reference/Result/get/test_get[6].pbz2 b/python/tests/reference/Result/get/test_get[6].pbz2 new file mode 100644 index 000000000..4459d2e59 Binary files /dev/null and b/python/tests/reference/Result/get/test_get[6].pbz2 differ diff --git a/python/tests/reference/Result/get/test_get[7].pbz2 b/python/tests/reference/Result/get/test_get[7].pbz2 new file mode 100644 index 000000000..c1e4c31e9 Binary files /dev/null and b/python/tests/reference/Result/get/test_get[7].pbz2 differ diff --git a/python/tests/reference/Result/place/test_place[0].pbz2 b/python/tests/reference/Result/place/test_place[0].pbz2 new file mode 100644 index 000000000..2d1d432b9 Binary files /dev/null and b/python/tests/reference/Result/place/test_place[0].pbz2 differ diff --git a/python/tests/reference/Result/place/test_place[1].pbz2 b/python/tests/reference/Result/place/test_place[1].pbz2 new file mode 100644 index 000000000..c1963d4ea Binary files /dev/null and b/python/tests/reference/Result/place/test_place[1].pbz2 differ diff --git a/python/tests/reference/Result/place/test_place[2].pbz2 b/python/tests/reference/Result/place/test_place[2].pbz2 new file mode 100644 index 000000000..cfc94f1ee Binary files /dev/null and b/python/tests/reference/Result/place/test_place[2].pbz2 differ diff --git a/python/tests/reference/Result/place/test_place[3].pbz2 b/python/tests/reference/Result/place/test_place[3].pbz2 new file mode 100644 index 000000000..96f782ae2 Binary files /dev/null and b/python/tests/reference/Result/place/test_place[3].pbz2 differ diff --git a/python/tests/reference/Result/place/test_place[4].pbz2 b/python/tests/reference/Result/place/test_place[4].pbz2 new file mode 100644 index 000000000..f2f744121 Binary files /dev/null and b/python/tests/reference/Result/place/test_place[4].pbz2 differ diff --git a/python/tests/reference/Result/place/test_place[5].pbz2 b/python/tests/reference/Result/place/test_place[5].pbz2 new file mode 100644 index 000000000..c1963d4ea Binary files /dev/null and b/python/tests/reference/Result/place/test_place[5].pbz2 differ diff --git a/python/tests/reference/Result/place/test_place[6].pbz2 b/python/tests/reference/Result/place/test_place[6].pbz2 new file mode 100644 index 000000000..2844addcb Binary files /dev/null and b/python/tests/reference/Result/place/test_place[6].pbz2 differ diff --git a/python/tests/reference/Result/place/test_place[7].pbz2 b/python/tests/reference/Result/place/test_place[7].pbz2 new file mode 100644 index 000000000..cd4f098a1 Binary files /dev/null and b/python/tests/reference/Result/place/test_place[7].pbz2 differ diff --git a/python/tests/reference/Result/save_VTK/test_vtk[0-0-0].md5 b/python/tests/reference/Result/save_VTK/test_vtk[0-0-0].md5 new file mode 100644 index 000000000..838037bb2 --- /dev/null +++ b/python/tests/reference/Result/save_VTK/test_vtk[0-0-0].md5 @@ -0,0 +1 @@ +3b83384def67552ab7dd211efc0d54fd \ No newline at end of file diff --git a/python/tests/reference/Result/save_VTK/test_vtk[0-0-1].md5 b/python/tests/reference/Result/save_VTK/test_vtk[0-0-1].md5 new file mode 100644 index 000000000..7ceffc337 --- /dev/null +++ b/python/tests/reference/Result/save_VTK/test_vtk[0-0-1].md5 @@ -0,0 +1 @@ +c32c86ed50dbb39a93ca2a2ebe47d9cb \ No newline at end of file diff --git a/python/tests/reference/Result/save_VTK/test_vtk[0-0-2].md5 b/python/tests/reference/Result/save_VTK/test_vtk[0-0-2].md5 new file mode 100644 index 000000000..f5b7daec3 --- /dev/null +++ b/python/tests/reference/Result/save_VTK/test_vtk[0-0-2].md5 @@ -0,0 +1 @@ +ead4f6fcaff174fddc041d701e54ac60 \ No newline at end of file diff --git a/python/tests/reference/Result/save_VTK/test_vtk[1-0-0].md5 b/python/tests/reference/Result/save_VTK/test_vtk[1-0-0].md5 new file mode 100644 index 000000000..df00a513f --- /dev/null +++ b/python/tests/reference/Result/save_VTK/test_vtk[1-0-0].md5 @@ -0,0 +1 @@ +bde8b728110c2c05a6a4740f7c5f9c06 \ No newline at end of file diff --git a/python/tests/reference/Result/save_VTK/test_vtk[1-0-1].md5 b/python/tests/reference/Result/save_VTK/test_vtk[1-0-1].md5 new file mode 100644 index 000000000..35c577900 --- /dev/null +++ b/python/tests/reference/Result/save_VTK/test_vtk[1-0-1].md5 @@ -0,0 +1 @@ +e09bfa9248283fc390003ad28d15d36e \ No newline at end of file diff --git a/python/tests/reference/Result/save_VTK/test_vtk[1-0-2].md5 b/python/tests/reference/Result/save_VTK/test_vtk[1-0-2].md5 new file mode 100644 index 000000000..d5874d88b --- /dev/null +++ b/python/tests/reference/Result/save_VTK/test_vtk[1-0-2].md5 @@ -0,0 +1 @@ +3f21254164f96de8ee4a28249ae72cc6 \ No newline at end of file diff --git a/python/tests/test_ConfigMaterial.py b/python/tests/test_ConfigMaterial.py index 98a557ae8..c968f429a 100644 --- a/python/tests/test_ConfigMaterial.py +++ b/python/tests/test_ConfigMaterial.py @@ -136,9 +136,9 @@ class TestConfigMaterial: def test_load_DREAM3D_reference(self,tmp_path,ref_path,update): cur = ConfigMaterial.load_DREAM3D(ref_path/'measured.dream3d') - ref = ConfigMaterial.load(ref_path/'measured.material_yaml') + ref = ConfigMaterial.load(ref_path/'measured.material.yaml') if update: - cur.save(ref_path/'measured.material_yaml') + cur.save(ref_path/'measured.material.yaml') for i,m in enumerate(ref['material']): assert Rotation(m['constituents'][0]['O']).isclose(Rotation(cur['material'][i]['constituents'][0]['O'])) assert cur.is_valid and cur['phase'] == ref['phase'] and cur['homogenization'] == ref['homogenization'] diff --git a/python/tests/test_Grid.py b/python/tests/test_Grid.py index 24d48bd56..9d82bde78 100644 --- a/python/tests/test_Grid.py +++ b/python/tests/test_Grid.py @@ -14,8 +14,7 @@ from damask import grid_filters def grid_equal(a,b): return np.all(a.material == b.material) and \ np.all(a.cells == b.cells) and \ - np.allclose(a.size, b.size) and \ - str(a.diff(b)) == str(b.diff(a)) + np.allclose(a.size, b.size) @pytest.fixture def default(): @@ -42,13 +41,9 @@ class TestGrid: def _patch_datetime_now(self, patch_datetime_now): print('patched datetime.datetime.now') - def test_diff_equal(self,default): - assert str(default.diff(default)) == '' - - def test_diff_not_equal(self,default): - new = Grid(default.material[1:,1:,1:]+1,default.size*.9,np.ones(3)-default.origin,comments=['modified']) - assert str(default.diff(new)) != '' + def test_equal(self,default): + assert default == default def test_repr(self,default): print(default) diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index 3ab5338ea..a659ce056 100644 --- a/python/tests/test_Result.py +++ b/python/tests/test_Result.py @@ -1,12 +1,14 @@ +import bz2 +import pickle import time import shutil import os import sys +import hashlib from datetime import datetime import pytest import numpy as np -import h5py from damask import Result from damask import Rotation @@ -21,8 +23,7 @@ def default(tmp_path,ref_path): fname = '12grains6x7x8_tensionY.hdf5' shutil.copy(ref_path/fname,tmp_path) f = Result(tmp_path/fname) - f.view('times',20.0) - return f + return f.view('times',20.0) @pytest.fixture def single_phase(tmp_path,ref_path): @@ -36,6 +37,17 @@ def ref_path(ref_path_base): """Directory containing reference results.""" return ref_path_base/'Result' +def dict_equal(d1, d2): + for k in d1: + if (k not in d2): + return False + else: + if type(d1[k]) is dict: + return dict_equal(d1[k],d2[k]) + else: + if not np.allclose(d1[k],d2[k]): + return False + return True class TestResult: @@ -44,51 +56,39 @@ class TestResult: def test_view_all(self,default): - default.view('increments',True) - a = default.get_dataset_location('F') - default.view('increments','*') - b = default.get_dataset_location('F') - default.view('increments',default.increments_in_range(0,np.iinfo(int).max)) - c = default.get_dataset_location('F') + a = default.view('increments',True).get('F') - default.view('times',True) - d = default.get_dataset_location('F') - default.view('times','*') - e = default.get_dataset_location('F') - default.view('times',default.times_in_range(0.0,np.inf)) - f = default.get_dataset_location('F') - assert a == b == c == d == e ==f + assert dict_equal(a,default.view('increments','*').get('F')) + assert dict_equal(a,default.view('increments',default.increments_in_range(0,np.iinfo(int).max)).get('F')) + + assert dict_equal(a,default.view('times',True).get('F')) + assert dict_equal(a,default.view('times','*').get('F')) + assert dict_equal(a,default.view('times',default.times_in_range(0.0,np.inf)).get('F')) @pytest.mark.parametrize('what',['increments','times','phases']) # ToDo: discuss homogenizations def test_view_none(self,default,what): - default.view(what,False) - a = default.get_dataset_location('F') - default.view(what,[]) - b = default.get_dataset_location('F') + a = default.view(what,False).get('F') + b = default.view(what,[]).get('F') - assert a == b == [] + assert a == b == {} @pytest.mark.parametrize('what',['increments','times','phases']) # ToDo: discuss homogenizations def test_view_more(self,default,what): - default.view(what,False) - default.view_more(what,'*') - a = default.get_dataset_location('F') + empty = default.view(what,False) - default.view(what,True) - b = default.get_dataset_location('F') + a = empty.view_more(what,'*').get('F') + b = empty.view_more(what,True).get('F') - assert a == b + assert dict_equal(a,b) @pytest.mark.parametrize('what',['increments','times','phases']) # ToDo: discuss homogenizations def test_view_less(self,default,what): - default.view(what,True) - default.view_less(what,'*') - a = default.get_dataset_location('F') + full = default.view(what,True) - default.view(what,False) - b = default.get_dataset_location('F') + a = full.view_less(what,'*').get('F') + b = full.view_less(what,True).get('F') - assert a == b == [] + assert a == b == {} def test_view_invalid(self,default): with pytest.raises(AttributeError): @@ -96,10 +96,8 @@ class TestResult: def test_add_absolute(self,default): default.add_absolute('F_e') - loc = {'F_e': default.get_dataset_location('F_e'), - '|F_e|': default.get_dataset_location('|F_e|')} - in_memory = np.abs(default.read_dataset(loc['F_e'],0)) - in_file = default.read_dataset(loc['|F_e|'],0) + in_memory = np.abs(default.place('F_e')) + in_file = default.place('|F_e|') assert np.allclose(in_memory,in_file) @pytest.mark.parametrize('mode',['direct','function']) @@ -115,77 +113,59 @@ class TestResult: default.enable_user_function(f.my_func) default.add_calculation('x','my_func(#F#)','-','my notes') - loc = {'F': default.get_dataset_location('F'), - 'x': default.get_dataset_location('x')} - in_memory = 2.0*np.abs(default.read_dataset(loc['F'],0))-1.0 - in_file = default.read_dataset(loc['x'],0) + in_memory = 2.0*np.abs(default.place('F'))-1.0 + in_file = default.place('x') assert np.allclose(in_memory,in_file) def test_add_stress_Cauchy(self,default): default.add_stress_Cauchy('P','F') - loc = {'F': default.get_dataset_location('F'), - 'P': default.get_dataset_location('P'), - 'sigma':default.get_dataset_location('sigma')} - in_memory = mechanics.stress_Cauchy(default.read_dataset(loc['P'],0), - default.read_dataset(loc['F'],0)) - in_file = default.read_dataset(loc['sigma'],0) + in_memory = mechanics.stress_Cauchy(default.place('P'), default.place('F')) + in_file = default.place('sigma') assert np.allclose(in_memory,in_file) def test_add_determinant(self,default): default.add_determinant('P') - loc = {'P': default.get_dataset_location('P'), - 'det(P)':default.get_dataset_location('det(P)')} - in_memory = np.linalg.det(default.read_dataset(loc['P'],0)).reshape(-1,1) - in_file = default.read_dataset(loc['det(P)'],0) + in_memory = np.linalg.det(default.place('P')) + in_file = default.place('det(P)') assert np.allclose(in_memory,in_file) def test_add_deviator(self,default): default.add_deviator('P') - loc = {'P' :default.get_dataset_location('P'), - 's_P':default.get_dataset_location('s_P')} - in_memory = tensor.deviatoric(default.read_dataset(loc['P'],0)) - in_file = default.read_dataset(loc['s_P'],0) + in_memory = tensor.deviatoric(default.place('P')) + in_file = default.place('s_P') assert np.allclose(in_memory,in_file) @pytest.mark.parametrize('eigenvalue,function',[('max',np.amax),('min',np.amin)]) def test_add_eigenvalue(self,default,eigenvalue,function): default.add_stress_Cauchy('P','F') default.add_eigenvalue('sigma',eigenvalue) - loc = {'sigma' :default.get_dataset_location('sigma'), - 'lambda':default.get_dataset_location(f'lambda_{eigenvalue}(sigma)')} - in_memory = function(tensor.eigenvalues(default.read_dataset(loc['sigma'],0)),axis=1,keepdims=True) - in_file = default.read_dataset(loc['lambda'],0) + in_memory = function(tensor.eigenvalues(default.place('sigma')),axis=1) + in_file = default.place(f'lambda_{eigenvalue}(sigma)') assert np.allclose(in_memory,in_file) @pytest.mark.parametrize('eigenvalue,idx',[('max',2),('mid',1),('min',0)]) def test_add_eigenvector(self,default,eigenvalue,idx): default.add_stress_Cauchy('P','F') default.add_eigenvector('sigma',eigenvalue) - loc = {'sigma' :default.get_dataset_location('sigma'), - 'v(sigma)':default.get_dataset_location(f'v_{eigenvalue}(sigma)')} - in_memory = tensor.eigenvectors(default.read_dataset(loc['sigma'],0))[:,idx] - in_file = default.read_dataset(loc['v(sigma)'],0) + in_memory = tensor.eigenvectors(default.place('sigma'))[:,idx] + in_file = default.place(f'v_{eigenvalue}(sigma)') assert np.allclose(in_memory,in_file) @pytest.mark.parametrize('d',[[1,0,0],[0,1,0],[0,0,1]]) def test_add_IPF_color(self,default,d): default.add_IPF_color(d,'O') - loc = {'O': default.get_dataset_location('O'), - 'color': default.get_dataset_location('IPFcolor_[{} {} {}]'.format(*d))} - qu = default.read_dataset(loc['O']).view(np.double).squeeze() - crystal_structure = default._get_attribute(default.get_dataset_location('O')[0],'lattice') + qu = default.place('O') + crystal_structure = qu.dtype.metadata['lattice'] c = Orientation(rotation=qu,lattice=crystal_structure) in_memory = np.uint8(c.IPF_color(np.array(d))*255) - in_file = default.read_dataset(loc['color']) + in_file = default.place('IPFcolor_({} {} {})'.format(*d)) assert np.allclose(in_memory,in_file) def test_add_maximum_shear(self,default): default.add_stress_Cauchy('P','F') default.add_maximum_shear('sigma') - loc = {'sigma' :default.get_dataset_location('sigma'), - 'max_shear(sigma)':default.get_dataset_location('max_shear(sigma)')} - in_memory = mechanics.maximum_shear(default.read_dataset(loc['sigma'],0)).reshape(-1,1) - in_file = default.read_dataset(loc['max_shear(sigma)'],0) + in_memory = mechanics.maximum_shear(default.place('sigma')) + in_file = default.place('max_shear(sigma)') assert np.allclose(in_memory,in_file) def test_add_Mises_strain(self,default): @@ -194,26 +174,22 @@ class TestResult: default.add_strain('F',t,m) label = f'epsilon_{t}^{m}(F)' default.add_equivalent_Mises(label) - loc = {label :default.get_dataset_location(label), - label+'_vM':default.get_dataset_location(label+'_vM')} - in_memory = mechanics.equivalent_strain_Mises(default.read_dataset(loc[label],0)).reshape(-1,1) - in_file = default.read_dataset(loc[label+'_vM'],0) + in_memory = mechanics.equivalent_strain_Mises(default.place(label)) + in_file = default.place(label+'_vM') assert np.allclose(in_memory,in_file) def test_add_Mises_stress(self,default): default.add_stress_Cauchy('P','F') default.add_equivalent_Mises('sigma') - loc = {'sigma' :default.get_dataset_location('sigma'), - 'sigma_vM':default.get_dataset_location('sigma_vM')} - in_memory = mechanics.equivalent_stress_Mises(default.read_dataset(loc['sigma'],0)).reshape(-1,1) - in_file = default.read_dataset(loc['sigma_vM'],0) + in_memory = mechanics.equivalent_stress_Mises(default.place('sigma')) + in_file = default.place('sigma_vM') assert np.allclose(in_memory,in_file) def test_add_Mises_invalid(self,default): default.add_stress_Cauchy('P','F') default.add_calculation('sigma_y','#sigma#',unit='y') default.add_equivalent_Mises('sigma_y') - assert default.get_dataset_location('sigma_y_vM') == [] + assert default.get('sigma_y_vM') == {} def test_add_Mises_stress_strain(self,default): default.add_stress_Cauchy('P','F') @@ -221,26 +197,20 @@ class TestResult: default.add_calculation('sigma_x','#sigma#',unit='x') default.add_equivalent_Mises('sigma_y',kind='strain') default.add_equivalent_Mises('sigma_x',kind='stress') - loc = {'y' :default.get_dataset_location('sigma_y_vM'), - 'x' :default.get_dataset_location('sigma_x_vM')} - assert not np.allclose(default.read_dataset(loc['y'],0),default.read_dataset(loc['x'],0)) + assert not np.allclose(default.place('sigma_y_vM'),default.place('sigma_x_vM')) - def test_add_norm(self,default): - default.add_norm('F',1) - loc = {'F': default.get_dataset_location('F'), - '|F|_1':default.get_dataset_location('|F|_1')} - in_memory = np.linalg.norm(default.read_dataset(loc['F'],0),ord=1,axis=(1,2),keepdims=True) - in_file = default.read_dataset(loc['|F|_1'],0) + @pytest.mark.parametrize('ord',[1,2]) + @pytest.mark.parametrize('dataset,axis',[('F',(1,2)),('xi_sl',(1,))]) + def test_add_norm(self,default,ord,dataset,axis): + default.add_norm(dataset,ord) + in_memory = np.linalg.norm(default.place(dataset),ord=ord,axis=axis,keepdims=True) + in_file = default.place(f'|{dataset}|_{ord}') assert np.allclose(in_memory,in_file) def test_add_stress_second_Piola_Kirchhoff(self,default): default.add_stress_second_Piola_Kirchhoff('P','F') - loc = {'F':default.get_dataset_location('F'), - 'P':default.get_dataset_location('P'), - 'S':default.get_dataset_location('S')} - in_memory = mechanics.stress_second_Piola_Kirchhoff(default.read_dataset(loc['P'],0), - default.read_dataset(loc['F'],0)) - in_file = default.read_dataset(loc['S'],0) + in_memory = mechanics.stress_second_Piola_Kirchhoff(default.place('P'),default.place('F')) + in_file = default.place('S') assert np.allclose(in_memory,in_file) @pytest.mark.skip(reason='requires rework of lattice.f90') @@ -248,30 +218,24 @@ class TestResult: def test_add_pole(self,default,polar): pole = np.array([1.,0.,0.]) default.add_pole('O',pole,polar) - loc = {'O': default.get_dataset_location('O'), - 'pole': default.get_dataset_location('p^{}_[1 0 0)'.format(u'rφ' if polar else 'xy'))} - rot = Rotation(default.read_dataset(loc['O']).view(np.double)) + rot = Rotation(default.place('O')) rotated_pole = rot * np.broadcast_to(pole,rot.shape+(3,)) xy = rotated_pole[:,0:2]/(1.+abs(pole[2])) in_memory = 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])]) - in_file = default.read_dataset(loc['pole']) + in_file = default.place('p^{}_[1 0 0)'.format(u'rφ' if polar else 'xy')) assert np.allclose(in_memory,in_file) def test_add_rotation(self,default): default.add_rotation('F') - loc = {'F': default.get_dataset_location('F'), - 'R(F)': default.get_dataset_location('R(F)')} - in_memory = mechanics.rotation(default.read_dataset(loc['F'],0)).as_matrix() - in_file = default.read_dataset(loc['R(F)'],0) + in_memory = mechanics.rotation(default.place('F')).as_matrix() + in_file = default.place('R(F)') assert np.allclose(in_memory,in_file) def test_add_spherical(self,default): default.add_spherical('P') - loc = {'P': default.get_dataset_location('P'), - 'p_P': default.get_dataset_location('p_P')} - in_memory = tensor.spherical(default.read_dataset(loc['P'],0),False).reshape(-1,1) - in_file = default.read_dataset(loc['p_P'],0) + in_memory = tensor.spherical(default.place('P'),False) + in_file = default.place('p_P') assert np.allclose(in_memory,in_file) def test_add_strain(self,default): @@ -279,26 +243,20 @@ class TestResult: m = np.random.random()*2.0 - 1.0 default.add_strain('F',t,m) label = f'epsilon_{t}^{m}(F)' - loc = {'F': default.get_dataset_location('F'), - label: default.get_dataset_location(label)} - in_memory = mechanics.strain(default.read_dataset(loc['F'],0),t,m) - in_file = default.read_dataset(loc[label],0) + in_memory = mechanics.strain(default.place('F'),t,m) + in_file = default.place(label) assert np.allclose(in_memory,in_file) def test_add_stretch_right(self,default): default.add_stretch_tensor('F','U') - loc = {'F': default.get_dataset_location('F'), - 'U(F)': default.get_dataset_location('U(F)')} - in_memory = mechanics.stretch_right(default.read_dataset(loc['F'],0)) - in_file = default.read_dataset(loc['U(F)'],0) + in_memory = mechanics.stretch_right(default.place('F')) + in_file = default.place('U(F)') assert np.allclose(in_memory,in_file) def test_add_stretch_left(self,default): default.add_stretch_tensor('F','V') - loc = {'F': default.get_dataset_location('F'), - 'V(F)': default.get_dataset_location('V(F)')} - in_memory = mechanics.stretch_left(default.read_dataset(loc['F'],0)) - in_file = default.read_dataset(loc['V(F)'],0) + in_memory = mechanics.stretch_left(default.place('F')) + in_file = default.place('V(F)') assert np.allclose(in_memory,in_file) def test_add_invalid(self,default): @@ -307,48 +265,40 @@ class TestResult: @pytest.mark.parametrize('overwrite',['off','on']) def test_add_overwrite(self,default,overwrite): - default.view('times',default.times_in_range(0,np.inf)[-1]) + last = default.view('times',default.times_in_range(0,np.inf)[-1]) - default.add_stress_Cauchy() - loc = default.get_dataset_location('sigma') - with h5py.File(default.fname,'r') as f: - # h5py3 compatibility - try: - created_first = f[loc[0]].attrs['created'].decode() - except AttributeError: - created_first = f[loc[0]].attrs['created'] + last.add_stress_Cauchy() + + created_first = last.place('sigma').dtype.metadata['created'] created_first = datetime.strptime(created_first,'%Y-%m-%d %H:%M:%S%z') if overwrite == 'on': - default.allow_modification() + last = last.allow_modification() else: - default.disallow_modification() + last = last.disallow_modification() time.sleep(2.) try: - default.add_calculation('sigma','#sigma#*0.0+311.','not the Cauchy stress') + last.add_calculation('sigma','#sigma#*0.0+311.','not the Cauchy stress') except ValueError: pass - with h5py.File(default.fname,'r') as f: - # h5py3 compatibility - try: - created_second = f[loc[0]].attrs['created'].decode() - except AttributeError: - created_second = f[loc[0]].attrs['created'] + + created_second = last.place('sigma').dtype.metadata['created'] created_second = datetime.strptime(created_second,'%Y-%m-%d %H:%M:%S%z') + if overwrite == 'on': - assert created_first < created_second and np.allclose(default.read_dataset(loc),311.) + assert created_first < created_second and np.allclose(last.place('sigma'),311.) else: - assert created_first == created_second and not np.allclose(default.read_dataset(loc),311.) + assert created_first == created_second and not np.allclose(last.place('sigma'),311.) @pytest.mark.parametrize('allowed',['off','on']) def test_rename(self,default,allowed): if allowed == 'on': - F = default.read_dataset(default.get_dataset_location('F')) - default.allow_modification() + F = default.place('F') + default = default.allow_modification() default.rename('F','new_name') - assert np.all(F == default.read_dataset(default.get_dataset_location('new_name'))) - default.disallow_modification() + assert np.all(F == default.place('new_name')) + default = default.disallow_modification() with pytest.raises(PermissionError): default.rename('P','another_new_name') @@ -363,10 +313,30 @@ class TestResult: b = default.coordinates0_node.reshape(tuple(default.cells+1)+(3,),order='F') assert np.allclose(a,b) - @pytest.mark.parametrize('output',['F',[],['F','P']]) - def test_vtk(self,tmp_path,default,output): + # need to wait for writing in parallel, output order might change if select more then one + @pytest.mark.parametrize('output',['F','*',['P']],ids=range(3)) + @pytest.mark.parametrize('fname',['12grains6x7x8_tensionY.hdf5'],ids=range(1)) + @pytest.mark.parametrize('inc',[4,0],ids=range(2)) + def test_vtk(self,request,tmp_path,ref_path,update,output,fname,inc): + result = Result(ref_path/fname).view('increments',inc) os.chdir(tmp_path) - default.save_VTK(output) + result.save_VTK(output) + fname = fname.split('.')[0]+f'_inc{(inc if type(inc) == int else inc[0]):0>2}.vtr' + last = '' + for i in range(10): + if os.path.isfile(tmp_path/fname): + with open(fname) as f: + cur = hashlib.md5(f.read().encode()).hexdigest() + if cur == last: + break + else: + last = cur + time.sleep(.5) + if update: + with open((ref_path/'save_VTK'/request.node.name).with_suffix('.md5'),'w') as f: + f.write(cur) + with open((ref_path/'save_VTK'/request.node.name).with_suffix('.md5')) as f: + assert cur == f.read() @pytest.mark.parametrize('mode',['point','cell']) def test_vtk_mode(self,tmp_path,single_phase,mode): @@ -387,3 +357,52 @@ class TestResult: def test_XDMF_invalid(self,default): with pytest.raises(TypeError): default.save_XDMF() + + @pytest.mark.parametrize('view,output,flatten,prune', + [({},['F','P','F','L_p','F_e','F_p'],True,True), + ({'increments':3},'F',True,True), + ({'increments':[1,8,3,4,5,6,7]},['F','P'],True,True), + ({'phases':['A','B']},['F','P'],True,True), + ({'phases':['A','C'],'homogenizations':False},['F','P','O'],True,True), + ({'phases':False,'homogenizations':False},['F','P','O'],True,True), + ({'phases':False},['Delta_V'],True,True), + ({},['u_p','u_n'],False,False)], + ids=list(range(8))) + def test_get(self,update,request,ref_path,view,output,flatten,prune): + result = Result(ref_path/'4grains2x4x3_compressionY.hdf5') + for key,value in view.items(): + result = result.view(key,value) + + fname = request.node.name + cur = result.get(output,flatten,prune) + if update: + with bz2.BZ2File((ref_path/'get'/fname).with_suffix('.pbz2'),'w') as f: + pickle.dump(cur,f) + + with bz2.BZ2File((ref_path/'get'/fname).with_suffix('.pbz2')) as f: + assert dict_equal(cur,pickle.load(f)) + + + @pytest.mark.parametrize('view,output,flatten,constituents,prune', + [({},['F','P','F','L_p','F_e','F_p'],True,True,None), + ({'increments':3},'F',True,True,[0,1,2,3,4,5,6,7]), + ({'increments':[1,8,3,4,5,6,7]},['F','P'],True,True,1), + ({'phases':['A','B']},['F','P'],True,True,[1,2]), + ({'phases':['A','C'],'homogenizations':False},['F','P','O'],True,True,[0,7]), + ({'phases':False,'homogenizations':False},['F','P','O'],True,True,[1,2,3,4]), + ({'phases':False},['Delta_V'],True,True,[1,2,4]), + ({},['u_p','u_n'],False,False,None)], + ids=list(range(8))) + def test_place(self,update,request,ref_path,view,output,flatten,prune,constituents): + result = Result(ref_path/'4grains2x4x3_compressionY.hdf5') + for key,value in view.items(): + result = result.view(key,value) + + fname = request.node.name + cur = result.place(output,flatten,prune,constituents) + if update: + with bz2.BZ2File((ref_path/'place'/fname).with_suffix('.pbz2'),'w') as f: + pickle.dump(cur,f) + + with bz2.BZ2File((ref_path/'place'/fname).with_suffix('.pbz2')) as f: + assert dict_equal(cur,pickle.load(f)) diff --git a/python/tests/test_VTK.py b/python/tests/test_VTK.py index a3cba354f..4328dc810 100644 --- a/python/tests/test_VTK.py +++ b/python/tests/test_VTK.py @@ -4,6 +4,7 @@ import time import pytest import numpy as np +import numpy.ma as ma from damask import VTK from damask import grid_filters @@ -134,6 +135,15 @@ class TestVTK: with pytest.raises(TypeError): default.add('invalid_type','valid') + def test_add_masked(self,default): + data = np.random.rand(5*6*7,3) + masked = ma.MaskedArray(data,mask=data<.4,fill_value=42.) + default.add(masked,'D') + result_masked = str(default) + default.add(np.where(masked.mask,masked.fill_value,masked),'D') + assert result_masked == str(default) + + def test_comments(self,tmp_path,default): default.add_comments(['this is a comment']) default.save(tmp_path/'with_comments',parallel=False) diff --git a/python/tests/test_util.py b/python/tests/test_util.py index b96908a9a..3ed3bf908 100644 --- a/python/tests/test_util.py +++ b/python/tests/test_util.py @@ -136,3 +136,25 @@ class TestUtil: else: with pytest.raises(ValueError): util.DREAM3D_cell_data_group(tmp_path/'cell_data_group.dream3d') + + + @pytest.mark.parametrize('full,reduced',[({}, {}), + ({'A':{}}, {}), + ({'A':{'B':{}}}, {}), + ({'A':{'B':'C'}},)*2, + ({'A':{'B':{},'C':'D'}}, {'A':{'C':'D'}})]) + def test_prune(self,full,reduced): + assert util.dict_prune(full) == reduced + + + @pytest.mark.parametrize('full,reduced',[({}, {}), + ({'A':{}}, {}), + ({'A':'F'}, 'F'), + ({'A':{'B':{}}}, {}), + ({'A':{'B':'C'}}, 'C'), + ({'A':1,'B':2},)*2, + ({'A':{'B':'C','D':'E'}}, {'B':'C','D':'E'}), + ({'B':'C','D':'E'},)*2, + ({'A':{'B':{},'C':'D'}}, {'B':{},'C':'D'})]) + def test_flatten(self,full,reduced): + assert util.dict_flatten(full) == reduced diff --git a/src/IO.f90 b/src/IO.f90 index b6dbf5664..9ea35503b 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -604,7 +604,7 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg) msg = 'read only the first document' case default msg = 'unknown warning number' - end select + end select !$OMP CRITICAL (write2out) write(IO_STDERR,'(/,a)') ' ┌'//IO_DIVIDER//'┐' @@ -658,7 +658,7 @@ subroutine selfTest if(any([1,1,1] /= IO_stringPos('a'))) error stop 'IO_stringPos' if(any([2,2,3,5,5] /= IO_stringPos(' aa b'))) error stop 'IO_stringPos' - str=' 1.0 xxx' + str = ' 1.0 xxx' chunkPos = IO_stringPos(str) if(dNeq(1.0_pReal,IO_floatValue(str,chunkPos,1))) error stop 'IO_floatValue' diff --git a/src/YAML_types.f90 b/src/YAML_types.f90 index 3d83831e6..2a6bd64f9 100644 --- a/src/YAML_types.f90 +++ b/src/YAML_types.f90 @@ -266,7 +266,7 @@ subroutine selfTest if(any(dNeq(l2%get_as1dFloat(1),[2.0_pReal,3.0_pReal]))) error stop 'byIndex_as1dFloat' call l2%append(l3) x = l2%as2dFloat() - if(x(2,1)/= 4.0_pReal) error stop 'byKey_as2dFloat' + if(dNeq(x(2,1),4.0_pReal)) error stop 'byKey_as2dFloat' if(any(dNeq(pack(l2%as2dFloat(),.true.),& [2.0_pReal,4.0_pReal,3.0_pReal,5.0_pReal]))) error stop 'byKey_as2dFloat' n => l2 diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index 585f242b5..1eac4f8fc 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -47,5 +47,8 @@ #include "homogenization_mechanical_isostrain.f90" #include "homogenization_mechanical_RGC.f90" #include "homogenization_thermal.f90" +#include "homogenization_thermal_pass.f90" +#include "homogenization_thermal_isotemperature.f90" #include "homogenization_damage.f90" +#include "homogenization_damage_pass.f90" #include "CPFEM.f90" diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 6f99365f7..002a15708 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -196,7 +196,7 @@ function grid_damage_spectral_solution(timeinc) result(solution) ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 - call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),ce) + call homogenization_set_phi(phi_current(i,j,k),ce) enddo; enddo; enddo call VecMin(solution_vec,devNull,phi_min,ierr); CHKERRQ(ierr) @@ -233,7 +233,7 @@ subroutine grid_damage_spectral_forward(cutBack) call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 - call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),ce) + call homogenization_set_phi(phi_current(i,j,k),ce) enddo; enddo; enddo else phi_lastInc = phi_current @@ -259,7 +259,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) PetscObject :: dummy PetscErrorCode :: ierr integer :: i, j, k, ce - real(pReal) :: phiDot, dPhiDot_dPhi, mobility + real(pReal) :: phiDot, mobility phi_current = x_scal !-------------------------------------------------------------------------------------------------- @@ -281,7 +281,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 - call damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi_current(i,j,k),ce) + call damage_nonlocal_getSourceAndItsTangent(phiDot, phi_current(i,j,k),ce) mobility = damage_nonlocal_getMobility(ce) scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + phiDot) & + mobility*(phi_lastInc(i,j,k) - phi_current(i,j,k)) & diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index dd431ad9b..14e2affb9 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -282,9 +282,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) ce = ce + 1 call thermal_conduction_getSource(Tdot,1,ce) scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + Tdot) & - + thermal_conduction_getMassDensity (ce)* & - thermal_conduction_getSpecificHeat(ce)*(T_lastInc(i,j,k) - & - T_current(i,j,k))& + + homogenization_thermal_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) & + mu_ref*T_current(i,j,k) enddo; enddo; enddo @@ -314,7 +312,7 @@ subroutine updateReference do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 K_ref = K_ref + thermal_conduction_getConductivity(ce) - mu_ref = mu_ref + thermal_conduction_getMassDensity(ce)* thermal_conduction_getSpecificHeat(ce) + mu_ref = mu_ref + homogenization_thermal_mu_T(ce) enddo; enddo; enddo K_ref = K_ref*wgt call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 7dfaf5b37..4ce75feca 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -39,8 +39,6 @@ module homogenization thermal_type !< thermal transport model integer(kind(DAMAGE_none_ID)), dimension(:), allocatable :: & damage_type !< nonlocal damage model - integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:), allocatable :: & - homogenization_type !< type of each homogenization type, private :: tNumerics_damage real(pReal) :: & @@ -117,6 +115,16 @@ module homogenization integer, intent(in) :: ho end subroutine mechanical_results + module subroutine damage_results(ho,group) + integer, intent(in) :: ho + character(len=*), intent(in) :: group + end subroutine damage_results + + module subroutine thermal_results(ho,group) + integer, intent(in) :: ho + character(len=*), intent(in) :: group + end subroutine thermal_results + module function mechanical_updateState(subdt,subF,ce) result(doneAndHappy) real(pReal), intent(in) :: & subdt !< current time step @@ -133,26 +141,16 @@ module homogenization real(pReal), dimension(3,3) :: K end function thermal_conduction_getConductivity - module function thermal_conduction_getSpecificHeat(ce) result(c_P) + module function homogenization_thermal_mu_T(ce) result(mu_T) integer, intent(in) :: ce - real(pReal) :: c_P - end function thermal_conduction_getSpecificHeat - - module function thermal_conduction_getMassDensity(ce) result(rho) - integer, intent(in) :: ce - real(pReal) :: rho - end function thermal_conduction_getMassDensity + real(pReal) :: mu_T + end function homogenization_thermal_mu_T module subroutine homogenization_thermal_setField(T,dot_T, ce) integer, intent(in) :: ce real(pReal), intent(in) :: T, dot_T end subroutine homogenization_thermal_setField - module subroutine thermal_conduction_results(ho,group) - integer, intent(in) :: ho - character(len=*), intent(in) :: group - end subroutine thermal_conduction_results - module function homogenization_thermal_T(ce) result(T) integer, intent(in) :: ce real(pReal) :: T @@ -170,37 +168,31 @@ module homogenization real(pReal) :: M end function damage_nonlocal_getMobility - module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ce) + module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, phi, ce) integer, intent(in) :: ce - real(pReal), intent(in) :: & + real(pReal), intent(in) :: & phi - real(pReal) :: & - phiDot, dPhiDot_dPhi + real(pReal), intent(out) :: & + phiDot end subroutine damage_nonlocal_getSourceAndItsTangent - module subroutine damage_nonlocal_putNonLocalDamage(phi,ce) + module subroutine homogenization_set_phi(phi,ce) integer, intent(in) :: ce real(pReal), intent(in) :: & phi - end subroutine damage_nonlocal_putNonLocalDamage - - module subroutine damage_nonlocal_results(ho,group) - integer, intent(in) :: ho - character(len=*), intent(in) :: group - end subroutine damage_nonlocal_results + end subroutine homogenization_set_phi end interface public :: & homogenization_init, & materialpoint_stressAndItsTangent, & - thermal_conduction_getSpecificHeat, & + homogenization_thermal_mu_T, & thermal_conduction_getConductivity, & - thermal_conduction_getMassDensity, & thermal_conduction_getSource, & damage_nonlocal_getMobility, & damage_nonlocal_getSourceAndItsTangent, & - damage_nonlocal_putNonLocalDamage, & + homogenization_set_phi, & homogenization_thermal_setfield, & homogenization_thermal_T, & homogenization_forward, & @@ -211,7 +203,6 @@ module homogenization DAMAGE_NONLOCAL_ID public :: & - damage_nonlocal_init, & damage_nonlocal_getDiffusion contains @@ -242,8 +233,6 @@ subroutine homogenization_init() call mechanical_init(num_homog) call thermal_init() call damage_init() - call damage_nonlocal_init() - end subroutine homogenization_init @@ -259,25 +248,25 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE NiterationMPstate, & ip, & !< integration point number el, & !< element number - myNgrains, co, ce, ho, me, ph + myNgrains, co, ce, ho, en, ph logical :: & converged logical, dimension(2) :: & doneAndHappy !$OMP PARALLEL - !$OMP DO PRIVATE(ce,me,ho,myNgrains,NiterationMPstate,converged,doneAndHappy) + !$OMP DO PRIVATE(ce,en,ho,myNgrains,NiterationMPstate,converged,doneAndHappy) do el = FEsolving_execElem(1),FEsolving_execElem(2) ho = material_homogenizationAt(el) myNgrains = homogenization_Nconstituents(ho) do ip = FEsolving_execIP(1),FEsolving_execIP(2) ce = (el-1)*discretization_nIPs + ip - me = material_homogenizationMemberAt2(ce) + en = material_homogenizationEntry(ce) call phase_restore(ce,.false.) ! wrong name (is more a forward function) - if(homogState(ho)%sizeState > 0) homogState(ho)%state(:,me) = homogState(ho)%state0(:,me) - if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%state(:,me) = damageState_h(ho)%state0(:,me) + if(homogState(ho)%sizeState > 0) homogState(ho)%state(:,en) = homogState(ho)%state0(:,en) + if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%state(:,en) = damageState_h(ho)%state0(:,en) call damage_partition(ce) doneAndHappy = [.false.,.true.] @@ -371,14 +360,14 @@ subroutine homogenization_results case(DAMAGE_NONLOCAL_ID) group = trim(group_base)//'/damage' call results_closeGroup(results_addGroup(group)) - call damage_nonlocal_results(ho,group) + call damage_results(ho,group) end select select case(thermal_type(ho)) case(THERMAL_CONDUCTION_ID) group = trim(group_base)//'/thermal' call results_closeGroup(results_addGroup(group)) - call thermal_conduction_results(ho,group) + call thermal_results(ho,group) end select enddo @@ -458,41 +447,6 @@ subroutine homogenization_restartRead(fileHandle) end subroutine homogenization_restartRead - -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -subroutine damage_nonlocal_init - - integer :: Ninstances,Nmaterialpoints,h - class(tNode), pointer :: & - num_generic, & - material_homogenization - - print'(/,a)', ' <<<+- damage_nonlocal init -+>>>'; flush(6) - -!------------------------------------------------------------------------------------ -! read numerics parameter - num_generic => config_numerics%get('generic',defaultVal= emptyDict) - num_damage%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal) - - Ninstances = count(damage_type == DAMAGE_nonlocal_ID) - - material_homogenization => config_material%get('homogenization') - do h = 1, material_homogenization%length - if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle - - Nmaterialpoints = count(material_homogenizationAt == h) - damageState_h(h)%sizeState = 1 - allocate(damageState_h(h)%state0 (1,Nmaterialpoints), source=1.0_pReal) - allocate(damageState_h(h)%state (1,Nmaterialpoints), source=1.0_pReal) - - enddo - -end subroutine damage_nonlocal_init - - !-------------------------------------------------------------------------------------------------- !> @brief returns homogenized non local damage diffusion tensor in reference configuration !-------------------------------------------------------------------------------------------------- @@ -505,12 +459,12 @@ function damage_nonlocal_getDiffusion(ce) ho, & co - ho = material_homogenizationAt2(ce) + ho = material_homogenizationID(ce) damage_nonlocal_getDiffusion = 0.0_pReal do co = 1, homogenization_Nconstituents(ho) damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + & - crystallite_push33ToRef(co,ce,lattice_D(1:3,1:3,material_phaseAt2(co,ce))) + crystallite_push33ToRef(co,ce,lattice_D(1:3,1:3,material_phaseID(co,ce))) enddo damage_nonlocal_getDiffusion = & @@ -521,14 +475,12 @@ end function damage_nonlocal_getDiffusion !-------------------------------------------------------------------------------------------------- !> @brief parses the homogenization part from the material configuration -! ToDo: This should be done in homogenization !-------------------------------------------------------------------------------------------------- subroutine material_parseHomogenization class(tNode), pointer :: & material_homogenization, & homog, & - homogMech, & homogThermal, & homogDamage @@ -536,23 +488,11 @@ subroutine material_parseHomogenization material_homogenization => config_material%get('homogenization') - allocate(homogenization_type(size(material_name_homogenization)), source=HOMOGENIZATION_undefined_ID) allocate(thermal_type(size(material_name_homogenization)), source=THERMAL_isothermal_ID) allocate(damage_type (size(material_name_homogenization)), source=DAMAGE_none_ID) do h=1, size(material_name_homogenization) homog => material_homogenization%get(h) - homogMech => homog%get('mechanical') - select case (homogMech%get_asString('type')) - case('pass') - homogenization_type(h) = HOMOGENIZATION_NONE_ID - case('isostrain') - homogenization_type(h) = HOMOGENIZATION_ISOSTRAIN_ID - case('RGC') - homogenization_type(h) = HOMOGENIZATION_RGC_ID - case default - call IO_error(500,ext_msg=homogMech%get_asString('type')) - end select if (homog%contains('thermal')) then homogThermal => homog%get('thermal') diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90 index 495cb0beb..a067dde43 100644 --- a/src/homogenization_damage.f90 +++ b/src/homogenization_damage.f90 @@ -1,10 +1,17 @@ !-------------------------------------------------------------------------------------------------- !> @author Martin Diehl, KU Leuven !-------------------------------------------------------------------------------------------------- -submodule(homogenization) homogenization_damage +submodule(homogenization) damage use lattice + interface + + module subroutine pass_init + end subroutine pass_init + + end interface + type :: tDataContainer real(pReal), dimension(:), allocatable :: phi end type tDataContainer @@ -30,19 +37,22 @@ module subroutine damage_init() class(tNode), pointer :: & configHomogenizations, & configHomogenization, & - configHomogenizationDamage + configHomogenizationDamage, & + num_generic, & + material_homogenization integer :: ho + integer :: Ninstances,Nmaterialpoints,h print'(/,a)', ' <<<+- homogenization:damage init -+>>>' - print'(/,a)', ' <<<+- homogenization:damage:isodamage init -+>>>' + print'(/,a)', ' <<<+- homogenization:damage:pass init -+>>>' configHomogenizations => config_material%get('homogenization') allocate(param(configHomogenizations%length)) allocate(current(configHomogenizations%length)) do ho = 1, configHomogenizations%length - allocate(current(ho)%phi(count(material_homogenizationAt2==ho)), source=1.0_pReal) + allocate(current(ho)%phi(count(material_homogenizationID==ho)), source=1.0_pReal) configHomogenization => configHomogenizations%get(ho) associate(prm => param(ho)) if (configHomogenization%contains('damage')) then @@ -58,6 +68,24 @@ module subroutine damage_init() end associate enddo +!------------------------------------------------------------------------------------ +! read numerics parameter + num_generic => config_numerics%get('generic',defaultVal= emptyDict) + num_damage%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal) + + Ninstances = count(damage_type == DAMAGE_nonlocal_ID) + + material_homogenization => config_material%get('homogenization') + do h = 1, material_homogenization%length + if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle + + Nmaterialpoints = count(material_homogenizationAt == h) + damageState_h(h)%sizeState = 1 + allocate(damageState_h(h)%state0 (1,Nmaterialpoints), source=1.0_pReal) + allocate(damageState_h(h)%state (1,Nmaterialpoints), source=1.0_pReal) + + enddo + end subroutine damage_init @@ -69,14 +97,10 @@ module subroutine damage_partition(ce) real(pReal) :: phi integer, intent(in) :: ce - integer :: co - - if(damageState_h(material_homogenizationAt2(ce))%sizeState < 1) return - phi = damagestate_h(material_homogenizationAt2(ce))%state(1,material_homogenizationMemberAt2(ce)) - do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) - call phase_damage_set_phi(phi,co,ce) - enddo + if(damageState_h(material_homogenizationID(ce))%sizeState < 1) return + phi = damagestate_h(material_homogenizationID(ce))%state(1,material_homogenizationEntry(ce)) + call phase_damage_set_phi(phi,1,ce) end subroutine damage_partition @@ -88,17 +112,9 @@ end subroutine damage_partition module function damage_nonlocal_getMobility(ce) result(M) integer, intent(in) :: ce - integer :: & - co real(pReal) :: M - M = 0.0_pReal - - do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) - M = M + lattice_M(material_phaseAt2(co,ce)) - enddo - - M = M/real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) + M = lattice_M(material_phaseID(1,ce)) end function damage_nonlocal_getMobility @@ -106,20 +122,15 @@ end function damage_nonlocal_getMobility !-------------------------------------------------------------------------------------------------- !> @brief calculates homogenized damage driving forces !-------------------------------------------------------------------------------------------------- -module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ce) +module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, phi, ce) integer, intent(in) :: ce - real(pReal), intent(in) :: & + real(pReal), intent(in) :: & phi - real(pReal) :: & - phiDot, dPhiDot_dPhi + real(pReal), intent(out) :: & + phiDot - phiDot = 0.0_pReal - dPhiDot_dPhi = 0.0_pReal - - call phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce) - phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) - dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) + phiDot = phase_damage_phi_dot(phi, 1, ce) end subroutine damage_nonlocal_getSourceAndItsTangent @@ -127,26 +138,27 @@ end subroutine damage_nonlocal_getSourceAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief updated nonlocal damage field with solution from damage phase field PDE !-------------------------------------------------------------------------------------------------- -module subroutine damage_nonlocal_putNonLocalDamage(phi,ce) +module subroutine homogenization_set_phi(phi,ce) integer, intent(in) :: ce real(pReal), intent(in) :: & phi integer :: & ho, & - me + en - ho = material_homogenizationAt2(ce) - me = material_homogenizationMemberAt2(ce) - damagestate_h(ho)%state(1,me) = phi + ho = material_homogenizationID(ce) + en = material_homogenizationEntry(ce) + damagestate_h(ho)%state(1,en) = phi + current(ho)%phi(en) = phi -end subroutine damage_nonlocal_putNonLocalDamage +end subroutine homogenization_set_phi !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -module subroutine damage_nonlocal_results(ho,group) +module subroutine damage_results(ho,group) integer, intent(in) :: ho character(len=*), intent(in) :: group @@ -163,6 +175,6 @@ module subroutine damage_nonlocal_results(ho,group) enddo outputsLoop end associate -end subroutine damage_nonlocal_results +end subroutine damage_results -end submodule homogenization_damage +end submodule damage diff --git a/src/homogenization_damage_pass.f90 b/src/homogenization_damage_pass.f90 new file mode 100644 index 000000000..cbb7ec79f --- /dev/null +++ b/src/homogenization_damage_pass.f90 @@ -0,0 +1,14 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, KU Leuven +!> @brief Dummy homogenization scheme for 1 constituent per material point +!-------------------------------------------------------------------------------------------------- +submodule(homogenization:damage) damage_pass + +contains + +module subroutine pass_init + + +end subroutine pass_init + +end submodule damage_pass diff --git a/src/homogenization_mechanical.f90 b/src/homogenization_mechanical.f90 index f7074bc40..7d1c64445 100644 --- a/src/homogenization_mechanical.f90 +++ b/src/homogenization_mechanical.f90 @@ -7,51 +7,51 @@ submodule(homogenization) mechanical interface - module subroutine mechanical_pass_init - end subroutine mechanical_pass_init + module subroutine pass_init + end subroutine pass_init - module subroutine mechanical_isostrain_init - end subroutine mechanical_isostrain_init + module subroutine isostrain_init + end subroutine isostrain_init - module subroutine mechanical_RGC_init(num_homogMech) + module subroutine RGC_init(num_homogMech) class(tNode), pointer, intent(in) :: & num_homogMech !< pointer to mechanical homogenization numerics data - end subroutine mechanical_RGC_init + end subroutine RGC_init - module subroutine mechanical_isostrain_partitionDeformation(F,avgF) + module subroutine isostrain_partitionDeformation(F,avgF) real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point - end subroutine mechanical_isostrain_partitionDeformation + end subroutine isostrain_partitionDeformation - module subroutine mechanical_RGC_partitionDeformation(F,avgF,ce) + module subroutine RGC_partitionDeformation(F,avgF,ce) real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point integer, intent(in) :: & ce - end subroutine mechanical_RGC_partitionDeformation + end subroutine RGC_partitionDeformation - module subroutine mechanical_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho) + module subroutine isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho) real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses integer, intent(in) :: ho - end subroutine mechanical_isostrain_averageStressAndItsTangent + end subroutine isostrain_averageStressAndItsTangent - module subroutine mechanical_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho) + module subroutine RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho) real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses integer, intent(in) :: ho - end subroutine mechanical_RGC_averageStressAndItsTangent + end subroutine RGC_averageStressAndItsTangent - module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) + module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) logical, dimension(2) :: doneAndHappy real(pReal), dimension(:,:,:), intent(in) :: & P,& !< partitioned stresses @@ -61,16 +61,19 @@ submodule(homogenization) mechanical real(pReal), intent(in) :: dt !< time increment integer, intent(in) :: & ce !< cell - end function mechanical_RGC_updateState + end function RGC_updateState - module subroutine mechanical_RGC_results(ho,group) + module subroutine RGC_results(ho,group) integer, intent(in) :: ho !< homogenization type character(len=*), intent(in) :: group !< group name in HDF5 file - end subroutine mechanical_RGC_results + end subroutine RGC_results end interface + integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:), allocatable :: & + homogenization_type !< type of each homogenization + contains !-------------------------------------------------------------------------------------------------- @@ -86,15 +89,17 @@ module subroutine mechanical_init(num_homog) print'(/,a)', ' <<<+- homogenization:mechanical init -+>>>' + call material_parseHomogenization2() + allocate(homogenization_dPdF(3,3,3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal) homogenization_F0 = spread(math_I3,3,discretization_nIPs*discretization_Nelems) ! initialize to identity homogenization_F = homogenization_F0 ! initialize to identity allocate(homogenization_P(3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal) num_homogMech => num_homog%get('mech',defaultVal=emptyDict) - if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call mechanical_pass_init - if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call mechanical_isostrain_init - if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call mechanical_RGC_init(num_homogMech) + if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call pass_init + if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call isostrain_init + if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call RGC_init(num_homogMech) end subroutine mechanical_init @@ -110,24 +115,24 @@ module subroutine mechanical_partition(subF,ce) ce integer :: co - real(pReal), dimension (3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) :: Fs + real(pReal), dimension (3,3,homogenization_Nconstituents(material_homogenizationID(ce))) :: Fs - chosenHomogenization: select case(homogenization_type(material_homogenizationAt2(ce))) + chosenHomogenization: select case(homogenization_type(material_homogenizationID(ce))) case (HOMOGENIZATION_NONE_ID) chosenHomogenization Fs(1:3,1:3,1) = subF case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization - call mechanical_isostrain_partitionDeformation(Fs,subF) + call isostrain_partitionDeformation(Fs,subF) case (HOMOGENIZATION_RGC_ID) chosenHomogenization - call mechanical_RGC_partitionDeformation(Fs,subF,ce) + call RGC_partitionDeformation(Fs,subF,ce) end select chosenHomogenization - do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce)) - call phase_mechanical_setF(Fs(1:3,1:3,co),co,ce) + do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) + call phase_set_F(Fs(1:3,1:3,co),co,ce) enddo @@ -143,37 +148,37 @@ module subroutine mechanical_homogenize(dt,ce) integer, intent(in) :: ce integer :: co - real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) - real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) + real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationID(ce))) + real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationID(ce))) - chosenHomogenization: select case(homogenization_type(material_homogenizationAt2(ce))) + chosenHomogenization: select case(homogenization_type(material_homogenizationID(ce))) case (HOMOGENIZATION_NONE_ID) chosenHomogenization - homogenization_P(1:3,1:3,ce) = phase_mechanical_getP(1,ce) + homogenization_P(1:3,1:3,ce) = phase_P(1,ce) homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(dt,1,ce) case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization - do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) + do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ce) - Ps(:,:,co) = phase_mechanical_getP(co,ce) + Ps(:,:,co) = phase_P(co,ce) enddo - call mechanical_isostrain_averageStressAndItsTangent(& + call isostrain_averageStressAndItsTangent(& homogenization_P(1:3,1:3,ce), & homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& Ps,dPdFs, & - material_homogenizationAt2(ce)) + material_homogenizationID(ce)) case (HOMOGENIZATION_RGC_ID) chosenHomogenization - do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) + do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ce) - Ps(:,:,co) = phase_mechanical_getP(co,ce) + Ps(:,:,co) = phase_P(co,ce) enddo - call mechanical_RGC_averageStressAndItsTangent(& + call RGC_averageStressAndItsTangent(& homogenization_P(1:3,1:3,ce), & homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& Ps,dPdFs, & - material_homogenizationAt2(ce)) + material_homogenizationID(ce)) end select chosenHomogenization @@ -195,18 +200,18 @@ module function mechanical_updateState(subdt,subF,ce) result(doneAndHappy) logical, dimension(2) :: doneAndHappy integer :: co - real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) - real(pReal) :: Fs(3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) - real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) + real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationID(ce))) + real(pReal) :: Fs(3,3,homogenization_Nconstituents(material_homogenizationID(ce))) + real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationID(ce))) - if (homogenization_type(material_homogenizationAt2(ce)) == HOMOGENIZATION_RGC_ID) then - do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) + if (homogenization_type(material_homogenizationID(ce)) == HOMOGENIZATION_RGC_ID) then + do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(subdt,co,ce) - Fs(:,:,co) = phase_mechanical_getF(co,ce) - Ps(:,:,co) = phase_mechanical_getP(co,ce) + Fs(:,:,co) = phase_F(co,ce) + Ps(:,:,co) = phase_P(co,ce) enddo - doneAndHappy = mechanical_RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ce) + doneAndHappy = RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ce) else doneAndHappy = .true. endif @@ -230,7 +235,7 @@ module subroutine mechanical_results(group_base,ho) select case(homogenization_type(ho)) case(HOMOGENIZATION_rgc_ID) - call mechanical_RGC_results(ho,group) + call RGC_results(ho,group) end select @@ -244,4 +249,38 @@ module subroutine mechanical_results(group_base,ho) end subroutine mechanical_results +!-------------------------------------------------------------------------------------------------- +!> @brief parses the homogenization part from the material configuration +!-------------------------------------------------------------------------------------------------- +subroutine material_parseHomogenization2() + + class(tNode), pointer :: & + material_homogenization, & + homog, & + homogMech + + integer :: h + + material_homogenization => config_material%get('homogenization') + + allocate(homogenization_type(size(material_name_homogenization)), source=HOMOGENIZATION_undefined_ID) + + do h=1, size(material_name_homogenization) + homog => material_homogenization%get(h) + homogMech => homog%get('mechanical') + select case (homogMech%get_asString('type')) + case('pass') + homogenization_type(h) = HOMOGENIZATION_NONE_ID + case('isostrain') + homogenization_type(h) = HOMOGENIZATION_ISOSTRAIN_ID + case('RGC') + homogenization_type(h) = HOMOGENIZATION_RGC_ID + case default + call IO_error(500,ext_msg=homogMech%get_asString('type')) + end select + enddo + +end subroutine material_parseHomogenization2 + + end submodule mechanical diff --git a/src/homogenization_mechanical_RGC.f90 b/src/homogenization_mechanical_RGC.f90 index 0de7a470b..772d1c4f5 100644 --- a/src/homogenization_mechanical_RGC.f90 +++ b/src/homogenization_mechanical_RGC.f90 @@ -71,7 +71,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all necessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_RGC_init(num_homogMech) +module subroutine RGC_init(num_homogMech) class(tNode), pointer, intent(in) :: & num_homogMech !< pointer to mechanical homogenization numerics data @@ -152,7 +152,7 @@ module subroutine mechanical_RGC_init(num_homogMech) prm%N_constituents = homogMech%get_as1dInt('cluster_size',requiredSize=3) if (homogenization_Nconstituents(ho) /= product(prm%N_constituents)) & - call IO_error(211,ext_msg='N_constituents (mechanical_RGC)') + call IO_error(211,ext_msg='N_constituents (RGC)') prm%xi_alpha = homogMech%get_asFloat('xi_alpha') prm%c_alpha = homogMech%get_asFloat('c_alpha') @@ -187,13 +187,13 @@ module subroutine mechanical_RGC_init(num_homogMech) enddo -end subroutine mechanical_RGC_init +end subroutine RGC_init !-------------------------------------------------------------------------------------------------- !> @brief partitions the deformation gradient onto the constituents !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_RGC_partitionDeformation(F,avgF,ce) +module subroutine RGC_partitionDeformation(F,avgF,ce) real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned F per grain @@ -204,12 +204,12 @@ module subroutine mechanical_RGC_partitionDeformation(F,avgF,ce) real(pReal), dimension(3) :: aVect,nVect integer, dimension(4) :: intFace integer, dimension(3) :: iGrain3 - integer :: iGrain,iFace,i,j,ho,me + integer :: iGrain,iFace,i,j,ho,en - associate(prm => param(material_homogenizationAt2(ce))) + associate(prm => param(material_homogenizationID(ce))) - ho = material_homogenizationAt2(ce) - me = material_homogenizationMemberAt2(ce) + ho = material_homogenizationID(ce) + en = material_homogenizationEntry(ce) !-------------------------------------------------------------------------------------------------- ! compute the deformation gradient of individual grains due to relaxations F = 0.0_pReal @@ -217,8 +217,8 @@ module subroutine mechanical_RGC_partitionDeformation(F,avgF,ce) iGrain3 = grain1to3(iGrain,prm%N_constituents) do iFace = 1,6 intFace = getInterface(iFace,iGrain3) ! identifying 6 interfaces of each grain - aVect = relaxationVector(intFace,ho,me) ! get the relaxation vectors for each interface from global relaxation vector array - nVect = interfaceNormal(intFace,ho,me) + aVect = relaxationVector(intFace,ho,en) ! get the relaxation vectors for each interface from global relaxation vector array + nVect = interfaceNormal(intFace,ho,en) forall (i=1:3,j=1:3) & F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! calculating deformation relaxations due to interface relaxation enddo @@ -227,14 +227,14 @@ module subroutine mechanical_RGC_partitionDeformation(F,avgF,ce) end associate -end subroutine mechanical_RGC_partitionDeformation +end subroutine RGC_partitionDeformation !-------------------------------------------------------------------------------------------------- !> @brief update the internal state of the homogenization scheme and tell whether "done" and ! "happy" with result !-------------------------------------------------------------------------------------------------- -module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) +module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) logical, dimension(2) :: doneAndHappy real(pReal), dimension(:,:,:), intent(in) :: & P,& !< partitioned stresses @@ -247,7 +247,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa integer, dimension(4) :: intFaceN,intFaceP,faceID integer, dimension(3) :: nGDim,iGr3N,iGr3P - integer :: ho,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ipert,iGrain,nGrain, me + integer :: ho,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ipert,nGrain, en real(pReal), dimension(3,3,size(P,3)) :: R,pF,pR,D,pD real(pReal), dimension(3,size(P,3)) :: NN,devNull real(pReal), dimension(3) :: normP,normN,mornP,mornN @@ -261,9 +261,9 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa return endif zeroTimeStep - ho = material_homogenizationAt2(ce) + ho = material_homogenizationID(ce) + en = material_homogenizationEntry(ce) - me = material_homogenizationMemberAt2(ce) associate(stt => state(ho), st0 => state0(ho), dst => dependentState(ho), prm => param(ho)) !-------------------------------------------------------------------------------------------------- @@ -278,16 +278,16 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa ! allocate the size of the global relaxation arrays/jacobian matrices depending on the size of the cluster allocate(resid(3*nIntFaceTot), source=0.0_pReal) allocate(tract(nIntFaceTot,3), source=0.0_pReal) - relax = stt%relaxationVector(:,me) - drelax = stt%relaxationVector(:,me) - st0%relaxationVector(:,me) + relax = stt%relaxationVector(:,en) + drelax = stt%relaxationVector(:,en) - st0%relaxationVector(:,en) !-------------------------------------------------------------------------------------------------- ! computing interface mismatch and stress penalty tensor for all interfaces of all grains - call stressPenalty(R,NN,avgF,F,ho,me) + call stressPenalty(R,NN,avgF,F,ho,en) !-------------------------------------------------------------------------------------------------- ! calculating volume discrepancy and stress penalty related to overall volume discrepancy - call volumePenalty(D,dst%volumeDiscrepancy(me),avgF,F,nGrain) + call volumePenalty(D,dst%volumeDiscrepancy(en),avgF,F,nGrain) !------------------------------------------------------------------------------------------------ ! computing the residual stress from the balance of traction at all (interior) interfaces @@ -299,7 +299,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate system (3-dimensional index) iGrN = grain3to1(iGr3N,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) intFaceN = getInterface(2*faceID(1),iGr3N) - normN = interfaceNormal(intFaceN,ho,me) + normN = interfaceNormal(intFaceN,ho,en) !-------------------------------------------------------------------------------------------------- ! identify the right/up/front grain (+|P) @@ -307,7 +307,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate system (3-dimensional index) iGrP = grain3to1(iGr3P,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) intFaceP = getInterface(2*faceID(1)-1,iGr3P) - normP = interfaceNormal(intFaceP,ho,me) + normP = interfaceNormal(intFaceP,ho,en) !-------------------------------------------------------------------------------------------------- ! compute the residual of traction at the interface (in local system, 4-dimensional index) @@ -335,9 +335,9 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa if (residMax < num%rtol*stresMax .or. residMax < num%atol) then doneAndHappy = .true. - dst%mismatch(1:3,me) = sum(NN,2)/real(nGrain,pReal) - dst%relaxationRate_avg(me) = sum(abs(drelax))/dt/real(3*nIntFaceTot,pReal) - dst%relaxationRate_max(me) = maxval(abs(drelax))/dt + dst%mismatch(1:3,en) = sum(NN,2)/real(nGrain,pReal) + dst%relaxationRate_avg(en) = sum(abs(drelax))/dt/real(3*nIntFaceTot,pReal) + dst%relaxationRate_max(en) = maxval(abs(drelax))/dt return @@ -363,10 +363,10 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate sytem iGrN = grain3to1(iGr3N,param(ho)%N_constituents) ! translate into global grain ID intFaceN = getInterface(2*faceID(1),iGr3N) ! identifying the connecting interface in local coordinate system - normN = interfaceNormal(intFaceN,ho,me) + normN = interfaceNormal(intFaceN,ho,en) do iFace = 1,6 intFaceN = getInterface(iFace,iGr3N) ! identifying all interfaces that influence relaxation of the above interface - mornN = interfaceNormal(intFaceN,ho,me) + mornN = interfaceNormal(intFaceN,ho,en) iMun = interface4to1(intFaceN,param(ho)%N_constituents) ! translate the interfaces ID into local 4-dimensional index if (iMun > 0) then ! get the corresponding tangent do i=1,3; do j=1,3; do k=1,3; do l=1,3 @@ -384,10 +384,10 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate sytem iGrP = grain3to1(iGr3P,param(ho)%N_constituents) ! translate into global grain ID intFaceP = getInterface(2*faceID(1)-1,iGr3P) ! identifying the connecting interface in local coordinate system - normP = interfaceNormal(intFaceP,ho,me) + normP = interfaceNormal(intFaceP,ho,en) do iFace = 1,6 intFaceP = getInterface(iFace,iGr3P) ! identifying all interfaces that influence relaxation of the above interface - mornP = interfaceNormal(intFaceP,ho,me) + mornP = interfaceNormal(intFaceP,ho,en) iMun = interface4to1(intFaceP,param(ho)%N_constituents) ! translate the interfaces ID into local 4-dimensional index if (iMun > 0) then ! get the corresponding tangent do i=1,3; do j=1,3; do k=1,3; do l=1,3 @@ -408,9 +408,9 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa do ipert = 1,3*nIntFaceTot p_relax = relax p_relax(ipert) = relax(ipert) + num%pPert ! perturb the relaxation vector - stt%relaxationVector(:,me) = p_relax - call grainDeformation(pF,avgF,ho,me) ! rain deformation from perturbed state - call stressPenalty(pR,DevNull, avgF,pF,ho,me) ! stress penalty due to interface mismatch from perturbed state + stt%relaxationVector(:,en) = p_relax + call grainDeformation(pF,avgF,ho,en) ! rain deformation from perturbed state + call stressPenalty(pR,DevNull, avgF,pF,ho,en) ! stress penalty due to interface mismatch from perturbed state call volumePenalty(pD,devNull(1,1), avgF,pF,nGrain) ! stress penalty due to volume discrepancy from perturbed state !-------------------------------------------------------------------------------------------------- @@ -424,7 +424,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa iGr3N = faceID(2:4) ! identify the grain ID in local coordinate system (3-dimensional index) iGrN = grain3to1(iGr3N,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) intFaceN = getInterface(2*faceID(1),iGr3N) ! identify the interface ID of the grain - normN = interfaceNormal(intFaceN,ho,me) + normN = interfaceNormal(intFaceN,ho,en) !-------------------------------------------------------------------------------------------------- ! identify the right/up/front grain (+|P) @@ -432,7 +432,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identify the grain ID in local coordinate system (3-dimensional index) iGrP = grain3to1(iGr3P,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) intFaceP = getInterface(2*faceID(1)-1,iGr3P) ! identify the interface ID of the grain - normP = interfaceNormal(intFaceP,ho,me) + normP = interfaceNormal(intFaceP,ho,en) !-------------------------------------------------------------------------------------------------- ! compute the residual stress (contribution of mismatch and volume penalties) from perturbed state @@ -472,7 +472,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa do i = 1,3*nIntFaceTot;do j = 1,3*nIntFaceTot drelax(i) = drelax(i) - jnverse(i,j)*resid(j) ! Calculate the correction for the state variable enddo; enddo - stt%relaxationVector(:,me) = relax + drelax ! Updateing the state variable for the next iteration + stt%relaxationVector(:,en) = relax + drelax ! Updateing the state variable for the next iteration if (any(abs(drelax) > num%maxdRelax)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large doneAndHappy = [.true.,.false.] !$OMP CRITICAL (write2out) @@ -488,14 +488,14 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa !------------------------------------------------------------------------------------------------ !> @brief calculate stress-like penalty due to deformation mismatch !------------------------------------------------------------------------------------------------ - subroutine stressPenalty(rPen,nMis,avgF,fDef,ho,me) + subroutine stressPenalty(rPen,nMis,avgF,fDef,ho,en) real(pReal), dimension (:,:,:), intent(out) :: rPen !< stress-like penalty real(pReal), dimension (:,:), intent(out) :: nMis !< total amount of mismatch real(pReal), dimension (:,:,:), intent(in) :: fDef !< deformation gradients real(pReal), dimension (3,3), intent(in) :: avgF !< initial effective stretch tensor - integer, intent(in) :: ho, me + integer, intent(in) :: ho, en integer, dimension (4) :: intFace integer, dimension (3) :: iGrain3,iGNghb3,nGDim @@ -515,7 +515,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa ! get the correction factor the modulus of penalty stress representing the evolution of area of ! the interfaces due to deformations - surfCorr = surfaceCorrection(avgF,ho,me) + surfCorr = surfaceCorrection(avgF,ho,en) associate(prm => param(ho)) @@ -527,7 +527,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa interfaceLoop: do iFace = 1,6 intFace = getInterface(iFace,iGrain3) ! get the 4-dimensional index of the interface in local numbering system of the grain - nVect = interfaceNormal(intFace,ho,me) + nVect = interfaceNormal(intFace,ho,en) iGNghb3 = iGrain3 ! identify the neighboring grain across the interface iGNghb3(abs(intFace(1))) = iGNghb3(abs(intFace(1))) & + int(real(intFace(1),pReal)/real(abs(intFace(1)),pReal)) @@ -611,14 +611,14 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa !> @brief compute the correction factor accouted for surface evolution (area change) due to ! deformation !-------------------------------------------------------------------------------------------------- - function surfaceCorrection(avgF,ho,me) + function surfaceCorrection(avgF,ho,en) real(pReal), dimension(3) :: surfaceCorrection real(pReal), dimension(3,3), intent(in) :: avgF !< average F integer, intent(in) :: & ho, & - me + en real(pReal), dimension(3,3) :: invC real(pReal), dimension(3) :: nVect real(pReal) :: detF @@ -629,7 +629,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa surfaceCorrection = 0.0_pReal do iBase = 1,3 - nVect = interfaceNormal([iBase,1,1,1],ho,me) + nVect = interfaceNormal([iBase,1,1,1],ho,en) do i = 1,3; do j = 1,3 surfaceCorrection(iBase) = surfaceCorrection(iBase) + invC(i,j)*nVect(i)*nVect(j) ! compute the component of (the inverse of) the stretch in the direction of the normal enddo; enddo @@ -651,7 +651,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa real(pReal), dimension(6,6) :: C - C = phase_homogenizedC(material_phaseAt2(grainID,ce),material_phaseMemberAt2(grainID,ce)) + C = phase_homogenizedC(material_phaseID(grainID,ce),material_phaseEntry(grainID,ce)) equivalentMu = lattice_equivalent_mu(C,'voigt') end function equivalentMu @@ -661,14 +661,14 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa !> @brief calculating the grain deformation gradient (the same with ! homogenization_RGC_partitionDeformation, but used only for perturbation scheme) !------------------------------------------------------------------------------------------------- - subroutine grainDeformation(F, avgF, ho, me) + subroutine grainDeformation(F, avgF, ho, en) real(pReal), dimension(:,:,:), intent(out) :: F !< partitioned F per grain real(pReal), dimension(:,:), intent(in) :: avgF !< averaged F integer, intent(in) :: & ho, & - me + en real(pReal), dimension(3) :: aVect,nVect integer, dimension(4) :: intFace @@ -685,8 +685,8 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa iGrain3 = grain1to3(iGrain,prm%N_constituents) do iFace = 1,6 intFace = getInterface(iFace,iGrain3) - aVect = relaxationVector(intFace,ho,me) - nVect = interfaceNormal(intFace,ho,me) + aVect = relaxationVector(intFace,ho,en) + nVect = interfaceNormal(intFace,ho,en) forall (i=1:3,j=1:3) & F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations enddo @@ -697,13 +697,13 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa end subroutine grainDeformation -end function mechanical_RGC_updateState +end function RGC_updateState !-------------------------------------------------------------------------------------------------- !> @brief derive average stress and stiffness from constituent quantities !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho) +module subroutine RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho) real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point @@ -715,13 +715,13 @@ module subroutine mechanical_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dP avgP = sum(P,3) /real(product(param(ho)%N_constituents),pReal) dAvgPdAvgF = sum(dPdF,5)/real(product(param(ho)%N_constituents),pReal) -end subroutine mechanical_RGC_averageStressAndItsTangent +end subroutine RGC_averageStressAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_RGC_results(ho,group) +module subroutine RGC_results(ho,group) integer, intent(in) :: ho character(len=*), intent(in) :: group @@ -747,17 +747,17 @@ module subroutine mechanical_RGC_results(ho,group) enddo outputsLoop end associate -end subroutine mechanical_RGC_results +end subroutine RGC_results !-------------------------------------------------------------------------------------------------- !> @brief collect relaxation vectors of an interface !-------------------------------------------------------------------------------------------------- -pure function relaxationVector(intFace,ho,me) +pure function relaxationVector(intFace,ho,en) real(pReal), dimension (3) :: relaxationVector - integer, intent(in) :: ho,me + integer, intent(in) :: ho,en integer, dimension(4), intent(in) :: intFace !< set of interface ID in 4D array (normal and position) integer :: iNum @@ -770,7 +770,7 @@ pure function relaxationVector(intFace,ho,me) iNum = interface4to1(intFace,prm%N_constituents) ! identify the position of the interface in global state array if (iNum > 0) then - relaxationVector = stt%relaxationVector((3*iNum-2):(3*iNum),me) + relaxationVector = stt%relaxationVector((3*iNum-2):(3*iNum),en) else relaxationVector = 0.0_pReal endif @@ -783,14 +783,14 @@ end function relaxationVector !-------------------------------------------------------------------------------------------------- !> @brief identify the normal of an interface !-------------------------------------------------------------------------------------------------- -pure function interfaceNormal(intFace,ho,me) +pure function interfaceNormal(intFace,ho,en) real(pReal), dimension(3) :: interfaceNormal integer, dimension(4), intent(in) :: intFace !< interface ID in 4D array (normal and position) integer, intent(in) :: & ho, & - me + en integer :: nPos associate (dst => dependentState(ho)) @@ -801,7 +801,7 @@ pure function interfaceNormal(intFace,ho,me) nPos = abs(intFace(1)) ! identify the position of the interface in global state array interfaceNormal(nPos) = real(intFace(1)/abs(intFace(1)),pReal) ! get the normal vector w.r.t. cluster axis - interfaceNormal = matmul(dst%orientation(1:3,1:3,me),interfaceNormal) ! map the normal vector into sample coordinate system (basis) + interfaceNormal = matmul(dst%orientation(1:3,1:3,en),interfaceNormal) ! map the normal vector into sample coordinate system (basis) end associate diff --git a/src/homogenization_mechanical_isostrain.f90 b/src/homogenization_mechanical_isostrain.f90 index e7d4cff4c..9a3704575 100644 --- a/src/homogenization_mechanical_isostrain.f90 +++ b/src/homogenization_mechanical_isostrain.f90 @@ -26,7 +26,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_isostrain_init +module subroutine isostrain_init integer :: & h, & @@ -56,7 +56,7 @@ module subroutine mechanical_isostrain_init case ('avg') prm%mapping = average_ID case default - call IO_error(211,ext_msg='sum'//' (mechanical_isostrain)') + call IO_error(211,ext_msg='sum'//' (isostrain)') end select Nmaterialpoints = count(material_homogenizationAt == h) @@ -68,13 +68,13 @@ module subroutine mechanical_isostrain_init enddo -end subroutine mechanical_isostrain_init +end subroutine isostrain_init !-------------------------------------------------------------------------------------------------- !> @brief partitions the deformation gradient onto the constituents !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_isostrain_partitionDeformation(F,avgF) +module subroutine isostrain_partitionDeformation(F,avgF) real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient @@ -82,13 +82,13 @@ module subroutine mechanical_isostrain_partitionDeformation(F,avgF) F = spread(avgF,3,size(F,3)) -end subroutine mechanical_isostrain_partitionDeformation +end subroutine isostrain_partitionDeformation !-------------------------------------------------------------------------------------------------- !> @brief derive average stress and stiffness from constituent quantities !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho) +module subroutine isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho) real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point @@ -110,6 +110,6 @@ module subroutine mechanical_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvg end associate -end subroutine mechanical_isostrain_averageStressAndItsTangent +end subroutine isostrain_averageStressAndItsTangent end submodule isostrain diff --git a/src/homogenization_mechanical_pass.f90 b/src/homogenization_mechanical_pass.f90 index 6217e6836..23fe74f44 100644 --- a/src/homogenization_mechanical_pass.f90 +++ b/src/homogenization_mechanical_pass.f90 @@ -4,14 +4,14 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief dummy homogenization homogenization scheme for 1 constituent per material point !-------------------------------------------------------------------------------------------------- -submodule(homogenization:mechanical) none +submodule(homogenization:mechanical) mechanical_pass contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all necessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_pass_init +module subroutine pass_init integer :: & Ninstances, & @@ -27,7 +27,7 @@ module subroutine mechanical_pass_init if(homogenization_type(h) /= HOMOGENIZATION_NONE_ID) cycle if(homogenization_Nconstituents(h) /= 1) & - call IO_error(211,ext_msg='N_constituents (mechanical_pass)') + call IO_error(211,ext_msg='N_constituents (pass)') Nmaterialpoints = count(material_homogenizationAt == h) homogState(h)%sizeState = 0 @@ -36,6 +36,6 @@ module subroutine mechanical_pass_init enddo -end subroutine mechanical_pass_init +end subroutine pass_init -end submodule none +end submodule mechanical_pass diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 index 7602e50d2..8dcddda7a 100644 --- a/src/homogenization_thermal.f90 +++ b/src/homogenization_thermal.f90 @@ -5,6 +5,16 @@ submodule(homogenization) thermal use lattice + interface + + module subroutine pass_init + end subroutine pass_init + + module subroutine isotemperature_init + end subroutine isotemperature_init + + end interface + type :: tDataContainer real(pReal), dimension(:), allocatable :: T, dot_T end type tDataContainer @@ -35,7 +45,7 @@ module subroutine thermal_init() print'(/,a)', ' <<<+- homogenization:thermal init -+>>>' - print'(/,a)', ' <<<+- homogenization:thermal:isotemperature init -+>>>' + print'(/,a)', ' <<<+- homogenization:thermal:pass init -+>>>' @@ -44,8 +54,8 @@ module subroutine thermal_init() allocate(current(configHomogenizations%length)) do ho = 1, configHomogenizations%length - allocate(current(ho)%T(count(material_homogenizationAt2==ho)), source=300.0_pReal) - allocate(current(ho)%dot_T(count(material_homogenizationAt2==ho)), source=0.0_pReal) + allocate(current(ho)%T(count(material_homogenizationID==ho)), source=300.0_pReal) + allocate(current(ho)%dot_T(count(material_homogenizationID==ho)), source=0.0_pReal) configHomogenization => configHomogenizations%get(ho) associate(prm => param(ho)) if (configHomogenization%contains('thermal')) then @@ -75,9 +85,9 @@ module subroutine thermal_partition(ce) integer :: co - T = current(material_homogenizationAt2(ce))%T(material_homogenizationMemberAt2(ce)) - dot_T = current(material_homogenizationAt2(ce))%dot_T(material_homogenizationMemberAt2(ce)) - do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) + T = current(material_homogenizationID(ce))%T(material_homogenizationEntry(ce)) + dot_T = current(material_homogenizationID(ce))%dot_T(material_homogenizationEntry(ce)) + do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) call phase_thermal_setField(T,dot_T,co,ce) enddo @@ -109,19 +119,29 @@ module function thermal_conduction_getConductivity(ce) result(K) K = 0.0_pReal - do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) - K = K + crystallite_push33ToRef(co,ce,lattice_K(:,:,material_phaseAt2(co,ce))) + do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) + K = K + crystallite_push33ToRef(co,ce,lattice_K(:,:,material_phaseID(co,ce))) enddo - K = K / real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) + K = K / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) end function thermal_conduction_getConductivity +module function homogenization_thermal_mu_T(ce) result(mu_T) + + integer, intent(in) :: ce + real(pReal) :: mu_T + + mu_T = c_P(ce) * rho(ce) + +end function homogenization_thermal_mu_T + + !-------------------------------------------------------------------------------------------------- !> @brief returns homogenized specific heat capacity !-------------------------------------------------------------------------------------------------- -module function thermal_conduction_getSpecificHeat(ce) result(c_P) +function c_P(ce) integer, intent(in) :: ce real(pReal) :: c_P @@ -129,21 +149,20 @@ module function thermal_conduction_getSpecificHeat(ce) result(c_P) integer :: co - c_P = 0.0_pReal - - do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) - c_P = c_P + lattice_c_p(material_phaseAt2(co,ce)) + c_P = lattice_c_p(material_phaseID(1,ce)) + do co = 2, homogenization_Nconstituents(material_homogenizationID(ce)) + c_P = c_P + lattice_c_p(material_phaseID(co,ce)) enddo - c_P = c_P / real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) + c_P = c_P / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) -end function thermal_conduction_getSpecificHeat +end function c_P !-------------------------------------------------------------------------------------------------- !> @brief returns homogenized mass density !-------------------------------------------------------------------------------------------------- -module function thermal_conduction_getMassDensity(ce) result(rho) +function rho(ce) integer, intent(in) :: ce real(pReal) :: rho @@ -151,15 +170,14 @@ module function thermal_conduction_getMassDensity(ce) result(rho) integer :: co - rho = 0.0_pReal - - do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) - rho = rho + lattice_rho(material_phaseAt2(co,ce)) + rho = lattice_rho(material_phaseID(1,ce)) + do co = 2, homogenization_Nconstituents(material_homogenizationID(ce)) + rho = rho + lattice_rho(material_phaseID(co,ce)) enddo - rho = rho / real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) + rho = rho / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) -end function thermal_conduction_getMassDensity +end function rho @@ -172,8 +190,8 @@ module subroutine homogenization_thermal_setField(T,dot_T, ce) real(pReal), intent(in) :: T, dot_T - current(material_homogenizationAt2(ce))%T(material_homogenizationMemberAt2(ce)) = T - current(material_homogenizationAt2(ce))%dot_T(material_homogenizationMemberAt2(ce)) = dot_T + current(material_homogenizationID(ce))%T(material_homogenizationEntry(ce)) = T + current(material_homogenizationID(ce))%dot_T(material_homogenizationEntry(ce)) = dot_T end subroutine homogenization_thermal_setField @@ -183,7 +201,7 @@ end subroutine homogenization_thermal_setField !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -module subroutine thermal_conduction_results(ho,group) +module subroutine thermal_results(ho,group) integer, intent(in) :: ho character(len=*), intent(in) :: group @@ -199,7 +217,7 @@ module subroutine thermal_conduction_results(ho,group) enddo outputsLoop end associate -end subroutine thermal_conduction_results +end subroutine thermal_results module function homogenization_thermal_T(ce) result(T) @@ -207,7 +225,7 @@ module function homogenization_thermal_T(ce) result(T) integer, intent(in) :: ce real(pReal) :: T - T = current(material_homogenizationAt2(ce))%T(material_homogenizationMemberAt2(ce)) + T = current(material_homogenizationID(ce))%T(material_homogenizationEntry(ce)) end function homogenization_thermal_T diff --git a/src/homogenization_thermal_isotemperature.f90 b/src/homogenization_thermal_isotemperature.f90 new file mode 100644 index 000000000..7358ecf08 --- /dev/null +++ b/src/homogenization_thermal_isotemperature.f90 @@ -0,0 +1,14 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, KU Leuven +!> @brief Dummy homogenization scheme for 1 constituent per material point +!-------------------------------------------------------------------------------------------------- +submodule(homogenization:thermal) isotemperature + +contains + +module subroutine isotemperature_init + + +end subroutine isotemperature_init + +end submodule isotemperature diff --git a/src/homogenization_thermal_pass.f90 b/src/homogenization_thermal_pass.f90 new file mode 100644 index 000000000..940a15024 --- /dev/null +++ b/src/homogenization_thermal_pass.f90 @@ -0,0 +1,14 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, KU Leuven +!> @brief Dummy homogenization scheme for 1 constituent per material point +!-------------------------------------------------------------------------------------------------- +submodule(homogenization:thermal) thermal_pass + +contains + +module subroutine pass_init + + +end subroutine pass_init + +end submodule thermal_pass diff --git a/src/material.f90 b/src/material.f90 index aecfaa1dd..6575872ed 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -28,14 +28,14 @@ module material integer, dimension(:), allocatable, public, protected :: & ! (elem) material_homogenizationAt, & !< homogenization ID of each element - material_homogenizationAt2, & !< per cell - material_homogenizationMemberAt2 !< cell + material_homogenizationID, & !< per cell + material_homogenizationEntry !< cell integer, dimension(:,:), allocatable :: & ! (ip,elem) material_homogenizationMemberAt !< position of the element within its homogenization instance integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem) material_phaseAt, & !< phase ID of each element - material_phaseAt2, & !< per constituent,cell - material_phaseMemberAt2 !< per constituent, cell + material_phaseID, & !< per constituent,cell + material_phaseEntry !< per constituent, cell integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,IP,elem) material_phaseMemberAt !< position of the element within its phase instance @@ -117,10 +117,10 @@ subroutine material_parseMaterial allocate(material_phaseMemberAt(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems),source=0) - allocate(material_homogenizationAt2(discretization_nIPs*discretization_Nelems),source=0) - allocate(material_homogenizationMemberAt2(discretization_nIPs*discretization_Nelems),source=0) - allocate(material_phaseAt2(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0) - allocate(material_phaseMemberAt2(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0) + allocate(material_homogenizationID(discretization_nIPs*discretization_Nelems),source=0) + allocate(material_homogenizationEntry(discretization_nIPs*discretization_Nelems),source=0) + allocate(material_phaseID(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0) + allocate(material_phaseEntry(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0) do el = 1, discretization_Nelems material => materials%get(discretization_materialAt(el)) @@ -131,8 +131,8 @@ subroutine material_parseMaterial ce = (el-1)*discretization_nIPs + ip counterHomogenization(material_homogenizationAt(el)) = counterHomogenization(material_homogenizationAt(el)) + 1 material_homogenizationMemberAt(ip,el) = counterHomogenization(material_homogenizationAt(el)) - material_homogenizationAt2(ce) = material_homogenizationAt(el) - material_homogenizationMemberAt2(ce) = material_homogenizationMemberAt(ip,el) + material_homogenizationID(ce) = material_homogenizationAt(el) + material_homogenizationEntry(ce) = material_homogenizationMemberAt(ip,el) enddo frac = 0.0_pReal @@ -146,8 +146,8 @@ subroutine material_parseMaterial counterPhase(material_phaseAt(co,el)) = counterPhase(material_phaseAt(co,el)) + 1 material_phaseMemberAt(co,ip,el) = counterPhase(material_phaseAt(co,el)) - material_phaseAt2(co,ce) = material_phaseAt(co,el) - material_phaseMemberAt2(co,ce) = material_phaseMemberAt(co,ip,el) + material_phaseID(co,ce) = material_phaseAt(co,el) + material_phaseEntry(co,ce) = material_phaseMemberAt(co,ip,el) enddo enddo diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index 9ec02f9fb..979845484 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/src/mesh/DAMASK_mesh.f90 @@ -85,8 +85,8 @@ program DAMASK_mesh stagItMax = num_mesh%get_asInt('maxStaggeredIter',defaultVal=10) maxCutBack = num_mesh%get_asInt('maxCutBack',defaultVal=3) - if (stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter') - if (maxCutBack < 0) call IO_error(301,ext_msg='maxCutBack') + if (stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter') + if (maxCutBack < 0) call IO_error(301,ext_msg='maxCutBack') ! reading basic information from load case file and allocate data structure containing load cases call DMGetDimension(geomMesh,dimPlex,ierr) !< dimension of mesh (2D or 3D) diff --git a/src/phase.f90 b/src/phase.f90 index 9470ab8af..1adba03d9 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -100,14 +100,12 @@ module phase integer, intent(in) :: ph end subroutine damage_results - module subroutine mechanical_windForward(ph,me) - integer, intent(in) :: ph, me - end subroutine mechanical_windForward - - module subroutine mechanical_forward() end subroutine mechanical_forward + module subroutine damage_forward() + end subroutine damage_forward + module subroutine thermal_forward() end subroutine thermal_forward @@ -147,20 +145,20 @@ module phase real(pReal), dimension(3,3) :: L_p end function mechanical_L_p - module function phase_mechanical_getF(co,ce) result(F) + module function phase_F(co,ce) result(F) integer, intent(in) :: co, ce real(pReal), dimension(3,3) :: F - end function phase_mechanical_getF + end function phase_F module function mechanical_F_e(ph,me) result(F_e) integer, intent(in) :: ph,me real(pReal), dimension(3,3) :: F_e end function mechanical_F_e - module function phase_mechanical_getP(co,ce) result(P) + module function phase_P(co,ce) result(P) integer, intent(in) :: co, ce real(pReal), dimension(3,3) :: P - end function phase_mechanical_getP + end function phase_P module function phase_damage_get_phi(co,ip,el) result(phi) integer, intent(in) :: co, ip, el @@ -183,10 +181,10 @@ module phase end function damage_phi - module subroutine phase_mechanical_setF(F,co,ce) + module subroutine phase_set_F(F,co,ce) real(pReal), dimension(3,3), intent(in) :: F integer, intent(in) :: co, ce - end subroutine phase_mechanical_setF + end subroutine phase_set_F module subroutine phase_thermal_setField(T,dot_T, co,ce) real(pReal), intent(in) :: T, dot_T @@ -229,14 +227,13 @@ module phase end function phase_homogenizedC - module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce) - integer, intent(in) :: ce + module function phase_damage_phi_dot(phi,co,ce) result(phi_dot) + integer, intent(in) :: ce,co real(pReal), intent(in) :: & phi !< damage parameter - real(pReal), intent(inout) :: & - phiDot, & - dPhiDot_dPhi - end subroutine phase_damage_getRateAndItsTangents + real(pReal) :: & + phi_dot + end function phase_damage_phi_dot module subroutine phase_thermal_getRate(TDot, ph,me) integer, intent(in) :: ph, me @@ -261,7 +258,7 @@ module phase end subroutine plastic_dependentState - module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) + module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) integer, intent(in) :: ph, me real(pReal), intent(in), dimension(3,3) :: & S @@ -269,9 +266,9 @@ module phase Ld !< damage velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) - end subroutine kinematics_cleavage_opening_LiAndItsTangent + end subroutine damage_anisobrittle_LiAndItsTangent - module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) + module subroutine damage_isoductile_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) integer, intent(in) :: ph, me real(pReal), intent(in), dimension(3,3) :: & S @@ -279,7 +276,7 @@ module phase Ld !< damage velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) - end subroutine kinematics_slipplane_opening_LiAndItsTangent + end subroutine damage_isoductile_LiAndItsTangent end interface @@ -303,7 +300,7 @@ module phase public :: & phase_init, & phase_homogenizedC, & - phase_damage_getRateAndItsTangents, & + phase_damage_phi_dot, & phase_thermal_getRate, & phase_results, & phase_allocateState, & @@ -323,21 +320,19 @@ module phase phase_thermal_setField, & phase_damage_set_phi, & phase_damage_get_phi, & - phase_mechanical_getP, & - phase_mechanical_setF, & - phase_mechanical_getF, & - phase_windForward + phase_P, & + phase_set_F, & + phase_F contains !-------------------------------------------------------------------------------------------------- -!> @brief Initialze constitutive models for individual physics +!> @brief Initialize constitutive models for individual physics !-------------------------------------------------------------------------------------------------- subroutine phase_init integer :: & - ph, & !< counter in phase loop - so !< counter in source loop + ph class (tNode), pointer :: & debug_constitutive, & materials, & @@ -382,12 +377,12 @@ end subroutine phase_init !> @brief Allocate the components of the state structure for a given phase !-------------------------------------------------------------------------------------------------- subroutine phase_allocateState(state, & - Nconstituents,sizeState,sizeDotState,sizeDeltaState) + NEntries,sizeState,sizeDotState,sizeDeltaState) class(tState), intent(out) :: & state integer, intent(in) :: & - Nconstituents, & + NEntries, & sizeState, & sizeDotState, & sizeDeltaState @@ -398,13 +393,13 @@ subroutine phase_allocateState(state, & state%sizeDeltaState = sizeDeltaState state%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition - allocate(state%atol (sizeState), source=0.0_pReal) - allocate(state%state0 (sizeState,Nconstituents), source=0.0_pReal) - allocate(state%state (sizeState,Nconstituents), source=0.0_pReal) + allocate(state%atol (sizeState), source=0.0_pReal) + allocate(state%state0 (sizeState,NEntries), source=0.0_pReal) + allocate(state%state (sizeState,NEntries), source=0.0_pReal) - allocate(state%dotState (sizeDotState,Nconstituents), source=0.0_pReal) + allocate(state%dotState (sizeDotState,NEntries), source=0.0_pReal) - allocate(state%deltaState (sizeDeltaState,Nconstituents), source=0.0_pReal) + allocate(state%deltaState (sizeDeltaState,NEntries), source=0.0_pReal) end subroutine phase_allocateState @@ -422,10 +417,10 @@ subroutine phase_restore(ce,includeL) co - do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce)) - if (damageState(material_phaseAt2(co,ce))%sizeState > 0) & - damageState(material_phaseAt2(co,ce))%state( :,material_phasememberAt2(co,ce)) = & - damageState(material_phaseAt2(co,ce))%state0(:,material_phasememberAt2(co,ce)) + do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) + if (damageState(material_phaseID(co,ce))%sizeState > 0) & + damageState(material_phaseID(co,ce))%state( :,material_phaseEntry(co,ce)) = & + damageState(material_phaseID(co,ce))%state0(:,material_phaseEntry(co,ce)) enddo call mechanical_restore(ce,includeL) @@ -435,21 +430,13 @@ end subroutine phase_restore !-------------------------------------------------------------------------------------------------- !> @brief Forward data after successful increment. -! ToDo: Any guessing for the current states possible? !-------------------------------------------------------------------------------------------------- subroutine phase_forward() - integer :: ph - - call mechanical_forward() + call damage_forward() call thermal_forward() - do ph = 1, size(damageState) - if (damageState(ph)%sizeState > 0) & - damageState(ph)%state0 = damageState(ph)%state - enddo - end subroutine phase_forward @@ -484,11 +471,9 @@ subroutine crystallite_init() integer :: & ph, & - me, & co, & !< counter in integration point component loop ip, & !< counter in integration point loop el, & !< counter in element loop - so, & cMax, & !< maximum number of integration point components iMax, & !< maximum number of integration points eMax !< maximum number of elements @@ -550,12 +535,10 @@ subroutine crystallite_init() flush(IO_STDOUT) - !$OMP PARALLEL DO PRIVATE(ph,me) + !$OMP PARALLEL DO do el = 1, size(material_phaseMemberAt,3) do ip = 1, size(material_phaseMemberAt,2) do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) - ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) call crystallite_orientations(co,ip,el) call plastic_dependentState(co,ip,el) ! update dependent state variables to be consistent with basic states enddo @@ -567,34 +550,6 @@ subroutine crystallite_init() end subroutine crystallite_init -!-------------------------------------------------------------------------------------------------- -!> @brief Wind homog inc forward. -!-------------------------------------------------------------------------------------------------- -subroutine phase_windForward(ip,el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - - integer :: & - co, & !< constituent number - so, ph, me - - - do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) - ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) - - call mechanical_windForward(ph,me) - - if(damageState(ph)%sizeState > 0) damageState(ph)%state0(:,me) = damageState(ph)%state(:,me) - - - enddo - -end subroutine phase_windForward - - !-------------------------------------------------------------------------------------------------- !> @brief calculates orientations !-------------------------------------------------------------------------------------------------- @@ -629,11 +584,11 @@ function crystallite_push33ToRef(co,ce, tensor33) real(pReal), dimension(3,3) :: crystallite_push33ToRef real(pReal), dimension(3,3) :: T - integer :: ph, me + integer :: ph, en - ph = material_phaseAt2(co,ce) - me = material_phaseMemberAt2(co,ce) - T = matmul(material_orientation0(co,ph,me)%asMatrix(),transpose(math_inv33(phase_mechanical_getF(co,ce)))) ! ToDo: initial orientation correct? + ph = material_phaseID(co,ce) + en = material_phaseEntry(co,ce) + T = matmul(material_orientation0(co,ph,en)%asMatrix(),transpose(math_inv33(phase_F(co,ce)))) ! ToDo: initial orientation correct? crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T)) @@ -658,8 +613,7 @@ end function converged !-------------------------------------------------------------------------------------------------- -!> @brief Write current restart information (Field and constitutive data) to file. -! ToDo: Merge data into one file for MPI +!> @brief Write restart data to file. !-------------------------------------------------------------------------------------------------- subroutine phase_restartWrite(fileHandle) @@ -687,8 +641,7 @@ end subroutine phase_restartWrite !-------------------------------------------------------------------------------------------------- -!> @brief Read data for restart -! ToDo: Merge data into one file for MPI +!> @brief Read restart data from file. !-------------------------------------------------------------------------------------------------- subroutine phase_restartRead(fileHandle) diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index c6b430b74..62ce70368 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -1,7 +1,7 @@ !---------------------------------------------------------------------------------------------------- !> @brief internal microstructure state for all damage sources and kinematics constitutive models !---------------------------------------------------------------------------------------------------- -submodule(phase) damagee +submodule(phase) damage enum, bind(c); enumerator :: & DAMAGE_UNDEFINED_ID, & DAMAGE_ISOBRITTLE_ID, & @@ -65,43 +65,6 @@ submodule(phase) damagee integer, intent(in) :: ph,me end subroutine isoductile_dotState - - module subroutine anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) - integer, intent(in) :: ph,me - real(pReal), intent(in) :: & - phi !< damage parameter - real(pReal), intent(out) :: & - localphiDot, & - dLocalphiDot_dPhi - end subroutine anisobrittle_getRateAndItsTangent - - module subroutine anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me) - integer, intent(in) :: ph,me - real(pReal), intent(in) :: & - phi !< damage parameter - real(pReal), intent(out) :: & - localphiDot, & - dLocalphiDot_dPhi - end subroutine anisoductile_getRateAndItsTangent - - module subroutine isobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me) - integer, intent(in) :: ph,me - real(pReal), intent(in) :: & - phi !< damage parameter - real(pReal), intent(out) :: & - localphiDot, & - dLocalphiDot_dPhi - end subroutine isobrittle_getRateAndItsTangent - - module subroutine isoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me) - integer, intent(in) :: ph,me - real(pReal), intent(in) :: & - phi !< damage parameter - real(pReal), intent(out) :: & - localphiDot, & - dLocalphiDot_dPhi - end subroutine isoductile_getRateAndItsTangent - module subroutine anisobrittle_results(phase,group) integer, intent(in) :: phase character(len=*), intent(in) :: group @@ -151,7 +114,7 @@ module subroutine damage_init do ph = 1,phases%length - Nmembers = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseID == ph) allocate(current(ph)%phi(Nmembers),source=1.0_pReal) allocate(current(ph)%d_phi_d_dot_phi(Nmembers),source=0.0_pReal) @@ -179,53 +142,30 @@ end subroutine damage_init !---------------------------------------------------------------------------------------------- !< @brief returns local part of nonlocal damage driving force !---------------------------------------------------------------------------------------------- -module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce) +module function phase_damage_phi_dot(phi,co,ce) result(phi_dot) - integer, intent(in) :: ce + integer, intent(in) :: ce,co real(pReal), intent(in) :: & phi !< damage parameter - real(pReal), intent(inout) :: & - phiDot, & - dPhiDot_dPhi - real(pReal) :: & - localphiDot, & - dLocalphiDot_dPhi + phi_dot + integer :: & ph, & - co, & - me + en - phiDot = 0.0_pReal - dPhiDot_dPhi = 0.0_pReal + ph = material_phaseID(co,ce) + en = material_phaseEntry(co,ce) - do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) - ph = material_phaseAt2(co,ce) - me = material_phasememberAt2(co,ce) + select case(phase_source(ph)) + case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ISODUCTILE_ID,DAMAGE_ANISOBRITTLE_ID,DAMAGE_ANISODUCTILE_ID) + phi_dot = 1.0_pReal & + - phi*damageState(ph)%state(1,en) + case default + phi_dot = 0.0_pReal + end select - select case(phase_source(ph)) - case (DAMAGE_ISOBRITTLE_ID) - call isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, ph, me) - - case (DAMAGE_ISODUCTILE_ID) - call isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, ph, me) - - case (DAMAGE_ANISOBRITTLE_ID) - call anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) - - case (DAMAGE_ANISODUCTILE_ID) - call anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) - - case default - localphiDot = 0.0_pReal - dLocalphiDot_dPhi = 0.0_pReal - - end select - phiDot = phiDot + localphiDot - dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi - enddo - -end subroutine phase_damage_getRateAndItsTangents +end function phase_damage_phi_dot @@ -477,7 +417,7 @@ module subroutine phase_damage_set_phi(phi,co,ce) integer, intent(in) :: ce, co - current(material_phaseAt2(co,ce))%phi(material_phaseMemberAt2(co,ce)) = phi + current(material_phaseID(co,ce))%phi(material_phaseEntry(co,ce)) = phi end subroutine phase_damage_set_phi @@ -503,4 +443,21 @@ module function damage_phi(ph,me) result(phi) end function damage_phi -end submodule damagee + +!-------------------------------------------------------------------------------------------------- +!> @brief Forward data after successful increment. +!-------------------------------------------------------------------------------------------------- +module subroutine damage_forward() + + integer :: ph + + + do ph = 1, size(damageState) + if (damageState(ph)%sizeState > 0) & + damageState(ph)%state0 = damageState(ph)%state + enddo + +end subroutine damage_forward + + +end submodule damage diff --git a/src/phase_damage_anisobrittle.f90 b/src/phase_damage_anisobrittle.f90 index 3c573f33a..ddcca9044 100644 --- a/src/phase_damage_anisobrittle.f90 +++ b/src/phase_damage_anisobrittle.f90 @@ -4,7 +4,7 @@ !> @brief material subroutine incorporating anisotropic brittle damage source mechanism !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule (phase:damagee) anisobrittle +submodule (phase:damage) anisobrittle type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & @@ -120,9 +120,6 @@ module subroutine anisobrittle_dotState(S, ph,me) S integer :: & - sourceOffset, & - damageOffset, & - homog, & i real(pReal) :: & traction_d, traction_t, traction_n, traction_crit @@ -148,29 +145,6 @@ module subroutine anisobrittle_dotState(S, ph,me) end subroutine anisobrittle_dotState -!-------------------------------------------------------------------------------------------------- -!> @brief returns local part of nonlocal damage driving force -!-------------------------------------------------------------------------------------------------- -module subroutine anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) - - integer, intent(in) :: & - ph, & - me - real(pReal), intent(in) :: & - phi - real(pReal), intent(out) :: & - localphiDot, & - dLocalphiDot_dPhi - - - dLocalphiDot_dPhi = -damageState(ph)%state(1,me) - - localphiDot = 1.0_pReal & - + dLocalphiDot_dPhi*phi - -end subroutine anisobrittle_getRateAndItsTangent - - !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- @@ -197,7 +171,7 @@ end subroutine anisobrittle_results !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) +module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) integer, intent(in) :: & ph,me @@ -253,6 +227,6 @@ module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, enddo end associate -end subroutine kinematics_cleavage_opening_LiAndItsTangent +end subroutine damage_anisobrittle_LiAndItsTangent end submodule anisobrittle diff --git a/src/phase_damage_anisoductile.f90 b/src/phase_damage_anisoductile.f90 index f047a5dac..049f148bc 100644 --- a/src/phase_damage_anisoductile.f90 +++ b/src/phase_damage_anisoductile.f90 @@ -4,7 +4,7 @@ !> @brief material subroutine incorporating anisotropic ductile damage source mechanism !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(phase:damagee) anisoductile +submodule(phase:damage) anisoductile type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & @@ -35,7 +35,7 @@ module function anisoductile_init() result(mySources) pl, & sources, & src - integer :: Ninstances,Nmembers,p + integer :: Ninstances,Nmembers,ph integer, dimension(:), allocatable :: N_sl character(len=pStringLen) :: extmsg = '' @@ -50,15 +50,15 @@ module function anisoductile_init() result(mySources) phases => config_material%get('phase') allocate(param(phases%length)) - do p = 1, phases%length - if(mySources(p)) then - phase => phases%get(p) + do ph = 1, phases%length + if(mySources(ph)) then + phase => phases%get(ph) mech => phase%get('mechanical') pl => mech%get('plastic') sources => phase%get('damage') - associate(prm => param(p)) + associate(prm => param(ph)) src => sources%get(1) N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) @@ -78,10 +78,10 @@ module function anisoductile_init() result(mySources) if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit' - Nmembers=count(material_phaseAt2==p) - call phase_allocateState(damageState(p),Nmembers,1,1,0) - damageState(p)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal) - if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' + Nmembers=count(material_phaseID==ph) + call phase_allocateState(damageState(ph),Nmembers,1,1,0) + damageState(ph)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' end associate @@ -113,29 +113,6 @@ module subroutine anisoductile_dotState(ph,me) end subroutine anisoductile_dotState -!-------------------------------------------------------------------------------------------------- -!> @brief returns local part of nonlocal damage driving force -!-------------------------------------------------------------------------------------------------- -module subroutine anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me) - - integer, intent(in) :: & - ph, & - me - real(pReal), intent(in) :: & - phi - real(pReal), intent(out) :: & - localphiDot, & - dLocalphiDot_dPhi - - - dLocalphiDot_dPhi = -damageState(ph)%state(1,me) - - localphiDot = 1.0_pReal & - + dLocalphiDot_dPhi*phi - -end subroutine anisoductile_getRateAndItsTangent - - !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_damage_isobrittle.f90 b/src/phase_damage_isobrittle.f90 index 3c7866e66..e88a87352 100644 --- a/src/phase_damage_isobrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -4,7 +4,7 @@ !> @brief material subroutine incoprorating isotropic brittle damage source mechanism !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(phase:damagee) isobrittle +submodule(phase:damage) isobrittle type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & @@ -31,7 +31,7 @@ module function isobrittle_init() result(mySources) phase, & sources, & src - integer :: Nmembers,p + integer :: Nmembers,ph character(len=pStringLen) :: extmsg = '' @@ -45,12 +45,12 @@ module function isobrittle_init() result(mySources) phases => config_material%get('phase') allocate(param(phases%length)) - do p = 1, phases%length - if(mySources(p)) then - phase => phases%get(p) + do ph = 1, phases%length + if(mySources(ph)) then + phase => phases%get(ph) sources => phase%get('damage') - associate(prm => param(p)) + associate(prm => param(ph)) src => sources%get(1) prm%W_crit = src%get_asFloat('W_crit') @@ -64,10 +64,10 @@ module function isobrittle_init() result(mySources) ! sanity checks if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit' - Nmembers = count(material_phaseAt2==p) - call phase_allocateState(damageState(p),Nmembers,1,1,1) - damageState(p)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal) - if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' + Nmembers = count(material_phaseID==ph) + call phase_allocateState(damageState(ph),Nmembers,1,1,1) + damageState(ph)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' end associate @@ -113,29 +113,6 @@ module subroutine isobrittle_deltaState(C, Fe, ph,me) end subroutine isobrittle_deltaState -!-------------------------------------------------------------------------------------------------- -!> @brief returns local part of nonlocal damage driving force -!-------------------------------------------------------------------------------------------------- -module subroutine isobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) - - integer, intent(in) :: & - ph, me - real(pReal), intent(in) :: & - phi - real(pReal), intent(out) :: & - localphiDot, & - dLocalphiDot_dPhi - - - associate(prm => param(ph)) - localphiDot = 1.0_pReal & - - phi*damageState(ph)%state(1,me) - dLocalphiDot_dPhi = - damageState(ph)%state(1,me) - end associate - -end subroutine isobrittle_getRateAndItsTangent - - !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_damage_isoductile.f90 b/src/phase_damage_isoductile.f90 index b5ca9a1c2..997b948fe 100644 --- a/src/phase_damage_isoductile.f90 +++ b/src/phase_damage_isoductile.f90 @@ -4,7 +4,7 @@ !> @brief material subroutine incorporating isotropic ductile damage source mechanism !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(phase:damagee) isoductile +submodule(phase:damage) isoductile type:: tParameters !< container type for internal constitutive parameters real(pReal) :: & @@ -33,7 +33,7 @@ module function isoductile_init() result(mySources) phase, & sources, & src - integer :: Ninstances,Nmembers,p + integer :: Ninstances,Nmembers,ph character(len=pStringLen) :: extmsg = '' @@ -47,12 +47,12 @@ module function isoductile_init() result(mySources) phases => config_material%get('phase') allocate(param(phases%length)) - do p = 1, phases%length - if(mySources(p)) then - phase => phases%get(p) + do ph = 1, phases%length + if(mySources(ph)) then + phase => phases%get(ph) sources => phase%get('damage') - associate(prm => param(p)) + associate(prm => param(ph)) src => sources%get(1) prm%q = src%get_asFloat('q') @@ -68,10 +68,10 @@ module function isoductile_init() result(mySources) if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit' - Nmembers=count(material_phaseAt2==p) - call phase_allocateState(damageState(p),Nmembers,1,1,0) - damageState(p)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal) - if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' + Nmembers=count(material_phaseID==ph) + call phase_allocateState(damageState(ph),Nmembers,1,1,0) + damageState(ph)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' end associate @@ -103,29 +103,6 @@ module subroutine isoductile_dotState(ph, me) end subroutine isoductile_dotState -!-------------------------------------------------------------------------------------------------- -!> @brief returns local part of nonlocal damage driving force -!-------------------------------------------------------------------------------------------------- -module subroutine isoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) - - integer, intent(in) :: & - ph, & - me - real(pReal), intent(in) :: & - phi - real(pReal), intent(out) :: & - localphiDot, & - dLocalphiDot_dPhi - - - dLocalphiDot_dPhi = -damageState(ph)%state(1,me) - - localphiDot = 1.0_pReal & - + dLocalphiDot_dPhi*phi - -end subroutine isoductile_getRateAndItsTangent - - !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index c507ce20c..9f2bc4df2 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -1060,26 +1060,6 @@ subroutine crystallite_results(group,ph) end subroutine crystallite_results -!-------------------------------------------------------------------------------------------------- -!> @brief Wind homog inc forward. -!-------------------------------------------------------------------------------------------------- -module subroutine mechanical_windForward(ph,me) - - integer, intent(in) :: ph, me - - - phase_mechanical_Fp0(ph)%data(1:3,1:3,me) = phase_mechanical_Fp(ph)%data(1:3,1:3,me) - phase_mechanical_Fi0(ph)%data(1:3,1:3,me) = phase_mechanical_Fi(ph)%data(1:3,1:3,me) - phase_mechanical_F0(ph)%data(1:3,1:3,me) = phase_mechanical_F(ph)%data(1:3,1:3,me) - phase_mechanical_Li0(ph)%data(1:3,1:3,me) = phase_mechanical_Li(ph)%data(1:3,1:3,me) - phase_mechanical_Lp0(ph)%data(1:3,1:3,me) = phase_mechanical_Lp(ph)%data(1:3,1:3,me) - phase_mechanical_S0(ph)%data(1:3,1:3,me) = phase_mechanical_S(ph)%data(1:3,1:3,me) - - plasticState(ph)%State0(:,me) = plasticState(ph)%state(:,me) - -end subroutine mechanical_windForward - - !-------------------------------------------------------------------------------------------------- !> @brief Forward data after successful increment. ! ToDo: Any guessing for the current states possible? @@ -1213,8 +1193,6 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) if (todo) then subF = subF0 & + subStep * (phase_mechanical_F(ph)%data(1:3,1:3,me) - phase_mechanical_F0(ph)%data(1:3,1:3,me)) - phase_mechanical_Fe(ph)%data(1:3,1:3,me) = matmul(subF,math_inv33(matmul(phase_mechanical_Fi(ph)%data(1:3,1:3,me), & - phase_mechanical_Fp(ph)%data(1:3,1:3,me)))) converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el) converged_ = converged_ .and. .not. integrateDamageState(subStep * dt,co,ip,el) endif @@ -1237,9 +1215,9 @@ module subroutine mechanical_restore(ce,includeL) co, ph, me - do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce)) - ph = material_phaseAt2(co,ce) - me = material_phaseMemberAt2(co,ce) + do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) + ph = material_phaseID(co,ce) + me = material_phaseEntry(co,ce) if (includeL) then phase_mechanical_Lp(ph)%data(1:3,1:3,me) = phase_mechanical_Lp0(ph)%data(1:3,1:3,me) phase_mechanical_Li(ph)%data(1:3,1:3,me) = phase_mechanical_Li0(ph)%data(1:3,1:3,me) @@ -1287,8 +1265,8 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF) logical :: error - ph = material_phaseAt2(co,ce) - me = material_phaseMemberAt2(co,ce) + ph = material_phaseID(co,ce) + me = material_phaseEntry(co,ce) call phase_hooke_SandItsTangents(devNull,dSdFe,dSdFi, & phase_mechanical_Fe(ph)%data(1:3,1:3,me), & @@ -1443,20 +1421,6 @@ module function mechanical_L_p(ph,me) result(L_p) end function mechanical_L_p -!---------------------------------------------------------------------------------------------- -!< @brief Get deformation gradient (for use by homogenization) -!---------------------------------------------------------------------------------------------- -module function phase_mechanical_getF(co,ce) result(F) - - integer, intent(in) :: co, ce - real(pReal), dimension(3,3) :: F - - - F = phase_mechanical_F(material_phaseAt2(co,ce))%data(1:3,1:3,material_phaseMemberAt2(co,ce)) - -end function phase_mechanical_getF - - !---------------------------------------------------------------------------------------------- !< @brief Get elastic deformation gradient (for use by non-mech physics) !---------------------------------------------------------------------------------------------- @@ -1471,31 +1435,46 @@ module function mechanical_F_e(ph,me) result(F_e) end function mechanical_F_e - !---------------------------------------------------------------------------------------------- !< @brief Get second Piola-Kichhoff stress (for use by homogenization) !---------------------------------------------------------------------------------------------- -module function phase_mechanical_getP(co,ce) result(P) +module function phase_P(co,ce) result(P) integer, intent(in) :: co, ce real(pReal), dimension(3,3) :: P - P = phase_mechanical_P(material_phaseAt2(co,ce))%data(1:3,1:3,material_phaseMemberAt2(co,ce)) + P = phase_mechanical_P(material_phaseID(co,ce))%data(1:3,1:3,material_phaseEntry(co,ce)) -end function phase_mechanical_getP +end function phase_P -! setter for homogenization -module subroutine phase_mechanical_setF(F,co,ce) +!---------------------------------------------------------------------------------------------- +!< @brief Get deformation gradient (for use by homogenization) +!---------------------------------------------------------------------------------------------- +module function phase_F(co,ce) result(F) + + integer, intent(in) :: co, ce + real(pReal), dimension(3,3) :: F + + + F = phase_mechanical_F(material_phaseID(co,ce))%data(1:3,1:3,material_phaseEntry(co,ce)) + +end function phase_F + + +!---------------------------------------------------------------------------------------------- +!< @brief Set deformation gradient (for use by homogenization) +!---------------------------------------------------------------------------------------------- +module subroutine phase_set_F(F,co,ce) real(pReal), dimension(3,3), intent(in) :: F integer, intent(in) :: co, ce - phase_mechanical_F(material_phaseAt2(co,ce))%data(1:3,1:3,material_phaseMemberAt2(co,ce)) = F + phase_mechanical_F(material_phaseID(co,ce))%data(1:3,1:3,material_phaseEntry(co,ce)) = F -end subroutine phase_mechanical_setF +end subroutine phase_set_F end submodule mechanical diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90 index 955f6ca5d..7c4e70c1c 100644 --- a/src/phase_mechanical_eigen.f90 +++ b/src/phase_mechanical_eigen.f90 @@ -9,13 +9,13 @@ submodule(phase:mechanical) eigen model_damage interface - module function kinematics_cleavage_opening_init() result(myKinematics) + module function damage_anisobrittle_init() result(myKinematics) logical, dimension(:), allocatable :: myKinematics - end function kinematics_cleavage_opening_init + end function damage_anisobrittle_init - module function kinematics_slipplane_opening_init() result(myKinematics) + module function damage_isoductile_init() result(myKinematics) logical, dimension(:), allocatable :: myKinematics - end function kinematics_slipplane_opening_init + end function damage_isoductile_init module function thermalexpansion_init(kinematics_length) result(myKinematics) integer, intent(in) :: kinematics_length @@ -46,7 +46,6 @@ module subroutine eigendeformation_init(phases) class(tNode), pointer :: & phase, & kinematics, & - damage, & mechanics print'(/,a)', ' <<<+- phase:mechanical:eigen init -+>>>' @@ -70,8 +69,8 @@ module subroutine eigendeformation_init(phases) allocate(model_damage(phases%length), source = KINEMATICS_UNDEFINED_ID) - where(kinematics_cleavage_opening_init()) model_damage = KINEMATICS_cleavage_opening_ID - where(kinematics_slipplane_opening_init()) model_damage = KINEMATICS_slipplane_opening_ID + where(damage_anisobrittle_init()) model_damage = KINEMATICS_cleavage_opening_ID + where(damage_isoductile_init()) model_damage = KINEMATICS_slipplane_opening_ID end subroutine eigendeformation_init @@ -198,12 +197,12 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & select case (model_damage(ph)) case (KINEMATICS_cleavage_opening_ID) - call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me) + call damage_anisobrittle_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me) Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS active = .true. case (KINEMATICS_slipplane_opening_ID) - call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me) + call damage_isoductile_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me) Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS active = .true. diff --git a/src/phase_mechanical_eigen_cleavageopening.f90 b/src/phase_mechanical_eigen_cleavageopening.f90 index bccf1cb5f..4243d09a4 100644 --- a/src/phase_mechanical_eigen_cleavageopening.f90 +++ b/src/phase_mechanical_eigen_cleavageopening.f90 @@ -13,7 +13,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function kinematics_cleavage_opening_init() result(myKinematics) +module function damage_anisobrittle_init() result(myKinematics) logical, dimension(:), allocatable :: myKinematics @@ -24,7 +24,7 @@ module function kinematics_cleavage_opening_init() result(myKinematics) print'(/,a)', ' <<<+- phase:mechanical:eigen:cleavageopening init -+>>>' print'(a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT) -end function kinematics_cleavage_opening_init +end function damage_anisobrittle_init end submodule cleavageopening diff --git a/src/phase_mechanical_eigen_slipplaneopening.f90 b/src/phase_mechanical_eigen_slipplaneopening.f90 index 18d99f750..2fd14700d 100644 --- a/src/phase_mechanical_eigen_slipplaneopening.f90 +++ b/src/phase_mechanical_eigen_slipplaneopening.f90 @@ -6,7 +6,7 @@ !-------------------------------------------------------------------------------------------------- submodule(phase:eigen) slipplaneopening - integer, dimension(:), allocatable :: kinematics_slipplane_opening_instance + integer, dimension(:), allocatable :: damage_isoductile_instance type :: tParameters !< container type for internal constitutive parameters integer :: & @@ -32,7 +32,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function kinematics_slipplane_opening_init() result(myKinematics) +module function damage_isoductile_init() result(myKinematics) logical, dimension(:), allocatable :: myKinematics @@ -107,13 +107,13 @@ module function kinematics_slipplane_opening_init() result(myKinematics) enddo -end function kinematics_slipplane_opening_init +end function damage_isoductile_init !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) +module subroutine damage_isoductile_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) integer, intent(in) :: & ph, me @@ -179,6 +179,6 @@ module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S end associate -end subroutine kinematics_slipplane_opening_LiAndItsTangent +end subroutine damage_isoductile_LiAndItsTangent end submodule slipplaneopening diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index b7ade48a8..2e7130756 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -220,7 +220,7 @@ module function plastic_dislotungsten_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nmembers = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseID == ph) sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl sizeState = sizeDotState diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 1e67db59e..43be614ac 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -406,7 +406,7 @@ module function plastic_dislotwin_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nmembers = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseID == ph) sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl & + size(['f_tw']) * prm%sum_N_tw & + size(['f_tr']) * prm%sum_N_tr diff --git a/src/phase_mechanical_plastic_isotropic.f90 b/src/phase_mechanical_plastic_isotropic.f90 index 5a37369f0..548e2fd75 100644 --- a/src/phase_mechanical_plastic_isotropic.f90 +++ b/src/phase_mechanical_plastic_isotropic.f90 @@ -119,7 +119,7 @@ module function plastic_isotropic_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nmembers = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseID == ph) sizeDotState = size(['xi ','gamma']) sizeState = sizeDotState diff --git a/src/phase_mechanical_plastic_kinehardening.f90 b/src/phase_mechanical_plastic_kinehardening.f90 index 0c02e11db..936a7f884 100644 --- a/src/phase_mechanical_plastic_kinehardening.f90 +++ b/src/phase_mechanical_plastic_kinehardening.f90 @@ -165,7 +165,7 @@ module function plastic_kinehardening_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nmembers = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseID == ph) sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%sum_N_sl !ToDo: adjust names like in material.yaml sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names like in material.yaml sizeState = sizeDotState + sizeDeltaState diff --git a/src/phase_mechanical_plastic_none.f90 b/src/phase_mechanical_plastic_none.f90 index 28e9fbc7c..544de31b1 100644 --- a/src/phase_mechanical_plastic_none.f90 +++ b/src/phase_mechanical_plastic_none.f90 @@ -30,7 +30,7 @@ module function plastic_none_init() result(myPlasticity) phases => config_material%get('phase') do ph = 1, phases%length if(.not. myPlasticity(ph)) cycle - call phase_allocateState(plasticState(ph),count(material_phaseAt2 == ph),0,0,0) + call phase_allocateState(plasticState(ph),count(material_phaseID == ph),0,0,0) enddo end function plastic_none_init diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index d86f7dc7e..ba7f76881 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -398,7 +398,7 @@ module function plastic_nonlocal_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nmembers = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseID == ph) sizeDotState = size([ 'rhoSglEdgePosMobile ','rhoSglEdgeNegMobile ', & 'rhoSglScrewPosMobile ','rhoSglScrewNegMobile ', & 'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', & @@ -527,7 +527,7 @@ module function plastic_nonlocal_init() result(myPlasticity) if(.not. myPlasticity(ph)) cycle phase => phases%get(ph) - Nmembers = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseID == ph) l = 0 do t = 1,4 do s = 1,param(ph)%sum_N_sl @@ -1824,9 +1824,9 @@ subroutine storeGeometry(ph) V = reshape(IPvolume,[product(shape(IPvolume))]) - do ce = 1, size(material_homogenizationMemberAt2,1) + do ce = 1, size(material_homogenizationEntry,1) do co = 1, homogenization_maxNconstituents - if (material_phaseAt2(co,ce) == ph) geom(ph)%V_0(material_phaseMemberAt2(co,ce)) = V(ce) + if (material_phaseID(co,ce) == ph) geom(ph)%V_0(material_phaseEntry(co,ce)) = V(ce) enddo enddo diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90 index 6be46f6a1..77c47907f 100644 --- a/src/phase_mechanical_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -223,7 +223,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nmembers = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseID == ph) sizeDotState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl & + size(['xi_tw ','gamma_tw']) * prm%sum_N_tw sizeState = sizeDotState diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index 21ec93c9c..1e157f338 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -89,7 +89,7 @@ module subroutine thermal_init(phases) allocate(thermal_Nsources(phases%length),source = 0) do ph = 1, phases%length - Nmembers = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseID == ph) allocate(current(ph)%T(Nmembers),source=300.0_pReal) allocate(current(ph)%dot_T(Nmembers),source=0.0_pReal) phase => phases%get(ph) @@ -268,8 +268,8 @@ module subroutine phase_thermal_setField(T,dot_T, co,ce) integer, intent(in) :: ce, co - current(material_phaseAt2(co,ce))%T(material_phaseMemberAt2(co,ce)) = T - current(material_phaseAt2(co,ce))%dot_T(material_phaseMemberAt2(co,ce)) = dot_T + current(material_phaseID(co,ce))%T(material_phaseEntry(co,ce)) = T + current(material_phaseID(co,ce))%dot_T(material_phaseEntry(co,ce)) = dot_T end subroutine phase_thermal_setField diff --git a/src/phase_thermal_dissipation.f90 b/src/phase_thermal_dissipation.f90 index d3e7094a1..b16ed3cf7 100644 --- a/src/phase_thermal_dissipation.f90 +++ b/src/phase_thermal_dissipation.f90 @@ -54,7 +54,7 @@ module function dissipation_init(source_length) result(mySources) src => sources%get(so) prm%kappa = src%get_asFloat('kappa') - Nmembers = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseID == ph) call phase_allocateState(thermalState(ph)%p(so),Nmembers,0,0,0) end associate diff --git a/src/phase_thermal_externalheat.f90 b/src/phase_thermal_externalheat.f90 index 760326044..d5bbc7c38 100644 --- a/src/phase_thermal_externalheat.f90 +++ b/src/phase_thermal_externalheat.f90 @@ -67,7 +67,7 @@ module function externalheat_init(source_length) result(mySources) prm%f_T = src%get_as1dFloat('f_T',requiredSize = size(prm%t_n)) - Nmembers = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseID == ph) call phase_allocateState(thermalState(ph)%p(so),Nmembers,1,1,0) end associate endif