223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge remote-tracking branch 'origin/development' into general-N_constituents

This commit is contained in:
Martin Diehl 2020-10-12 05:27:11 +02:00
parent 2b18161909 c8a9bc60d3
commit 72ba4645cb
24 changed files with 832 additions and 627 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitattributes vendored
View File

@ -8,6 +8,7 @@
*.jpg binary
*.hdf5 binary
*.pdf binary
*.dream3d binary
# ignore files from MSC.Marc in language statistics
installation/mods_MarcMentat/20*/* linguist-vendored

12
.gitlab-ci.yml
View File

@ -193,6 +193,8 @@ grid_mech_compile_Intel:
- module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel
- cp -r grid_mech_compile grid_mech_compile_Intel
- grid_mech_compile_Intel/test.py
- cd pytest
- pytest -k 'compile and grid'
except:
- master
- release
@ -203,6 +205,8 @@ Compile_FEM_Intel:
- module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel
- cp -r FEM_compile FEM_compile_Intel
- FEM_compile_Intel/test.py
- cd pytest
- pytest -k 'compile and mesh'
except:
- master
- release
@ -213,6 +217,8 @@ grid_mech_compile_GNU:
- module load $GNUCompiler $MPICH_GNU $PETSc_MPICH_GNU
- cp -r grid_mech_compile grid_mech_compile_GNU
- grid_mech_compile_GNU/test.py
- cd pytest
- pytest -k 'compile and grid'
except:
- master
- release
@ -223,6 +229,8 @@ Compile_FEM_GNU:
- module load $GNUCompiler $MPICH_GNU $PETSc_MPICH_GNU
- cp -r FEM_compile FEM_compile_GNU
- FEM_compile_GNU/test.py
- cd pytest
- pytest -k 'compile and mesh'
except:
- master
- release
@ -244,7 +252,7 @@ Pytest_grid:
script:
- module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel
- cd pytest
- pytest
- pytest -m 'not compile'
except:
- master
- release
@ -320,6 +328,8 @@ Marc_compileIfort:
script:
- module load $IntelMarc $HDF5Marc $MSC
- Marc_compileIfort/test.py
- cd pytest
- pytest -k 'compile and Marc'
except:
- master
- release

2
VERSION
View File

@ -1 +1 @@
v3.0.0-alpha-452-gd3f068cd7
v3.0.0-alpha-502-g8ab405e5d

1
examples/SpectralMethod/Polycrystal/material.yaml
View File

@ -1,3 +1,4 @@
---
homogenization:
SX:
N_constituents: 1

76
processing/post/addGaussian.py
View File

@ -1,76 +0,0 @@
#!/usr/bin/env python3
import os
import sys
from io import StringIO
from optparse import OptionParser
from scipy import ndimage
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog option [ASCIItable(s)]', description = """
Add column(s) containing Gaussian filtered values of requested column(s).
Operates on periodic and non-periodic ordered three-dimensional data sets.
For details see scipy.ndimage documentation.
""", version = scriptID)
parser.add_option('-p','--pos','--periodiccellcenter',
dest = 'pos',
type = 'string', metavar = 'string',
help = 'label of coordinates [%default]')
parser.add_option('-s','--scalar',
dest = 'labels',
action = 'extend', metavar = '<string LIST>',
help = 'label(s) of scalar field values')
parser.add_option('-o','--order',
dest = 'order',
type = int,
metavar = 'int',
help = 'order of the filter [%default]')
parser.add_option('--sigma',
dest = 'sigma',
type = float,
metavar = 'float',
help = 'standard deviation [%default]')
parser.add_option('--periodic',
dest = 'periodic',
action = 'store_true',
help = 'assume periodic grain structure')
parser.set_defaults(pos = 'pos',
order = 0,
sigma = 1,
)
(options,filenames) = parser.parse_args()
if filenames == []: filenames = [None]
if options.labels is None: parser.error('no data column specified.')
for name in filenames:
damask.util.report(scriptName,name)
table = damask.Table.load(StringIO(''.join(sys.stdin.read())) if name is None else name)
damask.grid_filters.coord0_check(table.get(options.pos))
for label in options.labels:
table = table.add('Gauss{}({})'.format(options.sigma,label),
ndimage.filters.gaussian_filter(table.get(label).reshape(-1),
options.sigma,options.order,
mode = 'wrap' if options.periodic else 'nearest'),
scriptID+' '+' '.join(sys.argv[1:]))
table.save((sys.stdout if name is None else name), legacy=True)

82
processing/pre/geom_fromDREAM3D.py
View File

@ -1,12 +1,8 @@
#!/usr/bin/env python3
import os
import sys
from optparse import OptionParser
import h5py
import numpy as np
import damask
@ -64,88 +60,12 @@ parser.set_defaults(pointwise = 'CellData',
if options.basegroup is None:
parser.error('No base group selected')
rootDir ='DataContainers'
if filenames == []: parser.error('no input file specified.')
for name in filenames:
damask.util.report(scriptName,name)
errors = []
inFile = h5py.File(name, 'r')
group_geom = os.path.join(rootDir,options.basegroup,'_SIMPL_GEOMETRY')
try:
size = inFile[os.path.join(group_geom,'DIMENSIONS')][...] \
* inFile[os.path.join(group_geom,'SPACING')][...]
grid = inFile[os.path.join(group_geom,'DIMENSIONS')][...]
origin = inFile[os.path.join(group_geom,'ORIGIN')][...]
except KeyError:
errors.append('Geometry data ({}) not found'.format(group_geom))
group_pointwise = os.path.join(rootDir,options.basegroup,options.pointwise)
if options.average is None:
label = 'Point'
dataset = os.path.join(group_pointwise,options.quaternion)
try:
quats = np.reshape(inFile[dataset][...],(np.product(grid),4))
rot = [damask.Rotation.from_quaternion(q,True,P=+1) for q in quats]
except KeyError:
errors.append('Pointwise orientation (quaternion) data ({}) not readable'.format(dataset))
dataset = os.path.join(group_pointwise,options.phase)
try:
phase = np.reshape(inFile[dataset][...],(np.product(grid)))
except KeyError:
errors.append('Pointwise phase data ({}) not readable'.format(dataset))
microstructure = np.arange(1,np.product(grid)+1,dtype=int).reshape(grid,order='F')
else:
label = 'Grain'
dataset = os.path.join(group_pointwise,options.microstructure)
try:
microstructure = np.transpose(inFile[dataset][...].reshape(grid[::-1]),(2,1,0)) # convert from C ordering
except KeyError:
errors.append('Link between pointwise and grain average data ({}) not readable'.format(dataset))
group_average = os.path.join(rootDir,options.basegroup,options.average)
dataset = os.path.join(group_average,options.quaternion)
try:
rot = [damask.Rotation.from_quaternion(q,True,P=+1) for q in inFile[dataset][...][1:]] # skip first entry (unindexed)
except KeyError:
errors.append('Average orientation data ({}) not readable'.format(dataset))
dataset = os.path.join(group_average,options.phase)
try:
phase = [i[0] for i in inFile[dataset][...]][1:] # skip first entry (unindexed)
except KeyError:
errors.append('Average phase data ({}) not readable'.format(dataset))
if errors != []:
damask.util.croak(errors)
continue
config_header = ['<texture>']
for i in range(np.nanmax(microstructure)):
config_header += ['[{}{}]'.format(label,i+1),
'(gauss)\tphi1 {:.2f}\tPhi {:.2f}\tphi2 {:.2f}'.format(*rot[i].as_Eulers(degrees = True)),
]
config_header += ['<microstructure>']
for i in range(np.nanmax(microstructure)):
config_header += ['[{}{}]'.format(label,i+1),
'(constituent)\tphase {}\ttexture {}\tfraction 1.0'.format(phase[i],i+1),
]
header = [scriptID + ' ' + ' '.join(sys.argv[1:])]\
+ config_header
geom = damask.Geom(microstructure,size,origin,comments=header)
geom = damask.Geom.load_DREAM3D(name,options.basegroup,options.pointwise)
damask.util.croak(geom)
geom.save_ASCII(os.path.splitext(name)[0]+'.geom',compress=False)

58
processing/pre/geom_fromTable.py
View File

@ -2,11 +2,8 @@
import os
import sys
from io import StringIO
from optparse import OptionParser
import numpy as np
import damask
@ -40,64 +37,21 @@ parser.add_option('-q', '--quaternion',
dest = 'quaternion',
type = 'string', metavar='string',
help = 'quaternion label')
parser.add_option('--axes',
dest = 'axes',
type = 'string', nargs = 3, metavar = ' '.join(['string']*3),
help = 'orientation coordinate frame in terms of position coordinate frame [+x +y +z]')
parser.set_defaults(pos = 'pos',
)
parser.set_defaults(pos= 'pos')
(options,filenames) = parser.parse_args()
if filenames == []: filenames = [None]
if np.sum([options.quaternion is not None,
options.microstructure is not None]) != 1:
parser.error('need either microstructure or quaternion (and optionally phase) as input.')
if options.microstructure is not None and options.phase is not None:
parser.error('need either microstructure or phase (and mandatory quaternion) as input.')
if options.axes is not None and not set(options.axes).issubset(set(['x','+x','-x','y','+y','-y','z','+z','-z'])):
parser.error('invalid axes {} {} {}.'.format(*options.axes))
for name in filenames:
damask.util.report(scriptName,name)
table = damask.Table.load(StringIO(''.join(sys.stdin.read())) if name is None else name)
table.sort_by(['{}_{}'.format(i,options.pos) for i in range(3,0,-1)]) # x fast, y slow
grid,size,origin = damask.grid_filters.cell_coord0_gridSizeOrigin(table.get(options.pos))
labels = []
for l in [options.quaternion,options.phase,options.microstructure]:
if l is not None: labels.append(l)
config_header = table.comments
if options.microstructure:
microstructure = table.get(options.microstructure).reshape(grid,order='F')
elif options.quaternion:
q = table.get(options.quaternion)
phase = table.get(options.phase).astype(int) if options.phase else \
np.ones((table.data.shape[0],1),dtype=int)
unique,unique_inverse = np.unique(np.hstack((q,phase)),return_inverse=True,axis=0)
microstructure = unique_inverse.reshape(grid,order='F') + 1
config_header = ['<texture>']
for i,data in enumerate(unique):
ori = damask.Rotation(data[0:4])
config_header += ['[Grain{}]'.format(i+1),
'(gauss)\tphi1 {:.2f}\tPhi {:.2f}\tphi2 {:.2f}'.format(*ori.as_Eulers(degrees = True)),
]
if options.axes is not None: config_header += ['axes\t{} {} {}'.format(*options.axes)]
config_header += ['<microstructure>']
for i,data in enumerate(unique):
config_header += ['[Grain{}]'.format(i+1),
'(constituent)\tphase {}\ttexture {}\tfraction 1.0'.format(int(data[4]),i+1),
]
header = [scriptID + ' ' + ' '.join(sys.argv[1:])]\
+ config_header
geom = damask.Geom(microstructure,size,origin,
comments=header)
t = damask.Table.load(name)
geom = damask.Geom.from_table(t,options.pos,labels)
damask.util.croak(geom)
geom.save_ASCII(sys.stdout if name is None else os.path.splitext(name)[0]+'.geom',compress=False)

23
python/damask/__init__.py
View File

@ -10,20 +10,21 @@ with open(_Path(__file__).parent/_Path('VERSION')) as _f:
# make classes directly accessible as damask.Class
from ._environment import Environment as _ # noqa
environment = _()
from ._table import Table # noqa
from ._vtk import VTK # noqa
from ._colormap import Colormap # noqa
from ._rotation import Rotation # noqa
from ._lattice import Symmetry, Lattice# noqa
from ._orientation import Orientation # noqa
from ._result import Result # noqa
from ._geom import Geom # noqa
from ._material import Material # noqa
from . import solver # noqa
from . import util # noqa
from . import seeds # noqa
from . import grid_filters # noqa
from . import mechanics # noqa
from . import solver # noqa
from . import grid_filters # noqa
from ._lattice import Symmetry, Lattice# noqa
from ._table import Table # noqa
from ._rotation import Rotation # noqa
from ._vtk import VTK # noqa
from ._colormap import Colormap # noqa
from ._orientation import Orientation # noqa
from ._config import Config # noqa
from ._configmaterial import ConfigMaterial # noqa
from ._geom import Geom # noqa
from ._result import Result # noqa

80
python/damask/_config.py Normal file
View File

@ -0,0 +1,80 @@
from io import StringIO
import abc
import yaml
class NiceDumper(yaml.SafeDumper):
"""Make YAML readable for humans."""
def write_line_break(self, data=None):
super().write_line_break(data)
if len(self.indents) == 1:
super().write_line_break()
def increase_indent(self, flow=False, indentless=False):
return super().increase_indent(flow, False)
class Config(dict):
"""YAML-based configuration."""
def __repr__(self):
"""Show as in file."""
output = StringIO()
self.save(output)
output.seek(0)
return ''.join(output.readlines())
@classmethod
def load(cls,fname):
"""
Load from yaml file.
Parameters
----------
fname : file, str, or pathlib.Path
Filename or file for writing.
"""
try:
fhandle = open(fname)
except TypeError:
fhandle = fname
return cls(yaml.safe_load(fhandle))
def save(self,fname,**kwargs):
"""
Save to yaml file.
Parameters
----------
fname : file, str, or pathlib.Path
Filename or file for writing.
**kwargs : dict
Keyword arguments parsed to yaml.dump.
"""
try:
fhandle = open(fname,'w')
except TypeError:
fhandle = fname
if 'width' not in kwargs:
kwargs['width'] = 256
if 'default_flow_style' not in kwargs:
kwargs['default_flow_style'] = None
fhandle.write(yaml.dump(dict(self),Dumper=NiceDumper,**kwargs))
@property
@abc.abstractmethod
def is_complete(self):
"""Check for completeness."""
pass
@property
@abc.abstractmethod
def is_valid(self):
"""Check for valid file layout."""
pass

299
python/damask/_configmaterial.py Normal file
View File

@ -0,0 +1,299 @@
import copy
import numpy as np
from . import grid_filters
from . import Config
from . import Lattice
from . import Rotation
class ConfigMaterial(Config):
"""Material configuration."""
def save(self,fname='material.yaml',**kwargs):
"""
Save to yaml file.
Parameters
----------
fname : file, str, or pathlib.Path, optional
Filename or file for writing. Defaults to 'material.yaml'.
**kwargs
Keyword arguments parsed to yaml.dump.
"""
super().save(fname,**kwargs)
@staticmethod
def from_table(table,coordinates=None,constituents={},**kwargs):
"""
Load from an ASCII table.
Parameters
----------
table : damask.Table
Table that contains material information.
coordinates : str, optional
Label of spatial coordiates. Used for sorting and performing a
sanity check. Default to None, in which case no sorting or checking is
peformed.
constituents : dict, optional
Entries for 'constituents'. The key is the name and the value specifies
the label of the data column in the table
**kwargs
Keyword arguments where the key is the name and the value specifies
the label of the data column in the table
Examples
--------
>>> import damask
>>> import damask.ConfigMaterial as cm
>>> t = damask.Table.load('small.txt')
>>> t
pos pos pos qu qu qu qu phase homog
0 0 0 0 0.19 0.8 0.24 -0.51 Aluminum SX
1 1 0 0 0.8 0.19 0.24 -0.51 Steel SX
>>> cm.from_table(t,'pos',{'O':'qu','phase':'phase'},homogenization='homog')
material:
- constituents:
- O: [0.19, 0.8, 0.24, -0.51]
fraction: 1.0
phase: Aluminum
homogenization: SX
- constituents:
- O: [0.8, 0.19, 0.24, -0.51]
fraction: 1.0
phase: Steel
homogenization: SX
"""
if coordinates is not None:
t = table.sort_by([f'{i}_{coordinates}' for i in range(3,0,-1)])
grid_filters.coord0_check(t.get(coordinates))
else:
t = table
constituents_ = {k:t.get(v) for k,v in constituents.items()}
kwargs_ = {k:t.get(v) for k,v in kwargs.items()}
_,idx = np.unique(np.hstack(list({**constituents_,**kwargs_}.values())),return_index=True,axis=0)
constituents_ = {k:v[idx].squeeze() for k,v in constituents_.items()}
kwargs_ = {k:v[idx].squeeze() for k,v in kwargs_.items()}
return ConfigMaterial().material_add(constituents_,**kwargs_)
@property
def is_complete(self):
"""Check for completeness."""
ok = True
for top_level in ['homogenization','phase','material']:
# ToDo: With python 3.8 as prerequisite we can shorten with :=
ok &= top_level in self
if top_level not in self: print(f'{top_level} entry missing')
if ok:
ok &= len(self['material']) > 0
if len(self['material']) < 1: print('Incomplete material definition')
if ok:
homogenization = set()
phase = set()
for i,v in enumerate(self['material']):
if 'homogenization' in v:
homogenization.add(v['homogenization'])
else:
print(f'No homogenization specified in material {i}')
ok = False
if 'constituents' in v:
for ii,vv in enumerate(v['constituents']):
if 'O' not in vv:
print('No orientation specified in constituent {ii} of material {i}')
ok = False
if 'phase' in vv:
phase.add(vv['phase'])
else:
print(f'No phase specified in constituent {ii} of material {i}')
ok = False
for k,v in self['phase'].items():
if 'lattice' not in v:
print(f'No lattice specified in phase {k}')
ok = False
for k,v in self['homogenization'].items():
if 'N_constituents' not in v:
print(f'No. of constituents not specified in homogenization {k}')
ok = False
if phase - set(self['phase']):
print(f'Phase(s) {phase-set(self["phase"])} missing')
ok = False
if homogenization - set(self['homogenization']):
print(f'Homogenization(s) {homogenization-set(self["homogenization"])} missing')
ok = False
return ok
@property
def is_valid(self):
"""Check for valid file layout."""
ok = True
if 'phase' in self:
for k,v in self['phase'].items():
if 'lattice' in v:
try:
Lattice(v['lattice'])
except KeyError:
s = v['lattice']
print(f"Invalid lattice: '{s}' in phase '{k}'")
ok = False
if 'material' in self:
for i,v in enumerate(self['material']):
if 'constituents' in v:
f = 0.0
for c in v['constituents']:
f+= float(c['fraction'])
if 'O' in c:
try:
Rotation.from_quaternion(c['O'])
except ValueError:
o = c['O']
print(f"Invalid orientation: '{o}' in material '{i}'")
ok = False
if not np.isclose(f,1.0):
print(f"Invalid total fraction '{f}' in material '{i}'")
ok = False
return ok
def material_rename_phase(self,mapping,ID=None,constituent=None):
"""
Change phase name in material.
Parameters
----------
mapping: dictionary
Mapping from old name to new name
ID: list of ints, optional
Limit renaming to selected material IDs.
constituent: list of ints, optional
Limit renaming to selected constituents.
"""
dup = copy.deepcopy(self)
for i,m in enumerate(dup['material']):
if ID and i not in ID: continue
for c in m['constituents']:
if constituent is not None and c not in constituent: continue
try:
c['phase'] = mapping[c['phase']]
except KeyError:
continue
return dup
def material_rename_homogenization(self,mapping,ID=None):
"""
Change homogenization name in material.
Parameters
----------
mapping: dictionary
Mapping from old name to new name
ID: list of ints, optional
Limit renaming to selected homogenization IDs.
"""
dup = copy.deepcopy(self)
for i,m in enumerate(dup['material']):
if ID and i not in ID: continue
try:
m['homogenization'] = mapping[m['homogenization']]
except KeyError:
continue
return dup
def material_add(self,constituents,**kwargs):
"""
Add material entries.
Parameters
----------
constituents : dict
Entries for 'constituents'. The key is the name and the value specifies
the label of the data column in the table
**kwargs
Keyword arguments where the key is the name and the value specifies
the label of the data column in the table
Examples
--------
>>> import damask
>>> m = damask.ConfigMaterial()
>>> O = damask.Rotation.from_random(3).as_quaternion()
>>> phase = ['Aluminum','Steel','Aluminum']
>>> m.material_add(constituents={'phase':phase,'O':O},homogenization='SX')
material:
- constituents:
- O: [0.577764, -0.146299, -0.617669, 0.513010]
fraction: 1.0
phase: Aluminum
homogenization: SX
- constituents:
- O: [0.184176, 0.340305, 0.737247, 0.553840]
fraction: 1.0
phase: Steel
homogenization: SX
- constituents:
- O: [0.0886257, -0.144848, 0.615674, -0.769487]
fraction: 1.0
phase: Aluminum
homogenization: SX
"""
c = [{'constituents':u} for u in ConfigMaterial._constituents(**constituents)]
for k,v in kwargs.items():
if hasattr(v,'__len__') and not isinstance(v,str):
if len(v) != len(c):
raise ValueError('len mismatch 1')
for i,vv in enumerate(v):
c[i][k] = [w.item() for w in vv] if isinstance(vv,np.ndarray) else vv.item()
else:
for i in range(len(c)):
c[i][k] = v
dup = copy.deepcopy(self)
if 'material' not in dup: dup['material'] = []
dup['material'] +=c
return dup
@staticmethod
def _constituents(N=1,**kwargs):
"""Construct list of constituents."""
for v in kwargs.values():
if hasattr(v,'__len__') and not isinstance(v,str): N_material = len(v)
if N == 1:
m = [[{'fraction':1.0}] for _ in range(N_material)]
for k,v in kwargs.items():
if hasattr(v,'__len__') and not isinstance(v,str):
if len(v) != N_material:
raise ValueError('len mismatch 2')
for i,vv in enumerate(np.array(v)):
m[i][0][k] = [w.item() for w in vv] if isinstance(vv,np.ndarray) else vv.item()
else:
for i in range(N_material):
m[i][0][k] = v
return m
else:
raise NotImplementedError

70
python/damask/_geom.py
View File

@ -1,15 +1,17 @@
import copy
import multiprocessing as mp
from functools import partial
from os import path
import numpy as np
import h5py
from scipy import ndimage,spatial
from . import environment
from . import Rotation
from . import VTK
from . import util
from . import grid_filters
from . import Rotation
class Geom:
@ -207,6 +209,68 @@ class Geom:
return Geom(material.reshape(grid,order='F'),size,origin,comments)
@staticmethod
def load_DREAM3D(fname,base_group,point_data=None,material='FeatureIds'):
"""
Load a DREAM.3D file.
Parameters
----------
fname : str
Filename of the DREAM.3D file
base_group : str
Name of the group (folder) below 'DataContainers'. For example
'SyntheticVolumeDataContainer'.
point_data : str, optional
Name of the group (folder) containing the point wise material data,
for example 'CellData'. Defaults to None, in which case points consecutively numbered.
material : str, optional
Name of the dataset containing the material ID. Defaults to
'FeatureIds'.
"""
root_dir ='DataContainers'
f = h5py.File(fname, 'r')
g = path.join(root_dir,base_group,'_SIMPL_GEOMETRY')
size = f[path.join(g,'DIMENSIONS')][()] * f[path.join(g,'SPACING')][()]
grid = f[path.join(g,'DIMENSIONS')][()]
origin = f[path.join(g,'ORIGIN')][()]
group_pointwise = path.join(root_dir,base_group,point_data)
ma = np.arange(1,np.product(grid)+1,dtype=int) if point_data is None else \
np.reshape(f[path.join(group_pointwise,material)],grid.prod())
return Geom(ma.reshape(grid,order='F'),size,origin,util.execution_stamp('Geom','load_DREAM3D'))
@staticmethod
def from_table(table,coordinates,labels):
"""
Load an ASCII table.
Parameters
----------
table : damask.Table
Table that contains material information.
coordinates : str
Label of the column containing the spatial coordinates.
labels : str or list of str
Label(s) of the columns containing the material definition.
Each unique combintation of values results in a material.
"""
t = table.sort_by([f'{i}_{coordinates}' for i in range(3,0,-1)])
grid,size,origin = grid_filters.cell_coord0_gridSizeOrigin(t.get(coordinates))
labels_ = [labels] if isinstance(labels,str) else labels
_,unique_inverse = np.unique(np.hstack([t.get(l) for l in labels_]),return_inverse=True,axis=0)
ma = unique_inverse.reshape(grid,order='F') + 1
return Geom(ma,size,origin,util.execution_stamp('Geom','from_table'))
@staticmethod
def _find_closest_seed(seeds, weights, point):
return np.argmin(np.sum((np.broadcast_to(point,(len(seeds),3))-seeds)**2,axis=1) - weights)
@ -399,7 +463,7 @@ class Geom:
def save(self,fname,compress=True):
"""
Generates vtk rectilinear grid.
Store as vtk rectilinear grid.
Parameters
----------
@ -536,7 +600,7 @@ class Geom:
coords_rot = R.broadcast_to(tuple(self.grid))@coords
with np.errstate(all='ignore'):
mask = np.where(np.sum(np.power(coords_rot/r,2.0**exponent),axis=-1) > 1.0,True,False)
mask = np.sum(np.power(coords_rot/r,2.0**np.array(exponent)),axis=-1) > 1.0
if periodic: # translate back to center
mask = np.roll(mask,((c-np.ones(3)*.5)*self.grid).astype(int),(0,1,2))

193
python/damask/_material.py
View File

@ -1,193 +0,0 @@
from io import StringIO
import copy
import yaml
import numpy as np
from . import Lattice
from . import Rotation
class NiceDumper(yaml.SafeDumper):
"""Make YAML readable for humans."""
def write_line_break(self, data=None):
super().write_line_break(data)
if len(self.indents) == 1:
super().write_line_break()
def increase_indent(self, flow=False, indentless=False):
return super().increase_indent(flow, False)
class Material(dict):
"""Material configuration."""
def __repr__(self):
"""Show as in file."""
output = StringIO()
self.save(output)
output.seek(0)
return ''.join(output.readlines())
@staticmethod
def load(fname):
"""Load from yaml file."""
try:
fhandle = open(fname)
except TypeError:
fhandle = fname
return Material(yaml.safe_load(fhandle))
def save(self,fname='material.yaml'):
"""
Save to yaml file.
Parameters
----------
fname : file, str, or pathlib.Path
Filename or file for reading.
"""
try:
fhandle = open(fname,'w')
except TypeError:
fhandle = fname
fhandle.write(yaml.dump(dict(self),width=256,default_flow_style=None,Dumper=NiceDumper))
@property
def is_complete(self):
"""Check for completeness."""
ok = True
for top_level in ['homogenization','phase','microstructure']:
# ToDo: With python 3.8 as prerequisite we can shorten with :=
ok &= top_level in self
if top_level not in self: print(f'{top_level} entry missing')
if ok:
ok &= len(self['microstructure']) > 0
if len(self['microstructure']) < 1: print('Incomplete microstructure definition')
if ok:
homogenization = set()
phase = set()
for i,v in enumerate(self['microstructure']):
if 'homogenization' in v:
homogenization.add(v['homogenization'])
else:
print(f'No homogenization specified in microstructure {i}')
ok = False
if 'constituents' in v:
for ii,vv in enumerate(v['constituents']):
if 'orientation' not in vv:
print('No orientation specified in constituent {ii} of microstructure {i}')
ok = False
if 'phase' in vv:
phase.add(vv['phase'])
else:
print(f'No phase specified in constituent {ii} of microstructure {i}')
ok = False
for k,v in self['phase'].items():
if 'lattice' not in v:
print(f'No lattice specified in phase {k}')
ok = False
#for k,v in self['homogenization'].items():
# if 'N_constituents' not in v:
# print(f'No. of constituents not specified in homogenization {k}'}
# ok = False
if phase - set(self['phase']):
print(f'Phase(s) {phase-set(self["phase"])} missing')
ok = False
if homogenization - set(self['homogenization']):
print(f'Homogenization(s) {homogenization-set(self["homogenization"])} missing')
ok = False
return ok
@property
def is_valid(self):
"""Check for valid file layout."""
ok = True
if 'phase' in self:
for k,v in self['phase'].items():
if 'lattice' in v:
try:
Lattice(v['lattice'])
except KeyError:
s = v['lattice']
print(f"Invalid lattice: '{s}' in phase '{k}'")
ok = False
if 'microstructure' in self:
for i,v in enumerate(self['microstructure']):
if 'constituents' in v:
f = 0.0
for c in v['constituents']:
f+= float(c['fraction'])
if 'orientation' in c:
try:
Rotation.from_quaternion(c['orientation'])
except ValueError:
o = c['orientation']
print(f"Invalid orientation: '{o}' in microstructure '{i}'")
ok = False
if not np.isclose(f,1.0):
print(f"Invalid total fraction '{f}' in microstructure '{i}'")
ok = False
return ok
def microstructure_rename_phase(self,mapping,ID=None,constituent=None):
"""
Change phase name in microstructure.
Parameters
----------
mapping: dictionary
Mapping from old name to new name
ID: list of ints, optional
Limit renaming to selected microstructure IDs.
constituent: list of ints, optional
Limit renaming to selected constituents.
"""
dup = copy.deepcopy(self)
for i,m in enumerate(dup['microstructure']):
if ID and i not in ID: continue
for c in m['constituents']:
if constituent is not None and c not in constituent: continue
try:
c['phase'] = mapping[c['phase']]
except KeyError:
continue
return dup
def microstructure_rename_homogenization(self,mapping,ID=None):
"""
Change homogenization name in microstructure.
Parameters
----------
mapping: dictionary
Mapping from old name to new name
ID: list of ints, optional
Limit renaming to selected homogenization IDs.
"""
dup = copy.deepcopy(self)
for i,m in enumerate(dup['microstructure']):
if ID and i not in ID: continue
try:
m['homogenization'] = mapping[m['homogenization']]
except KeyError:
continue
return dup

4
python/damask/_table.py
View File

@ -33,6 +33,10 @@ class Table:
"""Brief overview."""
return util.srepr(self.comments)+'\n'+self.data.__repr__()
def __len__(self):
"""Number of rows."""
return len(self.data)
def __copy__(self):
"""Copy Table."""
return copy.deepcopy(self)

24
python/tests/reference/Material/material.yaml → python/tests/reference/ConfigMaterial/material.yaml
View File

@ -1,33 +1,35 @@
homogenization:
SX:
N_constituents: 2
mech: {type: none}
Taylor:
mech: {type: isostrain, N_constituents: 2}
N_constituents: 2
mech: {type: isostrain}
microstructure:
material:
- constituents:
- fraction: 1.0
orientation: [1.0, 0.0, 0.0, 0.0]
O: [1.0, 0.0, 0.0, 0.0]
phase: Aluminum
homogenization: SX
- constituents:
- fraction: 1.0
orientation: [0.7936696712125002, -0.28765777461664166, -0.3436487135089419, 0.4113964260949434]
O: [0.7936696712125002, -0.28765777461664166, -0.3436487135089419, 0.4113964260949434]
phase: Aluminum
homogenization: SX
- constituents:
- fraction: 1.0
orientation: [0.3986143167493579, -0.7014883552495493, 0.2154871765709027, 0.5500781677772945]
O: [0.3986143167493579, -0.7014883552495493, 0.2154871765709027, 0.5500781677772945]
phase: Aluminum
homogenization: SX
- homogenization: Taylor
constituents:
- fraction: .5
orientation: [0.28645844315788244, -0.022571491243423537, -0.467933059311115, -0.8357456192708106]
- constituents:
- fraction: 0.5
O: [0.28645844315788244, -0.022571491243423537, -0.467933059311115, -0.8357456192708106]
phase: Aluminum
- fraction: .5
orientation: [0.3986143167493579, -0.7014883552495493, 0.2154871765709027, 0.5500781677772945]
- fraction: 0.5
O: [0.3986143167493579, -0.7014883552495493, 0.2154871765709027, 0.5500781677772945]
phase: Steel
homogenization: Taylor
phase:
Aluminum:

37
python/tests/test_Config.py Normal file
View File

@ -0,0 +1,37 @@
import pytest
from damask import Config
class TestConfig:
@pytest.mark.parametrize('flow_style',[None,True,False])
def test_load_save_str(self,tmp_path,flow_style):
config = Config()
config['A'] = 1
config['B'] = [2,3]
config.save(tmp_path/'config.yaml',default_flow_style=flow_style)
assert Config.load(tmp_path/'config.yaml') == config
def test_load_save_file(self,tmp_path):
config = Config()
config['A'] = 1
config['B'] = [2,3]
with open(tmp_path/'config.yaml','w') as f:
config.save(f)
with open(tmp_path/'config.yaml') as f:
assert Config.load(f) == config
def test_repr(self,tmp_path):
config = Config()
config['A'] = 1
config['B'] = [2,3]
with open(tmp_path/'config.yaml','w') as f:
f.write(config.__repr__())
assert Config.load(tmp_path/'config.yaml') == config
def test_abstract_is_valid(self):
assert Config().is_valid is None
def test_abstract_is_complete(self):
assert Config().is_complete is None

76
python/tests/test_ConfigMaterial.py Normal file
View File

@ -0,0 +1,76 @@
import os
import pytest
from damask import ConfigMaterial
@pytest.fixture
def reference_dir(reference_dir_base):
"""Directory containing reference results."""
return reference_dir_base/'ConfigMaterial'
class TestConfigMaterial:
@pytest.mark.parametrize('fname',[None,'test.yaml'])
def test_load_save(self,reference_dir,tmp_path,fname):
reference = ConfigMaterial.load(reference_dir/'material.yaml')
os.chdir(tmp_path)
if fname is None:
reference.save()
new = ConfigMaterial.load('material.yaml')
else:
reference.save(fname)
new = ConfigMaterial.load(fname)
assert reference == new
def test_valid_complete(self,reference_dir):
material_config = ConfigMaterial.load(reference_dir/'material.yaml')
assert material_config.is_valid and material_config.is_complete
def test_invalid_lattice(self,reference_dir):
material_config = ConfigMaterial.load(reference_dir/'material.yaml')
material_config['phase']['Aluminum']['lattice']='fxc'
assert not material_config.is_valid
def test_invalid_orientation(self,reference_dir):
material_config = ConfigMaterial.load(reference_dir/'material.yaml')
material_config['material'][0]['constituents'][0]['O']=[0,0,0,0]
assert not material_config.is_valid
def test_invalid_fraction(self,reference_dir):
material_config = ConfigMaterial.load(reference_dir/'material.yaml')
material_config['material'][0]['constituents'][0]['fraction']=.9
assert not material_config.is_valid
@pytest.mark.parametrize('item',['homogenization','phase','material'])
def test_incomplete_missing(self,reference_dir,item):
material_config = ConfigMaterial.load(reference_dir/'material.yaml')
del material_config[item]
assert not material_config.is_complete
@pytest.mark.parametrize('item',['O','phase'])
def test_incomplete_material_constituent(self,reference_dir,item):
material_config = ConfigMaterial.load(reference_dir/'material.yaml')
del material_config['material'][0]['constituents'][0][item]
assert not material_config.is_complete
def test_incomplete_material_homogenization(self,reference_dir):
material_config = ConfigMaterial.load(reference_dir/'material.yaml')
del material_config['material'][0]['homogenization']
assert not material_config.is_complete
def test_incomplete_phase_lattice(self,reference_dir):
material_config = ConfigMaterial.load(reference_dir/'material.yaml')
del material_config['phase']['Aluminum']['lattice']
assert not material_config.is_complete
def test_incomplete_wrong_phase(self,reference_dir):
material_config = ConfigMaterial.load(reference_dir/'material.yaml')
new = material_config.material_rename_phase({'Steel':'FeNbC'})
assert not new.is_complete
def test_incomplete_wrong_homogenization(self,reference_dir):
material_config = ConfigMaterial.load(reference_dir/'material.yaml')
new = material_config.material_rename_homogenization({'Taylor':'isostrain'})
assert not new.is_complete

61
python/tests/test_Material.py
View File

@ -1,61 +0,0 @@
import os
import pytest
from damask import Material
@pytest.fixture
def reference_dir(reference_dir_base):
"""Directory containing reference results."""
return reference_dir_base/'Material'
class TestMaterial:
@pytest.mark.parametrize('fname',[None,'test.yaml'])
def test_load_save(self,reference_dir,tmp_path,fname):
reference = Material.load(reference_dir/'material.yaml')
os.chdir(tmp_path)
if fname is None:
reference.save()
new = Material.load('material.yaml')
else:
reference.save(fname)
new = Material.load(fname)
assert reference == new
def test_valid_complete(self,reference_dir):
material_config = Material.load(reference_dir/'material.yaml')
assert material_config.is_valid and material_config.is_complete
def test_invalid_lattice(self,reference_dir):
material_config = Material.load(reference_dir/'material.yaml')
material_config['phase']['Aluminum']['lattice']='fxc'
assert not material_config.is_valid
def test_invalid_orientation(self,reference_dir):
material_config = Material.load(reference_dir/'material.yaml')
material_config['microstructure'][0]['constituents'][0]['orientation']=[0,0,0,0]
assert not material_config.is_valid
def test_invalid_fraction(self,reference_dir):
material_config = Material.load(reference_dir/'material.yaml')
material_config['microstructure'][0]['constituents'][0]['fraction']=.9
assert not material_config.is_valid
@pytest.mark.parametrize('item',['homogenization','phase','microstructure'])
def test_incomplete_missing(self,reference_dir,item):
material_config = Material.load(reference_dir/'material.yaml')
del material_config[item]
assert not material_config.is_complete
def test_incomplete_wrong_phase(self,reference_dir):
material_config = Material.load(reference_dir/'material.yaml')
new = material_config.microstructure_rename_phase({'Steel':'FeNbC'})
assert not new.is_complete
def test_incomplete_wrong_homogenization(self,reference_dir):
material_config = Material.load(reference_dir/'material.yaml')
new = material_config.microstructure_rename_homogenization({'Taylor':'isostrain'})
assert not new.is_complete

145
src/YAML_parse.f90
View File

@ -220,13 +220,26 @@ logical function isKey(line)
if(len(IO_rmComment(line)) == 0) then
isKey = .false.
else
isKey = IO_rmComment(line(len(IO_rmComment(line)):len(IO_rmComment(line)))) == ':' &
.and. .not. isFlow(line)
isKey = index(IO_rmComment(line),':',back=.false.) == len(IO_rmComment(line)) .and. &
index(IO_rmComment(line),':',back=.true.) == len(IO_rmComment(line)) .and. &
.not. isFlow(line)
endif
end function isKey
!--------------------------------------------------------------------------------------------------
! @brief check whether a string is a list in flow style
!--------------------------------------------------------------------------------------------------
logical function isFlowList(line)
character(len=*), intent(in) :: line
isFlowList = index(adjustl(line),'[') == 1
end function isFlowList
!--------------------------------------------------------------------------------------------------
! @brief skip empty lines
! @details update start position in the block by skipping empty lines if present.
@ -244,7 +257,6 @@ subroutine skip_empty_lines(blck,s_blck)
if(empty) s_blck = s_blck + index(blck(s_blck:),IO_EOL)
enddo
end subroutine skip_empty_lines
@ -273,6 +285,58 @@ subroutine skip_file_header(blck,s_blck)
end subroutine skip_file_header
!--------------------------------------------------------------------------------------------------
!> @brief check if a line in flow YAML starts and ends in the same line
!--------------------------------------------------------------------------------------------------
logical function flow_is_closed(str,e_char)
character(len=*), intent(in) :: str
character, intent(in) :: e_char !< end of list/dict ( '}' or ']')
integer :: N_sq, & !< number of open square brackets
N_cu, & !< number of open curly brackets
i
character(len=:), allocatable:: line
flow_is_closed = .false.
N_sq = 0
N_cu = 0
if(e_char == ']') line = str(index(str(:),'[')+1:)
if(e_char == '}') line = str(index(str(:),'{')+1:)
do i = 1, len_trim(line)
flow_is_closed = (N_sq==0 .and. N_cu==0 .and. scan(line(i:i),e_char) == 1)
N_sq = N_sq + merge(1,0,line(i:i) == '[')
N_cu = N_cu + merge(1,0,line(i:i) == '{')
N_sq = N_sq - merge(1,0,line(i:i) == ']')
N_cu = N_cu - merge(1,0,line(i:i) == '}')
enddo
end function flow_is_closed
!--------------------------------------------------------------------------------------------------
!> @brief return the flow YAML line without line break
!--------------------------------------------------------------------------------------------------
subroutine remove_line_break(blck,s_blck,e_char,flow_line)
character(len=*), intent(in) :: blck !< YAML in mixed style
integer, intent(inout) :: s_blck
character, intent(in) :: e_char !< end of list/dict ( '}' or ']')
character(len=:), allocatable, intent(out) :: flow_line
logical :: line_end
line_end =.false.
flow_line = ''
do while(.not.line_end)
flow_line = flow_line//IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))//' '
line_end = flow_is_closed(flow_line,e_char)
s_blck = s_blck + index(blck(s_blck:),IO_EOL)
enddo
end subroutine remove_line_break
!--------------------------------------------------------------------------------------------------
! @brief reads a line of YAML block which is already in flow style
! @details Dicts should be enlcosed within '{}' for it to be consistent with DAMASK YAML parser
@ -402,7 +466,7 @@ recursive subroutine lst(blck,flow,s_blck,s_flow,offset)
integer, intent(inout) :: s_blck, & !< start position in blck
s_flow, & !< start position in flow
offset !< stores leading '- ' in nested lists
character(len=:), allocatable :: line
character(len=:), allocatable :: line,flow_line
integer :: e_blck,indent
indent = indentDepth(blck(s_blck:),offset)
@ -437,8 +501,12 @@ recursive subroutine lst(blck,flow,s_blck,s_flow,offset)
s_blck = e_blck +2
offset = 0
elseif(isFlow(line)) then
call line_isFlow(flow,s_flow,line)
s_blck = e_blck +2
if(isFlowList(line)) then
call remove_line_break(blck,s_blck,']',flow_line)
else
call remove_line_break(blck,s_blck,'}',flow_line)
endif
call line_isFlow(flow,s_flow,flow_line)
offset = 0
endif
else ! list item in the same line
@ -448,8 +516,13 @@ recursive subroutine lst(blck,flow,s_blck,s_flow,offset)
s_blck = e_blck +2
offset = 0
elseif(isFlow(line)) then
call line_isFlow(flow,s_flow,line)
s_blck = e_blck +2
s_blck = s_blck + index(blck(s_blck:),'-')
if(isFlowList(line)) then
call remove_line_break(blck,s_blck,']',flow_line)
else
call remove_line_break(blck,s_blck,'}',flow_line)
endif
call line_isFlow(flow,s_flow,flow_line)
offset = 0
else ! non scalar list item
offset = offset + indentDepth(blck(s_blck:))+1 ! offset in spaces to be ignored
@ -486,8 +559,8 @@ recursive subroutine dct(blck,flow,s_blck,s_flow,offset)
s_flow, & !< start position in flow
offset
character(len=:), allocatable :: line
integer :: e_blck,indent
character(len=:), allocatable :: line,flow_line
integer :: e_blck,indent,col_pos
logical :: previous_isKey
previous_isKey = .false.
@ -521,12 +594,22 @@ recursive subroutine dct(blck,flow,s_blck,s_flow,offset)
endif
if(isKeyValue(line)) then
col_pos = index(line,':')
if(isFlow(line(col_pos+1:))) then
if(isFlowList(line(col_pos+1:))) then
call remove_line_break(blck,s_blck,']',flow_line)
else
call remove_line_break(blck,s_blck,'}',flow_line)
endif
call keyValue_toFlow(flow,s_flow,flow_line)
else
call keyValue_toFlow(flow,s_flow,line)
s_blck = e_blck + 2
endif
else
call line_toFlow(flow,s_flow,line)
s_blck = e_blck + 2
endif
s_blck = e_blck +2
end if
if(isScalar(line) .or. isKeyValue(line)) then
@ -559,7 +642,7 @@ recursive subroutine decide(blck,flow,s_blck,s_flow,offset)
s_flow, & !< start position in flow
offset
integer :: e_blck
character(len=:), allocatable :: line
character(len=:), allocatable :: line,flow_line
if(s_blck <= len(blck)) then
call skip_empty_lines(blck,s_blck)
@ -583,8 +666,12 @@ recursive subroutine decide(blck,flow,s_blck,s_flow,offset)
flow(s_flow:s_flow) = '}'
s_flow = s_flow + 1
elseif(isFlow(line)) then
if(isFlowList(line)) then
call remove_line_break(blck,s_blck,']',flow_line)
else
call remove_line_break(blck,s_blck,'}',flow_line)
endif
call line_isFlow(flow,s_flow,line)
s_blck = e_blck +2
else
line = line(indentDepth(line)+1:)
call line_toFlow(flow,s_flow,line)
@ -723,6 +810,36 @@ subroutine selfTest
if (.not. to_flow(flow_mixed_braces) == flow) error stop 'to_flow'
end block basic_flow
multi_line_flow1: block
character(len=*), parameter :: flow_multi = &
"%YAML 1.1"//IO_EOL//&
"---"//IO_EOL//&
"a: [b,"//IO_EOL//&
"c: "//IO_EOL//&
"d, e]"//IO_EOL
character(len=*), parameter :: flow = &
"{a: [b, {c: d}, e]}"
if( .not. to_flow(flow_multi) == flow) error stop 'to_flow'
end block multi_line_flow1
multi_line_flow2: block
character(len=*), parameter :: flow_multi = &
"%YAML 1.1"//IO_EOL//&
"---"//IO_EOL//&
"-"//IO_EOL//&
" a: {b:"//IO_EOL//&
"[c,"//IO_EOL//&
"d"//IO_EOL//&
"e, f]}"//IO_EOL
character(len=*), parameter :: flow = &
"[{a: {b: [c, d e, f]}}]"
if( .not. to_flow(flow_multi) == flow) error stop 'to_flow'
end block multi_line_flow2
basic_mixed: block
character(len=*), parameter :: block_flow = &
"%YAML 1.1"//IO_EOL//&

3
src/YAML_types.f90
View File

@ -952,6 +952,9 @@ function tNode_get_byKey_asIndex(self,key) result(keyIndex)
endif
enddo
if(keyIndex == -1) call IO_error(140,ext_msg=key)
end function tNode_get_byKey_asIndex

10
src/constitutive.f90
View File

@ -452,11 +452,11 @@ subroutine constitutive_init
PhaseLoop2:do p = 1,phases%length
!--------------------------------------------------------------------------------------------------
! partition and initialize state
plasticState(p)%partionedState0 = plasticState(p)%state0
plasticState(p)%state = plasticState(p)%partionedState0
plasticState(p)%partitionedState0 = plasticState(p)%state0
plasticState(p)%state = plasticState(p)%partitionedState0
forall(s = 1:phase_Nsources(p))
sourceState(p)%p(s)%partionedState0 = sourceState(p)%p(s)%state0
sourceState(p)%p(s)%state = sourceState(p)%p(s)%partionedState0
sourceState(p)%p(s)%partitionedState0 = sourceState(p)%p(s)%state0
sourceState(p)%p(s)%state = sourceState(p)%p(s)%partitionedState0
end forall
constitutive_source_maxSizeDotState = max(constitutive_source_maxSizeDotState, &
@ -922,7 +922,7 @@ subroutine constitutive_allocateState(state, &
allocate(state%atol (sizeState), source=0.0_pReal)
allocate(state%state0 (sizeState,NipcMyPhase), source=0.0_pReal)
allocate(state%partionedState0(sizeState,NipcMyPhase), source=0.0_pReal)
allocate(state%partitionedState0(sizeState,NipcMyPhase), source=0.0_pReal)
allocate(state%subState0 (sizeState,NipcMyPhase), source=0.0_pReal)
allocate(state%state (sizeState,NipcMyPhase), source=0.0_pReal)

170
src/crystallite.f90
View File

@ -44,30 +44,30 @@ module crystallite
!
crystallite_Fp, & !< current plastic def grad (end of converged time step)
crystallite_Fp0, & !< plastic def grad at start of FE inc
crystallite_partionedFp0,& !< plastic def grad at start of homog inc
crystallite_partitionedFp0,& !< plastic def grad at start of homog inc
crystallite_subFp0,& !< plastic def grad at start of crystallite inc
!
crystallite_Fi, & !< current intermediate def grad (end of converged time step)
crystallite_Fi0, & !< intermediate def grad at start of FE inc
crystallite_partionedFi0,& !< intermediate def grad at start of homog inc
crystallite_partitionedFi0,& !< intermediate def grad at start of homog inc
crystallite_subFi0,& !< intermediate def grad at start of crystallite inc
!
crystallite_Lp0, & !< plastic velocitiy grad at start of FE inc
crystallite_partionedLp0, & !< plastic velocity grad at start of homog inc
crystallite_partitionedLp0, & !< plastic velocity grad at start of homog inc
!
crystallite_Li, & !< current intermediate velocitiy grad (end of converged time step)
crystallite_Li0, & !< intermediate velocitiy grad at start of FE inc
crystallite_partionedLi0, & !< intermediate velocity grad at start of homog inc
crystallite_partitionedLi0, & !< intermediate velocity grad at start of homog inc
!
crystallite_S0, & !< 2nd Piola-Kirchhoff stress vector at start of FE inc
crystallite_partionedS0 !< 2nd Piola-Kirchhoff stress vector at start of homog inc
crystallite_partitionedS0 !< 2nd Piola-Kirchhoff stress vector at start of homog inc
real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: &
crystallite_P, & !< 1st Piola-Kirchhoff stress per grain
crystallite_Lp, & !< current plastic velocitiy grad (end of converged time step)
crystallite_S, & !< current 2nd Piola-Kirchhoff stress vector (end of converged time step)
crystallite_partionedF0 !< def grad at start of homog inc
crystallite_partitionedF0 !< def grad at start of homog inc
real(pReal), dimension(:,:,:,:,:), allocatable, public :: &
crystallite_partionedF !< def grad to be reached at end of homog inc
crystallite_partitionedF !< def grad to be reached at end of homog inc
logical, dimension(:,:,:), allocatable, public :: &
crystallite_requested !< used by upper level (homogenization) to request crystallite calculation
@ -166,20 +166,20 @@ subroutine crystallite_init
iMax = discretization_nIP
eMax = discretization_nElem
allocate(crystallite_partionedF(3,3,cMax,iMax,eMax),source=0.0_pReal)
allocate(crystallite_partitionedF(3,3,cMax,iMax,eMax),source=0.0_pReal)
allocate(crystallite_S0, &
crystallite_F0, crystallite_Fi0,crystallite_Fp0, &
crystallite_Li0,crystallite_Lp0, &
crystallite_partionedS0, &
crystallite_partionedF0,crystallite_partionedFp0,crystallite_partionedFi0, &
crystallite_partionedLp0,crystallite_partionedLi0, &
crystallite_partitionedS0, &
crystallite_partitionedF0,crystallite_partitionedFp0,crystallite_partitionedFi0, &
crystallite_partitionedLp0,crystallite_partitionedLi0, &
crystallite_S,crystallite_P, &
crystallite_Fe,crystallite_Fi,crystallite_Fp, &
crystallite_Li,crystallite_Lp, &
crystallite_subF,crystallite_subF0, &
crystallite_subFp0,crystallite_subFi0, &
source = crystallite_partionedF)
source = crystallite_partitionedF)
allocate(crystallite_dt(cMax,iMax,eMax),source=0.0_pReal)
allocate(crystallite_subdt,crystallite_subFrac,crystallite_subStep, &
@ -269,10 +269,10 @@ subroutine crystallite_init
!$OMP END PARALLEL DO
crystallite_partionedFp0 = crystallite_Fp0
crystallite_partionedFi0 = crystallite_Fi0
crystallite_partionedF0 = crystallite_F0
crystallite_partionedF = crystallite_F0
crystallite_partitionedFp0 = crystallite_Fp0
crystallite_partitionedFi0 = crystallite_Fi0
crystallite_partitionedF0 = crystallite_F0
crystallite_partitionedF = crystallite_F0
call crystallite_orientations()
@ -280,8 +280,8 @@ subroutine crystallite_init
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do c = 1,homogenization_Ngrains(material_homogenizationAt(e))
call constitutive_dependentState(crystallite_partionedF0(1:3,1:3,c,i,e), &
crystallite_partionedFp0(1:3,1:3,c,i,e), &
call constitutive_dependentState(crystallite_partitionedF0(1:3,1:3,c,i,e), &
crystallite_partitionedFp0(1:3,1:3,c,i,e), &
c,i,e) ! update dependent state variables to be consistent with basic states
enddo
enddo
@ -325,8 +325,8 @@ function crystallite_stress()
todo = .false.
subLp0 = crystallite_partionedLp0
subLi0 = crystallite_partionedLi0
subLp0 = crystallite_partitionedLp0
subLi0 = crystallite_partitionedLi0
@ -338,15 +338,15 @@ function crystallite_stress()
do i = FEsolving_execIP(1),FEsolving_execIP(2); do c = 1,homogenization_Ngrains(material_homogenizationAt(e))
homogenizationRequestsCalculation: if (crystallite_requested(c,i,e)) then
plasticState (material_phaseAt(c,e))%subState0( :,material_phaseMemberAt(c,i,e)) = &
plasticState (material_phaseAt(c,e))%partionedState0(:,material_phaseMemberAt(c,i,e))
plasticState (material_phaseAt(c,e))%partitionedState0(:,material_phaseMemberAt(c,i,e))
do s = 1, phase_Nsources(material_phaseAt(c,e))
sourceState(material_phaseAt(c,e))%p(s)%subState0( :,material_phaseMemberAt(c,i,e)) = &
sourceState(material_phaseAt(c,e))%p(s)%partionedState0(:,material_phaseMemberAt(c,i,e))
sourceState(material_phaseAt(c,e))%p(s)%partitionedState0(:,material_phaseMemberAt(c,i,e))
enddo
crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_partionedFp0(1:3,1:3,c,i,e)
crystallite_subFi0(1:3,1:3,c,i,e) = crystallite_partionedFi0(1:3,1:3,c,i,e)
crystallite_subF0(1:3,1:3,c,i,e) = crystallite_partionedF0(1:3,1:3,c,i,e)
crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_partitionedFp0(1:3,1:3,c,i,e)
crystallite_subFi0(1:3,1:3,c,i,e) = crystallite_partitionedFi0(1:3,1:3,c,i,e)
crystallite_subF0(1:3,1:3,c,i,e) = crystallite_partitionedF0(1:3,1:3,c,i,e)
crystallite_subFrac(c,i,e) = 0.0_pReal
crystallite_subStep(c,i,e) = 1.0_pReal/num%subStepSizeCryst
todo(c,i,e) = .true.
@ -426,8 +426,8 @@ function crystallite_stress()
! prepare for integration
if (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_subStep(c,i,e) *( crystallite_partitionedF (1:3,1:3,c,i,e) &
-crystallite_partitionedF0(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)))
@ -441,8 +441,6 @@ function crystallite_stress()
enddo elementLooping3
!$OMP END PARALLEL DO
call nonlocalConvergenceCheck
!--------------------------------------------------------------------------------------------------
! integrate --- requires fully defined state array (basic + dependent state)
where(.not. crystallite_converged .and. crystallite_subStep > num%subStepMinCryst) & ! do not try non-converged but fully cutbacked any further
@ -475,17 +473,17 @@ subroutine crystallite_initializeRestorationPoints(i,e)
s
do c = 1,homogenization_Ngrains(material_homogenizationAt(e))
crystallite_partionedFp0(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e)
crystallite_partionedLp0(1:3,1:3,c,i,e) = crystallite_Lp0(1:3,1:3,c,i,e)
crystallite_partionedFi0(1:3,1:3,c,i,e) = crystallite_Fi0(1:3,1:3,c,i,e)
crystallite_partionedLi0(1:3,1:3,c,i,e) = crystallite_Li0(1:3,1:3,c,i,e)
crystallite_partionedF0(1:3,1:3,c,i,e) = crystallite_F0(1:3,1:3,c,i,e)
crystallite_partionedS0(1:3,1:3,c,i,e) = crystallite_S0(1:3,1:3,c,i,e)
crystallite_partitionedFp0(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e)
crystallite_partitionedLp0(1:3,1:3,c,i,e) = crystallite_Lp0(1:3,1:3,c,i,e)
crystallite_partitionedFi0(1:3,1:3,c,i,e) = crystallite_Fi0(1:3,1:3,c,i,e)
crystallite_partitionedLi0(1:3,1:3,c,i,e) = crystallite_Li0(1:3,1:3,c,i,e)
crystallite_partitionedF0(1:3,1:3,c,i,e) = crystallite_F0(1:3,1:3,c,i,e)
crystallite_partitionedS0(1:3,1:3,c,i,e) = crystallite_S0(1:3,1:3,c,i,e)
plasticState(material_phaseAt(c,e))%partionedState0(:,material_phasememberAt(c,i,e)) = &
plasticState(material_phaseAt(c,e))%partitionedState0(:,material_phasememberAt(c,i,e)) = &
plasticState(material_phaseAt(c,e))%state0( :,material_phasememberAt(c,i,e))
do s = 1, phase_Nsources(material_phaseAt(c,e))
sourceState(material_phaseAt(c,e))%p(s)%partionedState0(:,material_phasememberAt(c,i,e)) = &
sourceState(material_phaseAt(c,e))%p(s)%partitionedState0(:,material_phasememberAt(c,i,e)) = &
sourceState(material_phaseAt(c,e))%p(s)%state0( :,material_phasememberAt(c,i,e))
enddo
enddo
@ -506,17 +504,17 @@ subroutine crystallite_windForward(i,e)
s
do c = 1,homogenization_Ngrains(material_homogenizationAt(e))
crystallite_partionedF0 (1:3,1:3,c,i,e) = crystallite_partionedF(1:3,1:3,c,i,e)
crystallite_partionedFp0(1:3,1:3,c,i,e) = crystallite_Fp (1:3,1:3,c,i,e)
crystallite_partionedLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e)
crystallite_partionedFi0(1:3,1:3,c,i,e) = crystallite_Fi (1:3,1:3,c,i,e)
crystallite_partionedLi0(1:3,1:3,c,i,e) = crystallite_Li (1:3,1:3,c,i,e)
crystallite_partionedS0 (1:3,1:3,c,i,e) = crystallite_S (1:3,1:3,c,i,e)
crystallite_partitionedF0 (1:3,1:3,c,i,e) = crystallite_partitionedF(1:3,1:3,c,i,e)
crystallite_partitionedFp0(1:3,1:3,c,i,e) = crystallite_Fp (1:3,1:3,c,i,e)
crystallite_partitionedLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e)
crystallite_partitionedFi0(1:3,1:3,c,i,e) = crystallite_Fi (1:3,1:3,c,i,e)
crystallite_partitionedLi0(1:3,1:3,c,i,e) = crystallite_Li (1:3,1:3,c,i,e)
crystallite_partitionedS0 (1:3,1:3,c,i,e) = crystallite_S (1:3,1:3,c,i,e)
plasticState (material_phaseAt(c,e))%partionedState0(:,material_phasememberAt(c,i,e)) = &
plasticState (material_phaseAt(c,e))%partitionedState0(:,material_phasememberAt(c,i,e)) = &
plasticState (material_phaseAt(c,e))%state (:,material_phasememberAt(c,i,e))
do s = 1, phase_Nsources(material_phaseAt(c,e))
sourceState(material_phaseAt(c,e))%p(s)%partionedState0(:,material_phasememberAt(c,i,e)) = &
sourceState(material_phaseAt(c,e))%p(s)%partitionedState0(:,material_phasememberAt(c,i,e)) = &
sourceState(material_phaseAt(c,e))%p(s)%state (:,material_phasememberAt(c,i,e))
enddo
enddo
@ -540,18 +538,18 @@ subroutine crystallite_restore(i,e,includeL)
do c = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (includeL) then
crystallite_Lp(1:3,1:3,c,i,e) = crystallite_partionedLp0(1:3,1:3,c,i,e)
crystallite_Li(1:3,1:3,c,i,e) = crystallite_partionedLi0(1:3,1:3,c,i,e)
crystallite_Lp(1:3,1:3,c,i,e) = crystallite_partitionedLp0(1:3,1:3,c,i,e)
crystallite_Li(1:3,1:3,c,i,e) = crystallite_partitionedLi0(1:3,1:3,c,i,e)
endif ! maybe protecting everything from overwriting makes more sense
crystallite_Fp(1:3,1:3,c,i,e) = crystallite_partionedFp0(1:3,1:3,c,i,e)
crystallite_Fi(1:3,1:3,c,i,e) = crystallite_partionedFi0(1:3,1:3,c,i,e)
crystallite_S (1:3,1:3,c,i,e) = crystallite_partionedS0 (1:3,1:3,c,i,e)
crystallite_Fp(1:3,1:3,c,i,e) = crystallite_partitionedFp0(1:3,1:3,c,i,e)
crystallite_Fi(1:3,1:3,c,i,e) = crystallite_partitionedFi0(1:3,1:3,c,i,e)
crystallite_S (1:3,1:3,c,i,e) = crystallite_partitionedS0 (1:3,1:3,c,i,e)
plasticState (material_phaseAt(c,e))%state( :,material_phasememberAt(c,i,e)) = &
plasticState (material_phaseAt(c,e))%partionedState0(:,material_phasememberAt(c,i,e))
plasticState (material_phaseAt(c,e))%partitionedState0(:,material_phasememberAt(c,i,e))
do s = 1, phase_Nsources(material_phaseAt(c,e))
sourceState(material_phaseAt(c,e))%p(s)%state( :,material_phasememberAt(c,i,e)) = &
sourceState(material_phaseAt(c,e))%p(s)%partionedState0(:,material_phasememberAt(c,i,e))
sourceState(material_phaseAt(c,e))%p(s)%partitionedState0(:,material_phasememberAt(c,i,e))
enddo
enddo
@ -758,7 +756,7 @@ subroutine crystallite_results
do o = 1, size(output_constituent(p)%label)
select case (output_constituent(p)%label(o))
case('F')
selected_tensors = select_tensors(crystallite_partionedF,p)
selected_tensors = select_tensors(crystallite_partitionedF,p)
call results_writeDataset(group,selected_tensors,output_constituent(p)%label(o),&
'deformation gradient','1')
case('Fe')
@ -943,7 +941,7 @@ function integrateStress(ipc,ip,el,timeFraction) result(broken)
F = crystallite_subF(1:3,1:3,ipc,ip,el)
endif
call constitutive_dependentState(crystallite_partionedF(1:3,1:3,ipc,ip,el), &
call constitutive_dependentState(crystallite_partitionedF(1:3,1:3,ipc,ip,el), &
crystallite_Fp(1:3,1:3,ipc,ip,el),ipc,ip,el)
Lpguess = crystallite_Lp(1:3,1:3,ipc,ip,el) ! take as first guess
@ -1120,9 +1118,9 @@ subroutine integrateStateFPI(g,i,e)
c = material_phaseMemberAt(g,i,e)
broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, &
crystallite_partitionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, &
crystallite_partitionedFp0, &
crystallite_subdt(g,i,e), g,i,e,p,c)
if(broken) return
@ -1152,9 +1150,9 @@ subroutine integrateStateFPI(g,i,e)
if(broken) exit iteration
broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, &
crystallite_partitionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, &
crystallite_partitionedFp0, &
crystallite_subdt(g,i,e), g,i,e,p,c)
if(broken) exit iteration
@ -1243,9 +1241,9 @@ subroutine integrateStateEuler(g,i,e)
c = material_phaseMemberAt(g,i,e)
broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, &
crystallite_partitionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, &
crystallite_partitionedFp0, &
crystallite_subdt(g,i,e), g,i,e,p,c)
if(broken) return
@ -1296,9 +1294,9 @@ subroutine integrateStateAdaptiveEuler(g,i,e)
c = material_phaseMemberAt(g,i,e)
broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, &
crystallite_partitionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, &
crystallite_partitionedFp0, &
crystallite_subdt(g,i,e), g,i,e,p,c)
if(broken) return
@ -1325,9 +1323,9 @@ subroutine integrateStateAdaptiveEuler(g,i,e)
if(broken) return
broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, &
crystallite_partitionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, &
crystallite_partitionedFp0, &
crystallite_subdt(g,i,e), g,i,e,p,c)
if(broken) return
@ -1434,9 +1432,9 @@ subroutine integrateStateRK(g,i,e,A,B,CC,DB)
c = material_phaseMemberAt(g,i,e)
broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, &
crystallite_partitionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, &
crystallite_partitionedFp0, &
crystallite_subdt(g,i,e), g,i,e,p,c)
if(broken) return
@ -1476,9 +1474,9 @@ subroutine integrateStateRK(g,i,e,A,B,CC,DB)
if(broken) exit
broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, &
crystallite_partitionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, &
crystallite_partitionedFp0, &
crystallite_subdt(g,i,e)*CC(stage), g,i,e,p,c)
if(broken) exit
@ -1526,38 +1524,6 @@ subroutine integrateStateRK(g,i,e,A,B,CC,DB)
end subroutine integrateStateRK
!--------------------------------------------------------------------------------------------------
!> @brief sets convergence flag for nonlocal calculations
!> @details one non-converged nonlocal sets all other nonlocals to non-converged to trigger cut back
!--------------------------------------------------------------------------------------------------
subroutine nonlocalConvergenceCheck
integer :: e,i,p
logical :: nonlocal_broken
nonlocal_broken = .false.
!$OMP PARALLEL DO PRIVATE(p)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
p = material_phaseAt(1,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
if(plasticState(p)%nonlocal .and. .not. crystallite_converged(1,i,e)) nonlocal_broken = .true.
enddo
enddo
!$OMP END PARALLEL DO
if(.not. nonlocal_broken) return
!$OMP PARALLEL DO PRIVATE(p)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
p = material_phaseAt(1,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
if(plasticState(p)%nonlocal) crystallite_converged(1,i,e) = .false.
enddo
enddo
!$OMP END PARALLEL DO
end subroutine nonlocalConvergenceCheck
!--------------------------------------------------------------------------------------------------
!> @brief determines whether a point is converged
!--------------------------------------------------------------------------------------------------
@ -1590,7 +1556,7 @@ subroutine crystallite_restartWrite
write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5'
fileHandle = HDF5_openFile(fileName,'a')
call HDF5_write(fileHandle,crystallite_partionedF,'F')
call HDF5_write(fileHandle,crystallite_partitionedF,'F')
call HDF5_write(fileHandle,crystallite_Fp, 'Fp')
call HDF5_write(fileHandle,crystallite_Fi, 'Fi')
call HDF5_write(fileHandle,crystallite_Lp, 'Lp')
@ -1665,7 +1631,7 @@ subroutine crystallite_forward
integer :: i, j
crystallite_F0 = crystallite_partionedF
crystallite_F0 = crystallite_partitionedF
crystallite_Fp0 = crystallite_Fp
crystallite_Lp0 = crystallite_Lp
crystallite_Fi0 = crystallite_Fi

10
src/homogenization.f90
View File

@ -404,16 +404,16 @@ subroutine partitionDeformation(subF,ip,el)
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization
crystallite_partionedF(1:3,1:3,1,ip,el) = subF
crystallite_partitionedF(1:3,1:3,1,ip,el) = subF
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
call mech_isostrain_partitionDeformation(&
crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
crystallite_partitionedF(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
subF)
case (HOMOGENIZATION_RGC_ID) chosenHomogenization
call mech_RGC_partitionDeformation(&
crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
crystallite_partitionedF(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
subF,&
ip, &
el)
@ -448,8 +448,8 @@ function updateState(subdt,subF,ip,el)
updateState = &
updateState .and. &
mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
crystallite_partionedF0(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el),&
crystallite_partitionedF(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
crystallite_partitionedF0(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el),&
subF,&
subdt, &
dPdFs, &

4
src/homogenization_mech_RGC.f90
View File

@ -212,7 +212,7 @@ end subroutine mech_RGC_init
!--------------------------------------------------------------------------------------------------
module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of)
real(pReal), dimension (:,:,:), intent(out) :: F !< partioned F per grain
real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned F per grain
real(pReal), dimension (3,3), intent(in) :: avgF !< averaged F
integer, intent(in) :: &
@ -867,7 +867,7 @@ module procedure mech_RGC_updateState
!--------------------------------------------------------------------------------------------------
subroutine grainDeformation(F, avgF, instance, of)
real(pReal), dimension(:,:,:), intent(out) :: F !< partioned F per grain
real(pReal), dimension(:,:,:), intent(out) :: F !< partitioned F per grain
real(pReal), dimension(:,:), intent(in) :: avgF !< averaged F
integer, intent(in) :: &

2
src/prec.f90
View File

@ -47,7 +47,7 @@ module prec
dotState, & !< rate of state change
deltaState !< increment of state change
real(pReal), allocatable, dimension(:,:) :: &
partionedState0, &
partitionedState0, &
subState0
end type