diff --git a/PRIVATE b/PRIVATE index ec615d249..6db5f4666 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit ec615d249d39e5d01446b01ab9a5b7e7601340ad +Subproject commit 6db5f4666fc651b4de3b44ceaed3f2b848170ac9 diff --git a/VERSION b/VERSION index fe11c9855..ff156acbf 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-1732-gab88ffd2 +v2.0.3-1747-ga2e8e5b1 diff --git a/processing/post/addCauchy.py b/processing/post/addCauchy.py index afc5a57be..8bbc7fb0e 100755 --- a/processing/post/addCauchy.py +++ b/processing/post/addCauchy.py @@ -42,8 +42,8 @@ for name in filenames: table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) table.add('Cauchy', - damask.mechanics.Cauchy(table.get(options.defgrad).reshape(-1,3,3), - table.get(options.stress ).reshape(-1,3,3)).reshape(-1,9), + damask.mechanics.Cauchy(table.get(options.stress ).reshape(-1,3,3), + table.get(options.defgrad).reshape(-1,3,3)).reshape(-1,9), scriptID+' '+' '.join(sys.argv[1:])) table.to_ASCII(sys.stdout if name is None else name) diff --git a/processing/post/addPK2.py b/processing/post/addPK2.py index 185160d79..2894a5a90 100755 --- a/processing/post/addPK2.py +++ b/processing/post/addPK2.py @@ -43,8 +43,8 @@ for name in filenames: table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) table.add('S', - damask.mechanics.PK2(table.get(options.defgrad).reshape(-1,3,3), - table.get(options.stress ).reshape(-1,3,3)).reshape(-1,9), + damask.mechanics.PK2(table.get(options.stress ).reshape(-1,3,3), + table.get(options.defgrad).reshape(-1,3,3)).reshape(-1,9), scriptID+' '+' '.join(sys.argv[1:])) table.to_ASCII(sys.stdout if name is None else name) diff --git a/processing/post/addSpectralDecomposition.py b/processing/post/addSpectralDecomposition.py index c711c6a45..5909ca147 100755 --- a/processing/post/addSpectralDecomposition.py +++ b/processing/post/addSpectralDecomposition.py @@ -2,6 +2,7 @@ import os import sys +from io import StringIO from optparse import OptionParser import numpy as np @@ -33,69 +34,27 @@ parser.add_option('--no-check', parser.set_defaults(rh = True, ) + (options,filenames) = parser.parse_args() - -if options.tensor is None: - parser.error('no data column specified.') - -# --- loop over input files ------------------------------------------------------------------------- - if filenames == []: filenames = [None] for name in filenames: - try: - table = damask.ASCIItable(name = name, - buffered = False) - except: continue - damask.util.report(scriptName,name) + damask.util.report(scriptName,name) -# ------------------------------------------ read header ------------------------------------------ + table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) - table.head_read() + for tensor in options.tensor: + t = table.get(tensor).reshape(-1,3,3) + (u,v) = np.linalg.eigh(damask.mechanics.symmetric(t)) + if options.rh: v[np.linalg.det(v) < 0.0,:,2] *= -1.0 -# ------------------------------------------ assemble header 1 ------------------------------------ + for i,o in enumerate(['Min','Mid','Max']): + table.add('eigval{}({})'.format(o,tensor),u[:,i], + scriptID+' '+' '.join(sys.argv[1:])) - items = { - 'tensor': {'dim': 9, 'shape': [3,3], 'labels':options.tensor, 'column': []}, - } - errors = [] - remarks = [] - - for type, data in items.items(): - for what in data['labels']: - dim = table.label_dimension(what) - if dim != data['dim']: remarks.append('column {} is not a {}...'.format(what,type)) - else: - items[type]['column'].append(table.label_index(what)) - for order in ['Min','Mid','Max']: - table.labels_append(['eigval{}({})'.format(order,what)]) # extend ASCII header with new labels - for order in ['Min','Mid','Max']: - table.labels_append(['{}_eigvec{}({})'.format(i+1,order,what) for i in range(3)]) # extend ASCII header with new labels - - if remarks != []: damask.util.croak(remarks) - if errors != []: - damask.util.croak(errors) - table.close(dismiss = True) - continue - -# ------------------------------------------ assemble header 2 ------------------------------------ - - table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) - table.head_write() - -# ------------------------------------------ process data ----------------------------------------- - - outputAlive = True - while outputAlive and table.data_read(): # read next data line of ASCII table - for type, data in items.items(): - for column in data['column']: - (u,v) = np.linalg.eigh(np.array(list(map(float,table.data[column:column+data['dim']]))).reshape(data['shape'])) - if options.rh and np.dot(np.cross(v[:,0], v[:,1]), v[:,2]) < 0.0 : v[:, 2] *= -1.0 # ensure right-handed eigenvector basis - table.data_append(list(u)) # vector of max,mid,min eigval - table.data_append(list(v.transpose().reshape(data['dim']))) # 3x3=9 combo vector of max,mid,min eigvec coordinates - outputAlive = table.data_write() # output processed line in accordance with column labeling - -# ------------------------------------------ output finalization ----------------------------------- - - table.close() # close input ASCII table (works for stdin) + for i,o in enumerate(['Min','Mid','Max']): + table.add('eigvec{}({})'.format(o,tensor),v[:,:,i], + scriptID+' '+' '.join(sys.argv[1:])) + + table.to_ASCII(sys.stdout if name is None else name) diff --git a/python/damask/__init__.py b/python/damask/__init__.py index 30fb37ced..2961ba6e3 100644 --- a/python/damask/__init__.py +++ b/python/damask/__init__.py @@ -15,6 +15,7 @@ from .config import Material # noqa from .colormaps import Colormap, Color # noqa from .orientation import Symmetry, Lattice, Rotation, Orientation # noqa from .dadf5 import DADF5 # noqa +from .dadf5 import DADF5 as Result # noqa from .geom import Geom # noqa from .solver import Solver # noqa diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 819b5603e..ccbd3be03 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -1,4 +1,4 @@ -from queue import Queue +from queue import Queue import re import glob import os @@ -11,998 +11,1165 @@ import numpy as np from . import util from . import version from . import mechanics +from . import Rotation +from . import Orientation +from . import Environment +from . import grid_filters -# ------------------------------------------------------------------ class DADF5(): - """ - Read and write to DADF5 files. - - DADF5 files contain DAMASK results. - """ - -# ------------------------------------------------------------------ - def __init__(self,fname): """ - Opens an existing DADF5 file. - - Parameters - ---------- - fname : str - name of the DADF5 file to be openend. - + Read and write to DADF5 files. + + DADF5 files contain DAMASK results. """ - with h5py.File(fname,'r') as f: - - try: - self.version_major = f.attrs['DADF5_version_major'] - self.version_minor = f.attrs['DADF5_version_minor'] - except KeyError: - self.version_major = f.attrs['DADF5-major'] - self.version_minor = f.attrs['DADF5-minor'] - if self.version_major != 0 or not 2 <= self.version_minor <= 6: - raise TypeError('Unsupported DADF5 version {}.{} '.format(f.attrs['DADF5_version_major'], - f.attrs['DADF5_version_minor'])) - - self.structured = 'grid' in f['geometry'].attrs.keys() - - if self.structured: - self.grid = f['geometry'].attrs['grid'] - self.size = f['geometry'].attrs['size'] - if self.version_major == 0 and self.version_minor >= 5: - self.origin = f['geometry'].attrs['origin'] + def __init__(self,fname): + """ + Opens an existing DADF5 file. - - r=re.compile('inc[0-9]+') - increments_unsorted = {int(i[3:]):i for i in f.keys() if r.match(i)} - self.increments = [increments_unsorted[i] for i in sorted(increments_unsorted)] - self.times = [round(f[i].attrs['time/s'],12) for i in self.increments] + Parameters + ---------- + fname : str + name of the DADF5 file to be openend. - self.Nmaterialpoints, self.Nconstituents = np.shape(f['mapping/cellResults/constituent']) - self.materialpoints = [m.decode() for m in np.unique(f['mapping/cellResults/materialpoint']['Name'])] - self.constituents = [c.decode() for c in np.unique(f['mapping/cellResults/constituent'] ['Name'])] + """ + with h5py.File(fname,'r') as f: - self.con_physics = [] - for c in self.constituents: - self.con_physics += f['/'.join([self.increments[0],'constituent',c])].keys() - self.con_physics = list(set(self.con_physics)) # make unique + try: + self.version_major = f.attrs['DADF5_version_major'] + self.version_minor = f.attrs['DADF5_version_minor'] + except KeyError: + self.version_major = f.attrs['DADF5-major'] + self.version_minor = f.attrs['DADF5-minor'] - self.mat_physics = [] - for m in self.materialpoints: - self.mat_physics += f['/'.join([self.increments[0],'materialpoint',m])].keys() - self.mat_physics = list(set(self.mat_physics)) # make unique + if self.version_major != 0 or not 2 <= self.version_minor <= 6: + raise TypeError('Unsupported DADF5 version {}.{} '.format(f.attrs['DADF5_version_major'], + f.attrs['DADF5_version_minor'])) - self.visible= {'increments': self.increments, - 'constituents': self.constituents, - 'materialpoints': self.materialpoints, - 'constituent': range(self.Nconstituents), # ToDo: stupid naming - 'con_physics': self.con_physics, - 'mat_physics': self.mat_physics} - - self.fname = fname + self.structured = 'grid' in f['geometry'].attrs.keys() + + if self.structured: + self.grid = f['geometry'].attrs['grid'] + self.size = f['geometry'].attrs['size'] + if self.version_major == 0 and self.version_minor >= 5: + self.origin = f['geometry'].attrs['origin'] - def __manage_visible(self,datasets,what,action): - """ - Manages the visibility of the groups. - - Parameters - ---------- - datasets : list of str or Boolean - name of datasets as list, supports ? and * wildcards. - True is equivalent to [*], False is equivalent to [] - what : str - attribute to change (must be in self.visible) - action : str - select from 'set', 'add', and 'del' + r=re.compile('inc[0-9]+') + increments_unsorted = {int(i[3:]):i for i in f.keys() if r.match(i)} + self.increments = [increments_unsorted[i] for i in sorted(increments_unsorted)] + self.times = [round(f[i].attrs['time/s'],12) for i in self.increments] - """ - # allow True/False and string arguments - if datasets is True: - datasets = ['*'] - elif datasets is False: - datasets = [] - choice = [datasets] if isinstance(datasets,str) else datasets - - valid = [e for e_ in [glob.fnmatch.filter(getattr(self,what),s) for s in choice] for e in e_] - existing = set(self.visible[what]) + self.Nmaterialpoints, self.Nconstituents = np.shape(f['mapping/cellResults/constituent']) + self.materialpoints = [m.decode() for m in np.unique(f['mapping/cellResults/materialpoint']['Name'])] + self.constituents = [c.decode() for c in np.unique(f['mapping/cellResults/constituent'] ['Name'])] - if action == 'set': - self.visible[what] = valid - elif action == 'add': - self.visible[what] = list(existing.union(valid)) - elif action == 'del': - self.visible[what] = list(existing.difference_update(valid)) - - - def __time_to_inc(self,start,end): - selected = [] - for i,time in enumerate(self.times): - if start <= time <= end: - selected.append(self.increments[i]) - return selected + self.con_physics = [] + for c in self.constituents: + self.con_physics += f['/'.join([self.increments[0],'constituent',c])].keys() + self.con_physics = list(set(self.con_physics)) # make unique + + self.mat_physics = [] + for m in self.materialpoints: + self.mat_physics += f['/'.join([self.increments[0],'materialpoint',m])].keys() + self.mat_physics = list(set(self.mat_physics)) # make unique + + self.selection= {'increments': self.increments, + 'constituents': self.constituents, + 'materialpoints': self.materialpoints, + 'con_physics': self.con_physics, + 'mat_physics': self.mat_physics} + + self.fname = fname - def set_by_time(self,start,end): - """ - Set active increments based on start and end time. + def __manage_visible(self,datasets,what,action): + """ + Manages the visibility of the groups. - Parameters - ---------- - start : float - start time (included) - end : float - end time (included) + Parameters + ---------- + datasets : list of str or Boolean + name of datasets as list, supports ? and * wildcards. + True is equivalent to [*], False is equivalent to [] + what : str + attribute to change (must be in self.selection) + action : str + select from 'set', 'add', and 'del' - """ - self.__manage_visible(self.__time_to_inc(start,end),'increments','set') + """ + # allow True/False and string arguments + if datasets is True: + datasets = ['*'] + elif datasets is False: + datasets = [] + choice = [datasets] if isinstance(datasets,str) else datasets + + valid = [e for e_ in [glob.fnmatch.filter(getattr(self,what),s) for s in choice] for e in e_] + existing = set(self.selection[what]) + + if action == 'set': + self.selection[what] = valid + elif action == 'add': + add=existing.union(valid) + add_sorted=sorted(add, key=lambda x: int("".join([i for i in x if i.isdigit()]))) + self.selection[what] = add_sorted + elif action == 'del': + diff=existing.difference(valid) + diff_sorted=sorted(diff, key=lambda x: int("".join([i for i in x if i.isdigit()]))) + self.selection[what] = diff_sorted + + def __time_to_inc(self,start,end): + selected = [] + for i,time in enumerate(self.times): + if start <= time <= end: + selected.append(self.increments[i]) + return selected - def add_by_time(self,start,end): - """ - Add to active increments based on start and end time. + def set_by_time(self,start,end): + """ + Set active increments based on start and end time. - Parameters - ---------- - start : float - start time (included) - end : float - end time (included) + Parameters + ---------- + start : float + start time (included) + end : float + end time (included) - """ - self.__manage_visible(self.__time_to_inc(start,end),'increments','add') + """ + self.__manage_visible(self.__time_to_inc(start,end),'increments','set') - def del_by_time(self,start,end): - """ - Delete from active increments based on start and end time. + def add_by_time(self,start,end): + """ + Add to active increments based on start and end time. - Parameters - ---------- - start : float - start time (included) - end : float - end time (included) + Parameters + ---------- + start : float + start time (included) + end : float + end time (included) - """ - self.__manage_visible(self.__time_to_inc(start,end),'increments','del') - - - def set_by_increment(self,start,end): - """ - Set active time increments based on start and end increment. - - Parameters - ---------- - start : int - start increment (included) - end : int - end increment (included) - - """ - if self.version_minor >= 4: - self.__manage_visible([ 'inc{}'.format(i) for i in range(start,end+1)],'increments','set') - else: - self.__manage_visible(['inc{:05d}'.format(i) for i in range(start,end+1)],'increments','set') + """ + self.__manage_visible(self.__time_to_inc(start,end),'increments','add') - def add_by_increment(self,start,end): - """ - Add to active time increments based on start and end increment. + def del_by_time(self,start,end): + """ + Delete from active increments based on start and end time. - Parameters - ---------- - start : int - start increment (included) - end : int - end increment (included) + Parameters + ---------- + start : float + start time (included) + end : float + end time (included) - """ - if self.version_minor >= 4: - self.__manage_visible([ 'inc{}'.format(i) for i in range(start,end+1)],'increments','add') - else: - self.__manage_visible(['inc{:05d}'.format(i) for i in range(start,end+1)],'increments','add') + """ + self.__manage_visible(self.__time_to_inc(start,end),'increments','del') - def del_by_increment(self,start,end): - """ - Delet from active time increments based on start and end increment. + def set_by_increment(self,start,end): + """ + Set active time increments based on start and end increment. - Parameters - ---------- - start : int - start increment (included) - end : int - end increment (included) + Parameters + ---------- + start : int + start increment (included) + end : int + end increment (included) - """ - if self.version_minor >= 4: - self.__manage_visible([ 'inc{}'.format(i) for i in range(start,end+1)],'increments','del') - else: - self.__manage_visible(['inc{:05d}'.format(i) for i in range(start,end+1)],'increments','del') + """ + if self.version_minor >= 4: + self.__manage_visible([ 'inc{}'.format(i) for i in range(start,end+1)],'increments','set') + else: + self.__manage_visible(['inc{:05d}'.format(i) for i in range(start,end+1)],'increments','set') - def iter_visible(self,what): - """ - Iterate over visible items by setting each one visible. + def add_by_increment(self,start,end): + """ + Add to active time increments based on start and end increment. - Parameters - ---------- - what : str - attribute to change (must be in self.visible) + Parameters + ---------- + start : int + start increment (included) + end : int + end increment (included) - """ - datasets = self.visible[what] - last_datasets = datasets.copy() - for dataset in datasets: - if last_datasets != self.visible[what]: + """ + if self.version_minor >= 4: + self.__manage_visible([ 'inc{}'.format(i) for i in range(start,end+1)],'increments','add') + else: + self.__manage_visible(['inc{:05d}'.format(i) for i in range(start,end+1)],'increments','add') + + + def del_by_increment(self,start,end): + """ + Delet from active time increments based on start and end increment. + + Parameters + ---------- + start : int + start increment (included) + end : int + end increment (included) + + """ + if self.version_minor >= 4: + self.__manage_visible([ 'inc{}'.format(i) for i in range(start,end+1)],'increments','del') + else: + self.__manage_visible(['inc{:05d}'.format(i) for i in range(start,end+1)],'increments','del') + + + def iter_visible(self,what): + """ + Iterate over visible items by setting each one visible. + + Parameters + ---------- + what : str + attribute to change (must be in self.selection) + + """ + datasets = self.selection[what] + last_datasets = datasets.copy() + for dataset in datasets: + if last_datasets != self.selection[what]: + self.__manage_visible(datasets,what,'set') + raise Exception + self.__manage_visible(dataset,what,'set') + last_datasets = self.selection[what] + yield dataset self.__manage_visible(datasets,what,'set') - raise Exception - self.__manage_visible(dataset,what,'set') - last_datasets = self.visible[what] - yield dataset - self.__manage_visible(datasets,what,'set') - - - def set_visible(self,what,datasets): - """ - Set active groups. - - Parameters - ---------- - datasets : list of str or Boolean - name of datasets as list, supports ? and * wildcards. - True is equivalent to [*], False is equivalent to [] - what : str - attribute to change (must be in self.visible) - """ - self.__manage_visible(datasets,what,'set') + def set_visible(self,what,datasets): + """ + Set active groups. - def add_visible(self,what,datasets): - """ - Add to active groups. - - Parameters - ---------- - datasets : list of str or Boolean - name of datasets as list, supports ? and * wildcards. - True is equivalent to [*], False is equivalent to [] - what : str - attribute to change (must be in self.visible) + Parameters + ---------- + datasets : list of str or Boolean + name of datasets as list, supports ? and * wildcards. + True is equivalent to [*], False is equivalent to [] + what : str + attribute to change (must be in self.selection) - """ - self.__manage_visible(datasets,what,'add') + """ + self.__manage_visible(datasets,what,'set') - def del_visible(self,what,datasets): - """ - Delete from active groupse. - - Parameters - ---------- - datasets : list of str or Boolean - name of datasets as list, supports ? and * wildcards. - True is equivalent to [*], False is equivalent to [] - what : str - attribute to change (must be in self.visible) + def add_visible(self,what,datasets): + """ + Add to active groups. - """ - self.__manage_visible(datasets,what,'del') + Parameters + ---------- + datasets : list of str or Boolean + name of datasets as list, supports ? and * wildcards. + True is equivalent to [*], False is equivalent to [] + what : str + attribute to change (must be in self.selection) + """ + self.__manage_visible(datasets,what,'add') - def groups_with_datasets(self,datasets): - """ - Get groups that contain all requested datasets. - - Only groups within inc?????/constituent/*_*/* inc?????/materialpoint/*_*/* - are considered as they contain the data. - Single strings will be treated as list with one entry. - - Wild card matching is allowed, but the number of arguments need to fit. - - Parameters - ---------- - datasets : iterable or str or boolean - - 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,str) else datasets - groups = [] - - with h5py.File(self.fname,'r') as f: - for i in self.iter_visible('increments'): - for o,p in zip(['constituents','materialpoints'],['con_physics','mat_physics']): - for oo in self.iter_visible(o): - for pp in self.iter_visible(p): - group = '/'.join([i,o[:-1],oo,pp]) # o[:-1]: plural/singular issue - if sets is True: - groups.append(group) - else: - match = [e for e_ in [glob.fnmatch.filter(f[group].keys(),s) for s in sets] for e in e_] - if len(set(match)) == len(sets) : groups.append(group) - return groups + def del_visible(self,what,datasets): + """ + Delete from active groupe. + Parameters + ---------- + datasets : list of str or Boolean + name of datasets as list, supports ? and * wildcards. + True is equivalent to [*], False is equivalent to [] + what : str + attribute to change (must be in self.selection) - def list_data(self): - """Return information on all active datasets in the file.""" - message = '' - with h5py.File(self.fname,'r') as f: - for i in self.iter_visible('increments'): - message+='\n{} ({}s)\n'.format(i,self.times[self.increments.index(i)]) - for o,p in zip(['constituents','materialpoints'],['con_physics','mat_physics']): - for oo in self.iter_visible(o): - message+=' {}\n'.format(oo) - for pp in self.iter_visible(p): - message+=' {}\n'.format(pp) - group = '/'.join([i,o[:-1],oo,pp]) # o[:-1]: plural/singular issue - for d in f[group].keys(): - try: - dataset = f['/'.join([group,d])] - message+=' {} / ({}): {}\n'.format(d,dataset.attrs['Unit'].decode(),dataset.attrs['Description'].decode()) - except KeyError: - pass - return message + """ + self.__manage_visible(datasets,what,'del') - 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.iter_visible('increments'): - k = '/'.join([i,'geometry',label]) - try: - f[k] - path.append(k) - except KeyError as e: - pass - for o,p in zip(['constituents','materialpoints'],['con_physics','mat_physics']): - for oo in self.iter_visible(o): - for pp in self.iter_visible(p): - k = '/'.join([i,o[:-1],oo,pp,label]) - try: - f[k] - path.append(k) - except KeyError as e: - pass - return path - - - def get_constituent_ID(self,c=0): - """Pointwise constituent ID.""" - with h5py.File(self.fname,'r') as f: - names = f['/mapping/cellResults/constituent']['Name'][:,c].astype('str') - return np.array([int(n.split('_')[0]) for n in names.tolist()],dtype=np.int32) + def groups_with_datasets(self,datasets): + """ + Get groups that contain all requested datasets. + Only groups within inc?????/constituent/*_*/* inc?????/materialpoint/*_*/* + are considered as they contain the data. + Single strings will be treated as list with one entry. - def get_crystal_structure(self): # ToDo: extension to multi constituents/phase - """Info about the crystal structure.""" - with h5py.File(self.fname,'r') as f: - return f[self.get_dataset_location('orientation')[0]].attrs['Lattice'].astype('str') # np.bytes_ to string + Wild card matching is allowed, but the number of arguments need to fit. + Parameters + ---------- + datasets : iterable or str or boolean - 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. - """ - with h5py.File(self.fname,'r') as f: - shape = (self.Nmaterialpoints,) + np.shape(f[path[0]])[1:] - if len(shape) == 1: shape = shape +(1,) - dataset = np.full(shape,np.nan,dtype=np.dtype(f[path[0]])) - for pa in path: - label = pa.split('/')[2] - - if (pa.split('/')[1] == 'geometry'): - dataset = np.array(f[pa]) - continue - - p = np.where(f['mapping/cellResults/constituent'][:,c]['Name'] == str.encode(label))[0] - if len(p)>0: - u = (f['mapping/cellResults/constituent']['Position'][p,c]) - a = np.array(f[pa]) - if len(a.shape) == 1: - a=a.reshape([a.shape[0],1]) - dataset[p,:] = a[u,:] + 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'] - p = np.where(f['mapping/cellResults/materialpoint']['Name'] == str.encode(label))[0] - if len(p)>0: - u = (f['mapping/cellResults/materialpoint']['Position'][p.tolist()]) - a = np.array(f[pa]) - if len(a.shape) == 1: - a=a.reshape([a.shape[0],1]) - dataset[p,:] = a[u,:] - - if plain and dataset.dtype.names is not None: - return dataset.view(('float64',len(dataset.dtype.names))) - else: - return dataset + """ + if datasets is False: return [] + sets = [datasets] if isinstance(datasets,str) else datasets + groups = [] - def cell_coordinates(self): - """Return initial coordinates of the cell centers.""" - if self.structured: - delta = self.size/self.grid*0.5 - z, y, x = np.meshgrid(np.linspace(delta[2],self.size[2]-delta[2],self.grid[2]), - np.linspace(delta[1],self.size[1]-delta[1],self.grid[1]), - np.linspace(delta[0],self.size[0]-delta[0],self.grid[0]), - ) - return np.concatenate((x[:,:,:,None],y[:,:,:,None],z[:,:,:,None]),axis = 3).reshape([np.product(self.grid),3]) - else: - with h5py.File(self.fname,'r') as f: - return f['geometry/x_c'][()] - - - def add_absolute(self,x): - """ - Add absolute value. - - Parameters - ---------- - x : str - Label of the dataset containing a scalar, vector, or tensor. - - """ - def __add_absolute(x): - - return { - 'data': np.abs(x['data']), - 'label': '|{}|'.format(x['label']), - 'meta': { - 'Unit': x['meta']['Unit'], - 'Description': 'Absolute value of {} ({})'.format(x['label'],x['meta']['Description']), - 'Creator': 'dadf5.py:add_abs v{}'.format(version) - } - } - - requested = [{'label':x,'arg':'x'}] - - self.__add_generic_pointwise(__add_absolute,requested) - - - def add_calculation(self,formula,label,unit='n/a',description=None,vectorized=True): - """ - Add result of a general formula. - - Parameters - ---------- - formula : str - Formula, refer to datasets by ‘#Label#‘. - label : str - Label of the dataset containing the result of the calculation. - unit : str, optional - Physical unit of the result. - description : str, optional - Human readable description of the result. - vectorized : bool, optional - Indicate whether the formula is written in vectorized form. Default is ‘True’. - - """ - if vectorized is False: - raise NotImplementedError - - def __add_calculation(**kwargs): - - formula = kwargs['formula'] - for d in re.findall(r'#(.*?)#',formula): - formula = formula.replace('#{}#'.format(d),"kwargs['{}']['data']".format(d)) - - return { - 'data': eval(formula), - 'label': kwargs['label'], - 'meta': { - 'Unit': kwargs['unit'], - 'Description': '{} (formula: {})'.format(kwargs['description'],kwargs['formula']), - 'Creator': 'dadf5.py:add_calculation v{}'.format(version) - } - } - - requested = [{'label':d,'arg':d} for d in set(re.findall(r'#(.*?)#',formula))] # datasets used in the formula - pass_through = {'formula':formula,'label':label,'unit':unit,'description':description} - - self.__add_generic_pointwise(__add_calculation,requested,pass_through) - - - def add_Cauchy(self,P='P',F='F'): - """ - Add Cauchy stress calculated from 1. Piola-Kirchhoff stress and deformation gradient. - - Parameters - ---------- - P : str, optional - Label of the dataset containing the 1. Piola-Kirchhoff stress. Default value is ‘P’. - F : str, optional - Label of the dataset containing the deformation gradient. Default value is ‘F’. - - """ - def __add_Cauchy(F,P): - - return { - 'data': mechanics.Cauchy(F['data'],P['data']), - 'label': 'sigma', - 'meta': { - 'Unit': P['meta']['Unit'], - 'Description': 'Cauchy stress calculated from {} ({}) '.format(P['label'],P['meta']['Description'])+\ - 'and deformation gradient {} ({})'.format(F['label'],F['meta']['Description']), - 'Creator': 'dadf5.py:add_Cauchy v{}'.format(version) - } - } - - requested = [{'label':F,'arg':'F'}, - {'label':P,'arg':'P'} ] - - self.__add_generic_pointwise(__add_Cauchy,requested) - - - def add_determinant(self,x): - """ - Add the determinant of a tensor. - - Parameters - ---------- - x : str - Label of the dataset containing a tensor. - - """ - def __add_determinant(x): - - return { - 'data': np.linalg.det(x['data']), - 'label': 'det({})'.format(x['label']), - 'meta': { - 'Unit': x['meta']['Unit'], - 'Description': 'Determinant of tensor {} ({})'.format(x['label'],x['meta']['Description']), - 'Creator': 'dadf5.py:add_determinant v{}'.format(version) - } - } - - requested = [{'label':x,'arg':'x'}] - - self.__add_generic_pointwise(__add_determinant,requested) - - - def add_deviator(self,x): - """ - Add the deviatoric part of a tensor. - - Parameters - ---------- - x : str - Label of the dataset containing a tensor. - - """ - def __add_deviator(x): - - if not np.all(np.array(x['data'].shape[1:]) == np.array([3,3])): - raise ValueError - - return { - 'data': mechanics.deviatoric_part(x['data']), - 'label': 's_{}'.format(x['label']), - 'meta': { - 'Unit': x['meta']['Unit'], - 'Description': 'Deviator of tensor {} ({})'.format(x['label'],x['meta']['Description']), - 'Creator': 'dadf5.py:add_deviator v{}'.format(version) - } - } - - requested = [{'label':x,'arg':'x'}] - - self.__add_generic_pointwise(__add_deviator,requested) - - - def add_maximum_shear(self,x): - """ - Add maximum shear components of symmetric tensor. - - Parameters - ---------- - x : str - Label of the dataset containing a symmetric tensor. - - """ - def __add_maximum_shear(x): - - return { - 'data': mechanics.maximum_shear(x['data']), - 'label': 'max_shear({})'.format(x['label']), - 'meta': { - 'Unit': x['meta']['Unit'], - 'Description': 'Maximum shear component of of {} ({})'.format(x['label'],x['meta']['Description']), - 'Creator': 'dadf5.py:add_maximum_shear v{}'.format(version) - } - } - - requested = [{'label':x,'arg':'x'}] - - self.__add_generic_pointwise(__add_maximum_shear,requested) - - - def add_Mises(self,x): - """ - Add the equivalent Mises stress or strain of a symmetric tensor. - - Parameters - ---------- - x : str - Label of the dataset containing a symmetric stress or strain tensor. - - """ - def __add_Mises(x): - - t = 'strain' if x['meta']['Unit'] == '1' else \ - 'stress' - return { - 'data': mechanics.Mises_strain(x['data']) if t=='strain' else mechanics.Mises_stress(x['data']), - 'label': '{}_vM'.format(x['label']), - 'meta': { - 'Unit': x['meta']['Unit'], - 'Description': 'Mises equivalent {} of {} ({})'.format(t,x['label'],x['meta']['Description']), - 'Creator': 'dadf5.py:add_Mises v{}'.format(version) - } - } - - requested = [{'label':x,'arg':'x'}] - - self.__add_generic_pointwise(__add_Mises,requested) - - - def add_norm(self,x,ord=None): - """ - Add the norm of vector or tensor. - - Parameters - ---------- - x : str - Label of the dataset containing a vector or tensor. - 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. - - """ - def __add_norm(x,ord): - - o = ord - if len(x['data'].shape) == 2: - axis = 1 - t = 'vector' - if o is None: o = 2 - elif len(x['data'].shape) == 3: - axis = (1,2) - t = 'tensor' - if o is None: o = 'fro' - else: - raise ValueError - - return { - 'data': np.linalg.norm(x['data'],ord=o,axis=axis,keepdims=True), - 'label': '|{}|_{}'.format(x['label'],o), - 'meta': { - 'Unit': x['meta']['Unit'], - 'Description': '{}-Norm of {} {} ({})'.format(ord,t,x['label'],x['meta']['Description']), - 'Creator': 'dadf5.py:add_norm v{}'.format(version) - } - } - - requested = [{'label':x,'arg':'x'}] - - self.__add_generic_pointwise(__add_norm,requested,{'ord':ord}) - - - def add_principal_components(self,x): - """ - Add principal components of symmetric tensor. - - The principal components are sorted in descending order, each repeated according to its multiplicity. - - Parameters - ---------- - x : str - Label of the dataset containing a symmetric tensor. - - """ - def __add_principal_components(x): - - return { - 'data': mechanics.principal_components(x['data']), - 'label': 'lambda_{}'.format(x['label']), - 'meta': { - 'Unit': x['meta']['Unit'], - 'Description': 'Pricipal components of {} ({})'.format(x['label'],x['meta']['Description']), - 'Creator': 'dadf5.py:add_principal_components v{}'.format(version) - } - } - - requested = [{'label':x,'arg':'x'}] - - self.__add_generic_pointwise(__add_principal_components,requested) - - - def add_spherical(self,x): - """ - Add the spherical (hydrostatic) part of a tensor. - - Parameters - ---------- - x : str - Label of the dataset containing a tensor. - - """ - def __add_spherical(x): - - if not np.all(np.array(x['data'].shape[1:]) == np.array([3,3])): - raise ValueError - - return { - 'data': mechanics.spherical_part(x['data']), - 'label': 'p_{}'.format(x['label']), - 'meta': { - 'Unit': x['meta']['Unit'], - 'Description': 'Spherical component of tensor {} ({})'.format(x['label'],x['meta']['Description']), - 'Creator': 'dadf5.py:add_spherical v{}'.format(version) - } - } - - requested = [{'label':x,'arg':'x'}] - - self.__add_generic_pointwise(__add_spherical,requested) - - - def add_strain_tensor(self,F='F',t='U',m=0): - """ - Add strain tensor calculated from a deformation gradient. - - For details refer to damask.mechanics.strain_tensor - - Parameters - ---------- - F : str, optional - Label of the dataset containing the deformation gradient. Default value is ‘F’. - t : {‘V’, ‘U’}, optional - Type of the polar decomposition, ‘V’ for right stretch tensor and ‘U’ for left stretch tensor. - Defaults value is ‘U’. - m : float, optional - Order of the strain calculation. Default value is ‘0.0’. - - """ - def __add_strain_tensor(F,t,m): - - return { - 'data': mechanics.strain_tensor(F['data'],t,m), - 'label': 'epsilon_{}^{}({})'.format(t,m,F['label']), - 'meta': { - 'Unit': F['meta']['Unit'], - 'Description': 'Strain tensor of {} ({})'.format(F['label'],F['meta']['Description']), - 'Creator': 'dadf5.py:add_strain_tensor v{}'.format(version) - } - } - - requested = [{'label':F,'arg':'F'}] - - self.__add_generic_pointwise(__add_strain_tensor,requested,{'t':t,'m':m}) - - - def __add_generic_pointwise(self,func,datasets_requested,extra_args={}): - """ - General function to add pointwise data. - - Parameters - ---------- - func : function - Function that calculates a new dataset from one or more datasets per HDF5 group. - datasets_requested : list of dictionaries - Details of the datasets to be used: label (in HDF5 file) and arg (argument to which the data is parsed in func). - extra_args : dictionary, optional - Any extra arguments parsed to func. - - """ - def job(args): - """Call function with input data + extra arguments, returns results + group.""" - args['results'].put({**args['func'](**args['in']),'group':args['group']}) - - - N_threads = 1 # ToDo: should be a parameter - - results = Queue(N_threads) - pool = util.ThreadPool(N_threads) - N_added = N_threads + 1 - - todo = [] - # ToDo: It would be more memory efficient to read only from file when required, i.e. do to it in pool.add_task - for group in self.groups_with_datasets([d['label'] for d in datasets_requested]): - with h5py.File(self.fname,'r') as f: - datasets_in = {} - for d in datasets_requested: - loc = f[group+'/'+d['label']] - data = loc[()] - meta = {k:loc.attrs[k].decode() for k in loc.attrs.keys()} - datasets_in[d['arg']] = {'data': data, 'meta' : meta, 'label' : d['label']} - - todo.append({'in':{**datasets_in,**extra_args},'func':func,'group':group,'results':results}) - - pool.map(job, todo[:N_added]) # initialize - - N_not_calculated = len(todo) - while N_not_calculated > 0: - result = results.get() - with h5py.File(self.fname,'a') as f: # write to file - dataset_out = f[result['group']].create_dataset(result['label'],data=result['data']) - for k in result['meta'].keys(): - dataset_out.attrs[k] = result['meta'][k].encode() - N_not_calculated-=1 - - if N_added < len(todo): # add more jobs - pool.add_task(job,todo[N_added]) - N_added +=1 - - pool.wait_completion() - - - def to_vtk(self,labels,mode='Cell'): - """ - Export to vtk cell/point data. - - Parameters - ---------- - labels : str or list of - Labels of the datasets to be exported. - mode : str, either 'Cell' or 'Point' - Export in cell format or point format. - Default value is 'Cell'. - - """ - if mode=='Cell': - - if self.structured: - - coordArray = [vtk.vtkDoubleArray(),vtk.vtkDoubleArray(),vtk.vtkDoubleArray()] - for dim in [0,1,2]: - for c in np.linspace(0,self.size[dim],1+self.grid[dim]): - coordArray[dim].InsertNextValue(c) - - vtk_geom = vtk.vtkRectilinearGrid() - vtk_geom.SetDimensions(*(self.grid+1)) - vtk_geom.SetXCoordinates(coordArray[0]) - vtk_geom.SetYCoordinates(coordArray[1]) - vtk_geom.SetZCoordinates(coordArray[2]) - - else: - - nodes = vtk.vtkPoints() with h5py.File(self.fname,'r') as f: - nodes.SetData(numpy_support.numpy_to_vtk(f['/geometry/x_n'][()],deep=True)) - - vtk_geom = vtk.vtkUnstructuredGrid() - vtk_geom.SetPoints(nodes) - vtk_geom.Allocate(f['/geometry/T_c'].shape[0]) + for i in self.iter_visible('increments'): + for o,p in zip(['constituents','materialpoints'],['con_physics','mat_physics']): + for oo in self.iter_visible(o): + for pp in self.iter_visible(p): + group = '/'.join([i,o[:-1],oo,pp]) # o[:-1]: plural/singular issue + if sets is True: + groups.append(group) + else: + match = [e for e_ in [glob.fnmatch.filter(f[group].keys(),s) for s in sets] for e in e_] + if len(set(match)) == len(sets) : groups.append(group) + return groups - if self.version_major == 0 and self.version_minor <= 5: - vtk_type = vtk.VTK_HEXAHEDRON - n_nodes = 8 - else: - if f['/geometry/T_c'].attrs['VTK_TYPE'] == b'TRIANGLE': - vtk_type = vtk.VTK_TRIANGLE - n_nodes = 3 - elif f['/geometry/T_c'].attrs['VTK_TYPE'] == b'QUAD': - vtk_type = vtk.VTK_QUAD - n_nodes = 4 - elif f['/geometry/T_c'].attrs['VTK_TYPE'] == b'TETRA': # not tested - vtk_type = vtk.VTK_TETRA - n_nodes = 4 - elif f['/geometry/T_c'].attrs['VTK_TYPE'] == b'HEXAHEDRON': - vtk_type = vtk.VTK_HEXAHEDRON - n_nodes = 8 - for i in f['/geometry/T_c']: - vtk_geom.InsertNextCell(vtk_type,n_nodes,i-1) + def list_data(self): + """Return information on all active datasets in the file.""" + message = '' + with h5py.File(self.fname,'r') as f: + for i in self.iter_visible('increments'): + message+='\n{} ({}s)\n'.format(i,self.times[self.increments.index(i)]) + for o,p in zip(['constituents','materialpoints'],['con_physics','mat_physics']): + for oo in self.iter_visible(o): + message+=' {}\n'.format(oo) + for pp in self.iter_visible(p): + message+=' {}\n'.format(pp) + group = '/'.join([i,o[:-1],oo,pp]) # o[:-1]: plural/singular issue + for d in f[group].keys(): + try: + dataset = f['/'.join([group,d])] + message+=' {} / ({}): {}\n'.\ + format(d,dataset.attrs['Unit'].decode(),dataset.attrs['Description'].decode()) + except KeyError: + pass + return message - elif mode == 'Point': - Points = vtk.vtkPoints() - Vertices = vtk.vtkCellArray() - for c in self.cell_coordinates(): - pointID = Points.InsertNextPoint(c) - Vertices.InsertNextCell(1) - Vertices.InsertCellPoint(pointID) - - vtk_geom = vtk.vtkPolyData() - vtk_geom.SetPoints(Points) - vtk_geom.SetVerts(Vertices) - vtk_geom.Modified() - - N_digits = int(np.floor(np.log10(int(self.increments[-1][3:]))))+1 - - for i,inc in enumerate(self.iter_visible('increments')): - vtk_data = [] - - materialpoints_backup = self.visible['materialpoints'].copy() - self.set_visible('materialpoints',False) - for label in (labels if isinstance(labels,list) else [labels]): - for p in self.iter_visible('con_physics'): - if p != 'generic': - for c in self.iter_visible('constituents'): - x = self.get_dataset_location(label) - if len(x) == 0: - continue - array = self.read_dataset(x,0) - shape = [array.shape[0],np.product(array.shape[1:])] - vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape), - deep=True,array_type= vtk.VTK_DOUBLE)) - vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) #ToDo: hard coded 1! - vtk_geom.GetCellData().AddArray(vtk_data[-1]) - else: - x = self.get_dataset_location(label) - if len(x) == 0: + 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.iter_visible('increments'): + k = '/'.join([i,'geometry',label]) + try: + f[k] + path.append(k) + except KeyError as e: + pass + for o,p in zip(['constituents','materialpoints'],['con_physics','mat_physics']): + for oo in self.iter_visible(o): + for pp in self.iter_visible(p): + k = '/'.join([i,o[:-1],oo,pp,label]) + try: + f[k] + path.append(k) + except KeyError as e: + pass + return path + + + def get_constituent_ID(self,c=0): + """Pointwise constituent ID.""" + with h5py.File(self.fname,'r') as f: + names = f['/mapping/cellResults/constituent']['Name'][:,c].astype('str') + return np.array([int(n.split('_')[0]) for n in names.tolist()],dtype=np.int32) + + + def get_crystal_structure(self): # ToDo: extension to multi constituents/phase + """Info about the crystal structure.""" + with h5py.File(self.fname,'r') as f: + return f[self.get_dataset_location('orientation')[0]].attrs['Lattice'].astype('str') # np.bytes_ to string + + + 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. + """ + with h5py.File(self.fname,'r') as f: + shape = (self.Nmaterialpoints,) + np.shape(f[path[0]])[1:] + if len(shape) == 1: shape = shape +(1,) + dataset = np.full(shape,np.nan,dtype=np.dtype(f[path[0]])) + for pa in path: + label = pa.split('/')[2] + + if (pa.split('/')[1] == 'geometry'): + dataset = np.array(f[pa]) continue - array = self.read_dataset(x,0) - shape = [array.shape[0],np.product(array.shape[1:])] - vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape), - deep=True,array_type= vtk.VTK_DOUBLE)) - ph_name = re.compile(r'(?<=(constituent\/))(.*?)(?=(generic))') # identify phase name - dset_name = '1_' + re.sub(ph_name,r'',x[0].split('/',1)[1]) # removing phase name - vtk_data[-1].SetName(dset_name) - vtk_geom.GetCellData().AddArray(vtk_data[-1]) - self.set_visible('materialpoints',materialpoints_backup) + p = np.where(f['mapping/cellResults/constituent'][:,c]['Name'] == str.encode(label))[0] + if len(p)>0: + u = (f['mapping/cellResults/constituent']['Position'][p,c]) + a = np.array(f[pa]) + if len(a.shape) == 1: + a=a.reshape([a.shape[0],1]) + dataset[p,:] = a[u,:] - constituents_backup = self.visible['constituents'].copy() - self.set_visible('constituents',False) - for label in (labels if isinstance(labels,list) else [labels]): - for p in self.iter_visible('mat_physics'): - if p != 'generic': - for m in self.iter_visible('materialpoints'): - x = self.get_dataset_location(label) - if len(x) == 0: - continue - array = self.read_dataset(x,0) - shape = [array.shape[0],np.product(array.shape[1:])] - vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape), - deep=True,array_type= vtk.VTK_DOUBLE)) - vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) #ToDo: why 1_? - vtk_geom.GetCellData().AddArray(vtk_data[-1]) + p = np.where(f['mapping/cellResults/materialpoint']['Name'] == str.encode(label))[0] + if len(p)>0: + u = (f['mapping/cellResults/materialpoint']['Position'][p.tolist()]) + a = np.array(f[pa]) + if len(a.shape) == 1: + a=a.reshape([a.shape[0],1]) + dataset[p,:] = a[u,:] + + if plain and dataset.dtype.names is not None: + return dataset.view(('float64',len(dataset.dtype.names))) + else: + return dataset + + + def cell_coordinates(self): + """Return initial coordinates of the cell centers.""" + if self.structured: + return grid_filters.cell_coord0(self.grid,self.size,self.origin) + else: + with h5py.File(self.fname,'r') as f: + return f['geometry/x_c'][()] + + + def add_absolute(self,x): + """ + Add absolute value. + + Parameters + ---------- + x : str + Label of scalar, vector, or tensor dataset to take absolute value of. + + """ + def _add_absolute(x): + + return { + 'data': np.abs(x['data']), + 'label': '|{}|'.format(x['label']), + 'meta': { + 'Unit': x['meta']['Unit'], + 'Description': 'Absolute value of {} ({})'.format(x['label'],x['meta']['Description']), + 'Creator': 'dadf5.py:add_abs v{}'.format(version) + } + } + + self.__add_generic_pointwise(_add_absolute,{'x':x}) + + + def add_calculation(self,label,formula,unit='n/a',description=None,vectorized=True): + """ + Add result of a general formula. + + Parameters + ---------- + label : str + Label of resulting dataset. + formula : str + Formula to calculate resulting dataset. Existing datasets are referenced by ‘#TheirLabel#‘. + unit : str, optional + Physical unit of the result. + description : str, optional + Human-readable description of the result. + vectorized : bool, optional + Indicate whether the formula can be used in vectorized form. Defaults to ‘True’. + + """ + if not vectorized: + raise NotImplementedError + + def _add_calculation(**kwargs): + + formula = kwargs['formula'] + for d in re.findall(r'#(.*?)#',formula): + formula = formula.replace('#{}#'.format(d),"kwargs['{}']['data']".format(d)) + + return { + 'data': eval(formula), + 'label': kwargs['label'], + 'meta': { + 'Unit': kwargs['unit'], + 'Description': '{} (formula: {})'.format(kwargs['description'],kwargs['formula']), + 'Creator': 'dadf5.py:add_calculation v{}'.format(version) + } + } + + dataset_mapping = {d:d for d in set(re.findall(r'#(.*?)#',formula))} # datasets used in the formula + args = {'formula':formula,'label':label,'unit':unit,'description':description} + + self.__add_generic_pointwise(_add_calculation,dataset_mapping,args) + + + def add_Cauchy(self,P='P',F='F'): + """ + Add Cauchy stress calculated from first Piola-Kirchhoff stress and deformation gradient. + + Parameters + ---------- + P : str, optional + 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’. + + """ + def _add_Cauchy(P,F): + + return { + 'data': mechanics.Cauchy(P['data'],F['data']), + 'label': 'sigma', + 'meta': { + 'Unit': P['meta']['Unit'], + 'Description': 'Cauchy stress calculated from {} ({}) '.format(P['label'], + P['meta']['Description'])+\ + 'and {} ({})'.format(F['label'],F['meta']['Description']), + 'Creator': 'dadf5.py:add_Cauchy v{}'.format(version) + } + } + + self.__add_generic_pointwise(_add_Cauchy,{'P':P,'F':F}) + + + def add_determinant(self,T): + """ + Add the determinant of a tensor. + + Parameters + ---------- + T : str + Label of tensor dataset. + + """ + def _add_determinant(T): + + return { + 'data': np.linalg.det(T['data']), + 'label': 'det({})'.format(T['label']), + 'meta': { + 'Unit': T['meta']['Unit'], + 'Description': 'Determinant of tensor {} ({})'.format(T['label'],T['meta']['Description']), + 'Creator': 'dadf5.py:add_determinant v{}'.format(version) + } + } + + self.__add_generic_pointwise(_add_determinant,{'T':T}) + + + def add_deviator(self,T): + """ + Add the deviatoric part of a tensor. + + Parameters + ---------- + T : str + Label of tensor dataset. + + """ + def _add_deviator(T): + + if not np.all(np.array(T['data'].shape[1:]) == np.array([3,3])): + raise ValueError + + return { + 'data': mechanics.deviatoric_part(T['data']), + 'label': 's_{}'.format(T['label']), + 'meta': { + 'Unit': T['meta']['Unit'], + 'Description': 'Deviator of tensor {} ({})'.format(T['label'],T['meta']['Description']), + 'Creator': 'dadf5.py:add_deviator v{}'.format(version) + } + } + + self.__add_generic_pointwise(_add_deviator,{'T':T}) + + + def add_eigenvalues(self,S): + """ + Add eigenvalues of symmetric tensor. + + Parameters + ---------- + S : str + Label of symmetric tensor dataset. + + """ + def _add_eigenvalue(S): + + return { + 'data': mechanics.eigenvalues(S['data']), + 'label': 'lambda({})'.format(S['label']), + 'meta' : { + 'Unit': S['meta']['Unit'], + 'Description': 'Eigenvalues of {} ({})'.format(S['label'],S['meta']['Description']), + 'Creator': 'dadf5.py:add_eigenvalues v{}'.format(version) + } + } + + self.__add_generic_pointwise(_add_eigenvalue,{'S':S}) + + + def add_eigenvectors(self,S): + """ + Add eigenvectors of symmetric tensor. + + Parameters + ---------- + S : str + Label of symmetric tensor dataset. + + """ + def _add_eigenvector(S): + + return { + 'data': mechanics.eigenvectors(S['data']), + 'label': 'v({})'.format(S['label']), + 'meta' : { + 'Unit': '1', + 'Description': 'Eigenvectors of {} ({})'.format(S['label'],S['meta']['Description']), + 'Creator': 'dadf5.py:add_eigenvectors v{}'.format(version) + } + } + + self.__add_generic_pointwise(_add_eigenvector,{'S':S}) + + + def add_IPFcolor(self,q,l): + """ + Add RGB color tuple of inverse pole figure (IPF) color. + + Parameters + ---------- + q : str + Label of the dataset containing the crystallographic orientation as quaternions. + l : numpy.array of shape (3) + Lab frame direction for inverse pole figure. + + """ + def _add_IPFcolor(q,l): + + d = np.array(l) + d_unit = d/np.linalg.norm(d) + m = util.scale_to_coprime(d) + colors = np.empty((len(q['data']),3),np.uint8) + + lattice = q['meta']['Lattice'] + + for i,q in enumerate(q['data']): + o = Orientation(np.array([q['w'],q['x'],q['y'],q['z']]),lattice).reduced() + colors[i] = np.uint8(o.IPFcolor(d_unit)*255) + + return { + 'data': colors, + 'label': 'IPFcolor_[{} {} {}]'.format(*m), + 'meta' : { + 'Unit': 'RGB (8bit)', + 'Lattice': lattice, + 'Description': 'Inverse Pole Figure (IPF) colors for direction/plane [{} {} {})'.format(*m), + 'Creator': 'dadf5.py:add_IPFcolor v{}'.format(version) + } + } + + self.__add_generic_pointwise(_add_IPFcolor,{'q':q},{'l':l}) + + + def add_maximum_shear(self,S): + """ + Add maximum shear components of symmetric tensor. + + Parameters + ---------- + S : str + Label of symmetric tensor dataset. + + """ + def _add_maximum_shear(S): + + return { + 'data': mechanics.maximum_shear(S['data']), + 'label': 'max_shear({})'.format(S['label']), + 'meta': { + 'Unit': S['meta']['Unit'], + 'Description': 'Maximum shear component of {} ({})'.format(S['label'],S['meta']['Description']), + 'Creator': 'dadf5.py:add_maximum_shear v{}'.format(version) + } + } + + self.__add_generic_pointwise(_add_maximum_shear,{'S':S}) + + + def add_Mises(self,S): + """ + Add the equivalent Mises stress or strain of a symmetric tensor. + + Parameters + ---------- + S : str + Label of symmetric tensorial stress or strain dataset. + + """ + def _add_Mises(S): + + t = 'strain' if S['meta']['Unit'] == '1' else \ + 'stress' + return { + 'data': mechanics.Mises_strain(S['data']) if t=='strain' else mechanics.Mises_stress(S['data']), + 'label': '{}_vM'.format(S['label']), + 'meta': { + 'Unit': S['meta']['Unit'], + 'Description': 'Mises equivalent {} of {} ({})'.format(t,S['label'],S['meta']['Description']), + 'Creator': 'dadf5.py:add_Mises v{}'.format(version) + } + } + + self.__add_generic_pointwise(_add_Mises,{'S':S}) + + + def add_norm(self,x,ord=None): + """ + Add the norm of vector or tensor. + + Parameters + ---------- + x : str + Label of vector or tensor dataset. + 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. + + """ + def _add_norm(x,ord): + + o = ord + if len(x['data'].shape) == 2: + axis = 1 + t = 'vector' + if o is None: o = 2 + elif len(x['data'].shape) == 3: + axis = (1,2) + t = 'tensor' + if o is None: o = 'fro' else: - x = self.get_dataset_location(label) - if len(x) == 0: - continue - array = self.read_dataset(x,0) - shape = [array.shape[0],np.product(array.shape[1:])] - vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape), - deep=True,array_type= vtk.VTK_DOUBLE)) - vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) - vtk_geom.GetCellData().AddArray(vtk_data[-1]) - self.set_visible('constituents',constituents_backup) - - if mode=='Cell': - writer = vtk.vtkXMLRectilinearGridWriter() if self.structured else \ - vtk.vtkXMLUnstructuredGridWriter() - x = self.get_dataset_location('u_n') - vtk_data.append(numpy_support.numpy_to_vtk(num_array=self.read_dataset(x,0), - deep=True,array_type=vtk.VTK_DOUBLE)) - vtk_data[-1].SetName('u') - vtk_geom.GetPointData().AddArray(vtk_data[-1]) - elif mode == 'Point': - writer = vtk.vtkXMLPolyDataWriter() + raise ValueError - - file_out = '{}_inc{}.{}'.format(os.path.splitext(os.path.basename(self.fname))[0], - inc[3:].zfill(N_digits), - writer.GetDefaultFileExtension()) - - writer.SetCompressorTypeToZLib() - writer.SetDataModeToBinary() - writer.SetFileName(file_out) - writer.SetInputData(vtk_geom) + return { + 'data': np.linalg.norm(x['data'],ord=o,axis=axis,keepdims=True), + 'label': '|{}|_{}'.format(x['label'],o), + 'meta': { + 'Unit': x['meta']['Unit'], + 'Description': '{}-norm of {} {} ({})'.format(ord,t,x['label'],x['meta']['Description']), + 'Creator': 'dadf5.py:add_norm v{}'.format(version) + } + } - writer.Write() + self.__add_generic_pointwise(_add_norm,{'x':x},{'ord':ord}) + + + def add_PK2(self,P='P',F='F'): + """ + Add 2. Piola-Kirchhoff calculated from first Piola-Kirchhoff stress and deformation gradient. + + Parameters + ---------- + P : str, optional + Label first Piola-Kirchhoff stress dataset. Defaults to ‘P’. + F : str, optional + Label of deformation gradient dataset. Defaults to ‘F’. + + """ + def _add_PK2(P,F): + + return { + 'data': mechanics.PK2(P['data'],F['data']), + 'label': 'S', + 'meta': { + 'Unit': P['meta']['Unit'], + 'Description': '2. Kirchhoff stress calculated from {} ({}) '.format(P['label'], + P['meta']['Description'])+\ + 'and {} ({})'.format(F['label'],F['meta']['Description']), + 'Creator': 'dadf5.py:add_PK2 v{}'.format(version) + } + } + + self.__add_generic_pointwise(_add_PK2,{'P':P,'F':F}) + + + def add_pole(self,q,p,polar=False): + """ + Add coordinates of stereographic projection of given pole in crystal frame. + + Parameters + ---------- + q : str + Label of the dataset containing the crystallographic orientation as quaternions. + p : numpy.array of shape (3) + Crystallographic direction or plane. + polar : bool, optional + Give pole in polar coordinates. Defaults to False. + + """ + def _add_pole(q,p,polar): + + pole = np.array(p) + unit_pole = pole/np.linalg.norm(pole) + m = util.scale_to_coprime(pole) + coords = np.empty((len(q['data']),2)) + + for i,q in enumerate(q['data']): + o = Rotation(np.array([q['w'],q['x'],q['y'],q['z']])) + rotatedPole = o*unit_pole # rotate pole according to crystal orientation + (x,y) = rotatedPole[0:2]/(1.+abs(unit_pole[2])) # stereographic projection + coords[i] = [np.sqrt(x*x+y*y),np.arctan2(y,x)] if polar else [x,y] + + return { + 'data': coords, + 'label': 'p^{}_[{} {} {})'.format(u'rφ' if polar else 'xy',*m), + 'meta' : { + 'Unit': '1', + 'Description': '{} coordinates of stereographic projection of pole (direction/plane) in crystal frame'\ + .format('Polar' if polar else 'Cartesian'), + 'Creator' : 'dadf5.py:add_pole v{}'.format(version) + } + } + + self.__add_generic_pointwise(_add_pole,{'q':q},{'p':p,'polar':polar}) + + + def add_rotational_part(self,F): + """ + Add rotational part of a deformation gradient. + + Parameters + ---------- + F : str, optional + Label of deformation gradient dataset. + + """ + def _add_rotational_part(F): + + return { + 'data': mechanics.rotational_part(F['data']), + 'label': 'R({})'.format(F['label']), + 'meta': { + 'Unit': F['meta']['Unit'], + 'Description': 'Rotational part of {} ({})'.format(F['label'],F['meta']['Description']), + 'Creator': 'dadf5.py:add_rotational_part v{}'.format(version) + } + } + + self.__add_generic_pointwise(_add_rotational_part,{'F':F}) + + + def add_spherical(self,T): + """ + Add the spherical (hydrostatic) part of a tensor. + + Parameters + ---------- + T : str + Label of tensor dataset. + + """ + def _add_spherical(T): + + if not np.all(np.array(T['data'].shape[1:]) == np.array([3,3])): + raise ValueError + + return { + 'data': mechanics.spherical_part(T['data']), + 'label': 'p_{}'.format(T['label']), + 'meta': { + 'Unit': T['meta']['Unit'], + 'Description': 'Spherical component of tensor {} ({})'.format(T['label'],T['meta']['Description']), + 'Creator': 'dadf5.py:add_spherical v{}'.format(version) + } + } + + self.__add_generic_pointwise(_add_spherical,{'T':T}) + + + def add_strain_tensor(self,F='F',t='V',m=0.0): + """ + Add strain tensor of a deformation gradient. + + For details refer to damask.mechanics.strain_tensor + + 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’. + m : float, optional + Order of the strain calculation. Defaults to ‘0.0’. + + """ + def _add_strain_tensor(F,t,m): + + return { + 'data': mechanics.strain_tensor(F['data'],t,m), + 'label': 'epsilon_{}^{}({})'.format(t,m,F['label']), + 'meta': { + 'Unit': F['meta']['Unit'], + 'Description': 'Strain tensor of {} ({})'.format(F['label'],F['meta']['Description']), + 'Creator': 'dadf5.py:add_strain_tensor v{}'.format(version) + } + } + + self.__add_generic_pointwise(_add_strain_tensor,{'F':F},{'t':t,'m':m}) + + + def add_stretch_tensor(self,F='F',t='V'): + """ + Add stretch tensor of a deformation gradient. + + 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’. + + """ + def _add_stretch_tensor(F,t): + + return { + 'data': mechanics.left_stretch(F['data']) if t == 'V' else mechanics.right_stretch(F['data']), + 'label': '{}({})'.format(t,F['label']), + 'meta': { + 'Unit': F['meta']['Unit'], + 'Description': '{} stretch tensor of {} ({})'.format('Left' if t == 'V' else 'Right', + F['label'],F['meta']['Description']), + 'Creator': 'dadf5.py:add_stretch_tensor v{}'.format(version) + } + } + + self.__add_generic_pointwise(_add_stretch_tensor,{'F':F},{'t':t}) + + + def __add_generic_pointwise(self,func,dataset_mapping,args={}): + """ + General function to add pointwise data. + + Parameters + ---------- + func : function + Function that calculates a new dataset from one or more datasets per HDF5 group. + dataset_mapping : dictionary + Mapping HDF5 data label to callback function argument + extra_args : dictionary, optional + Any extra arguments parsed to func. + + """ + def job(args): + """Call function with input data + extra arguments, returns results + group.""" + args['results'].put({**args['func'](**args['in']),'group':args['group']}) + + env = Environment() + N_threads = int(env.options['DAMASK_NUM_THREADS']) + N_threads //=N_threads # disable for the moment + + results = Queue(N_threads) + pool = util.ThreadPool(N_threads) + N_added = N_threads + 1 + + todo = [] + # ToDo: It would be more memory efficient to read only from file when required, i.e. do to it in pool.add_task + for group in self.groups_with_datasets(dataset_mapping.values()): + with h5py.File(self.fname,'r') as f: + datasets_in = {} + for arg,label in dataset_mapping.items(): + loc = f[group+'/'+label] + data = loc[()] + meta = {k:loc.attrs[k].decode() for k in loc.attrs.keys()} + datasets_in[arg] = {'data': data, 'meta': meta, 'label': label} + + todo.append({'in':{**datasets_in,**args},'func':func,'group':group,'results':results}) + + pool.map(job, todo[:N_added]) # initialize + + N_not_calculated = len(todo) + while N_not_calculated > 0: + result = results.get() + with h5py.File(self.fname,'a') as f: # write to file + dataset_out = f[result['group']].create_dataset(result['label'],data=result['data']) + for k in result['meta'].keys(): + dataset_out.attrs[k] = result['meta'][k].encode() + N_not_calculated-=1 + + if N_added < len(todo): # add more jobs + pool.add_task(job,todo[N_added]) + N_added +=1 + + pool.wait_completion() + + + def to_vtk(self,labels,mode='cell'): + """ + Export to vtk cell/point data. + + Parameters + ---------- + labels : str or list of + Labels of the datasets to be exported. + mode : str, either 'cell' or 'point' + Export in cell format or point format. + Defaults to 'cell'. + + """ + if mode.lower()=='cell': + + if self.structured: + + coordArray = [vtk.vtkDoubleArray(),vtk.vtkDoubleArray(),vtk.vtkDoubleArray()] + for dim in [0,1,2]: + for c in np.linspace(0,self.size[dim],1+self.grid[dim]): + coordArray[dim].InsertNextValue(c) + + vtk_geom = vtk.vtkRectilinearGrid() + vtk_geom.SetDimensions(*(self.grid+1)) + vtk_geom.SetXCoordinates(coordArray[0]) + vtk_geom.SetYCoordinates(coordArray[1]) + vtk_geom.SetZCoordinates(coordArray[2]) + + else: + + nodes = vtk.vtkPoints() + with h5py.File(self.fname,'r') as f: + nodes.SetData(numpy_support.numpy_to_vtk(f['/geometry/x_n'][()],deep=True)) + + vtk_geom = vtk.vtkUnstructuredGrid() + vtk_geom.SetPoints(nodes) + vtk_geom.Allocate(f['/geometry/T_c'].shape[0]) + + if self.version_major == 0 and self.version_minor <= 5: + vtk_type = vtk.VTK_HEXAHEDRON + n_nodes = 8 + else: + if f['/geometry/T_c'].attrs['VTK_TYPE'] == b'TRIANGLE': + vtk_type = vtk.VTK_TRIANGLE + n_nodes = 3 + elif f['/geometry/T_c'].attrs['VTK_TYPE'] == b'QUAD': + vtk_type = vtk.VTK_QUAD + n_nodes = 4 + elif f['/geometry/T_c'].attrs['VTK_TYPE'] == b'TETRA': # not tested + vtk_type = vtk.VTK_TETRA + n_nodes = 4 + elif f['/geometry/T_c'].attrs['VTK_TYPE'] == b'HEXAHEDRON': + vtk_type = vtk.VTK_HEXAHEDRON + n_nodes = 8 + + for i in f['/geometry/T_c']: + vtk_geom.InsertNextCell(vtk_type,n_nodes,i-1) + + elif mode.lower()=='point': + Points = vtk.vtkPoints() + Vertices = vtk.vtkCellArray() + for c in self.cell_coordinates(): + pointID = Points.InsertNextPoint(c) + Vertices.InsertNextCell(1) + Vertices.InsertCellPoint(pointID) + + vtk_geom = vtk.vtkPolyData() + vtk_geom.SetPoints(Points) + vtk_geom.SetVerts(Vertices) + vtk_geom.Modified() + + N_digits = int(np.floor(np.log10(int(self.increments[-1][3:]))))+1 + + for i,inc in enumerate(self.iter_visible('increments')): + vtk_data = [] + + materialpoints_backup = self.selection['materialpoints'].copy() + self.set_visible('materialpoints',False) + for label in (labels if isinstance(labels,list) else [labels]): + for p in self.iter_visible('con_physics'): + if p != 'generic': + for c in self.iter_visible('constituents'): + x = self.get_dataset_location(label) + if len(x) == 0: + continue + array = self.read_dataset(x,0) + shape = [array.shape[0],np.product(array.shape[1:])] + vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape), + deep=True,array_type= vtk.VTK_DOUBLE)) + vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) #ToDo: hard coded 1! + vtk_geom.GetCellData().AddArray(vtk_data[-1]) + + else: + x = self.get_dataset_location(label) + if len(x) == 0: + continue + array = self.read_dataset(x,0) + shape = [array.shape[0],np.product(array.shape[1:])] + vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape), + deep=True,array_type= vtk.VTK_DOUBLE)) + ph_name = re.compile(r'(?<=(constituent\/))(.*?)(?=(generic))') # identify phase name + dset_name = '1_' + re.sub(ph_name,r'',x[0].split('/',1)[1]) # removing phase name + vtk_data[-1].SetName(dset_name) + vtk_geom.GetCellData().AddArray(vtk_data[-1]) + + self.set_visible('materialpoints',materialpoints_backup) + + constituents_backup = self.selection['constituents'].copy() + self.set_visible('constituents',False) + for label in (labels if isinstance(labels,list) else [labels]): + for p in self.iter_visible('mat_physics'): + if p != 'generic': + for m in self.iter_visible('materialpoints'): + x = self.get_dataset_location(label) + if len(x) == 0: + continue + array = self.read_dataset(x,0) + shape = [array.shape[0],np.product(array.shape[1:])] + vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape), + deep=True,array_type= vtk.VTK_DOUBLE)) + vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) #ToDo: why 1_? + vtk_geom.GetCellData().AddArray(vtk_data[-1]) + else: + x = self.get_dataset_location(label) + if len(x) == 0: + continue + array = self.read_dataset(x,0) + shape = [array.shape[0],np.product(array.shape[1:])] + vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape), + deep=True,array_type= vtk.VTK_DOUBLE)) + vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) + vtk_geom.GetCellData().AddArray(vtk_data[-1]) + self.set_visible('constituents',constituents_backup) + + if mode.lower()=='cell': + writer = vtk.vtkXMLRectilinearGridWriter() if self.structured else \ + vtk.vtkXMLUnstructuredGridWriter() + x = self.get_dataset_location('u_n') + vtk_data.append(numpy_support.numpy_to_vtk(num_array=self.read_dataset(x,0), + deep=True,array_type=vtk.VTK_DOUBLE)) + vtk_data[-1].SetName('u') + vtk_geom.GetPointData().AddArray(vtk_data[-1]) + elif mode.lower()=='point': + writer = vtk.vtkXMLPolyDataWriter() + + + file_out = '{}_inc{}.{}'.format(os.path.splitext(os.path.basename(self.fname))[0], + inc[3:].zfill(N_digits), + writer.GetDefaultFileExtension()) + + writer.SetCompressorTypeToZLib() + writer.SetDataModeToBinary() + writer.SetFileName(file_out) + writer.SetInputData(vtk_geom) + + writer.Write() diff --git a/python/damask/environment.py b/python/damask/environment.py index 2fde2c329..4c9ab762c 100644 --- a/python/damask/environment.py +++ b/python/damask/environment.py @@ -8,7 +8,7 @@ class Environment(): def __init__(self): """Read and provide values of DAMASK configuration.""" self.options = {} - self.get_options() + self.__get_options() def relPath(self,relative = '.'): return os.path.join(self.rootDir(),relative) @@ -16,7 +16,7 @@ class Environment(): def rootDir(self): return os.path.normpath(os.path.join(os.path.realpath(__file__),'../../../')) - def get_options(self): + def __get_options(self): for item in ['DAMASK_NUM_THREADS', 'MSC_ROOT', 'MARC_VERSION', diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index 307f1d83d..9ffd76536 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -1,11 +1,11 @@ import numpy as np -def Cauchy(F,P): +def Cauchy(P,F): """ - Return Cauchy stress calculated from 1. Piola-Kirchhoff stress and deformation gradient. - + Return Cauchy stress calculated from first Piola-Kirchhoff stress and deformation gradient. + Resulting tensor is symmetrized as the Cauchy stress needs to be symmetric. - + Parameters ---------- F : numpy.array of shape (:,3,3) or (3,3) @@ -21,67 +21,10 @@ def Cauchy(F,P): return symmetric(sigma) -def PK2(F,P): - """ - Return 2. Piola-Kirchhoff stress calculated from 1. Piola-Kirchhoff stress and deformation gradient. - - Parameters - ---------- - F : numpy.array of shape (:,3,3) or (3,3) - Deformation gradient. - P : numpy.array of shape (:,3,3) or (3,3) - 1. Piola-Kirchhoff stress. - - """ - if np.shape(F) == np.shape(P) == (3,3): - S = np.dot(np.linalg.inv(F),P) - else: - S = np.einsum('ijk,ikl->ijl',np.linalg.inv(F),P) - return S - - -def strain_tensor(F,t,m): - """ - Return strain tensor calculated from deformation gradient. - - For details refer to https://en.wikipedia.org/wiki/Finite_strain_theory and - https://de.wikipedia.org/wiki/Verzerrungstensor - - Parameters - ---------- - F : numpy.array of shape (:,3,3) or (3,3) - Deformation gradient. - t : {‘V’, ‘U’} - Type of the polar decomposition, ‘V’ for left stretch tensor and ‘U’ for right stretch tensor. - m : float - Order of the strain. - - """ - F_ = F.reshape((1,3,3)) if F.shape == (3,3) else F - if t == 'V': - B = np.matmul(F_,transpose(F_)) - w,n = np.linalg.eigh(B) - elif t == 'U': - C = np.matmul(transpose(F_),F_) - w,n = np.linalg.eigh(C) - - if m > 0.0: - eps = 1.0/(2.0*abs(m)) * (+ np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n)) - - np.broadcast_to(np.eye(3),[F_.shape[0],3,3])) - elif m < 0.0: - eps = 1.0/(2.0*abs(m)) * (- np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n)) - + np.broadcast_to(np.eye(3),[F_.shape[0],3,3])) - else: - eps = np.matmul(n,np.einsum('ij,ikj->ijk',0.5*np.log(w),n)) - - return eps.reshape((3,3)) if np.shape(F) == (3,3) else \ - eps - - def deviatoric_part(x): """ Return deviatoric part of a tensor. - + Parameters ---------- x : numpy.array of shape (:,3,3) or (3,3) @@ -89,13 +32,151 @@ def deviatoric_part(x): """ return x - np.eye(3)*spherical_part(x) if np.shape(x) == (3,3) else \ - x - np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),[x.shape[0],3,3]),spherical_part(x)) + x - np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),[x.shape[0],3,3]),spherical_part(x)) + + +def eigenvalues(x): + """ + Return the eigenvalues, i.e. principal components, of a symmetric tensor. + + The eigenvalues are sorted in ascending order, each repeated according to + its multiplicity. + + Parameters + ---------- + x : numpy.array of shape (:,3,3) or (3,3) + Symmetric tensor of which the eigenvalues are computed. + + """ + return np.linalg.eigvalsh(symmetric(x)) + + +def eigenvectors(x,RHS=False): + """ + Return eigenvectors of a symmetric tensor. + + The eigenvalues are sorted in ascending order of their associated eigenvalues. + + Parameters + ---------- + x : numpy.array of shape (:,3,3) or (3,3) + Symmetric tensor of which the eigenvectors are computed. + RHS: bool, optional + Enforce right-handed coordinate system. Default is False. + + """ + (u,v) = np.linalg.eigh(symmetric(x)) + + if RHS: + if np.shape(x) == (3,3): + if np.linalg.det(v) < 0.0: v[:,2] *= -1.0 + else: + v[np.linalg.det(v) < 0.0,:,2] *= -1.0 + return v + + +def left_stretch(x): + """ + Return the left stretch of a tensor. + + Parameters + ---------- + x : numpy.array of shape (:,3,3) or (3,3) + Tensor of which the left stretch is computed. + + """ + return __polar_decomposition(x,'V')[0] + + +def maximum_shear(x): + """ + Return the maximum shear component of a symmetric tensor. + + Parameters + ---------- + x : numpy.array of shape (:,3,3) or (3,3) + Symmetric tensor of which the maximum shear is computed. + + """ + w = eigenvalues(x) + return (w[0] - w[2])*0.5 if np.shape(x) == (3,3) else \ + (w[:,0] - w[:,2])*0.5 + + +def Mises_strain(epsilon): + """ + Return the Mises equivalent of a strain tensor. + + Parameters + ---------- + epsilon : numpy.array of shape (:,3,3) or (3,3) + Symmetric strain tensor of which the von Mises equivalent is computed. + + """ + return __Mises(epsilon,2.0/3.0) + + +def Mises_stress(sigma): + """ + Return the Mises equivalent of a stress tensor. + + Parameters + ---------- + sigma : numpy.array of shape (:,3,3) or (3,3) + Symmetric stress tensor of which the von Mises equivalent is computed. + + """ + return __Mises(sigma,3.0/2.0) + + +def PK2(P,F): + """ + Calculate second Piola-Kirchhoff stress from first Piola-Kirchhoff stress and deformation gradient. + + Parameters + ---------- + P : numpy.array of shape (:,3,3) or (3,3) + 1. Piola-Kirchhoff stress. + F : numpy.array of shape (:,3,3) or (3,3) + Deformation gradient. + + """ + if np.shape(F) == np.shape(P) == (3,3): + S = np.dot(np.linalg.inv(F),P) + else: + S = np.einsum('ijk,ikl->ijl',np.linalg.inv(F),P) + return symmetric(S) + +def right_stretch(x): + """ + Return the right stretch of a tensor. + + Parameters + ---------- + x : numpy.array of shape (:,3,3) or (3,3) + Tensor of which the right stretch is computed. + + """ + return __polar_decomposition(x,'U')[0] + + +def rotational_part(x): + """ + Return the rotational part of a tensor. + + Parameters + ---------- + x : numpy.array of shape (:,3,3) or (3,3) + Tensor of which the rotational part is computed. + + """ + return __polar_decomposition(x,'R')[0] def spherical_part(x,tensor=False): """ Return spherical (hydrostatic) part of a tensor. - + Parameters ---------- x : numpy.array of shape (:,3,3) or (3,3) @@ -113,42 +194,50 @@ def spherical_part(x,tensor=False): return sph else: return np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),(x.shape[0],3,3)),sph) - - -def Mises_stress(sigma): + + +def strain_tensor(F,t,m): """ - Return the Mises equivalent of a stress tensor. - + Return strain tensor calculated from deformation gradient. + + For details refer to https://en.wikipedia.org/wiki/Finite_strain_theory and + https://de.wikipedia.org/wiki/Verzerrungstensor + Parameters ---------- - sigma : numpy.array of shape (:,3,3) or (3,3) - Symmetric stress tensor of which the von Mises equivalent is computed. + F : numpy.array of shape (:,3,3) or (3,3) + Deformation gradient. + t : {‘V’, ‘U’} + Type of the polar decomposition, ‘V’ for left stretch tensor and ‘U’ for right stretch tensor. + m : float + Order of the strain. """ - s = deviatoric_part(sigma) - return np.sqrt(3.0/2.0*(np.sum(s**2.0))) if np.shape(sigma) == (3,3) else \ - np.sqrt(3.0/2.0*np.einsum('ijk->i',s**2.0)) - - -def Mises_strain(epsilon): - """ - Return the Mises equivalent of a strain tensor. - - Parameters - ---------- - epsilon : numpy.array of shape (:,3,3) or (3,3) - Symmetric strain tensor of which the von Mises equivalent is computed. + F_ = F.reshape((1,3,3)) if F.shape == (3,3) else F + if t == 'V': + B = np.matmul(F_,transpose(F_)) + w,n = np.linalg.eigh(B) + elif t == 'U': + C = np.matmul(transpose(F_),F_) + w,n = np.linalg.eigh(C) - """ - s = deviatoric_part(epsilon) - return np.sqrt(2.0/3.0*(np.sum(s**2.0))) if np.shape(epsilon) == (3,3) else \ - np.sqrt(2.0/3.0*np.einsum('ijk->i',s**2.0)) + if m > 0.0: + eps = 1.0/(2.0*abs(m)) * (+ np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n)) + - np.broadcast_to(np.eye(3),[F_.shape[0],3,3])) + elif m < 0.0: + eps = 1.0/(2.0*abs(m)) * (- np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n)) + + np.broadcast_to(np.eye(3),[F_.shape[0],3,3])) + else: + eps = np.matmul(n,np.einsum('ij,ikj->ijk',0.5*np.log(w),n)) + + return eps.reshape((3,3)) if np.shape(F) == (3,3) else \ + eps def symmetric(x): """ Return the symmetrized tensor. - + Parameters ---------- x : numpy.array of shape (:,3,3) or (3,3) @@ -158,43 +247,10 @@ def symmetric(x): return (x+transpose(x))*0.5 -def maximum_shear(x): - """ - Return the maximum shear component of a symmetric tensor. - - Parameters - ---------- - x : numpy.array of shape (:,3,3) or (3,3) - Symmetric tensor of which the maximum shear is computed. - - """ - w = np.linalg.eigvalsh(symmetric(x)) # eigenvalues in ascending order - return (w[2] - w[0])*0.5 if np.shape(x) == (3,3) else \ - (w[:,2] - w[:,0])*0.5 - - -def principal_components(x): - """ - Return the principal components of a symmetric tensor. - - The principal components (eigenvalues) are sorted in descending order, each repeated according to - its multiplicity. - - Parameters - ---------- - x : numpy.array of shape (:,3,3) or (3,3) - Symmetric tensor of which the principal compontents are computed. - - """ - w = np.linalg.eigvalsh(symmetric(x)) # eigenvalues in ascending order - return w[::-1] if np.shape(x) == (3,3) else \ - w[:,::-1] - - def transpose(x): """ Return the transpose of a tensor. - + Parameters ---------- x : numpy.array of shape (:,3,3) or (3,3) @@ -205,62 +261,23 @@ def transpose(x): np.transpose(x,(0,2,1)) -def rotational_part(x): - """ - Return the rotational part of a tensor. - - Parameters - ---------- - x : numpy.array of shape (:,3,3) or (3,3) - Tensor of which the rotational part is computed. - - """ - return __polar_decomposition(x,'R')[0] - - -def left_stretch(x): - """ - Return the left stretch of a tensor. - - Parameters - ---------- - x : numpy.array of shape (:,3,3) or (3,3) - Tensor of which the left stretch is computed. - - """ - return __polar_decomposition(x,'V')[0] - - -def right_stretch(x): - """ - Return the right stretch of a tensor. - - Parameters - ---------- - x : numpy.array of shape (:,3,3) or (3,3) - Tensor of which the right stretch is computed. - - """ - return __polar_decomposition(x,'U')[0] - - def __polar_decomposition(x,requested): """ Singular value decomposition. - + Parameters ---------- x : numpy.array of shape (:,3,3) or (3,3) Tensor of which the singular values are computed. requested : iterable of str - Requested outputs: ‘R’ for the rotation tensor, + Requested outputs: ‘R’ for the rotation tensor, ‘V’ for left stretch tensor and ‘U’ for right stretch tensor. """ u, s, vh = np.linalg.svd(x) R = np.dot(u,vh) if np.shape(x) == (3,3) else \ np.einsum('ijk,ikl->ijl',u,vh) - + output = [] if 'R' in requested: output.append(R) @@ -268,5 +285,22 @@ def __polar_decomposition(x,requested): output.append(np.dot(x,R.T) if np.shape(x) == (3,3) else np.einsum('ijk,ilk->ijl',x,R)) if 'U' in requested: output.append(np.dot(R.T,x) if np.shape(x) == (3,3) else np.einsum('ikj,ikl->ijl',R,x)) - + return tuple(output) + + +def __Mises(x,s): + """ + Base equation for Mises equivalent of a stres or strain tensor. + + Parameters + ---------- + x : numpy.array of shape (:,3,3) or (3,3) + Symmetric tensor of which the von Mises equivalent is computed. + s : float + Scaling factor (2/3 for strain, 3/2 for stress). + + """ + d = deviatoric_part(x) + return np.sqrt(s*(np.sum(d**2.0))) if np.shape(x) == (3,3) else \ + np.sqrt(s*np.einsum('ijk->i',d**2.0)) diff --git a/python/damask/util.py b/python/damask/util.py index 0065daba5..1cb7b8602 100644 --- a/python/damask/util.py +++ b/python/damask/util.py @@ -3,14 +3,18 @@ import time import os import subprocess import shlex +from fractions import Fraction +from functools import reduce from optparse import Option from queue import Queue from threading import Thread +import numpy as np + class bcolors: """ ASCII Colors (Blender code). - + https://svn.blender.org/svnroot/bf-blender/trunk/blender/build_files/scons/tools/bcolors.py http://stackoverflow.com/questions/287871/print-in-terminal-with-colors-using-python """ @@ -36,7 +40,7 @@ class bcolors: self.BOLD = '' self.UNDERLINE = '' self.CROSSOUT = '' - + # ----------------------------- def srepr(arg,glue = '\n'): @@ -159,11 +163,29 @@ def progressBar(iteration, total, prefix='', bar_length=50): if iteration == total: sys.stderr.write('\n') sys.stderr.flush() - + + +def scale_to_coprime(v): + """Scale vector to co-prime (relatively prime) integers.""" + MAX_DENOMINATOR = 1000 + + def get_square_denominator(x): + """Denominator of the square of a number.""" + return Fraction(x ** 2).limit_denominator(MAX_DENOMINATOR).denominator + + def lcm(a, b): + """Least common multiple.""" + return a * b // np.gcd(a, b) + + denominators = [int(get_square_denominator(i)) for i in v] + s = reduce(lcm, denominators) ** 0.5 + m = (np.array(v)*s).astype(np.int) + return m//reduce(np.gcd,m) + class return_message(): """Object with formatted return message.""" - + def __init__(self,message): """ Sets return message. @@ -175,7 +197,7 @@ class return_message(): """ self.message = message - + def __repr__(self): """Return message suitable for interactive shells.""" return srepr(self.message) diff --git a/python/tests/test_DADF5.py b/python/tests/test_DADF5.py index 669ce8446..4cf3c4e2b 100644 --- a/python/tests/test_DADF5.py +++ b/python/tests/test_DADF5.py @@ -40,7 +40,7 @@ class TestDADF5: assert np.allclose(in_memory,in_file) def test_add_calculation(self,default): - default.add_calculation('2.0*np.abs(#F#)-1.0','x','-','test') + default.add_calculation('x','2.0*np.abs(#F#)-1.0','-','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 @@ -52,8 +52,8 @@ class TestDADF5: loc = {'F': default.get_dataset_location('F'), 'P': default.get_dataset_location('P'), 'sigma':default.get_dataset_location('sigma')} - in_memory = mechanics.Cauchy(default.read_dataset(loc['F'],0), - default.read_dataset(loc['P'],0)) + in_memory = mechanics.Cauchy(default.read_dataset(loc['P'],0), + default.read_dataset(loc['F'],0)) in_file = default.read_dataset(loc['sigma'],0) assert np.allclose(in_memory,in_file) @@ -73,6 +73,54 @@ class TestDADF5: in_file = default.read_dataset(loc['s_P'],0) assert np.allclose(in_memory,in_file) + def test_add_eigenvalues(self,default): + default.add_Cauchy('P','F') + default.add_eigenvalues('sigma') + loc = {'sigma' :default.get_dataset_location('sigma'), + 'lambda(sigma)':default.get_dataset_location('lambda(sigma)')} + in_memory = mechanics.eigenvalues(default.read_dataset(loc['sigma'],0)) + in_file = default.read_dataset(loc['lambda(sigma)'],0) + assert np.allclose(in_memory,in_file) + + def test_add_eigenvectors(self,default): + default.add_Cauchy('P','F') + default.add_eigenvectors('sigma') + loc = {'sigma' :default.get_dataset_location('sigma'), + 'v(sigma)':default.get_dataset_location('v(sigma)')} + in_memory = mechanics.eigenvectors(default.read_dataset(loc['sigma'],0)) + in_file = default.read_dataset(loc['v(sigma)'],0) + assert np.allclose(in_memory,in_file) + + def test_add_maximum_shear(self,default): + default.add_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) + assert np.allclose(in_memory,in_file) + + def test_add_Mises_strain(self,default): + t = ['V','U'][np.random.randint(0,2)] + m = np.random.random()*2.0 - 1.0 + default.add_strain_tensor('F',t,m) + label = 'epsilon_{}^{}(F)'.format(t,m) + default.add_Mises(label) + loc = {label :default.get_dataset_location(label), + label+'_vM':default.get_dataset_location(label+'_vM')} + in_memory = mechanics.Mises_strain(default.read_dataset(loc[label],0)).reshape(-1,1) + in_file = default.read_dataset(loc[label+'_vM'],0) + assert np.allclose(in_memory,in_file) + + def test_add_Mises_stress(self,default): + default.add_Cauchy('P','F') + default.add_Mises('sigma') + loc = {'sigma' :default.get_dataset_location('sigma'), + 'sigma_vM':default.get_dataset_location('sigma_vM')} + in_memory = mechanics.Mises_stress(default.read_dataset(loc['sigma'],0)).reshape(-1,1) + in_file = default.read_dataset(loc['sigma_vM'],0) + assert np.allclose(in_memory,in_file) + def test_add_norm(self,default): default.add_norm('F',1) loc = {'F': default.get_dataset_location('F'), @@ -81,6 +129,24 @@ class TestDADF5: in_file = default.read_dataset(loc['|F|_1'],0) assert np.allclose(in_memory,in_file) + def test_add_PK2(self,default): + default.add_PK2('P','F') + loc = {'F':default.get_dataset_location('F'), + 'P':default.get_dataset_location('P'), + 'S':default.get_dataset_location('S')} + in_memory = mechanics.PK2(default.read_dataset(loc['P'],0), + default.read_dataset(loc['F'],0)) + in_file = default.read_dataset(loc['S'],0) + assert np.allclose(in_memory,in_file) + + def test_add_rotational_part(self,default): + default.add_rotational_part('F') + loc = {'F': default.get_dataset_location('F'), + 'R(F)': default.get_dataset_location('R(F)')} + in_memory = mechanics.rotational_part(default.read_dataset(loc['F'],0)) + in_file = default.read_dataset(loc['R(F)'],0) + assert np.allclose(in_memory,in_file) + def test_add_spherical(self,default): default.add_spherical('P') loc = {'P': default.get_dataset_location('P'), @@ -88,3 +154,30 @@ class TestDADF5: in_memory = mechanics.spherical_part(default.read_dataset(loc['P'],0)).reshape(-1,1) in_file = default.read_dataset(loc['p_P'],0) assert np.allclose(in_memory,in_file) + + def test_add_strain(self,default): + t = ['V','U'][np.random.randint(0,2)] + m = np.random.random()*2.0 - 1.0 + default.add_strain_tensor('F',t,m) + label = 'epsilon_{}^{}(F)'.format(t,m) + loc = {'F': default.get_dataset_location('F'), + label: default.get_dataset_location(label)} + in_memory = mechanics.strain_tensor(default.read_dataset(loc['F'],0),t,m) + in_file = default.read_dataset(loc[label],0) + 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.right_stretch(default.read_dataset(loc['F'],0)) + in_file = default.read_dataset(loc['U(F)'],0) + 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.left_stretch(default.read_dataset(loc['F'],0)) + in_file = default.read_dataset(loc['V(F)'],0) + assert np.allclose(in_memory,in_file) diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py index 9e1d9bc0c..4248254ab 100644 --- a/python/tests/test_mechanics.py +++ b/python/tests/test_mechanics.py @@ -2,187 +2,224 @@ import numpy as np from damask import mechanics class TestMechanics: - + n = 1000 c = np.random.randint(n) - + def test_vectorize_Cauchy(self): - P = np.random.random((self.n,3,3)) - F = np.random.random((self.n,3,3)) - assert np.allclose(mechanics.Cauchy(F,P)[self.c], - mechanics.Cauchy(F[self.c],P[self.c])) - - - def test_vectorize_strain_tensor(self): - F = np.random.random((self.n,3,3)) - t = ['V','U'][np.random.randint(0,2)] - m = np.random.random()*10. -5.0 - assert np.allclose(mechanics.strain_tensor(F,t,m)[self.c], - mechanics.strain_tensor(F[self.c],t,m)) - + P = np.random.random((self.n,3,3)) + F = np.random.random((self.n,3,3)) + assert np.allclose(mechanics.Cauchy(P,F)[self.c], + mechanics.Cauchy(P[self.c],F[self.c])) def test_vectorize_deviatoric_part(self): - x = np.random.random((self.n,3,3)) - assert np.allclose(mechanics.deviatoric_part(x)[self.c], - mechanics.deviatoric_part(x[self.c])) + x = np.random.random((self.n,3,3)) + assert np.allclose(mechanics.deviatoric_part(x)[self.c], + mechanics.deviatoric_part(x[self.c])) + def test_vectorize_eigenvalues(self): + x = np.random.random((self.n,3,3)) + assert np.allclose(mechanics.eigenvalues(x)[self.c], + mechanics.eigenvalues(x[self.c])) - def test_vectorize_spherical_part(self): - x = np.random.random((self.n,3,3)) - assert np.allclose(mechanics.spherical_part(x,True)[self.c], - mechanics.spherical_part(x[self.c],True)) - - - def test_vectorize_Mises_stress(self): - sigma = np.random.random((self.n,3,3)) - assert np.allclose(mechanics.Mises_stress(sigma)[self.c], - mechanics.Mises_stress(sigma[self.c])) - - - def test_vectorize_Mises_strain(self): - epsilon = np.random.random((self.n,3,3)) - assert np.allclose(mechanics.Mises_strain(epsilon)[self.c], - mechanics.Mises_strain(epsilon[self.c])) - - - def test_vectorize_symmetric(self): - x = np.random.random((self.n,3,3)) - assert np.allclose(mechanics.symmetric(x)[self.c], - mechanics.symmetric(x[self.c])) - - - def test_vectorize_maximum_shear(self): - x = np.random.random((self.n,3,3)) - assert np.allclose(mechanics.maximum_shear(x)[self.c], - mechanics.maximum_shear(x[self.c])) - - - def test_vectorize_principal_components(self): - x = np.random.random((self.n,3,3)) - assert np.allclose(mechanics.principal_components(x)[self.c], - mechanics.principal_components(x[self.c])) - - - def test_vectorize_transpose(self): - x = np.random.random((self.n,3,3)) - assert np.allclose(mechanics.transpose(x)[self.c], - mechanics.transpose(x[self.c])) - - - def test_vectorize_rotational_part(self): - x = np.random.random((self.n,3,3)) - assert np.allclose(mechanics.rotational_part(x)[self.c], - mechanics.rotational_part(x[self.c])) - + def test_vectorize_eigenvectors(self): + x = np.random.random((self.n,3,3)) + assert np.allclose(mechanics.eigenvectors(x)[self.c], + mechanics.eigenvectors(x[self.c])) def test_vectorize_left_stretch(self): - x = np.random.random((self.n,3,3)) - assert np.allclose(mechanics.left_stretch(x)[self.c], - mechanics.left_stretch(x[self.c])) + x = np.random.random((self.n,3,3)) + assert np.allclose(mechanics.left_stretch(x)[self.c], + mechanics.left_stretch(x[self.c])) + def test_vectorize_maximum_shear(self): + x = np.random.random((self.n,3,3)) + assert np.allclose(mechanics.maximum_shear(x)[self.c], + mechanics.maximum_shear(x[self.c])) + + def test_vectorize_Mises_strain(self): + epsilon = np.random.random((self.n,3,3)) + assert np.allclose(mechanics.Mises_strain(epsilon)[self.c], + mechanics.Mises_strain(epsilon[self.c])) + + def test_vectorize_Mises_stress(self): + sigma = np.random.random((self.n,3,3)) + assert np.allclose(mechanics.Mises_stress(sigma)[self.c], + mechanics.Mises_stress(sigma[self.c])) + + def test_vectorize_PK2(self): + F = np.random.random((self.n,3,3)) + P = np.random.random((self.n,3,3)) + assert np.allclose(mechanics.PK2(P,F)[self.c], + mechanics.PK2(P[self.c],F[self.c])) def test_vectorize_right_stretch(self): - x = np.random.random((self.n,3,3)) - assert np.allclose(mechanics.right_stretch(x)[self.c], - mechanics.right_stretch(x[self.c])) + x = np.random.random((self.n,3,3)) + assert np.allclose(mechanics.right_stretch(x)[self.c], + mechanics.right_stretch(x[self.c])) + + def test_vectorize_rotational_part(self): + x = np.random.random((self.n,3,3)) + assert np.allclose(mechanics.rotational_part(x)[self.c], + mechanics.rotational_part(x[self.c])) + + def test_vectorize_spherical_part(self): + x = np.random.random((self.n,3,3)) + assert np.allclose(mechanics.spherical_part(x,True)[self.c], + mechanics.spherical_part(x[self.c],True)) + + def test_vectorize_strain_tensor(self): + F = np.random.random((self.n,3,3)) + t = ['V','U'][np.random.randint(0,2)] + m = np.random.random()*10. -5.0 + assert np.allclose(mechanics.strain_tensor(F,t,m)[self.c], + mechanics.strain_tensor(F[self.c],t,m)) + + def test_vectorize_symmetric(self): + x = np.random.random((self.n,3,3)) + assert np.allclose(mechanics.symmetric(x)[self.c], + mechanics.symmetric(x[self.c])) + + def test_vectorize_transpose(self): + x = np.random.random((self.n,3,3)) + assert np.allclose(mechanics.transpose(x)[self.c], + mechanics.transpose(x[self.c])) def test_Cauchy(self): - """Ensure Cauchy stress is symmetrized 1. Piola-Kirchhoff stress for no deformation.""" - P = np.random.random((self.n,3,3)) - assert np.allclose(mechanics.Cauchy(np.broadcast_to(np.eye(3),(self.n,3,3)),P), - mechanics.symmetric(P)) + """Ensure Cauchy stress is symmetrized 1. Piola-Kirchhoff stress for no deformation.""" + P = np.random.random((self.n,3,3)) + assert np.allclose(mechanics.Cauchy(P,np.broadcast_to(np.eye(3),(self.n,3,3))), + mechanics.symmetric(P)) + def test_polar_decomposition(self): - """F = RU = VR.""" - F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.random((self.n,3,3)) - R = mechanics.rotational_part(F) - V = mechanics.left_stretch(F) - U = mechanics.right_stretch(F) - assert np.allclose(np.matmul(R,U), - np.matmul(V,R)) - + """F = RU = VR.""" + F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.random((self.n,3,3)) + R = mechanics.rotational_part(F) + V = mechanics.left_stretch(F) + U = mechanics.right_stretch(F) + assert np.allclose(np.matmul(R,U), + np.matmul(V,R)) + + + def test_PK2(self): + """Ensure 2. Piola-Kirchhoff stress is symmetrized 1. Piola-Kirchhoff stress for no deformation.""" + P = np.random.random((self.n,3,3)) + assert np.allclose(mechanics.PK2(P,np.broadcast_to(np.eye(3),(self.n,3,3))), + mechanics.symmetric(P)) + def test_strain_tensor_no_rotation(self): - """Ensure that left and right stretch give same results for no rotation.""" - F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.random((self.n,3,3)) - m = np.random.random()*20.0-10.0 - assert np.allclose(mechanics.strain_tensor(F,'U',m), - mechanics.strain_tensor(F,'V',m)) - + """Ensure that left and right stretch give same results for no rotation.""" + F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.random((self.n,3,3)) + m = np.random.random()*20.0-10.0 + assert np.allclose(mechanics.strain_tensor(F,'U',m), + mechanics.strain_tensor(F,'V',m)) + def test_strain_tensor_rotation_equivalence(self): - """Ensure that left and right strain differ only by a rotation.""" - F = np.broadcast_to(np.eye(3),[self.n,3,3]) + (np.random.random((self.n,3,3))*0.5 - 0.25) - m = np.random.random()*5.0-2.5 - assert np.allclose(np.linalg.det(mechanics.strain_tensor(F,'U',m)), - np.linalg.det(mechanics.strain_tensor(F,'V',m))) + """Ensure that left and right strain differ only by a rotation.""" + F = np.broadcast_to(np.eye(3),[self.n,3,3]) + (np.random.random((self.n,3,3))*0.5 - 0.25) + m = np.random.random()*5.0-2.5 + assert np.allclose(np.linalg.det(mechanics.strain_tensor(F,'U',m)), + np.linalg.det(mechanics.strain_tensor(F,'V',m))) def test_strain_tensor_rotation(self): - """Ensure that pure rotation results in no strain.""" - F = mechanics.rotational_part(np.random.random((self.n,3,3))) - t = ['V','U'][np.random.randint(0,2)] - m = np.random.random()*2.0 - 1.0 - assert np.allclose(mechanics.strain_tensor(F,t,m), - 0.0) - - def test_rotation_determinant(self): - """ - Ensure that the determinant of the rotational part is +- 1. + """Ensure that pure rotation results in no strain.""" + F = mechanics.rotational_part(np.random.random((self.n,3,3))) + t = ['V','U'][np.random.randint(0,2)] + m = np.random.random()*2.0 - 1.0 + assert np.allclose(mechanics.strain_tensor(F,t,m), + 0.0) - Should be +1, but random F might contain a reflection. - """ - x = np.random.random((self.n,3,3)) - assert np.allclose(np.abs(np.linalg.det(mechanics.rotational_part(x))), - 1.0) + def test_rotation_determinant(self): + """ + Ensure that the determinant of the rotational part is +- 1. + + Should be +1, but random F might contain a reflection. + """ + x = np.random.random((self.n,3,3)) + assert np.allclose(np.abs(np.linalg.det(mechanics.rotational_part(x))), + 1.0) def test_spherical_deviatoric_part(self): - """Ensure that full tensor is sum of spherical and deviatoric part.""" - x = np.random.random((self.n,3,3)) - sph = mechanics.spherical_part(x,True) - assert np.allclose(sph + mechanics.deviatoric_part(x), - x) + """Ensure that full tensor is sum of spherical and deviatoric part.""" + x = np.random.random((self.n,3,3)) + sph = mechanics.spherical_part(x,True) + assert np.allclose(sph + mechanics.deviatoric_part(x), + x) def test_deviatoric_Mises(self): - """Ensure that Mises equivalent stress depends only on deviatoric part.""" - x = np.random.random((self.n,3,3)) - full = mechanics.Mises_stress(x) - dev = mechanics.Mises_stress(mechanics.deviatoric_part(x)) - assert np.allclose(full, - dev) + """Ensure that Mises equivalent stress depends only on deviatoric part.""" + x = np.random.random((self.n,3,3)) + full = mechanics.Mises_stress(x) + dev = mechanics.Mises_stress(mechanics.deviatoric_part(x)) + assert np.allclose(full, + dev) def test_spherical_mapping(self): - """Ensure that mapping to tensor is correct.""" - x = np.random.random((self.n,3,3)) - tensor = mechanics.spherical_part(x,True) - scalar = mechanics.spherical_part(x) - assert np.allclose(np.linalg.det(tensor), - scalar**3.0) + """Ensure that mapping to tensor is correct.""" + x = np.random.random((self.n,3,3)) + tensor = mechanics.spherical_part(x,True) + scalar = mechanics.spherical_part(x) + assert np.allclose(np.linalg.det(tensor), + scalar**3.0) def test_spherical_Mises(self): - """Ensure that Mises equivalent strrain of spherical strain is 0.""" - x = np.random.random((self.n,3,3)) - sph = mechanics.spherical_part(x,True) - assert np.allclose(mechanics.Mises_strain(sph), - 0.0) + """Ensure that Mises equivalent strrain of spherical strain is 0.""" + x = np.random.random((self.n,3,3)) + sph = mechanics.spherical_part(x,True) + assert np.allclose(mechanics.Mises_strain(sph), + 0.0) def test_symmetric(self): - """Ensure that a symmetric tensor is half of the sum of a tensor and its transpose.""" - x = np.random.random((self.n,3,3)) - assert np.allclose(mechanics.symmetric(x)*2.0, - mechanics.transpose(x)+x) + """Ensure that a symmetric tensor is half of the sum of a tensor and its transpose.""" + x = np.random.random((self.n,3,3)) + assert np.allclose(mechanics.symmetric(x)*2.0, + mechanics.transpose(x)+x) def test_transpose(self): - """Ensure that a symmetric tensor equals its transpose.""" - x = mechanics.symmetric(np.random.random((self.n,3,3))) - assert np.allclose(mechanics.transpose(x), - x) + """Ensure that a symmetric tensor equals its transpose.""" + x = mechanics.symmetric(np.random.random((self.n,3,3))) + assert np.allclose(mechanics.transpose(x), + x) def test_Mises(self): - """Ensure that equivalent stress is 3/2 of equivalent strain.""" - x = np.random.random((self.n,3,3)) - assert np.allclose(mechanics.Mises_stress(x)/mechanics.Mises_strain(x), - 1.5) + """Ensure that equivalent stress is 3/2 of equivalent strain.""" + x = np.random.random((self.n,3,3)) + assert np.allclose(mechanics.Mises_stress(x)/mechanics.Mises_strain(x), + 1.5) + + + def test_eigenvalues(self): + """Ensure that the characteristic polynomial can be solved.""" + A = mechanics.symmetric(np.random.random((self.n,3,3))) + lambd = mechanics.eigenvalues(A) + s = np.random.randint(self.n) + for i in range(3): + assert np.allclose(np.linalg.det(A[s]-lambd[s,i]*np.eye(3)),.0) + + def test_eigenvalues_and_vectors(self): + """Ensure that eigenvalues and -vectors are the solution to the characteristic polynomial.""" + A = mechanics.symmetric(np.random.random((self.n,3,3))) + lambd = mechanics.eigenvalues(A) + x = mechanics.eigenvectors(A) + s = np.random.randint(self.n) + for i in range(3): + assert np.allclose(np.dot(A[s]-lambd[s,i]*np.eye(3),x[s,:,i]),.0) + + def test_eigenvectors_RHS(self): + """Ensure that RHS coordinate system does only change sign of determinant.""" + A = mechanics.symmetric(np.random.random((self.n,3,3))) + LRHS = np.linalg.det(mechanics.eigenvectors(A,RHS=False)) + RHS = np.linalg.det(mechanics.eigenvectors(A,RHS=True)) + assert np.allclose(np.abs(LRHS),RHS) + + def test_spherical_no_shear(self): + """Ensure that sherical stress has max shear of 0.0.""" + A = mechanics.spherical_part(mechanics.symmetric(np.random.random((self.n,3,3))),True) + assert np.allclose(mechanics.maximum_shear(A),0.0) diff --git a/src/constitutive_plastic_disloUCLA.f90 b/src/constitutive_plastic_disloUCLA.f90 index 2b353ba26..00a1b8e21 100644 --- a/src/constitutive_plastic_disloUCLA.f90 +++ b/src/constitutive_plastic_disloUCLA.f90 @@ -6,21 +6,10 @@ !> @brief crystal plasticity model for bcc metals, especially Tungsten !-------------------------------------------------------------------------------------------------- submodule(constitutive) plastic_disloUCLA - + real(pReal), parameter :: & kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin - - enum, bind(c) - enumerator :: & - undefined_ID, & - rho_mob_ID, & - rho_dip_ID, & - dot_gamma_sl_ID, & - gamma_sl_ID, & - Lambda_sl_ID, & - tau_pass_ID - end enum - + type :: tParameters real(pReal) :: & aTol_rho, & @@ -28,7 +17,7 @@ submodule(constitutive) plastic_disloUCLA mu, & D_0, & !< prefactor for self-diffusion coefficient Q_cl !< activation energy for dislocation climb - real(pReal), dimension(:), allocatable :: & + real(pReal), allocatable, dimension(:) :: & rho_mob_0, & !< initial dislocation density rho_dip_0, & !< initial dipole density b_sl, & !< magnitude of burgers vector [m] @@ -46,30 +35,30 @@ submodule(constitutive) plastic_disloUCLA kink_height, & !< height of the kink pair w, & !< width of the kink pair omega !< attempt frequency for kink pair nucleation - real(pReal), dimension(:,:), allocatable :: & + real(pReal), allocatable, dimension(:,:) :: & h_sl_sl, & !< slip resistance from slip activity forestProjectionEdge - real(pReal), dimension(:,:,:), allocatable :: & + real(pReal), allocatable, dimension(:,:,:) :: & Schmid, & nonSchmid_pos, & nonSchmid_neg integer :: & sum_N_sl !< total number of active slip system - integer, dimension(:), allocatable :: & + integer, allocatable, dimension(:) :: & N_sl !< number of active slip systems for each family - integer(kind(undefined_ID)), dimension(:),allocatable :: & - outputID !< ID of each post result output + character(len=pStringLen), allocatable, dimension(:) :: & + output logical :: & dipoleFormation !< flag indicating consideration of dipole formation end type !< container type for internal constitutive parameters - + type :: tDisloUCLAState real(pReal), dimension(:,:), pointer :: & rho_mob, & rho_dip, & gamma_sl end type tDisloUCLAState - + type :: tDisloUCLAdependentState real(pReal), dimension(:,:), allocatable :: & Lambda_sl, & @@ -88,41 +77,36 @@ contains !-------------------------------------------------------------------------------------------------- -!> @brief module initialization +!> @brief Perform module initialization. !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- module subroutine plastic_disloUCLA_init - + integer :: & Ninstance, & p, i, & NipcMyPhase, & sizeState, sizeDotState, & startIndex, endIndex - - integer(kind(undefined_ID)) :: & - outputID - + character(len=pStringLen) :: & extmsg = '' - character(len=pStringLen), dimension(:), allocatable :: & - outputs - - write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_DISLOUCLA_label//' init -+>>>' - - write(6,'(/,a)') ' Cereceda et al., International Journal of Plasticity 78:242–256, 2016' - write(6,'(a)') ' https://dx.doi.org/10.1016/j.ijplas.2015.09.002' - + + write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_DISLOUCLA_label//' init -+>>>'; flush(6) + + write(6,'(/,a)') ' Cereceda et al., International Journal of Plasticity 78:242–256, 2016' + write(6,'(a)') ' https://dx.doi.org/10.1016/j.ijplas.2015.09.002' + Ninstance = count(phase_plasticity == PLASTICITY_DISLOUCLA_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - + allocate(param(Ninstance)) allocate(state(Ninstance)) allocate(dotState(Ninstance)) allocate(dependentState(Ninstance)) - - + + do p = 1, size(phase_plasticity) if (phase_plasticity(p) /= PLASTICITY_DISLOUCLA_ID) cycle associate(prm => param(phase_plasticityInstance(p)), & @@ -134,9 +118,9 @@ module subroutine plastic_disloUCLA_init !-------------------------------------------------------------------------------------------------- ! optional parameters that need to be defined prm%mu = lattice_mu(p) - + prm%aTol_rho = config%getFloat('atol_rho') - + ! sanity checks if (prm%aTol_rho <= 0.0_pReal) extmsg = trim(extmsg)//' atol_rho' @@ -147,7 +131,7 @@ module subroutine plastic_disloUCLA_init slipActive: if (prm%sum_N_sl > 0) then prm%Schmid = lattice_SchmidMatrix_slip(prm%N_sl,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) - + if(trim(config%getString('lattice_structure')) == 'bcc') then prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',& defaultVal = emptyRealArray) @@ -158,20 +142,20 @@ module subroutine plastic_disloUCLA_init prm%nonSchmid_pos = prm%Schmid prm%nonSchmid_neg = prm%Schmid endif - + prm%h_sl_sl = lattice_interaction_SlipBySlip(prm%N_sl, & config%getFloats('interaction_slipslip'), & config%getString('lattice_structure')) prm%forestProjectionEdge = lattice_forestProjection_edge(prm%N_sl,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) prm%forestProjectionEdge = transpose(prm%forestProjectionEdge) - + prm%rho_mob_0 = config%getFloats('rhoedge0', requiredSize=size(prm%N_sl)) prm%rho_dip_0 = config%getFloats('rhoedgedip0', requiredSize=size(prm%N_sl)) prm%v0 = config%getFloats('v0', requiredSize=size(prm%N_sl)) prm%b_sl = config%getFloats('slipburgers', requiredSize=size(prm%N_sl)) prm%delta_F = config%getFloats('qedge', requiredSize=size(prm%N_sl)) - + prm%i_sl = config%getFloats('clambdaslip', requiredSize=size(prm%N_sl)) prm%tau_0 = config%getFloats('tau_peierls', requiredSize=size(prm%N_sl)) prm%p = config%getFloats('p_slip', requiredSize=size(prm%N_sl), & @@ -182,14 +166,14 @@ module subroutine plastic_disloUCLA_init prm%w = config%getFloats('kink_width', requiredSize=size(prm%N_sl)) prm%omega = config%getFloats('omega', requiredSize=size(prm%N_sl)) prm%B = config%getFloats('friction_coeff', requiredSize=size(prm%N_sl)) - + prm%D = config%getFloat('grainsize') prm%D_0 = config%getFloat('d0') prm%Q_cl = config%getFloat('qsd') prm%atomicVolume = config%getFloat('catomicvolume') * prm%b_sl**3.0_pReal prm%D_a = config%getFloat('cedgedipmindistance') * prm%b_sl prm%dipoleformation = config%getFloat('dipoleformationfactor') > 0.0_pReal !should be on by default, ToDo: change to /key/-type key - + ! expand: family => system prm%rho_mob_0 = math_expand(prm%rho_mob_0, prm%N_sl) prm%rho_dip_0 = math_expand(prm%rho_dip_0, prm%N_sl) @@ -206,8 +190,8 @@ module subroutine plastic_disloUCLA_init prm%i_sl = math_expand(prm%i_sl, prm%N_sl) prm%atomicVolume = math_expand(prm%atomicVolume, prm%N_sl) prm%D_a = math_expand(prm%D_a, prm%N_sl) - - + + ! sanity checks if ( prm%D_0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0' if ( prm%Q_cl <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl' @@ -219,7 +203,7 @@ module subroutine plastic_disloUCLA_init if (any(prm%tau_0 < 0.0_pReal)) extmsg = trim(extmsg)//' tau_0' if (any(prm%D_a <= 0.0_pReal)) extmsg = trim(extmsg)//' cedgedipmindistance or slipb_sl' if (any(prm%atomicVolume <= 0.0_pReal)) extmsg = trim(extmsg)//' catomicvolume or slipb_sl' - + else slipActive allocate(prm%rho_mob_0(0)) allocate(prm%rho_dip_0(0)) @@ -232,39 +216,14 @@ module subroutine plastic_disloUCLA_init !-------------------------------------------------------------------------------------------------- ! output pararameters - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - do i=1, size(outputs) - outputID = undefined_ID - select case(trim(outputs(i))) - - case ('edge_density') - outputID = merge(rho_mob_ID,undefined_ID,prm%sum_N_sl>0) - case ('dipole_density') - outputID = merge(rho_dip_ID,undefined_ID,prm%sum_N_sl>0) - case ('shear_rate','shearrate','shear_rate_slip','shearrate_slip') - outputID = merge(dot_gamma_sl_ID,undefined_ID,prm%sum_N_sl>0) - case ('accumulated_shear','accumulatedshear','accumulated_shear_slip') - outputID = merge(gamma_sl_ID,undefined_ID,prm%sum_N_sl>0) - case ('mfp','mfp_slip') - outputID = merge(Lambda_sl_ID,undefined_ID,prm%sum_N_sl>0) - case ('threshold_stress','threshold_stress_slip') - outputID = merge(tau_pass_ID,undefined_ID,prm%sum_N_sl>0) - - end select - - if (outputID /= undefined_ID) then - prm%outputID = [prm%outputID, outputID] - endif - - enddo + prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) !-------------------------------------------------------------------------------------------------- ! allocate state arrays NipcMyPhase = count(material_phaseAt == p) * discretization_nIP sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl sizeState = sizeDotState - + call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) !-------------------------------------------------------------------------------------------------- @@ -275,14 +234,14 @@ module subroutine plastic_disloUCLA_init stt%rho_mob= spread(prm%rho_mob_0,2,NipcMyPhase) dot%rho_mob=>plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_rho - + startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl stt%rho_dip=>plasticState(p)%state(startIndex:endIndex,:) stt%rho_dip= spread(prm%rho_dip_0,2,NipcMyPhase) dot%rho_dip=>plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_rho - + startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl stt%gamma_sl=>plasticState(p)%state(startIndex:endIndex,:) @@ -290,21 +249,21 @@ module subroutine plastic_disloUCLA_init plasticState(p)%aTolState(startIndex:endIndex) = 1.0e6_pReal ! Don't use for convergence check ! global alias plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) - + allocate(dst%Lambda_sl(prm%sum_N_sl,NipcMyPhase), source=0.0_pReal) allocate(dst%threshold_stress(prm%sum_N_sl,NipcMyPhase), source=0.0_pReal) - + plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally - + end associate - + enddo end subroutine plastic_disloUCLA_init !-------------------------------------------------------------------------------------------------- -!> @brief calculates plastic velocity gradient and its tangent +!> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp, & Mp,T,instance,of) @@ -312,7 +271,7 @@ pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp, & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress real(pReal), intent(in) :: & @@ -320,18 +279,18 @@ pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp, & integer, intent(in) :: & instance, & of - + integer :: & i,k,l,m,n real(pReal), dimension(param(instance)%sum_N_sl) :: & dot_gamma_pos,dot_gamma_neg, & ddot_gamma_dtau_pos,ddot_gamma_dtau_neg - + Lp = 0.0_pReal dLp_dMp = 0.0_pReal - + associate(prm => param(instance)) - + call kinetics(Mp,T,instance,of,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg) do i = 1, prm%sum_N_sl Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%Schmid(1:3,1:3,i) @@ -340,17 +299,17 @@ pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp, & + ddot_gamma_dtau_pos(i) * prm%Schmid(k,l,i) * prm%nonSchmid_pos(m,n,i) & + ddot_gamma_dtau_neg(i) * prm%Schmid(k,l,i) * prm%nonSchmid_neg(m,n,i) enddo - + end associate end subroutine plastic_disloUCLA_LpAndItsTangent !-------------------------------------------------------------------------------------------------- -!> @brief calculates the rate of change of microstructure +!> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- module subroutine plastic_disloUCLA_dotState(Mp,T,instance,of) - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress real(pReal), intent(in) :: & @@ -358,7 +317,7 @@ module subroutine plastic_disloUCLA_dotState(Mp,T,instance,of) integer, intent(in) :: & instance, & of - + real(pReal) :: & VacancyDiffusion real(pReal), dimension(param(instance)%sum_N_sl) :: & @@ -369,16 +328,16 @@ module subroutine plastic_disloUCLA_dotState(Mp,T,instance,of) dot_rho_dip_formation, & dot_rho_dip_climb, & dip_distance - + associate(prm => param(instance), stt => state(instance),dot => dotState(instance), dst => dependentState(instance)) - + call kinetics(Mp,T,instance,of,& gdot_pos,gdot_neg, & tau_pos_out = tau_pos,tau_neg_out = tau_neg) - + dot%gamma_sl(:,of) = (gdot_pos+gdot_neg) ! ToDo: needs to be abs VacancyDiffusion = prm%D_0*exp(-prm%Q_cl/(kB*T)) - + where(dEq0(tau_pos)) ! ToDo: use avg of pos and neg dot_rho_dip_formation = 0.0_pReal dot_rho_dip_climb = 0.0_pReal @@ -393,73 +352,73 @@ module subroutine plastic_disloUCLA_dotState(Mp,T,instance,of) * (1.0_pReal/(dip_distance+prm%D_a)) dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,of))/(dip_distance-prm%D_a) ! ToDo: Discuss with Franz: Stress dependency? end where - + dot%rho_mob(:,of) = abs(dot%gamma_sl(:,of))/(prm%b_sl*dst%Lambda_sl(:,of)) & ! multiplication - dot_rho_dip_formation & - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of)) ! Spontaneous annihilation of 2 single edge dislocations dot%rho_dip(:,of) = dot_rho_dip_formation & - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_dip(:,of)*abs(dot%gamma_sl(:,of)) & ! Spontaneous annihilation of a single edge dislocation with a dipole constituent - dot_rho_dip_climb - + end associate end subroutine plastic_disloUCLA_dotState !-------------------------------------------------------------------------------------------------- -!> @brief calculates derived quantities from state +!> @brief Calculate derived quantities from state. !-------------------------------------------------------------------------------------------------- module subroutine plastic_disloUCLA_dependentState(instance,of) - + integer, intent(in) :: & instance, & of - + real(pReal), dimension(param(instance)%sum_N_sl) :: & dislocationSpacing - + associate(prm => param(instance), stt => state(instance),dst => dependentState(instance)) - + dislocationSpacing = sqrt(matmul(prm%forestProjectionEdge, & stt%rho_mob(:,of)+stt%rho_dip(:,of))) dst%threshold_stress(:,of) = prm%mu*prm%b_sl & * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,of)+stt%rho_dip(:,of))) - + dst%Lambda_sl(:,of) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl) - + end associate end subroutine plastic_disloUCLA_dependentState !-------------------------------------------------------------------------------------------------- -!> @brief writes results to HDF5 output file +!> @brief Write results to HDF5 output file. !-------------------------------------------------------------------------------------------------- module subroutine plastic_disloUCLA_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group - + integer :: o associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - case (rho_mob_ID) - call results_writeDataset(group,stt%rho_mob,'rho_mob',& - 'mobile dislocation density','1/m²') - case (rho_dip_ID) - call results_writeDataset(group,stt%rho_dip,'rho_dip',& - 'dislocation dipole density''1/m²') - case (dot_gamma_sl_ID) - call results_writeDataset(group,stt%gamma_sl,'dot_gamma_sl',& ! this is not dot!! - 'plastic shear','1') - case (Lambda_sl_ID) - call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',& - 'mean free path for slip','m') - case (tau_pass_ID) - call results_writeDataset(group,dst%threshold_stress,'tau_pass',& - 'threshold stress for slip','Pa') + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case('edge_density') ! ToDo: should be rho_mob + if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_mob,'rho_mob',& + 'mobile dislocation density','1/m²') + case('dipole_density') ! ToDo: should be rho_dip + if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_dip,'rho_dip',& + 'dislocation dipole density''1/m²') + case('shear_rate_slip') ! should be gamma + if(prm%sum_N_sl>0) call results_writeDataset(group,stt%gamma_sl,'dot_gamma_sl',& ! this is not dot!! + 'plastic shear','1') + case('mfp_slip') !ToDo: should be Lambda + if(prm%sum_N_sl>0) call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',& + 'mean free path for slip','m') + case('threshold_stress_slip') !ToDo: should be tau_pass + if(prm%sum_N_sl>0) call results_writeDataset(group,dst%threshold_stress,'tau_pass',& + 'threshold stress for slip','Pa') end select enddo outputsLoop end associate @@ -468,15 +427,15 @@ end subroutine plastic_disloUCLA_results !-------------------------------------------------------------------------------------------------- -!> @brief Shear rates on slip systems, their derivatives with respect to resolved stress and the -! resolved stresss +!> @brief Calculate shear rates on slip systems, their derivatives with respect to resolved +! stress, and the resolved stress. !> @details Derivatives and resolved stress are calculated only optionally. ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end !-------------------------------------------------------------------------------------------------- pure subroutine kinetics(Mp,T,instance,of, & dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg,tau_pos_out,tau_neg_out) - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress real(pReal), intent(in) :: & @@ -484,7 +443,7 @@ pure subroutine kinetics(Mp,T,instance,of, & integer, intent(in) :: & instance, & of - + real(pReal), intent(out), dimension(param(instance)%sum_N_sl) :: & dot_gamma_pos, & dot_gamma_neg @@ -501,82 +460,82 @@ pure subroutine kinetics(Mp,T,instance,of, & t_n, t_k, dtk,dtn, & needsGoodName ! ToDo: @Karo: any idea? integer :: j - + associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) - + do j = 1, prm%sum_N_sl tau_pos(j) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,j)) tau_neg(j) = math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,j)) enddo - - + + if (present(tau_pos_out)) tau_pos_out = tau_pos if (present(tau_neg_out)) tau_neg_out = tau_neg - + associate(BoltzmannRatio => prm%delta_F/(kB*T), & dot_gamma_0 => stt%rho_mob(:,of)*prm%b_sl*prm%v0, & effectiveLength => dst%Lambda_sl(:,of) - prm%w) - + significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check) StressRatio = (abs(tau_pos)-dst%threshold_stress(:,of))/prm%tau_0 StressRatio_p = StressRatio** prm%p StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) - + t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength) t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) - + vel = prm%kink_height/(t_n + t_k) - + dot_gamma_pos = dot_gamma_0 * sign(vel,tau_pos) * 0.5_pReal else where significantPositiveTau dot_gamma_pos = 0.0_pReal end where significantPositiveTau - + if (present(ddot_gamma_dtau_pos)) then significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check) dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_0 dtk = -1.0_pReal * t_k / tau_pos - + dvel = -1.0_pReal * prm%kink_height * (dtk + dtn) / (t_n + t_k)**2.0_pReal - + ddot_gamma_dtau_pos = dot_gamma_0 * dvel* 0.5_pReal else where significantPositiveTau2 ddot_gamma_dtau_pos = 0.0_pReal end where significantPositiveTau2 endif - + significantNegativeTau: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check) StressRatio = (abs(tau_neg)-dst%threshold_stress(:,of))/prm%tau_0 StressRatio_p = StressRatio** prm%p StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) - + t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength) t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) - + vel = prm%kink_height/(t_n + t_k) - + dot_gamma_neg = dot_gamma_0 * sign(vel,tau_neg) * 0.5_pReal else where significantNegativeTau dot_gamma_neg = 0.0_pReal end where significantNegativeTau - + if (present(ddot_gamma_dtau_neg)) then significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check) dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_0 dtk = -1.0_pReal * t_k / tau_neg - + dvel = -1.0_pReal * prm%kink_height * (dtk + dtn) / (t_n + t_k)**2.0_pReal - + ddot_gamma_dtau_neg = dot_gamma_0 * dvel * 0.5_pReal else where significantNegativeTau2 ddot_gamma_dtau_neg = 0.0_pReal end where significantNegativeTau2 end if - + end associate end associate diff --git a/src/constitutive_plastic_dislotwin.f90 b/src/constitutive_plastic_dislotwin.f90 index 94e0177a0..152c630ae 100644 --- a/src/constitutive_plastic_dislotwin.f90 +++ b/src/constitutive_plastic_dislotwin.f90 @@ -8,35 +8,17 @@ !> @details to be done !-------------------------------------------------------------------------------------------------- submodule(constitutive) plastic_dislotwin - - real(pReal), parameter :: & + + real(pReal), parameter :: & kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin - - enum, bind(c) - enumerator :: & - undefined_ID, & - rho_mob_ID, & - rho_dip_ID, & - dot_gamma_sl_ID, & - gamma_sl_ID, & - Lambda_sl_ID, & - resolved_stress_slip_ID, & - tau_pass_ID, & - edge_dipole_distance_ID, & - f_tw_ID, & - Lambda_tw_ID, & - resolved_stress_twin_ID, & - tau_hat_tw_ID, & - f_tr_ID - end enum - + type :: tParameters real(pReal) :: & mu, & nu, & D0, & !< prefactor for self-diffusion coefficient Qsd, & !< activation energy for dislocation climb - omega, & !< frequency factor for dislocation climb + omega, & !< frequency factor for dislocation climb D, & !< grain size p_sb, & !< p-exponent in shear band velocity q_sb, & !< q-exponent in shear band velocity @@ -58,8 +40,8 @@ submodule(constitutive) plastic_dislotwin aTol_f_tr, & !< absolute tolerance for integration of trans volume fraction gamma_fcc_hex, & !< Free energy difference between austensite and martensite i_tr, & !< - h !< Stack height of hex nucleus - real(pReal), dimension(:), allocatable :: & + h !< Stack height of hex nucleus + real(pReal), allocatable, dimension(:) :: & rho_mob_0, & !< initial unipolar dislocation density per slip system rho_dip_0, & !< initial dipole dislocation density per slip system b_sl, & !< absolute length of burgers vector [m] for each slip system @@ -79,40 +61,39 @@ submodule(constitutive) plastic_dislotwin s, & !< s-exponent in trans nucleation rate gamma_char, & !< characteristic shear for twins B !< drag coefficient - real(pReal), dimension(:,:), allocatable :: & - h_sl_sl, & !< + real(pReal), allocatable, dimension(:,:) :: & + h_sl_sl, & !< h_sl_tw, & !< - h_tw_tw, & !< - h_sl_tr, & !< - h_tr_tr !< - integer, dimension(:,:), allocatable :: & - fcc_twinNucleationSlipPair ! ToDo: Better name? Is also use for trans - real(pReal), dimension(:,:), allocatable :: & + h_tw_tw, & !< + h_sl_tr, & !< + h_tr_tr, & !< n0_sl, & !< slip system normal forestProjection, & C66 - real(pReal), dimension(:,:,:), allocatable :: & + real(pReal), allocatable, dimension(:,:,:) :: & P_tr, & P_sl, & P_tw, & C66_tw, & C66_tr - integer :: & + integer :: & sum_N_sl, & !< total number of active slip system sum_N_tw, & !< total number of active twin system - sum_N_tr !< total number of active transformation system - integer, dimension(:), allocatable :: & + sum_N_tr !< total number of active transformation system + integer, allocatable, dimension(:) :: & N_sl, & !< number of active slip systems for each family N_tw, & !< number of active twin systems for each family N_tr !< number of active transformation systems for each family - integer(kind(undefined_ID)), dimension(:), allocatable :: & - outputID !< ID of each post result output + integer, allocatable, dimension(:,:) :: & + fcc_twinNucleationSlipPair ! ToDo: Better name? Is also use for trans + character(len=pStringLen), allocatable, dimension(:) :: & + output logical :: & ExtendedDislocations, & !< consider split into partials for climb calculation fccTwinTransNucleation, & !< twinning and transformation models are for fcc dipoleFormation !< flag indicating consideration of dipole formation end type !< container type for internal constitutive parameters - + type :: tDislotwinState real(pReal), dimension(:,:), pointer :: & rho_mob, & @@ -121,7 +102,7 @@ submodule(constitutive) plastic_dislotwin f_tw, & f_tr end type tDislotwinState - + type :: tDislotwinMicrostructure real(pReal), dimension(:,:), allocatable :: & Lambda_sl, & !< mean free path between 2 obstacles seen by a moving dislocation @@ -135,7 +116,7 @@ submodule(constitutive) plastic_dislotwin tau_r_tw, & !< stress to bring partials close together (twin) tau_r_tr !< stress to bring partials close together (trans) end type tDislotwinMicrostructure - + !-------------------------------------------------------------------------------------------------- ! containers for parameters and state type(tParameters), allocatable, dimension(:) :: param @@ -148,7 +129,7 @@ contains !-------------------------------------------------------------------------------------------------- -!> @brief module initialization +!> @brief Perform module initialization. !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- module subroutine plastic_dislotwin_init @@ -156,39 +137,34 @@ module subroutine plastic_dislotwin_init integer :: & Ninstance, & p, i, & - NipcMyPhase, outputSize, & + NipcMyPhase, & sizeState, sizeDotState, & startIndex, endIndex - - integer(kind(undefined_ID)) :: & - outputID - + character(len=pStringLen) :: & extmsg = '' - character(len=pStringLen), dimension(:), allocatable :: & - outputs - - write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_DISLOTWIN_label//' init -+>>>' - - write(6,'(/,a)') ' Ma and Roters, Acta Materialia 52(12):3603–3612, 2004' - write(6,'(a)') ' https://doi.org/10.1016/j.actamat.2004.04.012' - - write(6,'(/,a)') ' Roters et al., Computational Materials Science 39:91–95, 2007' - write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2006.04.014' - - write(6,'(/,a)') ' Wong et al., Acta Materialia 118:140–151, 2016' - write(6,'(a,/)') ' https://doi.org/10.1016/j.actamat.2016.07.032' - + + write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_DISLOTWIN_label//' init -+>>>'; flush(6) + + write(6,'(/,a)') ' Ma and Roters, Acta Materialia 52(12):3603–3612, 2004' + write(6,'(a)') ' https://doi.org/10.1016/j.actamat.2004.04.012' + + write(6,'(/,a)') ' Roters et al., Computational Materials Science 39:91–95, 2007' + write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2006.04.014' + + write(6,'(/,a)') ' Wong et al., Acta Materialia 118:140–151, 2016' + write(6,'(a,/)') ' https://doi.org/10.1016/j.actamat.2016.07.032' + Ninstance = count(phase_plasticity == PLASTICITY_DISLOTWIN_ID) - + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - + allocate(param(Ninstance)) allocate(state(Ninstance)) allocate(dotState(Ninstance)) allocate(dependentState(Ninstance)) - + do p = 1, size(phase_plasticity) if (phase_plasticity(p) /= PLASTICITY_DISLOTWIN_ID) cycle associate(prm => param(phase_plasticityInstance(p)), & @@ -196,17 +172,16 @@ module subroutine plastic_dislotwin_init stt => state(phase_plasticityInstance(p)), & dst => dependentState(phase_plasticityInstance(p)), & config => config_phase(p)) - + prm%aTol_rho = config%getFloat('atol_rho', defaultVal=0.0_pReal) prm%aTol_f_tw = config%getFloat('atol_twinfrac', defaultVal=0.0_pReal) prm%aTol_f_tr = config%getFloat('atol_transfrac', defaultVal=0.0_pReal) - + ! This data is read in already in lattice prm%mu = lattice_mu(p) prm%nu = lattice_nu(p) prm%C66 = lattice_C66(1:6,1:6,p) - - + !-------------------------------------------------------------------------------------------------- ! slip related parameters prm%N_sl = config%getInts('nslip',defaultVal=emptyIntArray) @@ -220,14 +195,14 @@ module subroutine plastic_dislotwin_init prm%forestProjection = lattice_forestProjection_edge(prm%N_sl,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) prm%forestProjection = transpose(prm%forestProjection) - + prm%n0_sl = lattice_slip_normal(prm%N_sl,config%getString('lattice_structure'),& - config%getFloat('c/a',defaultVal=0.0_pReal)) + config%getFloat('c/a',defaultVal=0.0_pReal)) prm%fccTwinTransNucleation = merge(.true., .false., lattice_structure(p) == LATTICE_FCC_ID) & .and. (prm%N_sl(1) == 12) if(prm%fccTwinTransNucleation) & prm%fcc_twinNucleationSlipPair = lattice_fcc_twinNucleationSlipPair - + prm%rho_mob_0 = config%getFloats('rhoedge0', requiredSize=size(prm%N_sl)) prm%rho_dip_0 = config%getFloats('rhoedgedip0',requiredSize=size(prm%N_sl)) prm%v0 = config%getFloats('v0', requiredSize=size(prm%N_sl)) @@ -238,7 +213,7 @@ module subroutine plastic_dislotwin_init prm%q = config%getFloats('q_slip', requiredSize=size(prm%N_sl)) prm%B = config%getFloats('b', requiredSize=size(prm%N_sl), & defaultVal=[(0.0_pReal, i=1,size(prm%N_sl))]) - + prm%tau_0 = config%getFloat('solidsolutionstrength') prm%CEdgeDipMinDistance = config%getFloat('cedgedipmindistance') prm%D0 = config%getFloat('d0') @@ -249,15 +224,15 @@ module subroutine plastic_dislotwin_init prm%SFE_0K = config%getFloat('sfe_0k') prm%dSFE_dT = config%getFloat('dsfe_dt') endif - + ! multiplication factor according to crystal structure (nearest neighbors bcc vs fcc/hex) !@details: Refer: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 prm%omega = config%getFloat('omega', defaultVal = 1000.0_pReal) & * merge(12.0_pReal, & 8.0_pReal, & lattice_structure(p) == LATTICE_FCC_ID .or. lattice_structure(p) == LATTICE_HEX_ID) - - + + ! expand: family => system prm%rho_mob_0 = math_expand(prm%rho_mob_0, prm%N_sl) prm%rho_dip_0 = math_expand(prm%rho_dip_0, prm%N_sl) @@ -268,8 +243,8 @@ module subroutine plastic_dislotwin_init prm%p = math_expand(prm%p, prm%N_sl) prm%q = math_expand(prm%q, prm%N_sl) prm%B = math_expand(prm%B, prm%N_sl) - prm%atomicVolume = math_expand(prm%atomicVolume,prm%N_sl) - + prm%atomicVolume = math_expand(prm%atomicVolume,prm%N_sl) + ! sanity checks if ( prm%D0 <= 0.0_pReal) extmsg = trim(extmsg)//' D0' if ( prm%Qsd <= 0.0_pReal) extmsg = trim(extmsg)//' Qsd' @@ -282,11 +257,11 @@ module subroutine plastic_dislotwin_init if (any(prm%B < 0.0_pReal)) extmsg = trim(extmsg)//' B' if (any(prm%p<=0.0_pReal .or. prm%p>1.0_pReal)) extmsg = trim(extmsg)//' p' if (any(prm%q< 1.0_pReal .or. prm%q>2.0_pReal)) extmsg = trim(extmsg)//' q' - + else slipActive allocate(prm%b_sl(0)) endif slipActive - + !-------------------------------------------------------------------------------------------------- ! twin related parameters prm%N_tw = config%getInts('ntwin', defaultVal=emptyIntArray) @@ -297,31 +272,31 @@ module subroutine plastic_dislotwin_init prm%h_tw_tw = lattice_interaction_TwinByTwin(prm%N_tw,& config%getFloats('interaction_twintwin'), & config%getString('lattice_structure')) - + prm%b_tw = config%getFloats('twinburgers', requiredSize=size(prm%N_tw)) prm%t_tw = config%getFloats('twinsize', requiredSize=size(prm%N_tw)) prm%r = config%getFloats('r_twin', requiredSize=size(prm%N_tw)) - + prm%xc_twin = config%getFloat('xc_twin') prm%L_tw = config%getFloat('l0_twin') prm%i_tw = config%getFloat('cmfptwin') - + prm%gamma_char= lattice_characteristicShear_Twin(prm%N_tw,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) - + prm%C66_tw = lattice_C66_twin(prm%N_tw,prm%C66,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) - + if (.not. prm%fccTwinTransNucleation) then - prm%dot_N_0_tw = config%getFloats('ndot0_twin') + prm%dot_N_0_tw = config%getFloats('ndot0_twin') prm%dot_N_0_tw = math_expand(prm%dot_N_0_tw,prm%N_tw) endif - + ! expand: family => system prm%b_tw = math_expand(prm%b_tw,prm%N_tw) prm%t_tw = math_expand(prm%t_tw,prm%N_tw) prm%r = math_expand(prm%r,prm%N_tw) - + else allocate(prm%gamma_char(0)) allocate(prm%t_tw (0)) @@ -329,7 +304,7 @@ module subroutine plastic_dislotwin_init allocate(prm%r (0)) allocate(prm%h_tw_tw (0,0)) endif - + !-------------------------------------------------------------------------------------------------- ! transformation related parameters prm%N_tr = config%getInts('ntrans', defaultVal=emptyIntArray) @@ -337,29 +312,29 @@ module subroutine plastic_dislotwin_init if (prm%sum_N_tr > 0) then prm%b_tr = config%getFloats('transburgers') prm%b_tr = math_expand(prm%b_tr,prm%N_tr) - + prm%h = config%getFloat('transstackheight', defaultVal=0.0_pReal) ! ToDo: How to handle that??? prm%i_tr = config%getFloat('cmfptrans', defaultVal=0.0_pReal) ! ToDo: How to handle that??? prm%gamma_fcc_hex = config%getFloat('deltag') prm%xc_trans = config%getFloat('xc_trans', defaultVal=0.0_pReal) ! ToDo: How to handle that??? prm%L_tr = config%getFloat('l0_trans') - + prm%h_tr_tr = lattice_interaction_TransByTrans(prm%N_tr,& config%getFloats('interaction_transtrans'), & config%getString('lattice_structure')) - + prm%C66_tr = lattice_C66_trans(prm%N_tr,prm%C66, & config%getString('trans_lattice_structure'), & 0.0_pReal, & config%getFloat('a_bcc', defaultVal=0.0_pReal), & config%getFloat('a_fcc', defaultVal=0.0_pReal)) - + prm%P_tr = lattice_SchmidMatrix_trans(prm%N_tr, & config%getString('trans_lattice_structure'), & 0.0_pReal, & config%getFloat('a_bcc', defaultVal=0.0_pReal), & config%getFloat('a_fcc', defaultVal=0.0_pReal)) - + if (lattice_structure(p) /= LATTICE_fcc_ID) then prm%dot_N_0_tr = config%getFloats('ndot0_trans') prm%dot_N_0_tr = math_expand(prm%dot_N_0_tr,prm%N_tr) @@ -374,59 +349,57 @@ module subroutine plastic_dislotwin_init allocate(prm%s (0)) allocate(prm%h_tr_tr(0,0)) endif - + if (sum(prm%N_tw) > 0 .or. prm%sum_N_tr > 0) then prm%SFE_0K = config%getFloat('sfe_0k') prm%dSFE_dT = config%getFloat('dsfe_dt') prm%V_cs = config%getFloat('vcrossslip') endif - + if (prm%sum_N_sl > 0 .and. prm%sum_N_tw > 0) then prm%h_sl_tw = lattice_interaction_SlipByTwin(prm%N_sl,prm%N_tw,& config%getFloats('interaction_sliptwin'), & config%getString('lattice_structure')) if (prm%fccTwinTransNucleation .and. prm%sum_N_tw > 12) write(6,*) 'mist' ! ToDo: implement better test. The model will fail also if N_tw is [6,6] - endif - - if (prm%sum_N_sl > 0 .and. prm%sum_N_tr > 0) then + endif + + if (prm%sum_N_sl > 0 .and. prm%sum_N_tr > 0) then prm%h_sl_tr = lattice_interaction_SlipByTrans(prm%N_sl,prm%N_tr,& config%getFloats('interaction_sliptrans'), & config%getString('lattice_structure')) if (prm%fccTwinTransNucleation .and. prm%sum_N_tr > 12) write(6,*) 'mist' ! ToDo: implement better test. The model will fail also if N_tr is [6,6] - endif - + endif + !-------------------------------------------------------------------------------------------------- ! shearband related parameters prm%sbVelocity = config%getFloat('shearbandvelocity',defaultVal=0.0_pReal) - if (prm%sbVelocity > 0.0_pReal) then + if (prm%sbVelocity > 0.0_pReal) then prm%sbResistance = config%getFloat('shearbandresistance') prm%sbQedge = config%getFloat('qedgepersbsystem') prm%p_sb = config%getFloat('p_shearband') prm%q_sb = config%getFloat('q_shearband') - + ! sanity checks if (prm%sbResistance < 0.0_pReal) extmsg = trim(extmsg)//' shearbandresistance' if (prm%sbQedge < 0.0_pReal) extmsg = trim(extmsg)//' qedgepersbsystem' if (prm%p_sb <= 0.0_pReal) extmsg = trim(extmsg)//' p_shearband' if (prm%q_sb <= 0.0_pReal) extmsg = trim(extmsg)//' q_shearband' endif - - - + prm%D = config%getFloat('grainsize') - + if (config%keyExists('dipoleformationfactor')) call IO_error(1,ext_msg='use /nodipoleformation/') prm%dipoleformation = .not. config%keyExists('/nodipoleformation/') - - + + !if (Ndot0PerTwinFamily(f,p) < 0.0_pReal) & ! call IO_error(211,el=p,ext_msg='dot_N_0_tw ('//PLASTICITY_DISLOTWIN_label//')') - + if (any(prm%atomicVolume <= 0.0_pReal)) & call IO_error(211,el=p,ext_msg='cAtomicVolume ('//PLASTICITY_DISLOTWIN_label//')') if (prm%sum_N_tw > 0) then if (prm%aTol_rho <= 0.0_pReal) & - call IO_error(211,el=p,ext_msg='aTol_rho ('//PLASTICITY_DISLOTWIN_label//')') + call IO_error(211,el=p,ext_msg='aTol_rho ('//PLASTICITY_DISLOTWIN_label//')') if (prm%aTol_f_tw <= 0.0_pReal) & call IO_error(211,el=p,ext_msg='aTol_f_tw ('//PLASTICITY_DISLOTWIN_label//')') endif @@ -434,50 +407,9 @@ module subroutine plastic_dislotwin_init if (prm%aTol_f_tr <= 0.0_pReal) & call IO_error(211,el=p,ext_msg='aTol_f_tr ('//PLASTICITY_DISLOTWIN_label//')') endif - - outputs = config%getStrings('(output)', defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - do i= 1, size(outputs) - outputID = undefined_ID - select case(outputs(i)) - case ('rho_mob') - outputID = merge(rho_mob_ID,undefined_ID,prm%sum_N_sl > 0) - outputSize = prm%sum_N_sl - case ('rho_dip') - outputID = merge(rho_dip_ID,undefined_ID,prm%sum_N_sl > 0) - outputSize = prm%sum_N_sl - case ('gamma_sl') - outputID = merge(gamma_sl_ID,undefined_ID,prm%sum_N_sl > 0) - outputSize = prm%sum_N_sl - case ('lambda_sl') - outputID = merge(Lambda_sl_ID,undefined_ID,prm%sum_N_sl > 0) - outputSize = prm%sum_N_sl - case ('tau_pass') - outputID= merge(tau_pass_ID,undefined_ID,prm%sum_N_sl > 0) - outputSize = prm%sum_N_sl - - case ('f_tw') - outputID = merge(f_tw_ID,undefined_ID,prm%sum_N_tw >0) - outputSize = prm%sum_N_tw - case ('lambda_tw') - outputID = merge(Lambda_tw_ID,undefined_ID,prm%sum_N_tw >0) - outputSize = prm%sum_N_tw - case ('tau_hat_tw') - outputID = merge(tau_hat_tw_ID,undefined_ID,prm%sum_N_tw >0) - outputSize = prm%sum_N_tw - - case ('f_tr') - outputID = f_tr_ID - outputSize = prm%sum_N_tr - - end select - - if (outputID /= undefined_ID) then - prm%outputID = [prm%outputID, outputID] - endif - - enddo - + + prm%output = config%getStrings('(output)', defaultVal=emptyStringArray) + !-------------------------------------------------------------------------------------------------- ! allocate state arrays NipcMyPhase = count(material_phaseAt == p) * discretization_nIP @@ -485,9 +417,9 @@ module subroutine plastic_dislotwin_init + size(['f_tw']) * prm%sum_N_tw & + size(['f_tr']) * prm%sum_N_tr sizeState = sizeDotState - + call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) - + !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and aTolState startIndex = 1 @@ -496,14 +428,14 @@ module subroutine plastic_dislotwin_init stt%rho_mob= spread(prm%rho_mob_0,2,NipcMyPhase) dot%rho_mob=>plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_rho - + startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl stt%rho_dip=>plasticState(p)%state(startIndex:endIndex,:) stt%rho_dip= spread(prm%rho_dip_0,2,NipcMyPhase) dot%rho_dip=>plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_rho - + startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl stt%gamma_sl=>plasticState(p)%state(startIndex:endIndex,:) @@ -511,66 +443,66 @@ module subroutine plastic_dislotwin_init plasticState(p)%aTolState(startIndex:endIndex) = 1.0e6_pReal !ToDo: better make optional parameter ! global alias plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) - + startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tw stt%f_tw=>plasticState(p)%state(startIndex:endIndex,:) dot%f_tw=>plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_f_tw - + startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tr stt%f_tr=>plasticState(p)%state(startIndex:endIndex,:) dot%f_tr=>plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_f_tr - + allocate(dst%Lambda_sl (prm%sum_N_sl,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_pass (prm%sum_N_sl,NipcMyPhase),source=0.0_pReal) - + allocate(dst%Lambda_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_hat_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_r_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) allocate(dst%V_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) - + allocate(dst%Lambda_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_hat_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_r_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) allocate(dst%V_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) - - + + plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally - + end associate - + enddo end subroutine plastic_dislotwin_init !-------------------------------------------------------------------------------------------------- -!> @brief returns the homogenized elasticity matrix +!> @brief Return the homogenized elasticity matrix. !-------------------------------------------------------------------------------------------------- module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC) - + real(pReal), dimension(6,6) :: & homogenizedC integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element - + integer :: i, & of real(pReal) :: f_unrotated - + of = material_phasememberAt(ipc,ip,el) associate(prm => param(phase_plasticityInstance(material_phaseAt(ipc,el))),& stt => state(phase_plasticityInstance(material_phaseAT(ipc,el)))) - + f_unrotated = 1.0_pReal & - sum(stt%f_tw(1:prm%sum_N_tw,of)) & - sum(stt%f_tr(1:prm%sum_N_tr,of)) - + homogenizedC = f_unrotated * prm%C66 do i=1,prm%sum_N_tw homogenizedC = homogenizedC & @@ -580,23 +512,23 @@ module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC) homogenizedC = homogenizedC & + stt%f_tr(i,of)*prm%C66_tr(1:6,1:6,i) enddo - + end associate - + end function plastic_dislotwin_homogenizedC !-------------------------------------------------------------------------------------------------- -!> @brief calculates plastic velocity gradient and its tangent +!> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) - + real(pReal), dimension(3,3), intent(out) :: Lp real(pReal), dimension(3,3,3,3), intent(out) :: dLp_dMp real(pReal), dimension(3,3), intent(in) :: Mp integer, intent(in) :: instance,of real(pReal), intent(in) :: T - + integer :: i,k,l,m,n real(pReal) :: & f_unrotated,StressRatio_p,& @@ -632,16 +564,16 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) 0, 1,-1, & 0, 1, 1 & ],pReal),[ 3,6]) - + associate(prm => param(instance), stt => state(instance)) - + f_unrotated = 1.0_pReal & - sum(stt%f_tw(1:prm%sum_N_tw,of)) & - sum(stt%f_tr(1:prm%sum_N_tr,of)) - + Lp = 0.0_pReal - dLp_dMp = 0.0_pReal - + dLp_dMp = 0.0_pReal + call kinetics_slip(Mp,T,instance,of,dot_gamma_sl,ddot_gamma_dtau_slip) slipContribution: do i = 1, prm%sum_N_sl Lp = Lp + dot_gamma_sl(i)*prm%P_sl(1:3,1:3,i) @@ -649,37 +581,37 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) enddo slipContribution - + !ToDo: Why do this before shear banding? Lp = Lp * f_unrotated dLp_dMp = dLp_dMp * f_unrotated - + shearBandingContribution: if(dNeq0(prm%sbVelocity)) then - + BoltzmannRatio = prm%sbQedge/(kB*T) call math_eigenValuesVectorsSym(Mp,eigValues,eigVectors,error) - + do i = 1,6 P_sb = 0.5_pReal * math_outer(matmul(eigVectors,sb_sComposition(1:3,i)),& matmul(eigVectors,sb_mComposition(1:3,i))) tau = math_mul33xx33(Mp,P_sb) - + significantShearBandStress: if (abs(tau) > tol_math_check) then StressRatio_p = (abs(tau)/prm%sbResistance)**prm%p_sb dot_gamma_sb = sign(prm%sbVelocity*exp(-BoltzmannRatio*(1-StressRatio_p)**prm%q_sb), tau) ddot_gamma_dtau = abs(dot_gamma_sb)*BoltzmannRatio* prm%p_sb*prm%q_sb/ prm%sbResistance & * (abs(tau)/prm%sbResistance)**(prm%p_sb-1.0_pReal) & * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal) - + Lp = Lp + dot_gamma_sb * P_sb forall (k=1:3,l=1:3,m=1:3,n=1:3) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + ddot_gamma_dtau * P_sb(k,l) * P_sb(m,n) endif significantShearBandStress enddo - + endif shearBandingContribution - + call kinetics_twin(Mp,T,dot_gamma_sl,instance,of,dot_gamma_twin,ddot_gamma_dtau_twin) twinContibution: do i = 1, prm%sum_N_tw Lp = Lp + dot_gamma_twin(i)*prm%P_tw(1:3,1:3,i) * f_unrotated @@ -687,7 +619,7 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + ddot_gamma_dtau_twin(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) * f_unrotated enddo twinContibution - + call kinetics_trans(Mp,T,dot_gamma_sl,instance,of,dot_gamma_tr,ddot_gamma_dtau_trans) transContibution: do i = 1, prm%sum_N_tr Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i) * f_unrotated @@ -695,15 +627,15 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + ddot_gamma_dtau_trans(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i) * f_unrotated enddo transContibution - - + + end associate - + end subroutine plastic_dislotwin_LpAndItsTangent !-------------------------------------------------------------------------------------------------- -!> @brief calculates the rate of change of microstructure +!> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) @@ -720,10 +652,10 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) f_unrotated, & VacancyDiffusion, & rho_dip_distance, & - v_cl, & !< climb velocity + v_cl, & !< climb velocity Gamma, & !< stacking fault energy tau, & - sigma_cl, & !< climb stress + sigma_cl, & !< climb stress b_d !< ratio of burgers vector to stacking fault width real(pReal), dimension(param(instance)%sum_N_sl) :: & dot_rho_dip_formation, & @@ -745,9 +677,9 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) call kinetics_slip(Mp,T,instance,of,dot_gamma_sl) dot%gamma_sl(:,of) = abs(dot_gamma_sl) - + rho_dip_distance_min = prm%CEdgeDipMinDistance*prm%b_sl - + slipState: do i = 1, prm%sum_N_sl tau = math_mul33xx33(Mp,prm%P_sl(1:3,1:3,i)) @@ -771,15 +703,15 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) else !@details: Refer: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i))) - if (prm%ExtendedDislocations) then + if (prm%ExtendedDislocations) then Gamma = prm%SFE_0K + prm%dSFE_dT * T b_d = 24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu)* Gamma/(prm%mu*prm%b_sl(i)) else - b_d = 1.0_pReal + b_d = 1.0_pReal endif v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Qsd/(kB*T)) & * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal) - + dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,of) & / (rho_dip_distance-rho_dip_distance_min(i)) endif @@ -801,12 +733,12 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) dot%f_tr(:,of) = f_unrotated*dot_gamma_tr end associate - + end subroutine plastic_dislotwin_dotState !-------------------------------------------------------------------------------------------------- -!> @brief calculates derived quantities from state +!> @brief Calculate derived quantities from state. !-------------------------------------------------------------------------------------------------- module subroutine plastic_dislotwin_dependentState(T,instance,of) @@ -815,7 +747,7 @@ module subroutine plastic_dislotwin_dependentState(T,instance,of) of real(pReal), intent(in) :: & T - + real(pReal) :: & sumf_twin,Gamma,sumf_trans real(pReal), dimension(param(instance)%sum_N_sl) :: & @@ -830,35 +762,35 @@ module subroutine plastic_dislotwin_dependentState(T,instance,of) f_over_t_tr real(pReal), dimension(:), allocatable :: & x0 - - + + associate(prm => param(instance),& stt => state(instance),& dst => dependentState(instance)) - + sumf_twin = sum(stt%f_tw(1:prm%sum_N_tw,of)) sumf_trans = sum(stt%f_tr(1:prm%sum_N_tr,of)) - + Gamma = prm%SFE_0K + prm%dSFE_dT * T - + !* rescaled volume fraction for topology f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,of)/prm%t_tw ! this is per system ... f_over_t_tr = sumf_trans/prm%t_tr ! but this not ! ToDo ...Physically correct, but naming could be adjusted - + inv_lambda_sl_sl = sqrt(matmul(prm%forestProjection, & stt%rho_mob(:,of)+stt%rho_dip(:,of)))/prm%CLambdaSlip - + if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) & inv_lambda_sl_tw = matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_twin) - + inv_lambda_tw_tw = matmul(prm%h_tw_tw,f_over_t_tw)/(1.0_pReal-sumf_twin) - + if (prm%sum_N_tr > 0 .and. prm%sum_N_sl > 0) & inv_lambda_sl_tr = matmul(prm%h_sl_tr,f_over_t_tr)/(1.0_pReal-sumf_trans) - + inv_lambda_tr_tr = matmul(prm%h_tr_tr,f_over_t_tr)/(1.0_pReal-sumf_trans) - + if ((prm%sum_N_tw > 0) .or. (prm%sum_N_tr > 0)) then ! ToDo: better logic needed here dst%Lambda_sl(:,of) = prm%D & / (1.0_pReal+prm%D*(inv_lambda_sl_sl + inv_lambda_sl_tw + inv_lambda_sl_tr)) @@ -866,13 +798,13 @@ module subroutine plastic_dislotwin_dependentState(T,instance,of) dst%Lambda_sl(:,of) = prm%D & / (1.0_pReal+prm%D*inv_lambda_sl_sl) !!!!!! correct? endif - + dst%Lambda_tw(:,of) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw) dst%Lambda_tr(:,of) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr) - + !* threshold stress for dislocation motion dst%tau_pass(:,of) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,of)+stt%rho_dip(:,of))) - + !* threshold stress for growing twin/martensite if(prm%sum_N_tw == prm%sum_N_sl) & dst%tau_hat_tw(:,of) = Gamma/(3.0_pReal*prm%b_tw) & @@ -880,66 +812,67 @@ module subroutine plastic_dislotwin_dependentState(T,instance,of) if(prm%sum_N_tr == prm%sum_N_sl) & dst%tau_hat_tr(:,of) = Gamma/(3.0_pReal*prm%b_tr) & + 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip burgers here correct? - + prm%h*prm%gamma_fcc_hex/ (3.0_pReal*prm%b_tr) - + + prm%h*prm%gamma_fcc_hex/ (3.0_pReal*prm%b_tr) + dst%V_tw(:,of) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,of)**2.0_pReal dst%V_tr(:,of) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,of)**2.0_pReal - - + + x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the burgers vector for slip and is the same for twin and trans dst%tau_r_tw(:,of) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%xc_twin)+cos(pi/3.0_pReal)/x0) - + x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the burgers vector for slip dst%tau_r_tr(:,of) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%xc_trans)+cos(pi/3.0_pReal)/x0) - + end associate end subroutine plastic_dislotwin_dependentState !-------------------------------------------------------------------------------------------------- -!> @brief writes results to HDF5 output file +!> @brief Write results to HDF5 output file. !-------------------------------------------------------------------------------------------------- module subroutine plastic_dislotwin_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group + integer :: o - associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) + associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) - case (rho_mob_ID) - call results_writeDataset(group,stt%rho_mob,'rho_mob',& - 'mobile dislocation density','1/m²') - case (rho_dip_ID) - call results_writeDataset(group,stt%rho_dip,'rho_dip',& - 'dislocation dipole density''1/m²') - case (gamma_sl_ID) - call results_writeDataset(group,stt%gamma_sl,'gamma_sl',& - 'plastic shear','1') - case (Lambda_sl_ID) - call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',& - 'mean free path for slip','m') - case (tau_pass_ID) - call results_writeDataset(group,dst%tau_pass,'tau_pass',& - 'passing stress for slip','Pa') + case('rho_mob') + if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_mob,'rho_mob',& + 'mobile dislocation density','1/m²') + case('rho_dip') + if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_dip,'rho_dip',& + 'dislocation dipole density''1/m²') + case('gamma_sl') + if(prm%sum_N_sl>0) call results_writeDataset(group,stt%gamma_sl,'gamma_sl',& + 'plastic shear','1') + case('lambda_sl') + if(prm%sum_N_sl>0) call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',& + 'mean free path for slip','m') + case('tau_pass') + if(prm%sum_N_sl>0) call results_writeDataset(group,dst%tau_pass,'tau_pass',& + 'passing stress for slip','Pa') - case (f_tw_ID) - call results_writeDataset(group,stt%f_tw,'f_tw',& - 'twinned volume fraction','m³/m³') - case (Lambda_tw_ID) - call results_writeDataset(group,dst%Lambda_tw,'Lambda_tw',& - 'mean free path for twinning','m') - case (tau_hat_tw_ID) - call results_writeDataset(group,dst%tau_hat_tw,'tau_hat_tw',& - 'threshold stress for twinning','Pa') + case('f_tw') + if(prm%sum_N_tw>0) call results_writeDataset(group,stt%f_tw,'f_tw',& + 'twinned volume fraction','m³/m³') + case('lambda_tw') + if(prm%sum_N_tw>0) call results_writeDataset(group,dst%Lambda_tw,'Lambda_tw',& + 'mean free path for twinning','m') + case('tau_hat_tw') + if(prm%sum_N_tw>0) call results_writeDataset(group,dst%tau_hat_tw,'tau_hat_tw',& + 'threshold stress for twinning','Pa') + + case('f_tr') + if(prm%sum_N_tr>0) call results_writeDataset(group,stt%f_tr,'f_tr',& + 'martensite volume fraction','m³/m³') - case (f_tr_ID) - call results_writeDataset(group,stt%f_tr,'f_tr',& - 'martensite volume fraction','m³/m³') - end select enddo outputsLoop end associate @@ -948,8 +881,8 @@ end subroutine plastic_dislotwin_results !-------------------------------------------------------------------------------------------------- -!> @brief Shear rates on slip systems, their derivatives with respect to resolved stress and the -! resolved stresss +!> @brief Calculate shear rates on slip systems, their derivatives with respect to resolved +! stress, and the resolved stress. !> @details Derivatives and resolved stress are calculated only optionally. ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end @@ -964,7 +897,7 @@ pure subroutine kinetics_slip(Mp,T,instance,of, & integer, intent(in) :: & instance, & of - + real(pReal), dimension(param(instance)%sum_N_sl), intent(out) :: & dot_gamma_sl real(pReal), dimension(param(instance)%sum_N_sl), optional, intent(out) :: & @@ -972,7 +905,7 @@ pure subroutine kinetics_slip(Mp,T,instance,of, & tau_slip real(pReal), dimension(param(instance)%sum_N_sl) :: & ddot_gamma_dtau - + real(pReal), dimension(param(instance)%sum_N_sl) :: & tau, & stressRatio, & @@ -984,25 +917,25 @@ pure subroutine kinetics_slip(Mp,T,instance,of, & dV_run_inverse_dTau, & dV_dTau, & tau_eff !< effective resolved stress - integer :: i - + integer :: i + associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) - + do i = 1, prm%sum_N_sl tau(i) = math_mul33xx33(Mp,prm%P_sl(1:3,1:3,i)) enddo - + tau_eff = abs(tau)-dst%tau_pass(:,of) - + significantStress: where(tau_eff > tol_math_check) stressRatio = tau_eff/prm%tau_0 StressRatio_p = stressRatio** prm%p BoltzmannRatio = prm%Delta_F/(kB*T) v_wait_inverse = prm%v0**(-1.0_pReal) * exp(BoltzmannRatio*(1.0_pReal-StressRatio_p)** prm%q) v_run_inverse = prm%B/(tau_eff*prm%b_sl) - + dot_gamma_sl = sign(stt%rho_mob(:,of)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau) - + dV_wait_inverse_dTau = -1.0_pReal * v_wait_inverse * prm%p * prm%q * BoltzmannRatio & * (stressRatio**(prm%p-1.0_pReal)) & * (1.0_pReal-StressRatio_p)**(prm%q-1.0_pReal) & @@ -1015,17 +948,21 @@ pure subroutine kinetics_slip(Mp,T,instance,of, & dot_gamma_sl = 0.0_pReal ddot_gamma_dtau = 0.0_pReal end where significantStress - + end associate - + if(present(ddot_gamma_dtau_slip)) ddot_gamma_dtau_slip = ddot_gamma_dtau if(present(tau_slip)) tau_slip = tau - + end subroutine kinetics_slip !-------------------------------------------------------------------------------------------------- -!> @brief calculates shear rates on twin systems +!> @brief Calculate shear rates on twin systems and their derivatives with respect to resolved +! stress. +!> @details Derivatives are calculated only optionally. +! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to +! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,& dot_gamma_twin,ddot_gamma_dtau_twin) @@ -1039,22 +976,22 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,& of real(pReal), dimension(param(instance)%sum_N_sl), intent(in) :: & dot_gamma_sl - + real(pReal), dimension(param(instance)%sum_N_tw), intent(out) :: & dot_gamma_twin real(pReal), dimension(param(instance)%sum_N_tw), optional, intent(out) :: & ddot_gamma_dtau_twin - + real, dimension(param(instance)%sum_N_tw) :: & tau, & Ndot0, & stressRatio_r, & ddot_gamma_dtau - + integer :: i,s1,s2 - + associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) - + do i = 1, prm%sum_N_tw tau(i) = math_mul33xx33(Mp,prm%P_tw(1:3,1:3,i)) isFCC: if (prm%fccTwinTransNucleation) then @@ -1072,7 +1009,7 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,& Ndot0=prm%dot_N_0_tw(i) endif isFCC enddo - + significantStress: where(tau > tol_math_check) StressRatio_r = (dst%tau_hat_tw(:,of)/tau)**prm%r dot_gamma_twin = prm%gamma_char * dst%V_tw(:,of) * Ndot0*exp(-StressRatio_r) @@ -1081,16 +1018,20 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,& dot_gamma_twin = 0.0_pReal ddot_gamma_dtau = 0.0_pReal end where significantStress - + end associate - + if(present(ddot_gamma_dtau_twin)) ddot_gamma_dtau_twin = ddot_gamma_dtau end subroutine kinetics_twin !-------------------------------------------------------------------------------------------------- -!> @brief calculates shear rates on twin systems +!> @brief Calculate shear rates on transformation systems and their derivatives with respect to +! resolved stress. +!> @details Derivatives are calculated only optionally. +! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to +! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,& dot_gamma_tr,ddot_gamma_dtau_trans) @@ -1104,21 +1045,21 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,& of real(pReal), dimension(param(instance)%sum_N_sl), intent(in) :: & dot_gamma_sl - + real(pReal), dimension(param(instance)%sum_N_tr), intent(out) :: & dot_gamma_tr real(pReal), dimension(param(instance)%sum_N_tr), optional, intent(out) :: & ddot_gamma_dtau_trans - + real, dimension(param(instance)%sum_N_tr) :: & tau, & Ndot0, & stressRatio_s, & ddot_gamma_dtau - + integer :: i,s1,s2 associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) - + do i = 1, prm%sum_N_tr tau(i) = math_mul33xx33(Mp,prm%P_tr(1:3,1:3,i)) isFCC: if (prm%fccTwinTransNucleation) then @@ -1136,7 +1077,7 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,& Ndot0=prm%dot_N_0_tr(i) endif isFCC enddo - + significantStress: where(tau > tol_math_check) StressRatio_s = (dst%tau_hat_tr(:,of)/tau)**prm%s dot_gamma_tr = dst%V_tr(:,of) * Ndot0*exp(-StressRatio_s) @@ -1145,9 +1086,9 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,& dot_gamma_tr = 0.0_pReal ddot_gamma_dtau = 0.0_pReal end where significantStress - + end associate - + if(present(ddot_gamma_dtau_trans)) ddot_gamma_dtau_trans = ddot_gamma_dtau end subroutine kinetics_trans diff --git a/src/constitutive_plastic_isotropic.f90 b/src/constitutive_plastic_isotropic.f90 index ef8a823d4..beed2112a 100644 --- a/src/constitutive_plastic_isotropic.f90 +++ b/src/constitutive_plastic_isotropic.f90 @@ -8,14 +8,7 @@ !! untextured polycrystal !-------------------------------------------------------------------------------------------------- submodule(constitutive) plastic_isotropic - - enum, bind(c) - enumerator :: & - undefined_ID, & - xi_ID, & - dot_gamma_ID - end enum - + type :: tParameters real(pReal) :: & M, & !< Taylor factor @@ -34,12 +27,12 @@ submodule(constitutive) plastic_isotropic aTol_gamma integer :: & of_debug = 0 - integer(kind(undefined_ID)), allocatable, dimension(:) :: & - outputID logical :: & dilatation + character(len=pStringLen), allocatable, dimension(:) :: & + output end type tParameters - + type :: tIsotropicState real(pReal), pointer, dimension(:) :: & xi, & @@ -56,50 +49,45 @@ submodule(constitutive) plastic_isotropic contains !-------------------------------------------------------------------------------------------------- -!> @brief module initialization +!> @brief Perform module initialization. !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- module subroutine plastic_isotropic_init integer :: & Ninstance, & - p, i, & + p, & NipcMyPhase, & sizeState, sizeDotState - - integer(kind(undefined_ID)) :: & - outputID - + character(len=pStringLen) :: & extmsg = '' - character(len=pStringLen), dimension(:), allocatable :: & - outputs - - write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_ISOTROPIC_label//' init -+>>>' - - write(6,'(/,a)') ' Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018' - write(6,'(a)') ' https://doi.org/10.1016/j.scriptamat.2017.09.047' - + + write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_ISOTROPIC_label//' init -+>>>'; flush(6) + + write(6,'(/,a)') ' Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018' + write(6,'(a)') ' https://doi.org/10.1016/j.scriptamat.2017.09.047' + Ninstance = count(phase_plasticity == PLASTICITY_ISOTROPIC_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - + allocate(param(Ninstance)) allocate(state(Ninstance)) allocate(dotState(Ninstance)) - + do p = 1, size(phase_plasticity) if (phase_plasticity(p) /= PLASTICITY_ISOTROPIC_ID) cycle associate(prm => param(phase_plasticityInstance(p)), & dot => dotState(phase_plasticityInstance(p)), & stt => state(phase_plasticityInstance(p)), & config => config_phase(p)) - + #ifdef DEBUG if (p==material_phaseAt(debug_g,debug_e)) & prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e) #endif - + prm%xi_0 = config%getFloat('tau0') prm%xi_inf = config%getFloat('tausat') prm%dot_gamma_0 = config%getFloat('gdot0') @@ -114,9 +102,9 @@ module subroutine plastic_isotropic_init prm%a = config%getFloat('a') prm%aTol_xi = config%getFloat('atol_flowstress',defaultVal=1.0_pReal) prm%aTol_gamma = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal) - + prm%dilatation = config%keyExists('/dilatation/') - + !-------------------------------------------------------------------------------------------------- ! sanity checks extmsg = '' @@ -128,79 +116,62 @@ module subroutine plastic_isotropic_init if (prm%M <= 0.0_pReal) extmsg = trim(extmsg)//' m' if (prm%aTol_xi <= 0.0_pReal) extmsg = trim(extmsg)//' atol_xi' if (prm%aTol_gamma <= 0.0_pReal) extmsg = trim(extmsg)//' atol_shear' - + !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range if (extmsg /= '') & call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_ISOTROPIC_label//')') - + !-------------------------------------------------------------------------------------------------- ! output pararameters - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - do i=1, size(outputs) - outputID = undefined_ID - select case(outputs(i)) - - case ('flowstress') - outputID = xi_ID - case ('strainrate') - outputID = dot_gamma_ID - - end select - - if (outputID /= undefined_ID) then - prm%outputID = [prm%outputID, outputID] - endif - - enddo - + prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) + !-------------------------------------------------------------------------------------------------- ! allocate state arrays NipcMyPhase = count(material_phaseAt == p) * discretization_nIP sizeDotState = size(['xi ','accumulated_shear']) sizeState = sizeDotState - + call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) - + !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and aTolState stt%xi => plasticState(p)%state (1,:) stt%xi = prm%xi_0 dot%xi => plasticState(p)%dotState(1,:) plasticState(p)%aTolState(1) = prm%aTol_xi - + stt%gamma => plasticState(p)%state (2,:) dot%gamma => plasticState(p)%dotState(2,:) plasticState(p)%aTolState(2) = prm%aTol_gamma ! global alias plasticState(p)%slipRate => plasticState(p)%dotState(2:2,:) - + plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally - + end associate - + enddo end subroutine plastic_isotropic_init !-------------------------------------------------------------------------------------------------- -!> @brief calculates plastic velocity gradient and its tangent +!> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) - + real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & instance, & of - + real(pReal), dimension(3,3) :: & Mp_dev !< deviatoric part of the Mandel stress real(pReal) :: & @@ -209,16 +180,16 @@ module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) squarenorm_Mp_dev !< square of the norm of the deviatoric part of the Mandel stress integer :: & k, l, m, n - + associate(prm => param(instance), stt => state(instance)) - + Mp_dev = math_deviatoric33(Mp) squarenorm_Mp_dev = math_mul33xx33(Mp_dev,Mp_dev) norm_Mp_dev = sqrt(squarenorm_Mp_dev) - + if (norm_Mp_dev > 0.0_pReal) then dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(of))) **prm%n - + Lp = dot_gamma/prm%M * Mp_dev/norm_Mp_dev #ifdef DEBUG if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 & @@ -240,42 +211,42 @@ module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) Lp = 0.0_pReal dLp_dMp = 0.0_pReal end if - + end associate end subroutine plastic_isotropic_LpAndItsTangent !-------------------------------------------------------------------------------------------------- -!> @brief calculates plastic velocity gradient and its tangent +!> @brief Calculate inelastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) - + real(pReal), dimension(3,3), intent(out) :: & Li !< inleastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLi_dMi !< derivative of Li with respect to Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & - Mi !< Mandel stress + Mi !< Mandel stress integer, intent(in) :: & instance, & of - + real(pReal) :: & tr !< trace of spherical part of Mandel stress (= 3 x pressure) integer :: & k, l, m, n - + associate(prm => param(instance), stt => state(instance)) - + tr=math_trace33(math_spherical33(Mi)) if (prm%dilatation .and. abs(tr) > 0.0_pReal) then ! no stress or J2 plasticity --> Li and its derivative are zero Li = math_I3 & * prm%dot_gamma_0/prm%M * (3.0_pReal*prm%M*stt%xi(of))**(-prm%n) & * tr * abs(tr)**(prm%n-1.0_pReal) - + #ifdef DEBUG if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 & .and. (of == prm%of_debug .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0)) then @@ -292,38 +263,38 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) Li = 0.0_pReal dLi_dMi = 0.0_pReal endif - + end associate end subroutine plastic_isotropic_LiAndItsTangent !-------------------------------------------------------------------------------------------------- -!> @brief calculates the rate of change of microstructure +!> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- module subroutine plastic_isotropic_dotState(Mp,instance,of) - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & instance, & of - + real(pReal) :: & dot_gamma, & !< strainrate xi_inf_star, & !< saturation xi norm_Mp !< norm of the (deviatoric) Mandel stress - + associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) - + if (prm%dilatation) then norm_Mp = sqrt(math_mul33xx33(Mp,Mp)) else norm_Mp = sqrt(math_mul33xx33(math_deviatoric33(Mp),math_deviatoric33(Mp))) endif - + dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp /(prm%M*stt%xi(of))) **prm%n - + if (dot_gamma > 1e-12_pReal) then if (dEq0(prm%c_1)) then xi_inf_star = prm%xi_inf @@ -339,28 +310,28 @@ module subroutine plastic_isotropic_dotState(Mp,instance,of) else dot%xi(of) = 0.0_pReal endif - + dot%gamma(of) = dot_gamma ! ToDo: not really used - + end associate end subroutine plastic_isotropic_dotState !-------------------------------------------------------------------------------------------------- -!> @brief writes results to HDF5 output file +!> @brief Write results to HDF5 output file. !-------------------------------------------------------------------------------------------------- module subroutine plastic_isotropic_results(instance,group) - integer, intent(in) :: instance + integer, intent(in) :: instance character(len=*), intent(in) :: group - + integer :: o associate(prm => param(instance), stt => state(instance)) - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - case (xi_ID) + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case ('flowstress') ! ToDo: should be 'xi' call results_writeDataset(group,stt%xi,'xi','resistance against plastic flow','Pa') end select enddo outputsLoop diff --git a/src/constitutive_plastic_kinehardening.f90 b/src/constitutive_plastic_kinehardening.f90 index 11ce7c494..65751989c 100644 --- a/src/constitutive_plastic_kinehardening.f90 +++ b/src/constitutive_plastic_kinehardening.f90 @@ -7,26 +7,13 @@ !-------------------------------------------------------------------------------------------------- submodule(constitutive) plastic_kinehardening - enum, bind(c) - enumerator :: & - undefined_ID, & - crss_ID, & !< critical resolved stress - crss_back_ID, & !< critical resolved back stress - sense_ID, & !< sense of acting shear stress (-1 or +1) - chi0_ID, & !< backstress at last switch of stress sense (positive?) - gamma0_ID, & !< accumulated shear at last switch of stress sense (at current switch?) - accshear_ID, & - shearrate_ID, & - resolvedstress_ID - end enum - type :: tParameters real(pReal) :: & gdot0, & !< reference shear strain rate for slip n, & !< stress exponent for slip aTolResistance, & aTolShear - real(pReal), allocatable, dimension(:) :: & + real(pReal), allocatable, dimension(:) :: & crss0, & !< initial critical shear stress for slip theta0, & !< initial hardening rate of forward stress for each slip theta1, & !< asymptotic hardening rate of forward stress for each slip @@ -35,21 +22,21 @@ submodule(constitutive) plastic_kinehardening tau1, & tau1_b, & nonSchmidCoeff - real(pReal), allocatable, dimension(:,:) :: & + real(pReal), allocatable, dimension(:,:) :: & interaction_slipslip !< slip resistance from slip activity - real(pReal), allocatable, dimension(:,:,:) :: & + real(pReal), allocatable, dimension(:,:,:) :: & Schmid, & nonSchmid_pos, & nonSchmid_neg integer :: & totalNslip, & !< total number of active slip system of_debug = 0 - integer, allocatable, dimension(:) :: & + integer, allocatable, dimension(:) :: & Nslip !< number of active slip systems for each family - integer(kind(undefined_ID)), allocatable, dimension(:) :: & - outputID !< ID of each post result output + character(len=pStringLen), allocatable, dimension(:) :: & + output end type tParameters - + type :: tKinehardeningState real(pReal), pointer, dimension(:,:) :: & !< vectors along NipcMyInstance crss, & !< critical resolved stress @@ -59,7 +46,7 @@ submodule(constitutive) plastic_kinehardening gamma0, & !< accumulated shear at last switch of stress sense accshear !< accumulated (absolute) shear end type tKinehardeningState - + !-------------------------------------------------------------------------------------------------- ! containers for parameters and state type(tParameters), allocatable, dimension(:) :: param @@ -72,37 +59,32 @@ contains !-------------------------------------------------------------------------------------------------- -!> @brief module initialization +!> @brief Perform module initialization. !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- module subroutine plastic_kinehardening_init integer :: & Ninstance, & - p, i, o, & + p, o, & NipcMyPhase, & sizeState, sizeDeltaState, sizeDotState, & startIndex, endIndex - - integer(kind(undefined_ID)) :: & - outputID - + character(len=pStringLen) :: & extmsg = '' - character(len=pStringLen), dimension(:), allocatable :: & - outputs - - write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_KINEHARDENING_label//' init -+>>>' - + + write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_KINEHARDENING_label//' init -+>>>'; flush(6) + Ninstance = count(phase_plasticity == PLASTICITY_KINEHARDENING_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - + allocate(param(Ninstance)) allocate(state(Ninstance)) allocate(dotState(Ninstance)) allocate(deltaState(Ninstance)) - + do p = 1, size(phase_plasticityInstance) if (phase_plasticity(p) /= PLASTICITY_KINEHARDENING_ID) cycle associate(prm => param(phase_plasticityInstance(p)), & @@ -110,22 +92,22 @@ module subroutine plastic_kinehardening_init dlt => deltaState(phase_plasticityInstance(p)), & stt => state(phase_plasticityInstance(p)),& config => config_phase(p)) - + #ifdef DEBUG if (p==material_phaseAt(debug_g,debug_e)) then - prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e) + prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e) endif #endif - + !-------------------------------------------------------------------------------------------------- ! optional parameters that need to be defined prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal) prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal) - + ! sanity checks if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//' aTolresistance' if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//' aTolShear' - + !-------------------------------------------------------------------------------------------------- ! slip related parameters prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) @@ -133,7 +115,7 @@ module subroutine plastic_kinehardening_init slipActive: if (prm%totalNslip > 0) then prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) - + if(trim(config%getString('lattice_structure')) == 'bcc') then prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',& defaultVal = emptyRealArray) @@ -146,7 +128,7 @@ module subroutine plastic_kinehardening_init prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, & config%getFloats('interaction_slipslip'), & config%getString('lattice_structure')) - + prm%crss0 = config%getFloats('crss0', requiredSize=size(prm%Nslip)) prm%tau1 = config%getFloats('tau1', requiredSize=size(prm%Nslip)) prm%tau1_b = config%getFloats('tau1_b', requiredSize=size(prm%Nslip)) @@ -154,10 +136,10 @@ module subroutine plastic_kinehardening_init prm%theta1 = config%getFloats('theta1', requiredSize=size(prm%Nslip)) prm%theta0_b = config%getFloats('theta0_b', requiredSize=size(prm%Nslip)) prm%theta1_b = config%getFloats('theta1_b', requiredSize=size(prm%Nslip)) - - prm%gdot0 = config%getFloat('gdot0') - prm%n = config%getFloat('n_slip') - + + prm%gdot0 = config%getFloat('gdot0') + prm%n = config%getFloat('n_slip') + ! expand: family => system prm%crss0 = math_expand(prm%crss0, prm%Nslip) prm%tau1 = math_expand(prm%tau1, prm%Nslip) @@ -166,9 +148,9 @@ module subroutine plastic_kinehardening_init prm%theta1 = math_expand(prm%theta1, prm%Nslip) prm%theta0_b = math_expand(prm%theta0_b,prm%Nslip) prm%theta1_b = math_expand(prm%theta1_b,prm%Nslip) - - - + + + !-------------------------------------------------------------------------------------------------- ! sanity checks if ( prm%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0' @@ -176,58 +158,29 @@ module subroutine plastic_kinehardening_init if (any(prm%crss0 <= 0.0_pReal)) extmsg = trim(extmsg)//' crss0' if (any(prm%tau1 <= 0.0_pReal)) extmsg = trim(extmsg)//' tau1' if (any(prm%tau1_b <= 0.0_pReal)) extmsg = trim(extmsg)//' tau1_b' - + !ToDo: Any sensible checks for theta? - + endif slipActive - + !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range if (extmsg /= '') & call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_KINEHARDENING_label//')') - + !-------------------------------------------------------------------------------------------------- ! output pararameters - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - do i=1, size(outputs) - outputID = undefined_ID - select case(outputs(i)) - - case ('resistance') - outputID = merge(crss_ID,undefined_ID,prm%totalNslip>0) - case ('accumulatedshear') - outputID = merge(accshear_ID,undefined_ID,prm%totalNslip>0) - case ('shearrate') - outputID = merge(shearrate_ID,undefined_ID,prm%totalNslip>0) - case ('resolvedstress') - outputID = merge(resolvedstress_ID,undefined_ID,prm%totalNslip>0) - case ('backstress') - outputID = merge(crss_back_ID,undefined_ID,prm%totalNslip>0) - case ('sense') - outputID = merge(sense_ID,undefined_ID,prm%totalNslip>0) - case ('chi0') - outputID = merge(chi0_ID,undefined_ID,prm%totalNslip>0) - case ('gamma0') - outputID = merge(gamma0_ID,undefined_ID,prm%totalNslip>0) - - end select - - if (outputID /= undefined_ID) then - prm%outputID = [prm%outputID , outputID] - endif - - enddo - + prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) + !-------------------------------------------------------------------------------------------------- ! allocate state arrays NipcMyPhase = count(material_phaseAt == p) * discretization_nIP sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%totalNslip sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%totalNslip sizeState = sizeDotState + sizeDeltaState - + call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) - + !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and aTolState startIndex = 1 @@ -236,13 +189,13 @@ module subroutine plastic_kinehardening_init stt%crss = spread(prm%crss0, 2, NipcMyPhase) dot%crss => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance - + startIndex = endIndex + 1 endIndex = endIndex + prm%totalNslip stt%crss_back => plasticState(p)%state (startIndex:endIndex,:) dot%crss_back => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance - + startIndex = endIndex + 1 endIndex = endIndex + prm%totalNslip stt%accshear => plasticState(p)%state (startIndex:endIndex,:) @@ -250,34 +203,34 @@ module subroutine plastic_kinehardening_init plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear ! global alias plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) - + o = plasticState(p)%offsetDeltaState startIndex = endIndex + 1 endIndex = endIndex + prm%totalNslip stt%sense => plasticState(p)%state (startIndex :endIndex ,:) dlt%sense => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) - + startIndex = endIndex + 1 endIndex = endIndex + prm%totalNslip stt%chi0 => plasticState(p)%state (startIndex :endIndex ,:) dlt%chi0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) - + startIndex = endIndex + 1 endIndex = endIndex + prm%totalNslip stt%gamma0 => plasticState(p)%state (startIndex :endIndex ,:) dlt%gamma0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) - + plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally - + end associate - + enddo end subroutine plastic_kinehardening_init !-------------------------------------------------------------------------------------------------- -!> @brief calculates plastic velocity gradient and its tangent +!> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) @@ -285,26 +238,26 @@ pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,insta Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & instance, & of - + integer :: & i,k,l,m,n real(pReal), dimension(param(instance)%totalNslip) :: & gdot_pos,gdot_neg, & dgdot_dtau_pos,dgdot_dtau_neg - + Lp = 0.0_pReal dLp_dMp = 0.0_pReal - + associate(prm => param(instance)) - + call kinetics(Mp,instance,of,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) - + do i = 1, prm%totalNslip Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%Schmid(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -312,14 +265,14 @@ pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,insta + dgdot_dtau_pos(i) * prm%Schmid(k,l,i) * prm%nonSchmid_pos(m,n,i) & + dgdot_dtau_neg(i) * prm%Schmid(k,l,i) * prm%nonSchmid_neg(m,n,i) enddo - + end associate end subroutine plastic_kinehardening_LpAndItsTangent !-------------------------------------------------------------------------------------------------- -!> @brief calculates the rate of change of microstructure +!> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- module subroutine plastic_kinehardening_dotState(Mp,instance,of) @@ -328,40 +281,40 @@ module subroutine plastic_kinehardening_dotState(Mp,instance,of) integer, intent(in) :: & instance, & of - + real(pReal) :: & sumGamma real(pReal), dimension(param(instance)%totalNslip) :: & gdot_pos,gdot_neg - - + + associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) - + call kinetics(Mp,instance,of,gdot_pos,gdot_neg) dot%accshear(:,of) = abs(gdot_pos+gdot_neg) sumGamma = sum(stt%accshear(:,of)) - - + + dot%crss(:,of) = matmul(prm%interaction_SlipSlip,dot%accshear(:,of)) & * ( prm%theta1 & + (prm%theta0 - prm%theta1 + prm%theta0*prm%theta1*sumGamma/prm%tau1) & * exp(-sumGamma*prm%theta0/prm%tau1) & ) - + dot%crss_back(:,of) = stt%sense(:,of)*dot%accshear(:,of) * & ( prm%theta1_b + & (prm%theta0_b - prm%theta1_b & + prm%theta0_b*prm%theta1_b/(prm%tau1_b+stt%chi0(:,of))*(stt%accshear(:,of)-stt%gamma0(:,of))& ) *exp(-(stt%accshear(:,of)-stt%gamma0(:,of)) *prm%theta0_b/(prm%tau1_b+stt%chi0(:,of))) & ) - + end associate end subroutine plastic_kinehardening_dotState !-------------------------------------------------------------------------------------------------- -!> @brief calculates (instantaneous) incremental change of microstructure +!> @brief Calculate (instantaneous) incremental change of microstructure. !-------------------------------------------------------------------------------------------------- module subroutine plastic_kinehardening_deltaState(Mp,instance,of) @@ -370,18 +323,18 @@ module subroutine plastic_kinehardening_deltaState(Mp,instance,of) integer, intent(in) :: & instance, & of - + real(pReal), dimension(param(instance)%totalNslip) :: & gdot_pos,gdot_neg, & sense - + associate(prm => param(instance), stt => state(instance), dlt => deltaState(instance)) - + call kinetics(Mp,instance,of,gdot_pos,gdot_neg) sense = merge(state(instance)%sense(:,of), & ! keep existing... sign(1.0_pReal,gdot_pos+gdot_neg), & ! ...or have a defined dEq0(gdot_pos+gdot_neg,1e-10_pReal)) ! current sense of shear direction - + #ifdef DEBUG if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 & .and. (of == prm%of_debug & @@ -390,7 +343,7 @@ module subroutine plastic_kinehardening_deltaState(Mp,instance,of) write(6,*) sense,state(instance)%sense(:,of) endif #endif - + !-------------------------------------------------------------------------------------------------- ! switch in sense of shear? where(dNeq(sense,stt%sense(:,of),0.1_pReal)) @@ -402,45 +355,43 @@ module subroutine plastic_kinehardening_deltaState(Mp,instance,of) dlt%chi0 (:,of) = 0.0_pReal dlt%gamma0(:,of) = 0.0_pReal end where - + end associate end subroutine plastic_kinehardening_deltaState !-------------------------------------------------------------------------------------------------- -!> @brief writes results to HDF5 output file +!> @brief Write results to HDF5 output file. !-------------------------------------------------------------------------------------------------- module subroutine plastic_kinehardening_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group + integer :: o associate(prm => param(instance), stt => state(instance)) - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - case (crss_ID) - call results_writeDataset(group,stt%crss,'xi_sl', & - 'resistance against plastic slip','Pa') - - case(crss_back_ID) - call results_writeDataset(group,stt%crss_back,'tau_back', & - 'back stress against plastic slip','Pa') - - case (sense_ID) - call results_writeDataset(group,stt%sense,'sense_of_shear','tbd','1') - - case (chi0_ID) - call results_writeDataset(group,stt%chi0,'chi0','tbd','Pa') - - case (gamma0_ID) - call results_writeDataset(group,stt%gamma0,'gamma0','tbd','1') - - case (accshear_ID) - call results_writeDataset(group,stt%accshear,'gamma_sl', & - 'plastic shear','1') - + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case('resistance') + if(prm%totalNslip>0) call results_writeDataset(group,stt%crss,'xi_sl', & + 'resistance against plastic slip','Pa') + case('backstress') ! ToDo: should be 'tau_back' + if(prm%totalNslip>0) call results_writeDataset(group,stt%crss_back,'tau_back', & + 'back stress against plastic slip','Pa') + case ('sense') + if(prm%totalNslip>0) call results_writeDataset(group,stt%sense,'sense_of_shear', & + 'tbd','1') + case ('chi0') + if(prm%totalNslip>0) call results_writeDataset(group,stt%chi0,'chi0', & + 'tbd','Pa') + case ('gamma0') + if(prm%totalNslip>0) call results_writeDataset(group,stt%gamma0,'gamma0', & + 'tbd','1') + case ('accumulatedshear') + if(prm%totalNslip>0) call results_writeDataset(group,stt%accshear,'gamma_sl', & + 'plastic shear','1') end select enddo outputsLoop end associate @@ -449,10 +400,11 @@ end subroutine plastic_kinehardening_results !-------------------------------------------------------------------------------------------------- -!> @brief calculates shear rates on slip systems and derivatives with respect to resolved stress -!> @details: Shear rates are calculated only optionally. +!> @brief Calculate shear rates on slip systems and their derivatives with respect to resolved +! stress. +!> @details: Derivatives are calculated only optionally. ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to -! have the optional arguments at the end +! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- pure subroutine kinetics(Mp,instance,of, & gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) @@ -462,44 +414,44 @@ pure subroutine kinetics(Mp,instance,of, & integer, intent(in) :: & instance, & of - + real(pReal), intent(out), dimension(param(instance)%totalNslip) :: & gdot_pos, & gdot_neg real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: & dgdot_dtau_pos, & dgdot_dtau_neg - + real(pReal), dimension(param(instance)%totalNslip) :: & tau_pos, & tau_neg integer :: i logical :: nonSchmidActive - + associate(prm => param(instance), stt => state(instance)) - + nonSchmidActive = size(prm%nonSchmidCoeff) > 0 - + do i = 1, prm%totalNslip tau_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,of) tau_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)) - stt%crss_back(i,of), & 0.0_pReal, nonSchmidActive) enddo - + where(dNeq0(tau_pos)) gdot_pos = prm%gdot0 * merge(0.5_pReal,1.0_pReal, nonSchmidActive) & ! 1/2 if non-Schmid active * sign(abs(tau_pos/stt%crss(:,of))**prm%n, tau_pos) else where gdot_pos = 0.0_pReal end where - + where(dNeq0(tau_neg)) gdot_neg = prm%gdot0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2 * sign(abs(tau_neg/stt%crss(:,of))**prm%n, tau_neg) else where gdot_neg = 0.0_pReal end where - + if (present(dgdot_dtau_pos)) then where(dNeq0(gdot_pos)) dgdot_dtau_pos = gdot_pos*prm%n/tau_pos diff --git a/src/constitutive_plastic_none.f90 b/src/constitutive_plastic_none.f90 index 578fb4149..4840c9c09 100644 --- a/src/constitutive_plastic_none.f90 +++ b/src/constitutive_plastic_none.f90 @@ -18,19 +18,19 @@ module subroutine plastic_none_init Ninstance, & p, & NipcMyPhase - - write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_NONE_label//' init -+>>>' - + + write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_NONE_label//' init -+>>>'; flush(6) + Ninstance = count(phase_plasticity == PLASTICITY_NONE_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - + do p = 1, size(phase_plasticity) if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle - + NipcMyPhase = count(material_phaseAt == p) * discretization_nIP call material_allocatePlasticState(p,NipcMyPhase,0,0,0) - + enddo end subroutine plastic_none_init diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index 3739ab669..47db19385 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -48,32 +48,6 @@ submodule(constitutive) plastic_nonlocal real(pReal), dimension(:,:,:,:,:,:), allocatable :: & compatibility !< slip system compatibility between me and my neighbors - enum, bind(c) - enumerator :: & - undefined_ID, & - rho_sgl_mob_edg_pos_ID, & - rho_sgl_mob_edg_neg_ID, & - rho_sgl_mob_scr_pos_ID, & - rho_sgl_mob_scr_neg_ID, & - rho_sgl_imm_edg_pos_ID, & - rho_sgl_imm_edg_neg_ID, & - rho_sgl_imm_scr_pos_ID, & - rho_sgl_imm_scr_neg_ID, & - rho_dip_edg_ID, & - rho_dip_scr_ID, & - rho_forest_ID, & - resolvedstress_back_ID, & - tau_pass_ID, & - rho_dot_sgl_ID, & - rho_dot_sgl_mobile_ID, & - rho_dot_dip_ID, & - v_edg_pos_ID, & - v_edg_neg_ID, & - v_scr_pos_ID, & - v_scr_neg_ID, & - gamma_ID - end enum - type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & atomicVolume, & !< atomic volume @@ -135,14 +109,12 @@ submodule(constitutive) plastic_nonlocal integer, dimension(:) ,allocatable :: & Nslip,& colinearSystem !< colinear system to the active slip system (only valid for fcc!) - + character(len=pStringLen), allocatable, dimension(:) :: & + output logical :: & shortRangeStressCorrection, & !< flag indicating the use of the short range stress correction by a excess density gradient term probabilisticMultiplication - integer(kind(undefined_ID)), dimension(:), allocatable :: & - outputID !< ID of each post result output - end type tParameters type :: tNonlocalMicrostructure @@ -198,22 +170,19 @@ module subroutine plastic_nonlocal_init integer :: & sizeState, sizeDotState,sizeDependentState, sizeDeltaState, & maxNinstances, & - p, i, & + p, & l, & s1, s2, & s, & t, & c - integer(kind(undefined_ID)) :: & - outputID character(len=pStringLen) :: & extmsg = '', & structure - character(len=pStringLen), dimension(:), allocatable :: outputs integer :: NofMyPhase - write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONLOCAL_label//' init -+>>>' + write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONLOCAL_label//' init -+>>>'; flush(6) write(6,'(/,a)') ' Reuber et al., Acta Materialia 71:333–348, 2014' write(6,'(a)') ' https://doi.org/10.1016/j.actamat.2014.03.012' @@ -407,60 +376,7 @@ module subroutine plastic_nonlocal_init endif slipActive - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - do i=1, size(outputs) - outputID = undefined_ID - select case(trim(outputs(i))) - case ('rho_sgl_mob_edg_pos') - outputID = merge(rho_sgl_mob_edg_pos_ID,undefined_ID,prm%totalNslip>0) - case ('rho_sgl_mob_edg_neg') - outputID = merge(rho_sgl_mob_edg_neg_ID,undefined_ID,prm%totalNslip>0) - case ('rho_sgl_mob_scr_pos') - outputID = merge(rho_sgl_mob_scr_pos_ID,undefined_ID,prm%totalNslip>0) - case ('rho_sgl_mob_scr_neg') - outputID = merge(rho_sgl_mob_scr_neg_ID,undefined_ID,prm%totalNslip>0) - case ('rho_sgl_imm_edg_pos') - outputID = merge(rho_sgl_imm_edg_pos_ID,undefined_ID,prm%totalNslip>0) - case ('rho_sgl_imm_edg_neg') - outputID = merge(rho_sgl_imm_edg_neg_ID,undefined_ID,prm%totalNslip>0) - case ('rho_sgl_imm_scr_pos') - outputID = merge(rho_sgl_imm_scr_pos_ID,undefined_ID,prm%totalNslip>0) - case ('rho_sgl_imm_scr_neg') - outputID = merge(rho_sgl_imm_scr_neg_ID,undefined_ID,prm%totalNslip>0) - case ('rho_dip_edg') - outputID = merge(rho_dip_edg_ID,undefined_ID,prm%totalNslip>0) - case ('rho_dip_scr') - outputID = merge(rho_dip_scr_ID,undefined_ID,prm%totalNslip>0) - case ('rho_forest') - outputID = merge(rho_forest_ID,undefined_ID,prm%totalNslip>0) - case ('resolvedstress_back') - outputID = merge(resolvedstress_back_ID,undefined_ID,prm%totalNslip>0) - case ('tau_pass') - outputID = merge(tau_pass_ID,undefined_ID,prm%totalNslip>0) - case ('rho_dot_sgl') - outputID = merge(rho_dot_sgl_ID,undefined_ID,prm%totalNslip>0) - case ('rho_dot_sgl_mobile') - outputID = merge(rho_dot_sgl_mobile_ID,undefined_ID,prm%totalNslip>0) - case ('rho_dot_dip') - outputID = merge(rho_dot_dip_ID,undefined_ID,prm%totalNslip>0) - case ('v_edg_pos') - outputID = merge(v_edg_pos_ID,undefined_ID,prm%totalNslip>0) - case ('v_edg_neg') - outputID = merge(v_edg_neg_ID,undefined_ID,prm%totalNslip>0) - case ('v_scr_pos') - outputID = merge(v_scr_pos_ID,undefined_ID,prm%totalNslip>0) - case ('v_scr_neg') - outputID = merge(v_scr_neg_ID,undefined_ID,prm%totalNslip>0) - case ('gamma') - outputID = merge(gamma_ID,undefined_ID,prm%totalNslip>0) - end select - - if (outputID /= undefined_ID) then - prm%outputID = [prm%outputID , outputID] - endif - - enddo + prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) !-------------------------------------------------------------------------------------------------- ! allocate state arrays @@ -561,7 +477,7 @@ module subroutine plastic_nonlocal_init stt%v_scr_neg => plasticState(p)%state (15*prm%totalNslip + 1:16*prm%totalNslip ,1:NofMyPhase) allocate(dst%tau_pass(prm%totalNslip,NofMyPhase),source=0.0_pReal) - allocate(dst%tau_Back(prm%totalNslip,NofMyPhase), source=0.0_pReal) + allocate(dst%tau_back(prm%totalNslip,NofMyPhase),source=0.0_pReal) end associate if (NofMyPhase > 0) call stateInit(p,NofMyPhase) @@ -1968,62 +1884,63 @@ module subroutine plastic_nonlocal_results(instance,group) integer, intent(in) :: instance character(len=*),intent(in) :: group + integer :: o associate(prm => param(instance),dst => microstructure(instance),stt=>state(instance)) - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - case (rho_sgl_mob_edg_pos_ID) - call results_writeDataset(group,stt%rho_sgl_mob_edg_pos, 'rho_sgl_mob_edg_pos', & - 'positive mobile edge density','1/m²') - case (rho_sgl_imm_edg_pos_ID) - call results_writeDataset(group,stt%rho_sgl_imm_edg_pos, 'rho_sgl_imm_edg_pos',& - 'positive immobile edge density','1/m²') - case (rho_sgl_mob_edg_neg_ID) - call results_writeDataset(group,stt%rho_sgl_mob_edg_neg, 'rho_sgl_mob_edg_neg',& - 'negative mobile edge density','1/m²') - case (rho_sgl_imm_edg_neg_ID) - call results_writeDataset(group,stt%rho_sgl_imm_edg_neg, 'rho_sgl_imm_edg_neg',& - 'negative immobile edge density','1/m²') - case (rho_dip_edg_ID) - call results_writeDataset(group,stt%rho_dip_edg, 'rho_dip_edg',& - 'edge dipole density','1/m²') - case (rho_sgl_mob_scr_pos_ID) - call results_writeDataset(group,stt%rho_sgl_mob_scr_pos, 'rho_sgl_mob_scr_pos',& - 'positive mobile screw density','1/m²') - case (rho_sgl_imm_scr_pos_ID) - call results_writeDataset(group,stt%rho_sgl_imm_scr_pos, 'rho_sgl_imm_scr_pos',& - 'positive immobile screw density','1/m²') - case (rho_sgl_mob_scr_neg_ID) - call results_writeDataset(group,stt%rho_sgl_mob_scr_neg, 'rho_sgl_mob_scr_neg',& - 'negative mobile screw density','1/m²') - case (rho_sgl_imm_scr_neg_ID) - call results_writeDataset(group,stt%rho_sgl_imm_scr_neg, 'rho_sgl_imm_scr_neg',& - 'negative immobile screw density','1/m²') - case (rho_dip_scr_ID) - call results_writeDataset(group,stt%rho_dip_scr, 'rho_dip_scr',& - 'screw dipole density','1/m²') - case (rho_forest_ID) - call results_writeDataset(group,stt%rho_forest, 'rho_forest',& - 'forest density','1/m²') - case (v_edg_pos_ID) - call results_writeDataset(group,stt%v_edg_pos, 'v_edg_pos',& - 'positive edge velocity','m/s') - case (v_edg_neg_ID) - call results_writeDataset(group,stt%v_edg_neg, 'v_edg_neg',& - 'negative edge velocity','m/s') - case (v_scr_pos_ID) - call results_writeDataset(group,stt%v_scr_pos, 'v_scr_pos',& - 'positive srew velocity','m/s') - case (v_scr_neg_ID) - call results_writeDataset(group,stt%v_scr_neg, 'v_scr_neg',& - 'negative screw velocity','m/s') - case(gamma_ID) - call results_writeDataset(group,stt%gamma,'gamma',& - 'plastic shear','1') - case (tau_pass_ID) - call results_writeDataset(group,dst%tau_pass,'tau_pass',& - 'passing stress for slip','Pa') + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case('rho_sgl_mob_edg_pos') + if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_mob_edg_pos, 'rho_sgl_mob_edg_pos', & + 'positive mobile edge density','1/m²') + case('rho_sgl_imm_edg_pos') + if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_imm_edg_pos, 'rho_sgl_imm_edg_pos',& + 'positive immobile edge density','1/m²') + case('rho_sgl_mob_edg_neg') + if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_mob_edg_neg, 'rho_sgl_mob_edg_neg',& + 'negative mobile edge density','1/m²') + case('rho_sgl_imm_edg_neg') + if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_imm_edg_neg, 'rho_sgl_imm_edg_neg',& + 'negative immobile edge density','1/m²') + case('rho_dip_edg') + if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_dip_edg, 'rho_dip_edg',& + 'edge dipole density','1/m²') + case('rho_sgl_mob_scr_pos') + if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_mob_scr_pos, 'rho_sgl_mob_scr_pos',& + 'positive mobile screw density','1/m²') + case('rho_sgl_imm_scr_pos') + if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_imm_scr_pos, 'rho_sgl_imm_scr_pos',& + 'positive immobile screw density','1/m²') + case('rho_sgl_mob_scr_neg') + if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_mob_scr_neg, 'rho_sgl_mob_scr_neg',& + 'negative mobile screw density','1/m²') + case('rho_sgl_imm_scr_neg') + if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_imm_scr_neg, 'rho_sgl_imm_scr_neg',& + 'negative immobile screw density','1/m²') + case('rho_dip_scr') + if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_dip_scr, 'rho_dip_scr',& + 'screw dipole density','1/m²') + case('rho_forest') + if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_forest, 'rho_forest',& + 'forest density','1/m²') + case('v_edg_pos') + if(prm%totalNslip>0) call results_writeDataset(group,stt%v_edg_pos, 'v_edg_pos',& + 'positive edge velocity','m/s') + case('v_edg_neg') + if(prm%totalNslip>0) call results_writeDataset(group,stt%v_edg_neg, 'v_edg_neg',& + 'negative edge velocity','m/s') + case('v_scr_pos') + if(prm%totalNslip>0) call results_writeDataset(group,stt%v_scr_pos, 'v_scr_pos',& + 'positive srew velocity','m/s') + case('v_scr_neg') + if(prm%totalNslip>0) call results_writeDataset(group,stt%v_scr_neg, 'v_scr_neg',& + 'negative screw velocity','m/s') + case('gamma') + if(prm%totalNslip>0) call results_writeDataset(group,stt%gamma,'gamma',& + 'plastic shear','1') + case('tau_pass') + if(prm%totalNslip>0) call results_writeDataset(group,dst%tau_pass,'tau_pass',& + 'passing stress for slip','Pa') end select enddo outputsLoop end associate diff --git a/src/constitutive_plastic_phenopowerlaw.f90 b/src/constitutive_plastic_phenopowerlaw.f90 index 165ad0081..6dfa3fb16 100644 --- a/src/constitutive_plastic_phenopowerlaw.f90 +++ b/src/constitutive_plastic_phenopowerlaw.f90 @@ -6,19 +6,6 @@ !-------------------------------------------------------------------------------------------------- submodule(constitutive) plastic_phenopowerlaw - enum, bind(c) - enumerator :: & - undefined_ID, & - resistance_slip_ID, & - accumulatedshear_slip_ID, & - shearrate_slip_ID, & - resolvedstress_slip_ID, & - resistance_twin_ID, & - accumulatedshear_twin_ID, & - shearrate_twin_ID, & - resolvedstress_twin_ID - end enum - type :: tParameters real(pReal) :: & gdot0_slip, & !< reference shear strain rate for slip @@ -37,19 +24,19 @@ submodule(constitutive) plastic_phenopowerlaw aTolResistance, & !< absolute tolerance for integration of xi aTolShear, & !< absolute tolerance for integration of gamma aTolTwinfrac !< absolute tolerance for integration of f - real(pReal), allocatable, dimension(:) :: & + real(pReal), allocatable, dimension(:) :: & xi_slip_0, & !< initial critical shear stress for slip xi_twin_0, & !< initial critical shear stress for twin xi_slip_sat, & !< maximum critical shear stress for slip nonSchmidCoeff, & H_int, & !< per family hardening activity (optional) gamma_twin_char !< characteristic shear for twins - real(pReal), allocatable, dimension(:,:) :: & + real(pReal), allocatable, dimension(:,:) :: & interaction_SlipSlip, & !< slip resistance from slip activity interaction_SlipTwin, & !< slip resistance from twin activity interaction_TwinSlip, & !< twin resistance from slip activity interaction_TwinTwin !< twin resistance from twin activity - real(pReal), allocatable, dimension(:,:,:) :: & + real(pReal), allocatable, dimension(:,:,:) :: & Schmid_slip, & Schmid_twin, & nonSchmid_pos, & @@ -57,13 +44,13 @@ submodule(constitutive) plastic_phenopowerlaw integer :: & totalNslip, & !< total number of active slip system totalNtwin !< total number of active twin systems - integer, allocatable, dimension(:) :: & + integer, allocatable, dimension(:) :: & Nslip, & !< number of active slip systems for each family Ntwin !< number of active twin systems for each family - integer(kind(undefined_ID)), allocatable, dimension(:) :: & - outputID !< ID of each post result output + character(len=pStringLen), allocatable, dimension(:) :: & + output end type tParameters - + type :: tPhenopowerlawState real(pReal), pointer, dimension(:,:) :: & xi_slip, & @@ -71,7 +58,7 @@ submodule(constitutive) plastic_phenopowerlaw gamma_slip, & gamma_twin end type tPhenopowerlawState - + !-------------------------------------------------------------------------------------------------- ! containers for parameters and state type(tParameters), allocatable, dimension(:) :: param @@ -83,7 +70,7 @@ contains !-------------------------------------------------------------------------------------------------- -!> @brief module initialization +!> @brief Perform module initialization. !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- module subroutine plastic_phenopowerlaw_init @@ -91,51 +78,46 @@ module subroutine plastic_phenopowerlaw_init integer :: & Ninstance, & p, i, & - NipcMyPhase, outputSize, & + NipcMyPhase, & sizeState, sizeDotState, & startIndex, endIndex - - integer(kind(undefined_ID)) :: & - outputID - + character(len=pStringLen) :: & extmsg = '' - character(len=pStringLen), dimension(:), allocatable :: & - outputs - - write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>' - + + write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>'; flush(6) + Ninstance = count(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - + allocate(param(Ninstance)) allocate(state(Ninstance)) allocate(dotState(Ninstance)) - + do p = 1, size(phase_plasticity) if (phase_plasticity(p) /= PLASTICITY_PHENOPOWERLAW_ID) cycle associate(prm => param(phase_plasticityInstance(p)), & dot => dotState(phase_plasticityInstance(p)), & stt => state(phase_plasticityInstance(p)), & config => config_phase(p)) - + !-------------------------------------------------------------------------------------------------- ! optional parameters that need to be defined prm%c_1 = config%getFloat('twin_c',defaultVal=0.0_pReal) prm%c_2 = config%getFloat('twin_b',defaultVal=1.0_pReal) prm%c_3 = config%getFloat('twin_e',defaultVal=0.0_pReal) prm%c_4 = config%getFloat('twin_d',defaultVal=0.0_pReal) - + prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal) prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal) prm%aTolTwinfrac = config%getFloat('atol_twinfrac', defaultVal=1.0e-6_pReal) - + ! sanity checks if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//' aTolresistance' if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//' aTolShear' if (prm%aTolTwinfrac <= 0.0_pReal) extmsg = trim(extmsg)//' atoltwinfrac' - + !-------------------------------------------------------------------------------------------------- ! slip related parameters prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) @@ -143,7 +125,7 @@ module subroutine plastic_phenopowerlaw_init slipActive: if (prm%totalNslip > 0) then prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) - + if(trim(config%getString('lattice_structure')) == 'bcc') then prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',& defaultVal = emptyRealArray) @@ -157,22 +139,22 @@ module subroutine plastic_phenopowerlaw_init prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, & config%getFloats('interaction_slipslip'), & config%getString('lattice_structure')) - + prm%xi_slip_0 = config%getFloats('tau0_slip', requiredSize=size(prm%Nslip)) prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(prm%Nslip)) prm%H_int = config%getFloats('h_int', requiredSize=size(prm%Nslip), & defaultVal=[(0.0_pReal,i=1,size(prm%Nslip))]) - + prm%gdot0_slip = config%getFloat('gdot0_slip') prm%n_slip = config%getFloat('n_slip') prm%a_slip = config%getFloat('a_slip') prm%h0_SlipSlip = config%getFloat('h0_slipslip') - + ! expand: family => system prm%xi_slip_0 = math_expand(prm%xi_slip_0, prm%Nslip) prm%xi_slip_sat = math_expand(prm%xi_slip_sat,prm%Nslip) prm%H_int = math_expand(prm%H_int, prm%Nslip) - + ! sanity checks if ( prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_slip' if ( prm%a_slip <= 0.0_pReal) extmsg = trim(extmsg)//' a_slip' @@ -183,7 +165,7 @@ module subroutine plastic_phenopowerlaw_init allocate(prm%interaction_SlipSlip(0,0)) allocate(prm%xi_slip_0(0)) endif slipActive - + !-------------------------------------------------------------------------------------------------- ! twin related parameters prm%Ntwin = config%getInts('ntwin', defaultVal=emptyIntArray) @@ -196,17 +178,17 @@ module subroutine plastic_phenopowerlaw_init config%getString('lattice_structure')) prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,config%getString('lattice_structure'),& config%getFloat('c/a')) - + prm%xi_twin_0 = config%getFloats('tau0_twin',requiredSize=size(prm%Ntwin)) - + prm%gdot0_twin = config%getFloat('gdot0_twin') prm%n_twin = config%getFloat('n_twin') prm%spr = config%getFloat('s_pr') prm%h0_TwinTwin = config%getFloat('h0_twintwin') - + ! expand: family => system prm%xi_twin_0 = math_expand(prm%xi_twin_0, prm%Ntwin) - + ! sanity checks if (prm%gdot0_twin <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_twin' if (prm%n_twin <= 0.0_pReal) extmsg = trim(extmsg)//' n_twin' @@ -215,7 +197,7 @@ module subroutine plastic_phenopowerlaw_init allocate(prm%xi_twin_0(0)) allocate(prm%gamma_twin_char(0)) endif twinActive - + !-------------------------------------------------------------------------------------------------- ! slip-twin related parameters slipAndTwinActive: if (prm%totalNslip > 0 .and. prm%totalNtwin > 0) then @@ -231,63 +213,25 @@ module subroutine plastic_phenopowerlaw_init allocate(prm%interaction_TwinSlip(prm%TotalNtwin,prm%TotalNslip)) ! at least one dimension is 0 prm%h0_TwinSlip = 0.0_pReal endif slipAndTwinActive - + !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range if (extmsg /= '') & call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_PHENOPOWERLAW_label//')') - + !-------------------------------------------------------------------------------------------------- ! output pararameters - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - do i=1, size(outputs) - outputID = undefined_ID - select case(outputs(i)) - - case ('resistance_slip') - outputID = merge(resistance_slip_ID,undefined_ID,prm%totalNslip>0) - outputSize = prm%totalNslip - case ('accumulatedshear_slip') - outputID = merge(accumulatedshear_slip_ID,undefined_ID,prm%totalNslip>0) - outputSize = prm%totalNslip - case ('shearrate_slip') - outputID = merge(shearrate_slip_ID,undefined_ID,prm%totalNslip>0) - outputSize = prm%totalNslip - case ('resolvedstress_slip') - outputID = merge(resolvedstress_slip_ID,undefined_ID,prm%totalNslip>0) - outputSize = prm%totalNslip - - case ('resistance_twin') - outputID = merge(resistance_twin_ID,undefined_ID,prm%totalNtwin>0) - outputSize = prm%totalNtwin - case ('accumulatedshear_twin') - outputID = merge(accumulatedshear_twin_ID,undefined_ID,prm%totalNtwin>0) - outputSize = prm%totalNtwin - case ('shearrate_twin') - outputID = merge(shearrate_twin_ID,undefined_ID,prm%totalNtwin>0) - outputSize = prm%totalNtwin - case ('resolvedstress_twin') - outputID = merge(resolvedstress_twin_ID,undefined_ID,prm%totalNtwin>0) - outputSize = prm%totalNtwin - - end select - - if (outputID /= undefined_ID) then - prm%outputID = [prm%outputID, outputID] - endif - - enddo - + prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) + !-------------------------------------------------------------------------------------------------- ! allocate state arrays NipcMyPhase = count(material_phaseAt == p) * discretization_nIP sizeDotState = size(['tau_slip ','gamma_slip']) * prm%totalNslip & + size(['tau_twin ','gamma_twin']) * prm%totalNtwin sizeState = sizeDotState - + call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) - + !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and aTolState startIndex = 1 @@ -296,14 +240,14 @@ module subroutine plastic_phenopowerlaw_init stt%xi_slip = spread(prm%xi_slip_0, 2, NipcMyPhase) dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance - + startIndex = endIndex + 1 endIndex = endIndex + prm%totalNtwin stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:) stt%xi_twin = spread(prm%xi_twin_0, 2, NipcMyPhase) dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance - + startIndex = endIndex + 1 endIndex = endIndex + prm%totalNslip stt%gamma_slip => plasticState(p)%state (startIndex:endIndex,:) @@ -317,18 +261,18 @@ module subroutine plastic_phenopowerlaw_init stt%gamma_twin => plasticState(p)%state (startIndex:endIndex,:) dot%gamma_twin => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear - + plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally - + end associate - + enddo end subroutine plastic_phenopowerlaw_init !-------------------------------------------------------------------------------------------------- -!> @brief calculates plastic velocity gradient and its tangent +!> @brief Calculate plastic velocity gradient and its tangent. !> @details asummes that deformation by dislocation glide affects twinned and untwinned volume ! equally (Taylor assumption). Twinning happens only in untwinned volume !-------------------------------------------------------------------------------------------------- @@ -338,13 +282,13 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & instance, & of - + integer :: & i,k,l,m,n real(pReal), dimension(param(instance)%totalNslip) :: & @@ -352,12 +296,12 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta dgdot_dtauslip_pos,dgdot_dtauslip_neg real(pReal), dimension(param(instance)%totalNtwin) :: & gdot_twin,dgdot_dtautwin - + Lp = 0.0_pReal dLp_dMp = 0.0_pReal - + associate(prm => param(instance)) - + call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) slipSystems: do i = 1, prm%totalNslip Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i) @@ -366,7 +310,7 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta + dgdot_dtauslip_pos(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_pos(m,n,i) & + dgdot_dtauslip_neg(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_neg(m,n,i) enddo slipSystems - + call kinetics_twin(Mp,instance,of,gdot_twin,dgdot_dtautwin) twinSystems: do i = 1, prm%totalNtwin Lp = Lp + gdot_twin(i)*prm%Schmid_twin(1:3,1:3,i) @@ -374,14 +318,14 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + dgdot_dtautwin(i)*prm%Schmid_twin(k,l,i)*prm%Schmid_twin(m,n,i) enddo twinSystems - + end associate end subroutine plastic_phenopowerlaw_LpAndItsTangent !-------------------------------------------------------------------------------------------------- -!> @brief calculates the rate of change of microstructure +!> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) @@ -390,7 +334,7 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) integer, intent(in) :: & instance, & of - + real(pReal) :: & c_SlipSlip,c_TwinSlip,c_TwinTwin, & xi_slip_sat_offset,& @@ -398,9 +342,9 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) real(pReal), dimension(param(instance)%totalNslip) :: & left_SlipSlip,right_SlipSlip, & gdot_slip_pos,gdot_slip_neg - + associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) - + sumGamma = sum(stt%gamma_slip(:,of)) sumF = sum(stt%gamma_twin(:,of)/prm%gamma_twin_char) @@ -409,26 +353,26 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) c_SlipSlip = prm%h0_slipslip * (1.0_pReal + prm%c_1*sumF** prm%c_2) c_TwinSlip = prm%h0_TwinSlip * sumGamma**prm%c_3 c_TwinTwin = prm%h0_TwinTwin * sumF**prm%c_4 - + !-------------------------------------------------------------------------------------------------- ! calculate left and right vectors left_SlipSlip = 1.0_pReal + prm%H_int xi_slip_sat_offset = prm%spr*sqrt(sumF) right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) **prm%a_slip & * sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) - + !-------------------------------------------------------------------------------------------------- ! shear rates call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg) dot%gamma_slip(:,of) = abs(gdot_slip_pos+gdot_slip_neg) call kinetics_twin(Mp,instance,of,dot%gamma_twin(:,of)) - + !-------------------------------------------------------------------------------------------------- ! hardening dot%xi_slip(:,of) = c_SlipSlip * left_SlipSlip * & matmul(prm%interaction_SlipSlip,dot%gamma_slip(:,of)*right_SlipSlip) & + matmul(prm%interaction_SlipTwin,dot%gamma_twin(:,of)) - + dot%xi_twin(:,of) = c_TwinSlip * matmul(prm%interaction_TwinSlip,dot%gamma_slip(:,of)) & + c_TwinTwin * matmul(prm%interaction_TwinTwin,dot%gamma_twin(:,of)) end associate @@ -437,33 +381,33 @@ end subroutine plastic_phenopowerlaw_dotState !-------------------------------------------------------------------------------------------------- -!> @brief writes results to HDF5 output file +!> @brief Write results to HDF5 output file. !-------------------------------------------------------------------------------------------------- module subroutine plastic_phenopowerlaw_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group - + integer :: o associate(prm => param(instance), stt => state(instance)) - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + + case('resistance_slip') + if(prm%totalNslip>0) call results_writeDataset(group,stt%xi_slip, 'xi_sl', & + 'resistance against plastic slip','Pa') + case('accumulatedshear_slip') + if(prm%totalNslip>0) call results_writeDataset(group,stt%gamma_slip,'gamma_sl', & + 'plastic shear','1') + + case('resistance_twin') + if(prm%totalNtwin>0) call results_writeDataset(group,stt%xi_twin, 'xi_tw', & + 'resistance against twinning','Pa') + case('accumulatedshear_twin') + if(prm%totalNtwin>0) call results_writeDataset(group,stt%gamma_twin,'gamma_tw', & + 'twinning shear','1') - case (resistance_slip_ID) - call results_writeDataset(group,stt%xi_slip, 'xi_sl', & - 'resistance against plastic slip','Pa') - case (accumulatedshear_slip_ID) - call results_writeDataset(group,stt%gamma_slip,'gamma_sl', & - 'plastic shear','1') - - case (resistance_twin_ID) - call results_writeDataset(group,stt%xi_twin, 'xi_tw', & - 'resistance against twinning','Pa') - case (accumulatedshear_twin_ID) - call results_writeDataset(group,stt%gamma_twin,'gamma_tw', & - 'twinning shear','1') - end select enddo outputsLoop end associate @@ -472,10 +416,11 @@ end subroutine plastic_phenopowerlaw_results !-------------------------------------------------------------------------------------------------- -!> @brief Shear rates on slip systems and their derivatives with respect to resolved stress +!> @brief Calculate shear rates on slip systems and their derivatives with respect to resolved +! stress. !> @details Derivatives are calculated only optionally. ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to -! have the optional arguments at the end +! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- pure subroutine kinetics_slip(Mp,instance,of, & gdot_slip_pos,gdot_slip_neg,dgdot_dtau_slip_pos,dgdot_dtau_slip_neg) @@ -485,44 +430,44 @@ pure subroutine kinetics_slip(Mp,instance,of, & integer, intent(in) :: & instance, & of - + real(pReal), intent(out), dimension(param(instance)%totalNslip) :: & gdot_slip_pos, & gdot_slip_neg real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: & dgdot_dtau_slip_pos, & dgdot_dtau_slip_neg - + real(pReal), dimension(param(instance)%totalNslip) :: & tau_slip_pos, & tau_slip_neg integer :: i logical :: nonSchmidActive - + associate(prm => param(instance), stt => state(instance)) - + nonSchmidActive = size(prm%nonSchmidCoeff) > 0 - + do i = 1, prm%totalNslip tau_slip_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i)) tau_slip_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)), & 0.0_pReal, nonSchmidActive) enddo - + where(dNeq0(tau_slip_pos)) gdot_slip_pos = prm%gdot0_slip * merge(0.5_pReal,1.0_pReal, nonSchmidActive) & ! 1/2 if non-Schmid active * sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_slip, tau_slip_pos) else where gdot_slip_pos = 0.0_pReal end where - + where(dNeq0(tau_slip_neg)) gdot_slip_neg = prm%gdot0_slip * 0.5_pReal & ! only used if non-Schmid active, always 1/2 * sign(abs(tau_slip_neg/stt%xi_slip(:,of))**prm%n_slip, tau_slip_neg) else where gdot_slip_neg = 0.0_pReal end where - + if (present(dgdot_dtau_slip_pos)) then where(dNeq0(gdot_slip_pos)) dgdot_dtau_slip_pos = gdot_slip_pos*prm%n_slip/tau_slip_pos @@ -543,9 +488,9 @@ end subroutine kinetics_slip !-------------------------------------------------------------------------------------------------- -!> @brief Shear rates on twin systems and their derivatives with respect to resolved stress. -! twinning is assumed to take place only in untwinned volume. -!> @details Derivates are calculated only optionally. +!> @brief Calculate shear rates on twin systems and their derivatives with respect to resolved +! stress. Twinning is assumed to take place only in untwinned volume. +!> @details Derivatives are calculated only optionally. ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- @@ -557,29 +502,29 @@ pure subroutine kinetics_twin(Mp,instance,of,& integer, intent(in) :: & instance, & of - + real(pReal), dimension(param(instance)%totalNtwin), intent(out) :: & gdot_twin real(pReal), dimension(param(instance)%totalNtwin), intent(out), optional :: & dgdot_dtau_twin - + real(pReal), dimension(param(instance)%totalNtwin) :: & tau_twin integer :: i - + associate(prm => param(instance), stt => state(instance)) - + do i = 1, prm%totalNtwin tau_twin(i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i)) enddo - + where(tau_twin > 0.0_pReal) gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)) & ! only twin in untwinned volume fraction * prm%gdot0_twin*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_twin else where gdot_twin = 0.0_pReal end where - + if (present(dgdot_dtau_twin)) then where(dNeq0(gdot_twin)) dgdot_dtau_twin = gdot_twin*prm%n_twin/tau_twin @@ -587,7 +532,7 @@ pure subroutine kinetics_twin(Mp,instance,of,& dgdot_dtau_twin = 0.0_pReal end where endif - + end associate end subroutine kinetics_twin diff --git a/src/crystallite.f90 b/src/crystallite.f90 index a1610b254..1a8252cce 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -57,9 +57,7 @@ module crystallite crystallite_Li0, & !< intermediate velocitiy grad at start of FE inc crystallite_partionedLi0 !< intermediate velocity grad at start of homog inc real(pReal), dimension(:,:,:,:,:), allocatable :: & - crystallite_invFp, & !< inverse of current plastic def grad (end of converged time step) crystallite_subFp0,& !< plastic def grad at start of crystallite inc - crystallite_invFi, & !< inverse of current intermediate def grad (end of converged time step) crystallite_subFi0,& !< intermediate def grad at start of crystallite inc crystallite_subF, & !< def grad to be reached at end of crystallite inc crystallite_subF0, & !< def grad at start of crystallite inc @@ -145,12 +143,10 @@ subroutine crystallite_init allocate(crystallite_partionedFp0(3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_subFp0(3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_Fp(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_invFp(3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_Fi0(3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_partionedFi0(3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_subFi0(3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_Fi(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_invFi(3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_Fe(3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_Lp0(3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_partionedLp0(3,3,cMax,iMax,eMax), source=0.0_pReal) @@ -408,9 +404,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) else crystallite_subStep(c,i,e) = num%subStepSizeCryst * crystallite_subStep(c,i,e) crystallite_Fp (1:3,1:3,c,i,e) = crystallite_subFp0(1:3,1:3,c,i,e) - crystallite_invFp(1:3,1:3,c,i,e) = math_inv33(crystallite_Fp (1:3,1:3,c,i,e)) crystallite_Fi (1:3,1:3,c,i,e) = crystallite_subFi0(1:3,1:3,c,i,e) - crystallite_invFi(1:3,1:3,c,i,e) = math_inv33(crystallite_Fi (1:3,1:3,c,i,e)) crystallite_S (1:3,1:3,c,i,e) = crystallite_S0 (1:3,1:3,c,i,e) if (crystallite_subStep(c,i,e) < 1.0_pReal) then ! actual (not initial) cutback crystallite_Lp (1:3,1:3,c,i,e) = crystallite_subLp0(1:3,1:3,c,i,e) @@ -431,11 +425,11 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) ! prepare for integration if (crystallite_todo(c,i,e)) then crystallite_subF(1:3,1:3,c,i,e) = crystallite_subF0(1:3,1:3,c,i,e) & - + crystallite_subStep(c,i,e) * (crystallite_partionedF (1:3,1:3,c,i,e) & - - crystallite_partionedF0(1:3,1:3,c,i,e)) - crystallite_Fe(1:3,1:3,c,i,e) = matmul(matmul(crystallite_subF (1:3,1:3,c,i,e), & - crystallite_invFp(1:3,1:3,c,i,e)), & - crystallite_invFi(1:3,1:3,c,i,e)) + + crystallite_subStep(c,i,e) *( crystallite_partionedF (1:3,1:3,c,i,e) & + -crystallite_partionedF0(1:3,1:3,c,i,e)) + crystallite_Fe(1:3,1:3,c,i,e) = matmul(matmul(crystallite_subF(1:3,1:3,c,i,e), & + math_inv33(crystallite_Fp(1:3,1:3,c,i,e))), & + math_inv33(crystallite_Fi(1:3,1:3,c,i,e))) crystallite_subdt(c,i,e) = crystallite_subStep(c,i,e) * crystallite_dt(c,i,e) crystallite_converged(c,i,e) = .false. endif @@ -477,7 +471,9 @@ subroutine crystallite_stressTangent o, & p - real(pReal), dimension(3,3) :: temp_33_1, devNull,invSubFi0, temp_33_2, temp_33_3, temp_33_4 + real(pReal), dimension(3,3) :: devNull, & + invSubFp0,invSubFi0,invFp,invFi, & + temp_33_1, temp_33_2, temp_33_3, temp_33_4 real(pReal), dimension(3,3,3,3) :: dSdFe, & dSdF, & dSdFi, & @@ -493,7 +489,8 @@ subroutine crystallite_stressTangent real(pReal), dimension(9,9):: temp_99 logical :: error - !$OMP PARALLEL DO PRIVATE(dSdF,dSdFe,dSdFi,dLpdS,dLpdFi,dFpinvdF,dLidS,dLidFi,dFidS,invSubFi0,o,p, & + !$OMP PARALLEL DO PRIVATE(dSdF,dSdFe,dSdFi,dLpdS,dLpdFi,dFpinvdF,dLidS,dLidFi,dFidS,o,p, & + !$OMP invSubFp0,invSubFi0,invFp,invFi, & !$OMP rhs_3333,lhs_3333,temp_99,temp_33_1,temp_33_2,temp_33_3,temp_33_4,temp_3333,error) elementLooping: do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) @@ -507,16 +504,20 @@ subroutine crystallite_stressTangent crystallite_Fi(1:3,1:3,c,i,e), & c,i,e) + invFp = math_inv33(crystallite_Fp(1:3,1:3,c,i,e)) + invFi = math_inv33(crystallite_Fi(1:3,1:3,c,i,e)) + invSubFp0 = math_inv33(crystallite_subFp0(1:3,1:3,c,i,e)) + invSubFi0 = math_inv33(crystallite_subFi0(1:3,1:3,c,i,e)) + if (sum(abs(dLidS)) < tol_math_check) then dFidS = 0.0_pReal else - invSubFi0 = math_inv33(crystallite_subFi0(1:3,1:3,c,i,e)) lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal do o=1,3; do p=1,3 lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) & + crystallite_subdt(c,i,e)*matmul(invSubFi0,dLidFi(1:3,1:3,o,p)) lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) & - + crystallite_invFi(1:3,1:3,c,i,e)*crystallite_invFi(p,o,c,i,e) + + invFi*invFi(p,o) rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) & - crystallite_subdt(c,i,e)*matmul(invSubFi0,dLidS(1:3,1:3,o,p)) enddo; enddo @@ -538,18 +539,13 @@ subroutine crystallite_stressTangent !-------------------------------------------------------------------------------------------------- ! calculate dSdF - temp_33_1 = transpose(matmul(crystallite_invFp(1:3,1:3,c,i,e), & - crystallite_invFi(1:3,1:3,c,i,e))) - temp_33_2 = matmul(crystallite_subF(1:3,1:3,c,i,e), & - math_inv33(crystallite_subFp0(1:3,1:3,c,i,e))) - temp_33_3 = matmul(matmul(crystallite_subF(1:3,1:3,c,i,e), & - crystallite_invFp (1:3,1:3,c,i,e)), & - math_inv33(crystallite_subFi0(1:3,1:3,c,i,e))) + temp_33_1 = transpose(matmul(invFp,invFi)) + temp_33_2 = matmul(crystallite_subF(1:3,1:3,c,i,e),invSubFp0) + temp_33_3 = matmul(matmul(crystallite_subF(1:3,1:3,c,i,e),invFp), invSubFi0) do o=1,3; do p=1,3 rhs_3333(p,o,1:3,1:3) = matmul(dSdFe(p,o,1:3,1:3),temp_33_1) - temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), & - crystallite_invFi(1:3,1:3,c,i,e)) & + temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), invFi) & + matmul(temp_33_3,dLidS(1:3,1:3,p,o)) enddo; enddo lhs_3333 = crystallite_subdt(c,i,e)*math_mul3333xx3333(dSdFe,temp_3333) & @@ -569,15 +565,14 @@ subroutine crystallite_stressTangent temp_3333 = math_mul3333xx3333(dLpdS,dSdF) do o=1,3; do p=1,3 dFpinvdF(1:3,1:3,p,o) = -crystallite_subdt(c,i,e) & - * matmul(math_inv33(crystallite_subFp0(1:3,1:3,c,i,e)), & - matmul(temp_3333(1:3,1:3,p,o),crystallite_invFi(1:3,1:3,c,i,e))) + * matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi)) enddo; enddo !-------------------------------------------------------------------------------------------------- ! assemble dPdF - temp_33_1 = matmul(crystallite_S(1:3,1:3,c,i,e),transpose(crystallite_invFp(1:3,1:3,c,i,e))) - temp_33_2 = matmul(crystallite_invFp(1:3,1:3,c,i,e),temp_33_1) - temp_33_3 = matmul(crystallite_subF(1:3,1:3,c,i,e),crystallite_invFp(1:3,1:3,c,i,e)) + temp_33_1 = matmul(crystallite_S(1:3,1:3,c,i,e),transpose(invFp)) + temp_33_2 = matmul(invFp,temp_33_1) + temp_33_3 = matmul(crystallite_subF(1:3,1:3,c,i,e),invFp) temp_33_4 = matmul(temp_33_3,crystallite_S(1:3,1:3,c,i,e)) crystallite_dPdF(1:3,1:3,1:3,1:3,c,i,e) = 0.0_pReal @@ -589,7 +584,7 @@ subroutine crystallite_stressTangent + matmul(matmul(crystallite_subF(1:3,1:3,c,i,e), & dFpinvdF(1:3,1:3,p,o)),temp_33_1) & + matmul(matmul(temp_33_3,dSdF(1:3,1:3,p,o)), & - transpose(crystallite_invFp(1:3,1:3,c,i,e))) & + transpose(invFp)) & + matmul(temp_33_4,transpose(dFpinvdF(1:3,1:3,p,o))) enddo; enddo @@ -793,28 +788,28 @@ logical function integrateStress(ipc,ip,el,timeFraction) real(pReal), dimension(3,3):: F, & ! deformation gradient at end of timestep Fp_new, & ! plastic deformation gradient at end of timestep - Fe_new, & ! elastic deformation gradient at end of timestep invFp_new, & ! inverse of Fp_new - Fi_new, & ! gradient of intermediate deformation stages - invFi_new, & invFp_current, & ! inverse of Fp_current - invFi_current, & ! inverse of Fp_current Lpguess, & ! current guess for plastic velocity gradient Lpguess_old, & ! known last good guess for plastic velocity gradient Lp_constitutive, & ! plastic velocity gradient resulting from constitutive law residuumLp, & ! current residuum of plastic velocity gradient residuumLp_old, & ! last residuum of plastic velocity gradient deltaLp, & ! direction of next guess + Fi_new, & ! gradient of intermediate deformation stages + invFi_new, & + invFi_current, & ! inverse of Fi_current Liguess, & ! current guess for intermediate velocity gradient Liguess_old, & ! known last good guess for intermediate velocity gradient Li_constitutive, & ! intermediate velocity gradient resulting from constitutive law residuumLi, & ! current residuum of intermediate velocity gradient residuumLi_old, & ! last residuum of intermediate velocity gradient deltaLi, & ! direction of next guess + Fe, & ! elastic deformation gradient + Fe_new, & S, & ! 2nd Piola-Kirchhoff Stress in plastic (lattice) configuration A, & B, & - Fe, & ! elastic deformation gradient temp_33 real(pReal), dimension(9) :: temp_9 ! needed for matrix inversion by LAPACK integer, dimension(9) :: devNull_9 ! needed for matrix inversion by LAPACK @@ -993,16 +988,13 @@ logical function integrateStress(ipc,ip,el,timeFraction) Fe_new = matmul(matmul(F,invFp_new),invFi_new) integrateStress = .true. - crystallite_P (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) ! ToDo: We propably do not need to store P! + crystallite_P (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) crystallite_S (1:3,1:3,ipc,ip,el) = S crystallite_Lp (1:3,1:3,ipc,ip,el) = Lpguess crystallite_Li (1:3,1:3,ipc,ip,el) = Liguess crystallite_Fp (1:3,1:3,ipc,ip,el) = Fp_new crystallite_Fi (1:3,1:3,ipc,ip,el) = Fi_new crystallite_Fe (1:3,1:3,ipc,ip,el) = Fe_new - crystallite_invFp(1:3,1:3,ipc,ip,el) = invFp_new - crystallite_invFi(1:3,1:3,ipc,ip,el) = invFi_new - end function integrateStress diff --git a/src/homogenization_mech_RGC.f90 b/src/homogenization_mech_RGC.f90 index 9cb7efdc8..1f1d9d6ec 100644 --- a/src/homogenization_mech_RGC.f90 +++ b/src/homogenization_mech_RGC.f90 @@ -9,17 +9,6 @@ submodule(homogenization) homogenization_mech_RGC use rotations - enum, bind(c) - enumerator :: & - undefined_ID, & - constitutivework_ID, & - penaltyenergy_ID, & - volumediscrepancy_ID, & - averagerelaxrate_ID,& - maximumrelaxrate_ID,& - magnitudemismatch_ID - end enum - type :: tParameters integer, dimension(:), allocatable :: & Nconstituents @@ -31,8 +20,8 @@ submodule(homogenization) homogenization_mech_RGC angles integer :: & of_debug = 0 - integer(kind(undefined_ID)), dimension(:), allocatable :: & - outputID + character(len=pStringLen), allocatable, dimension(:) :: & + output end type tParameters type :: tRGCstate @@ -71,23 +60,18 @@ module subroutine mech_RGC_init integer :: & Ninstance, & - h, i, & + h, & NofMyHomog, & sizeState, nIntFaceTot - integer(kind(undefined_ID)) :: & - outputID - - character(len=pStringLen), dimension(:), allocatable :: & - outputs - write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_RGC_label//' init -+>>>' + write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_RGC_label//' init -+>>>'; flush(6) - write(6,'(/,a)') ' Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009' - write(6,'(a)') ' https://doi.org/10.1007/s12289-009-0619-1' + write(6,'(/,a)') ' Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009' + write(6,'(a)') ' https://doi.org/10.1007/s12289-009-0619-1' - write(6,'(/,a)') ' Tjahjanto et al., Modelling and Simulation in Materials Science and Engineering 18:015006, 2010' - write(6,'(a)') ' https://doi.org/10.1088/0965-0393/18/1/015006' + write(6,'(/,a)') ' Tjahjanto et al., Modelling and Simulation in Materials Science and Engineering 18:015006, 2010' + write(6,'(a)') ' https://doi.org/10.1088/0965-0393/18/1/015006' Ninstance = count(homogenization_type == HOMOGENIZATION_RGC_ID) if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0) & @@ -123,34 +107,8 @@ module subroutine mech_RGC_init prm%dAlpha = config%getFloats('grainsize', requiredSize=3) prm%angles = config%getFloats('clusterorientation',requiredSize=3) - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - - do i=1, size(outputs) - outputID = undefined_ID - select case(outputs(i)) - - case('constitutivework') - outputID = constitutivework_ID - case('penaltyenergy') - outputID = penaltyenergy_ID - case('volumediscrepancy') - outputID = volumediscrepancy_ID - case('averagerelaxrate') - outputID = averagerelaxrate_ID - case('maximumrelaxrate') - outputID = maximumrelaxrate_ID - case('magnitudemismatch') - outputID = magnitudemismatch_ID - - end select - - if (outputID /= undefined_ID) then - prm%outputID = [prm%outputID , outputID] - endif - - enddo - + prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) + NofMyHomog = count(material_homogenizationAt == h) nIntFaceTot = 3*( (prm%Nconstituents(1)-1)*prm%Nconstituents(2)*prm%Nconstituents(3) & + prm%Nconstituents(1)*(prm%Nconstituents(2)-1)*prm%Nconstituents(3) & @@ -934,26 +892,24 @@ module subroutine mech_RGC_results(instance,group) integer :: o associate(stt => state(instance), dst => dependentState(instance), prm => param(instance)) - - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - - case (constitutivework_ID) + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case('constitutivework') call results_writeDataset(group,stt%work,'W',& 'work density','J/m³') - case (magnitudemismatch_ID) + case('magnitudemismatch') call results_writeDataset(group,dst%mismatch,'N',& 'average mismatch tensor','1') - case (penaltyenergy_ID) + case('penaltyenergy') call results_writeDataset(group,stt%penaltyEnergy,'R',& 'mismatch penalty density','J/m³') - case (volumediscrepancy_ID) + case('volumediscrepancy') call results_writeDataset(group,dst%volumeDiscrepancy,'Delta_V',& 'volume discrepancy','m³') - case (maximumrelaxrate_ID) + case('maximumrelaxrate') call results_writeDataset(group,dst%relaxationrate_max,'max_alpha_dot',& 'maximum relaxation rate','m/s') - case (averagerelaxrate_ID) + case('averagerelaxrate') call results_writeDataset(group,dst%relaxationrate_avg,'avg_alpha_dot',& 'average relaxation rate','m/s') end select diff --git a/src/thermal_adiabatic.f90 b/src/thermal_adiabatic.f90 index 5fa6b3281..2bef1c3ee 100644 --- a/src/thermal_adiabatic.f90 +++ b/src/thermal_adiabatic.f90 @@ -16,14 +16,9 @@ module thermal_adiabatic implicit none private - enum, bind(c) - enumerator :: undefined_ID, & - temperature_ID - end enum - type :: tParameters - integer(kind(undefined_ID)), dimension(:), allocatable :: & - outputID + character(len=pStringLen), allocatable, dimension(:) :: & + output end type tParameters type(tparameters), dimension(:), allocatable :: & @@ -46,10 +41,9 @@ contains !-------------------------------------------------------------------------------------------------- subroutine thermal_adiabatic_init - integer :: maxNinstance,o,h,NofMyHomog - character(len=pStringLen), dimension(:), allocatable :: outputs + integer :: maxNinstance,h,NofMyHomog - write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_ADIABATIC_label//' init -+>>>'; flush(6) + write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_ADIABATIC_label//' init -+>>>'; flush(6) maxNinstance = count(thermal_type == THERMAL_adiabatic_ID) if (maxNinstance == 0) return @@ -60,15 +54,7 @@ subroutine thermal_adiabatic_init if (thermal_type(h) /= THERMAL_adiabatic_ID) cycle associate(prm => param(thermal_typeInstance(h)),config => config_homogenization(h)) - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - - do o=1, size(outputs) - select case(outputs(o)) - case('temperature') - prm%outputID = [prm%outputID, temperature_ID] - end select - enddo + prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) NofMyHomog=count(material_homogenizationAt==h) thermalState(h)%sizeState = 1 @@ -76,7 +62,6 @@ subroutine thermal_adiabatic_init allocate(thermalState(h)%subState0(1,NofMyHomog), source=thermal_initialT(h)) allocate(thermalState(h)%state (1,NofMyHomog), source=thermal_initialT(h)) - nullify(thermalMapping(h)%p) thermalMapping(h)%p => material_homogenizationMemberAt deallocate(temperature(h)%p) temperature(h)%p => thermalState(h)%state(1,:) @@ -246,14 +231,13 @@ subroutine thermal_adiabatic_results(homog,group) integer, intent(in) :: homog character(len=*), intent(in) :: group + integer :: o associate(prm => param(damage_typeInstance(homog))) - - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - - case (temperature_ID) + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case('temperature') ! ToDo: should be 'T' call results_writeDataset(group,temperature(homog)%p,'T',& 'temperature','K') end select diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index 9ec8082ec..dad3fecd8 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -15,15 +15,9 @@ module thermal_conduction implicit none private - enum, bind(c) - enumerator :: & - undefined_ID, & - temperature_ID - end enum - type :: tParameters - integer(kind(undefined_ID)), dimension(:), allocatable :: & - outputID + character(len=pStringLen), allocatable, dimension(:) :: & + output end type tParameters type(tparameters), dimension(:), allocatable :: & @@ -48,10 +42,9 @@ contains subroutine thermal_conduction_init - integer :: maxNinstance,o,NofMyHomog,h - character(len=pStringLen), dimension(:), allocatable :: outputs + integer :: maxNinstance,NofMyHomog,h - write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_CONDUCTION_label//' init -+>>>'; flush(6) + write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_CONDUCTION_label//' init -+>>>'; flush(6) maxNinstance = count(thermal_type == THERMAL_conduction_ID) if (maxNinstance == 0) return @@ -62,15 +55,7 @@ subroutine thermal_conduction_init if (thermal_type(h) /= THERMAL_conduction_ID) cycle associate(prm => param(thermal_typeInstance(h)),config => config_homogenization(h)) - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - - do o=1, size(outputs) - select case(outputs(o)) - case('temperature') - prm%outputID = [prm%outputID, temperature_ID] - end select - enddo + prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) NofMyHomog=count(material_homogenizationAt==h) thermalState(h)%sizeState = 0 @@ -78,7 +63,6 @@ subroutine thermal_conduction_init allocate(thermalState(h)%subState0(0,NofMyHomog)) allocate(thermalState(h)%state (0,NofMyHomog)) - nullify(thermalMapping(h)%p) thermalMapping(h)%p => material_homogenizationMemberAt deallocate(temperature (h)%p) allocate (temperature (h)%p(NofMyHomog), source=thermal_initialT(h)) @@ -259,14 +243,13 @@ subroutine thermal_conduction_results(homog,group) integer, intent(in) :: homog character(len=*), intent(in) :: group + integer :: o associate(prm => param(damage_typeInstance(homog))) - - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - - case (temperature_ID) + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case('temperature') ! ToDo: should be 'T' call results_writeDataset(group,temperature(homog)%p,'T',& 'temperature','K') end select