Merge branch 'development' into allclose-rotation-orientation

This commit is contained in:
Martin Diehl 2021-04-10 08:05:25 +02:00
commit 63b343ad87
80 changed files with 1993 additions and 1382 deletions

3
.gitattributes vendored
View File

@ -9,10 +9,11 @@
*.hdf5 binary *.hdf5 binary
*.pdf binary *.pdf binary
*.dream3d binary *.dream3d binary
*.pbz2 binary
# ignore files from MSC.Marc in language statistics # ignore files from MSC.Marc in language statistics
installation/mods_MarcMentat/20*/* linguist-vendored installation/mods_MarcMentat/20*/* linguist-vendored
src/marc/include/* linguist-vendored src/Marc/include/* linguist-vendored
# ignore reference files for tests in language statistics # ignore reference files for tests in language statistics
python/tests/reference/* linguist-vendored python/tests/reference/* linguist-vendored

@ -1 +1 @@
Subproject commit f1ac733d5eb90de1cdcaf79a261157c9034f8136 Subproject commit 1298124143e7e2901d0b9c2e79ab6388cb78a1e3

View File

@ -1 +1 @@
v3.0.0-alpha2-753-g565dab120 v3.0.0-alpha2-829-g73b07eda4

View File

@ -1,5 +1,3 @@
import os
import numpy as np import numpy as np
import h5py import h5py
@ -159,18 +157,18 @@ class ConfigMaterial(Config):
f = h5py.File(fname,'r') f = h5py.File(fname,'r')
if grain_data is None: if grain_data is None:
phase = f[os.path.join(b,c,phases)][()].flatten() phase = f['/'.join([b,c,phases])][()].flatten()
O = Rotation.from_Euler_angles(f[os.path.join(b,c,Euler_angles)]).as_quaternion().reshape(-1,4) # noqa O = Rotation.from_Euler_angles(f['/'.join([b,c,Euler_angles])]).as_quaternion().reshape(-1,4) # noqa
_,idx = np.unique(np.hstack([O,phase.reshape(-1,1)]),return_index=True,axis=0) _,idx = np.unique(np.hstack([O,phase.reshape(-1,1)]),return_index=True,axis=0)
idx = np.sort(idx) idx = np.sort(idx)
else: else:
phase = f[os.path.join(b,grain_data,phases)][()] phase = f['/'.join([b,grain_data,phases])][()]
O = Rotation.from_Euler_angles(f[os.path.join(b,grain_data,Euler_angles)]).as_quaternion() # noqa O = Rotation.from_Euler_angles(f['/'.join([b,grain_data,Euler_angles])]).as_quaternion() # noqa
idx = np.arange(phase.size) idx = np.arange(phase.size)
if cell_ensemble_data is not None and phase_names is not None: if cell_ensemble_data is not None and phase_names is not None:
try: try:
names = np.array([s.decode() for s in f[os.path.join(b,cell_ensemble_data,phase_names)]]) names = np.array([s.decode() for s in f['/'.join([b,cell_ensemble_data,phase_names])]])
phase = names[phase] phase = names[phase]
except KeyError: except KeyError:
pass pass

View File

@ -69,9 +69,9 @@ class Grid:
copy = __copy__ copy = __copy__
def diff(self,other): def __eq__(self,other):
""" """
Report property differences of self relative to other. Test equality of other.
Parameters Parameters
---------- ----------
@ -79,28 +79,10 @@ class Grid:
Grid to compare self against. Grid to compare self against.
""" """
message = [] return (np.allclose(other.size,self.size)
if np.any(other.cells != self.cells): and np.allclose(other.origin,self.origin)
message.append(util.deemph(f'cells a b c: {util.srepr(other.cells," x ")}')) and np.all(other.cells == self.cells)
message.append(util.emph( f'cells a b c: {util.srepr( self.cells," x ")}')) and np.all(other.material == self.material))
if not np.allclose(other.size,self.size):
message.append(util.deemph(f'size x y z: {util.srepr(other.size," x ")}'))
message.append(util.emph( f'size x y z: {util.srepr( self.size," x ")}'))
if not np.allclose(other.origin,self.origin):
message.append(util.deemph(f'origin x y z: {util.srepr(other.origin," ")}'))
message.append(util.emph( f'origin x y z: {util.srepr( self.origin," ")}'))
if other.N_materials != self.N_materials:
message.append(util.deemph(f'# materials: {other.N_materials}'))
message.append(util.emph( f'# materials: { self.N_materials}'))
if np.nanmax(other.material) != np.nanmax(self.material):
message.append(util.deemph(f'max material: {np.nanmax(other.material)}'))
message.append(util.emph( f'max material: {np.nanmax( self.material)}'))
return util.return_message(message)
@property @property
@ -270,7 +252,7 @@ class Grid:
""" """
Load DREAM.3D (HDF5) file. Load DREAM.3D (HDF5) file.
Data in DREAM.3D files can be stored per cell ('CellData') and/or Data in DREAM.3D files can be stored per cell ('CellData') and/or
per grain ('Grain Data'). Per default, cell-wise data is assumed. per grain ('Grain Data'). Per default, cell-wise data is assumed.
damask.ConfigMaterial.load_DREAM3D gives the corresponding material definition. damask.ConfigMaterial.load_DREAM3D gives the corresponding material definition.
@ -305,18 +287,18 @@ class Grid:
c = util.DREAM3D_cell_data_group(fname) if cell_data is None else cell_data c = util.DREAM3D_cell_data_group(fname) if cell_data is None else cell_data
f = h5py.File(fname, 'r') f = h5py.File(fname, 'r')
cells = f[os.path.join(b,'_SIMPL_GEOMETRY','DIMENSIONS')][()] cells = f['/'.join([b,'_SIMPL_GEOMETRY','DIMENSIONS'])][()]
size = f[os.path.join(b,'_SIMPL_GEOMETRY','SPACING')] * cells size = f['/'.join([b,'_SIMPL_GEOMETRY','SPACING'])] * cells
origin = f[os.path.join(b,'_SIMPL_GEOMETRY','ORIGIN')][()] origin = f['/'.join([b,'_SIMPL_GEOMETRY','ORIGIN'])][()]
if feature_IDs is None: if feature_IDs is None:
phase = f[os.path.join(b,c,phases)][()].reshape(-1,1) phase = f['/'.join([b,c,phases])][()].reshape(-1,1)
O = Rotation.from_Euler_angles(f[os.path.join(b,c,Euler_angles)]).as_quaternion().reshape(-1,4) # noqa O = Rotation.from_Euler_angles(f['/'.join([b,c,Euler_angles])]).as_quaternion().reshape(-1,4) # noqa
unique,unique_inverse = np.unique(np.hstack([O,phase]),return_inverse=True,axis=0) unique,unique_inverse = np.unique(np.hstack([O,phase]),return_inverse=True,axis=0)
ma = np.arange(cells.prod()) if len(unique) == cells.prod() else \ ma = np.arange(cells.prod()) if len(unique) == cells.prod() else \
np.arange(unique.size)[np.argsort(pd.unique(unique_inverse))][unique_inverse] np.arange(unique.size)[np.argsort(pd.unique(unique_inverse))][unique_inverse]
else: else:
ma = f[os.path.join(b,c,feature_IDs)][()].flatten() ma = f['/'.join([b,c,feature_IDs])][()].flatten()
return Grid(ma.reshape(cells,order='F'),size,origin,util.execution_stamp('Grid','load_DREAM3D')) return Grid(ma.reshape(cells,order='F'),size,origin,util.execution_stamp('Grid','load_DREAM3D'))

File diff suppressed because it is too large Load Diff

View File

@ -230,7 +230,6 @@ class Table:
f = fname f = fname
f.seek(0) f.seek(0)
f.seek(0)
comments = [] comments = []
line = f.readline().strip() line = f.readline().strip()
while line.startswith('#'): while line.startswith('#'):
@ -515,7 +514,7 @@ class Table:
""" """
if set(self.shapes) & set(other.shapes) or self.data.shape[0] != other.data.shape[0]: if set(self.shapes) & set(other.shapes) or self.data.shape[0] != other.data.shape[0]:
raise KeyError('Dublicated keys or row count mismatch') raise KeyError('Duplicated keys or row count mismatch')
else: else:
dup = self.copy() dup = self.copy()
dup.data = dup.data.join(other.data) dup.data = dup.data.join(other.data)

View File

@ -2,8 +2,8 @@ import os
import multiprocessing as mp import multiprocessing as mp
from pathlib import Path from pathlib import Path
import pandas as pd
import numpy as np import numpy as np
import numpy.ma as ma
import vtk import vtk
from vtk.util.numpy_support import numpy_to_vtk as np_to_vtk from vtk.util.numpy_support import numpy_to_vtk as np_to_vtk
from vtk.util.numpy_support import numpy_to_vtkIdTypeArray as np_to_vtkIdTypeArray from vtk.util.numpy_support import numpy_to_vtkIdTypeArray as np_to_vtkIdTypeArray
@ -224,14 +224,14 @@ class VTK:
# Check https://blog.kitware.com/ghost-and-blanking-visibility-changes/ for missing data # Check https://blog.kitware.com/ghost-and-blanking-visibility-changes/ for missing data
# Needs support for pd.DataFrame and/or table # Needs support for damask.Table
def add(self,data,label=None): def add(self,data,label=None):
""" """
Add data to either cells or points. Add data to either cells or points.
Parameters Parameters
---------- ----------
data : numpy.ndarray data : numpy.ndarray or numpy.ma.MaskedArray
Data to add. First dimension needs to match either Data to add. First dimension needs to match either
number of cells or number of points. number of cells or number of points.
label : str label : str
@ -246,8 +246,10 @@ class VTK:
raise ValueError('No label defined for numpy.ndarray') raise ValueError('No label defined for numpy.ndarray')
N_data = data.shape[0] N_data = data.shape[0]
d = np_to_vtk((data.astype(np.single) if data.dtype in [np.double, np.longdouble] else data_ = np.where(data.mask,data.fill_value,data) if isinstance(data,ma.MaskedArray) else\
data).reshape(N_data,-1),deep=True) # avoid large files data
d = np_to_vtk((data_.astype(np.single) if data_.dtype in [np.double, np.longdouble] else
data_).reshape(N_data,-1),deep=True) # avoid large files
d.SetName(label) d.SetName(label)
if N_data == N_points: if N_data == N_points:
@ -256,8 +258,6 @@ class VTK:
self.vtk_data.GetCellData().AddArray(d) self.vtk_data.GetCellData().AddArray(d)
else: else:
raise ValueError(f'Cell / point count ({N_cells} / {N_points}) differs from data ({N_data}).') raise ValueError(f'Cell / point count ({N_cells} / {N_points}) differs from data ({N_data}).')
elif isinstance(data,pd.DataFrame):
raise NotImplementedError('pd.DataFrame')
elif isinstance(data,Table): elif isinstance(data,Table):
raise NotImplementedError('damask.Table') raise NotImplementedError('damask.Table')
else: else:

View File

@ -17,15 +17,16 @@ __all__=[
'srepr', 'srepr',
'emph','deemph','warn','strikeout', 'emph','deemph','warn','strikeout',
'execute', 'execute',
'natural_sort',
'show_progress', 'show_progress',
'scale_to_coprime', 'scale_to_coprime',
'project_stereographic', 'project_stereographic',
'hybrid_IA', 'hybrid_IA',
'return_message',
'execution_stamp', 'execution_stamp',
'shapeshifter', 'shapeblender', 'shapeshifter', 'shapeblender',
'extend_docstring', 'extended_docstring', 'extend_docstring', 'extended_docstring',
'DREAM3D_base_group', 'DREAM3D_cell_data_group' 'DREAM3D_base_group', 'DREAM3D_cell_data_group',
'dict_prune', 'dict_flatten'
] ]
# https://svn.blender.org/svnroot/bf-blender/trunk/blender/build_files/scons/tools/bcolors.py # https://svn.blender.org/svnroot/bf-blender/trunk/blender/build_files/scons/tools/bcolors.py
@ -58,11 +59,12 @@ def srepr(arg,glue = '\n'):
Glue used for joining operation. Defaults to \n. Glue used for joining operation. Defaults to \n.
""" """
if (not hasattr(arg, "strip") and if (not hasattr(arg, 'strip') and
(hasattr(arg, "__getitem__") or (hasattr(arg, '__getitem__') or
hasattr(arg, "__iter__"))): hasattr(arg, '__iter__'))):
return glue.join(str(x) for x in arg) return glue.join(str(x) for x in arg)
return arg if isinstance(arg,str) else repr(arg) else:
return arg if isinstance(arg,str) else repr(arg)
def emph(what): def emph(what):
@ -112,6 +114,11 @@ def execute(cmd,wd='./',env=None):
return process.stdout, process.stderr return process.stdout, process.stderr
def natural_sort(key):
convert = lambda text: int(text) if text.isdigit() else text
return [ convert(c) for c in re.split('([0-9]+)', key) ]
def show_progress(iterable,N_iter=None,prefix='',bar_length=50): def show_progress(iterable,N_iter=None,prefix='',bar_length=50):
""" """
Decorate a loop with a status bar. Decorate a loop with a status bar.
@ -130,11 +137,15 @@ def show_progress(iterable,N_iter=None,prefix='',bar_length=50):
Character length of bar. Defaults to 50. Character length of bar. Defaults to 50.
""" """
status = _ProgressBar(N_iter if N_iter else len(iterable),prefix,bar_length) if N_iter == 1 or (hasattr(iterable,'__len__') and len(iterable) == 1):
for item in iterable:
yield item
else:
status = _ProgressBar(N_iter if N_iter is not None else len(iterable),prefix,bar_length)
for i,item in enumerate(iterable): for i,item in enumerate(iterable):
yield item yield item
status.update(i) status.update(i)
def scale_to_coprime(v): def scale_to_coprime(v):
@ -388,7 +399,7 @@ def DREAM3D_cell_data_group(fname):
""" """
base_group = DREAM3D_base_group(fname) base_group = DREAM3D_base_group(fname)
with h5py.File(fname,'r') as f: with h5py.File(fname,'r') as f:
cells = tuple(f[os.path.join(base_group,'_SIMPL_GEOMETRY','DIMENSIONS')][()][::-1]) cells = tuple(f['/'.join([base_group,'_SIMPL_GEOMETRY','DIMENSIONS'])][()][::-1])
cell_data_group = f[base_group].visititems(lambda path,obj: path.split('/')[0] \ cell_data_group = f[base_group].visititems(lambda path,obj: path.split('/')[0] \
if isinstance(obj,h5py._hl.dataset.Dataset) and np.shape(obj)[:-1] == cells \ if isinstance(obj,h5py._hl.dataset.Dataset) and np.shape(obj)[:-1] == cells \
else None) else None)
@ -399,29 +410,59 @@ def DREAM3D_cell_data_group(fname):
return cell_data_group return cell_data_group
def dict_prune(d):
"""
Recursively remove empty dictionaries.
Parameters
----------
d : dict
Dictionary to prune.
Returns
-------
pruned : dict
Pruned dictionary.
"""
# https://stackoverflow.com/questions/48151953
new = {}
for k,v in d.items():
if isinstance(v, dict):
v = dict_prune(v)
if not isinstance(v,dict) or v != {}:
new[k] = v
return new
def dict_flatten(d):
"""
Recursively remove keys of single-entry dictionaries.
Parameters
----------
d : dict
Dictionary to flatten.
Returns
-------
flattened : dict
Flattened dictionary.
"""
if isinstance(d,dict) and len(d) == 1:
entry = d[list(d.keys())[0]]
new = dict_flatten(entry.copy()) if isinstance(entry,dict) else entry
else:
new = {k: (dict_flatten(v) if isinstance(v, dict) else v) for k,v in d.items()}
return new
#################################################################################################### ####################################################################################################
# Classes # Classes
#################################################################################################### ####################################################################################################
class return_message:
"""Object with formatted return message."""
def __init__(self,message):
"""
Set return message.
Parameters
----------
message : str or list of str
message for output to screen
"""
self.message = message
def __repr__(self):
"""Return message suitable for interactive shells."""
return srepr(self.message)
class _ProgressBar: class _ProgressBar:
""" """
Report progress of an interation as a status bar. Report progress of an interation as a status bar.

View File

@ -0,0 +1,678 @@
material:
- constituents:
- phase: A
O: [0.44698096036898677, -0.035696623617731814, 0.050331402131209624, -0.8924127532086364]
v: 0.125
- phase: A
O: [0.6226205515354057, 0.4257792272120092, -0.6365250917893261, 0.16090837766592916]
v: 0.125
- phase: A
O: [0.4319357791504236, 0.5804552577310466, 0.4870786336250474, 0.48913963356904305]
v: 0.125
- phase: A
O: [0.6712122425649047, -0.5315470950319967, 0.11982471915778993, -0.5025672570639611]
v: 0.125
- phase: A
O: [0.03979536866120591, -0.01919843067074171, -0.5618628311727827, 0.8260495795286167]
v: 0.125
- phase: A
O: [0.6951890945163803, -0.46325436790720337, 0.03393768848173515, -0.5485943371753935]
v: 0.125
- phase: A
O: [0.3489469292795834, -0.5959323325634578, -0.6067706771532161, 0.3936115355256417]
v: 0.125
- phase: A
O: [0.39500485674202723, -0.07573446369604235, -0.015761384830863343, -0.9154163167144754]
v: 0.125
homogenization: RGC
- constituents:
- phase: B
O: [0.9221278604171202, 0.32701451341562027, 0.15952215192749652, -0.1315081750406033]
v: 0.125
- phase: B
O: [0.8117843306311869, 0.4102779561929707, -0.18584663447455768, 0.37167085930736465]
v: 0.125
- phase: B
O: [0.6548302235705483, 0.6510678661349768, -0.3832730015357094, 0.02024396894884073]
v: 0.125
- phase: B
O: [0.04255709858343937, -0.3005730931619659, -0.3060125731417352, -0.9023308783957146]
v: 0.125
- phase: B
O: [0.6960740097259569, 0.16329320418261428, 0.6991228405747408, 0.006599715032396515]
v: 0.125
- phase: B
O: [0.48818278044780217, -0.664294296073913, 0.2679908412891581, 0.4985695238008905]
v: 0.125
- phase: B
O: [0.545426202615872, 0.179298582176746, -0.7847831868864313, -0.23340442478628137]
v: 0.125
- phase: B
O: [0.8566603715426528, 0.26967937521434965, -0.03159131558814077, 0.43864339866435076]
v: 0.125
homogenization: RGC
- constituents:
- phase: C
O: [0.140304599816322, -0.017971728003341014, -0.9895572782263894, 0.027713342853867628]
v: 0.125
- phase: C
O: [0.8096406385129331, 0.2395956967578664, -0.25531574170054405, -0.47105181307727034]
v: 0.125
- phase: C
O: [0.15570676199790476, -0.7501738841096782, 0.2328253832435299, -0.5989882209070811]
v: 0.125
- phase: C
O: [0.4595649941522109, 0.3331746723324176, 0.23957383093695328, 0.7876541331042811]
v: 0.125
- phase: C
O: [0.2969510261483044, 0.0027640644222379847, -0.8133130954773664, 0.5003341450894221]
v: 0.125
- phase: C
O: [0.114170281047612, -0.9068850613039258, 0.3701597984948465, -0.16585040273553361]
v: 0.125
- phase: C
O: [0.29658840046433166, -0.5331902828876833, 0.7365409179409275, 0.2919776004129383]
v: 0.125
- phase: C
O: [0.931315147629592, -0.2017743466108175, 0.30303719188324046, 0.010376376100021569]
v: 0.125
homogenization: RGC
- constituents:
- phase: A
O: [0.20901169650798976, 0.5778057521557359, -0.3628020593274987, 0.7005920990464581]
v: 0.125
- phase: B
O: [0.6025573623737024, 0.7345291803280383, 0.3086529661391562, 0.04609614722911203]
v: 0.125
- phase: C
O: [0.6402291970447783, -0.5402560454335561, -0.5285008279381848, -0.13753856002062462]
v: 0.125
- phase: A
O: [0.12506307996119548, 0.5336909245447153, 0.8344001590799459, 0.05752910234470594]
v: 0.125
- phase: B
O: [0.5504712846116686, -0.6241210958640075, 0.42061504390637533, 0.3612993320712448]
v: 0.125
- phase: C
O: [0.3135636745691277, 0.6574027353927014, 0.6044800436379738, -0.32265049563317466]
v: 0.125
- phase: A
O: [0.1619452502554878, -0.4269267792831907, -0.8080415555188633, 0.3722581169097929]
v: 0.125
- phase: A
O: [0.17044595054343245, 0.8646695166660998, -0.30004868613437863, -0.365055599656809]
v: 0.125
homogenization: RGC
- constituents:
- phase: B
O: [0.43255174377606787, 0.2026672781261794, 0.35000490906160386, 0.8058048938583007]
v: 0.125
- phase: B
O: [0.19591912697365416, -0.8181375036268687, 0.3260880890153691, -0.4311998133665892]
v: 0.125
- phase: A
O: [0.1373248436058589, -0.43524771427953446, -0.8885320171891298, -0.047033700395389226]
v: 0.125
- phase: A
O: [0.250519065956112, -0.7636612069861932, 0.5672425862778043, -0.17971534951065127]
v: 0.125
- phase: B
O: [0.600521695167639, -0.2273230963128902, 0.7563911883237204, -0.12478090295368073]
v: 0.125
- phase: A
O: [0.4281846857945264, 0.5284504248063085, 0.21428541450090705, 0.7010561921167584]
v: 0.125
- phase: A
O: [0.34007599923699156, 0.5500172687924186, 0.6138467460818403, -0.45279298923219496]
v: 0.125
- phase: A
O: [0.8581382626679844, -0.20034169756111123, 0.3134907520728513, -0.35381559424127107]
v: 0.125
homogenization: RGC
- constituents:
- phase: C
O: [0.18637512817696594, 0.9354389080016575, 0.2967480269088709, -0.04646471262558266]
v: 0.125
- phase: C
O: [0.33425800495572805, 0.24948308477875414, 0.7297193930640264, 0.5417927499686225]
v: 0.125
- phase: A
O: [0.11918214701666907, -0.4185869380028542, 0.618905484805619, -0.6538628235673086]
v: 0.125
- phase: A
O: [0.40732205687168244, -0.6166051531766635, -0.6737097202702335, -0.0014282419993970551]
v: 0.125
- phase: A
O: [0.68397144415648, 0.6848389352552103, 0.24111983008651763, 0.07099242125789083]
v: 0.125
- phase: C
O: [0.05719438496521803, 0.2575447319897859, 0.5020591650211859, -0.8236116245968056]
v: 0.125
- phase: A
O: [0.10565206937493385, -0.1109750475948143, 0.8781814355546643, 0.45312199824690885]
v: 0.125
- phase: A
O: [0.22821192086563635, -0.9595048060506683, -0.07745760541549483, -0.1458429487626446]
v: 0.125
homogenization: RGC
- constituents:
- phase: A
O: [0.5745354940858776, 0.6366614824775295, 0.3307586622977252, -0.39391601907009377]
v: 0.125
- phase: B
O: [0.3505113084912295, -0.6046500894621007, 0.6274019946374426, -0.34337563841687757]
v: 0.125
- phase: A
O: [0.030949493058962222, 0.4416700940450431, 0.8819578143863066, 0.16161704906526728]
v: 0.125
- phase: A
O: [0.13563849174261922, 0.2687388407244776, 0.809782239640578, -0.5036212459840637]
v: 0.125
- phase: A
O: [0.2305505587522992, -0.6837828098041563, 0.3813348695443346, 0.577815910255975]
v: 0.125
- phase: A
O: [0.10699977862890839, 0.8478562366641894, 0.26664808851893185, -0.4456339823355066]
v: 0.125
- phase: A
O: [0.9203073703838264, 0.35592932457392046, 0.15147737103605818, -0.058337517855697775]
v: 0.125
- phase: B
O: [0.8544895693446002, 0.43485422381677186, -0.2696115267291114, 0.08977195867747481]
v: 0.125
homogenization: RGC
- constituents:
- phase: A
O: [0.22196756947117874, 0.940293527559, -0.2573416285961528, -0.018808676859083388]
v: 0.125
- phase: A
O: [0.28669138183611736, -0.5847463393505065, 0.6896415010701208, 0.3166612862331461]
v: 0.125
- phase: A
O: [0.18214681844078473, 0.07037956730772883, 0.6195004937983107, -0.7603212421214639]
v: 0.125
- phase: A
O: [0.5744613003376511, 0.17468957401795682, -0.7944671360487727, 0.09110289173380101]
v: 0.125
- phase: C
O: [0.4329446256398051, -0.7351014072423584, 0.5116672951703971, -0.10188940697110307]
v: 0.125
- phase: C
O: [0.2946245183318559, -0.8241105780568372, 0.40376354741611825, -0.2664829189845001]
v: 0.125
- phase: A
O: [0.22212591383189398, -0.06959072743195248, -0.8166434230649855, -0.5281199945320578]
v: 0.125
- phase: A
O: [0.038576201076459135, 0.9876698316278196, 0.06851090552108202, -0.13537516843004982]
v: 0.125
homogenization: RGC
- constituents:
- phase: C
O: [0.004892880719601024, -0.02167422880443418, 0.3339942463386166, 0.9423131809205983]
v: 0.125
- phase: A
O: [0.3719521264362446, -0.4055315244649798, -0.39177426332273346, 0.7373660725193388]
v: 0.125
- phase: A
O: [0.010947854228877327, 0.16142752376310676, -0.19859031041582603, -0.9666349816080735]
v: 0.125
- phase: B
O: [0.9101129345691923, -0.17592965401443023, -0.30944622338368544, 0.21209959453471436]
v: 0.125
- phase: A
O: [0.5280075626141939, -0.4828698476824812, -0.1580687189774305, 0.6804843893155447]
v: 0.125
- phase: A
O: [0.8213575048478461, 0.18229514326534527, -0.35486977770450967, -0.4076858727549186]
v: 0.125
- phase: A
O: [0.16635046542653203, -0.8410358537737095, -0.47641083611371204, 0.19498443669415633]
v: 0.125
- phase: B
O: [0.5406080589325533, 0.5405688205487984, -0.4980750162732086, -0.4092060056158767]
v: 0.125
homogenization: RGC
- constituents:
- phase: A
O: [0.41778322805525603, 0.23584193505422144, 0.8768801833720732, -0.03028035724639885]
v: 0.125
- phase: A
O: [0.8727428843707065, 0.0764211249039828, 0.3927224252414984, -0.2797298092108618]
v: 0.125
- phase: A
O: [0.8766166377340643, -0.2644877070285106, 0.022928146762306964, -0.40132757613285325]
v: 0.125
- phase: A
O: [0.6643649540361558, -0.04353445577200276, -0.2630693201131682, -0.6982252443333509]
v: 0.125
- phase: A
O: [0.43011162904952843, -0.7133156350005553, 0.22368449997894274, -0.5061126711408104]
v: 0.125
- phase: A
O: [0.42119993016818996, -0.2666762728079305, 0.7798784815692744, 0.37850223028772895]
v: 0.125
- phase: A
O: [0.17194138099812353, 0.48707029434265314, 0.28918455812138427, -0.8059596647559721]
v: 0.125
- phase: A
O: [0.7411029340880614, -0.6424443807902435, 0.06269338980849715, 0.18466509565000905]
v: 0.125
homogenization: RGC
- constituents:
- phase: A
O: [0.2238963963710994, -0.0650586111271435, -0.4249089470195672, -0.8746943280672199]
v: 0.125
- phase: A
O: [0.11834025356863993, -0.8856449810417544, 0.417012782917471, 0.16651994122112387]
v: 0.125
- phase: A
O: [0.11716309944373668, 0.0038756181092136927, -0.7000966032560261, -0.7043596622623864]
v: 0.125
- phase: A
O: [0.45702665319142677, -0.16783958203387148, -0.15365667123574342, 0.8598523945190182]
v: 0.125
- phase: A
O: [0.5197700015250588, 0.41758011415389684, -0.06820554688133636, 0.7421684425738383]
v: 0.125
- phase: A
O: [0.675852597832269, -0.5552996378928718, 0.4213174293100185, 0.23949363648960756]
v: 0.125
- phase: A
O: [0.5936583232008587, 0.6945289973655595, -0.40553899121228437, -0.027154994370431035]
v: 0.125
- phase: A
O: [0.1865986177378815, 0.3352356341235895, 0.35399611880081633, -0.8529271793922534]
v: 0.125
homogenization: RGC
- constituents:
- phase: A
O: [0.04264635665959766, -0.25662858077651657, 0.9118921459098731, -0.3174520027030545]
v: 0.125
- phase: A
O: [0.8821155978581948, 0.3761807268996467, 0.1631361963837017, 0.23183337584133834]
v: 0.125
- phase: A
O: [0.162154410316598, -0.8783087060961443, 0.0990284993753847, -0.4387175860642615]
v: 0.125
- phase: A
O: [0.5420568096555354, 0.5088379086577298, 0.34515629818391275, -0.572822422433749]
v: 0.125
- phase: A
O: [0.4535105167294993, -0.351824950997862, 0.5869707747562448, -0.5709752399650516]
v: 0.125
- phase: A
O: [0.8098666527705685, 0.36698878904202453, -0.4428476393801139, 0.11541751055678014]
v: 0.125
- phase: A
O: [0.2576165318782538, 0.537377864604212, -0.6660199792713344, 0.44863809506978913]
v: 0.125
- phase: A
O: [0.665871020774723, -0.24779933206030424, -0.6964987066828571, 0.10050286718299539]
v: 0.125
homogenization: RGC
- constituents:
- phase: A
O: [0.5998961636312351, -0.5496518322385892, 0.10155928244475286, 0.5724449041843198]
v: 0.125
- phase: A
O: [0.7039639404137236, 0.03679879173607231, 0.58207433274894, 0.4053024681380869]
v: 0.125
- phase: A
O: [0.40054268633071116, 0.37951399049616763, 0.5867236594596857, -0.5926972539795397]
v: 0.125
- phase: A
O: [0.9355573844663049, 0.24834752099784146, 0.19019564545030732, 0.16395580391231754]
v: 0.125
- phase: A
O: [0.7776851244299001, 0.5132289287723679, -0.3500697782737797, 0.0961928492714775]
v: 0.125
- phase: A
O: [0.14184146561701433, -0.12106060201428649, 0.06736696083733706, 0.9801464287845448]
v: 0.125
- phase: A
O: [0.35378845612452303, 0.6454299208817863, -0.6758747727133072, 0.03804257027716481]
v: 0.125
- phase: A
O: [0.5450137677982273, -0.6067354536251314, -0.5773779505037945, 0.03829862264786783]
v: 0.125
homogenization: RGC
- constituents:
- phase: A
O: [0.8960165930555501, 0.24799888375444512, -0.3529216677489853, 0.10534284531447227]
v: 0.125
- phase: A
O: [0.018356029815615817, -0.2889337455527527, 0.9520448827103584, 0.09894891689798985]
v: 0.125
- phase: A
O: [0.853080700130541, -0.3850382673796441, 0.3507108308253732, 0.03163486778610022]
v: 0.125
- phase: A
O: [0.3115000788768896, -0.36279737244767013, 0.6484519555216468, -0.5923308440263011]
v: 0.125
- phase: A
O: [0.44453200197416054, 0.5366087739174971, -0.4155131463908387, -0.5846290688564767]
v: 0.125
- phase: A
O: [0.6294708415113138, -0.5730458394479108, 0.3971797612768184, -0.34297691294104227]
v: 0.125
- phase: A
O: [0.7012827551365688, -0.6926674077585473, -0.1570646347205316, -0.06119689614043776]
v: 0.125
- phase: A
O: [0.344127722728084, -0.02877253391263471, -0.12761292039202274, 0.9297651285627186]
v: 0.125
homogenization: RGC
- constituents:
- phase: A
O: [0.6978484686818835, 0.22305447962726752, 0.026182556001199425, -0.6801240237175888]
v: 0.125
- phase: A
O: [0.9628977633574072, 0.16069766981498335, 0.2167078702347067, -0.006469560701840579]
v: 0.125
- phase: A
O: [0.5032342094162214, -0.32097872787787546, 0.8022984104181043, -0.006726616067118365]
v: 0.125
- phase: A
O: [0.8004358150893949, 0.0875805156227694, -0.3568768378095974, -0.47357267851983204]
v: 0.125
- phase: A
O: [0.5880748138588325, 0.45495016865260346, 0.5556078663040557, -0.37214010298397254]
v: 0.125
- phase: A
O: [0.0242946564098417, 0.6005127022483121, 0.7620365834257689, -0.24102802664656872]
v: 0.125
- phase: A
O: [0.34550850696510604, 0.7299878005482168, -0.29949045445085315, -0.5079834154363129]
v: 0.125
- phase: A
O: [0.10265979489126274, 0.34339772969920185, 0.0799881993576351, 0.9301294822302113]
v: 0.125
homogenization: RGC
- constituents:
- phase: A
O: [0.5849010547063512, -0.46890126598541493, -0.141130615667434, -0.6466100125129551]
v: 0.125
- phase: A
O: [0.3232591731614976, -0.2858051445366468, 0.8855713416105193, -0.17199513144701628]
v: 0.125
- phase: A
O: [0.7530648506011821, -0.3316323026169325, -0.30797159948874325, -0.477563441396382]
v: 0.125
- phase: A
O: [0.5841938587340145, -0.5195581188838846, 0.45610755437383976, 0.4251385601923404]
v: 0.125
- phase: A
O: [0.5699287377636442, -0.7161652717835199, -0.40080048103879345, -0.04058955236816678]
v: 0.125
- phase: A
O: [0.05536551119448681, -0.40478746018638956, -0.6031583160532334, -0.6850414717532457]
v: 0.125
- phase: A
O: [0.2227977302479141, -0.48882932798587353, -0.15904170725500447, -0.8283192590122909]
v: 0.125
- phase: A
O: [0.05275211816268626, -0.223660120175789, 0.6943402501868806, 0.6819713935662708]
v: 0.125
homogenization: RGC
- constituents:
- phase: A
O: [0.7102673449281939, 0.43290229875199215, 0.4268837400347354, 0.35480441225815024]
v: 0.125
- phase: A
O: [0.2541648349387051, 0.3725842249682419, -0.13402399303768983, -0.8823937903655196]
v: 0.125
- phase: A
O: [0.4091596391203177, 0.10748924891263324, 0.2807193957630742, -0.8615283349522196]
v: 0.125
- phase: A
O: [0.16364908456142396, -0.9217193970202214, 0.11860582627366922, 0.33103623404821947]
v: 0.125
- phase: A
O: [0.80868984786943, -0.5530964468638372, 0.18400280484186798, -0.07904440669549186]
v: 0.125
- phase: A
O: [0.785539568992473, -0.42139749663876175, 0.03451711851851934, 0.45184101618034067]
v: 0.125
- phase: A
O: [0.07920711616878613, -0.6177377682037178, -0.7622215143727179, 0.1764784562213612]
v: 0.125
- phase: A
O: [0.6079704911212058, -0.06435941038671157, 0.3854845277687917, 0.6911088388028229]
v: 0.125
homogenization: RGC
- constituents:
- phase: A
O: [0.6650964387751211, -0.34961073112696767, -0.4316460762174112, -0.4990999185490129]
v: 0.125
- phase: A
O: [0.1643442326674122, 0.21205555067771598, 0.18237328729946675, -0.9459193415378059]
v: 0.125
- phase: A
O: [0.5079103592861302, -0.4522879915621313, 0.7008603696203147, 0.21507529359320504]
v: 0.125
- phase: A
O: [0.4588629674996609, 0.7103482146539822, -0.14094577428597305, -0.5147664321867083]
v: 0.125
- phase: A
O: [0.044305288641285946, -0.3875768281048735, -0.7143901476515312, -0.5809199261972351]
v: 0.125
- phase: A
O: [0.8455643302581954, -0.4342346154649184, -0.29235576797838947, -0.10483018199359342]
v: 0.125
- phase: A
O: [0.5790795198800359, 0.5645942182636775, -0.26425181035873874, 0.5254248367567554]
v: 0.125
- phase: A
O: [0.6475654656469717, 0.04725837905512413, 0.4718561325845725, -0.5964707901086468]
v: 0.125
homogenization: RGC
- constituents:
- phase: A
O: [0.9334013387182871, -0.3485154208652402, -0.04006434439320742, 0.07545721043330728]
v: 0.125
- phase: A
O: [0.38741862953083567, -0.9054265066714982, 0.12281627778824236, 0.12257980428821813]
v: 0.125
- phase: A
O: [0.8578234268453049, 0.25286192172706, 0.21735917871049656, -0.3910943675459601]
v: 0.125
- phase: A
O: [0.4841852989554279, -0.5324884938481191, 0.6128392956763806, -0.32626461326610684]
v: 0.125
- phase: A
O: [0.24479559088002964, -0.39790158622306016, -0.6885075670842895, 0.5547132380199182]
v: 0.125
- phase: A
O: [0.006342119419357721, 0.8016474931043842, -0.09875135732320732, 0.589550034982232]
v: 0.125
- phase: A
O: [0.5815233878895042, -0.41292157227416415, -0.6523344380499, -0.2564880219859532]
v: 0.125
- phase: A
O: [0.08441325811399321, -0.9287379499833768, -0.04047714959468877, 0.3587224867163255]
v: 0.125
homogenization: RGC
- constituents:
- phase: A
O: [0.9272909114072972, 0.22505658735105372, -0.12087397257258231, 0.27362489080097285]
v: 0.125
- phase: A
O: [0.22477798467053944, -0.16730490431874273, 0.7619888558500476, 0.5838295214860949]
v: 0.125
- phase: A
O: [0.6944164222400087, -0.3089622138004769, 0.47430811205294454, -0.44425217816873547]
v: 0.125
- phase: A
O: [0.9319464859005648, -0.33314612927635295, 0.09454266632036588, 0.10747598899664906]
v: 0.125
- phase: A
O: [0.7631364561921935, 0.6348137033139695, -0.07779160762534851, 0.09264327875398014]
v: 0.125
- phase: A
O: [0.30213352773158436, 0.7607414448583625, -0.5489414690920735, -0.1692662075144206]
v: 0.125
- phase: A
O: [0.44393342791682283, -0.20509632166957376, 0.6929260547202925, 0.529822699688679]
v: 0.125
- phase: A
O: [0.5553300723266665, 0.8100017529476545, -0.18277171291698383, 0.045827633026130514]
v: 0.125
homogenization: RGC
- constituents:
- phase: A
O: [0.3369089353527777, 0.11156804371056815, -0.8851691518465128, -0.30086627182417736]
v: 0.125
- phase: A
O: [0.1904884889827469, -0.2955994647758964, -0.5609071888421532, 0.7494786304455029]
v: 0.125
- phase: A
O: [0.14357497570029057, -0.8835693969815234, -0.07717841484904021, -0.4390157620766682]
v: 0.125
- phase: A
O: [0.2844319189561915, -0.26616939939342976, -0.09666838434612286, 0.9159189689996324]
v: 0.125
- phase: A
O: [0.3986740568181236, -0.009203742364828175, 0.3059441820684439, 0.8645070531841439]
v: 0.125
- phase: A
O: [0.6689177202684424, -0.4911686233022882, 0.011604901743798118, -0.5578241597938561]
v: 0.125
- phase: A
O: [0.24260666144446064, -0.7362708394584294, 0.555535585659367, 0.3007116091075576]
v: 0.125
- phase: A
O: [0.8264968442974631, 0.08346800006449286, -0.5223744566842243, 0.19251230177687406]
v: 0.125
homogenization: RGC
- constituents:
- phase: A
O: [0.5964950823065459, -0.03914167887903209, -0.764986315653944, 0.23970290490697482]
v: 0.125
- phase: A
O: [0.8080255887350464, 0.369717513486498, 0.4541656797460721, 0.06432063052809132]
v: 0.125
- phase: A
O: [0.9384315445539018, 0.21987525747326847, 0.26439717033746396, 0.033094465621676714]
v: 0.125
- phase: A
O: [0.3398146205195601, 0.19417091506154938, 0.7316894460339152, 0.5580808489707294]
v: 0.125
- phase: A
O: [0.3246633791904426, 0.47411988229633434, 0.8166638950637336, 0.05351737963768399]
v: 0.125
- phase: A
O: [0.3795740054462351, -0.8471170612518093, 0.01606588144238899, 0.37156176657331014]
v: 0.125
- phase: A
O: [0.3324564774831089, 0.5840719896629913, 0.33174639891221797, 0.6620248698345199]
v: 0.125
- phase: A
O: [0.5142455012854315, 0.2115524033587584, -0.25751707686577474, -0.7902417985422785]
v: 0.125
homogenization: RGC
- constituents:
- phase: A
O: [0.5917389153130538, 0.014487837662305996, 0.1546249763526937, -0.7910286185416618]
v: 0.125
- phase: A
O: [0.13783983489759696, 0.19317471156196087, 0.9651184311351051, -0.11058989380440512]
v: 0.125
- phase: A
O: [0.8586698401906224, 0.4868991116453938, -0.04218722674591906, 0.15438781857849304]
v: 0.125
- phase: A
O: [0.9890970857322728, 0.05214494425397627, -0.1362936500688951, -0.019796482909146578]
v: 0.125
- phase: A
O: [0.9020138640071081, 0.25707037984798153, 0.30534577613336306, -0.16446813047303388]
v: 0.125
- phase: A
O: [0.3175523590958934, -0.563302190026314, -0.6464135983869942, -0.40496987760149405]
v: 0.125
- phase: A
O: [0.752248294923397, -0.30191235437425856, 0.44098110520789696, -0.3853661867764942]
v: 0.125
- phase: A
O: [0.07541991372724928, 0.21981465995106916, -0.8076581447316371, 0.5419240473835979]
v: 0.125
homogenization: RGC
- constituents:
- phase: A
O: [0.3735354245257942, -0.19362809638145217, 0.7860424876438061, 0.45289806196843757]
v: 0.125
- phase: A
O: [0.0320727653734856, -0.8866094001869422, 0.46133859224884993, 0.007862094078374177]
v: 0.125
- phase: A
O: [0.4437387774582536, 0.2752317472571834, 0.5589728539449029, -0.6441216742466462]
v: 0.125
- phase: A
O: [0.4189270336917096, -0.7997419507282093, 0.08375652034255085, 0.4217793238031132]
v: 0.125
- phase: A
O: [0.07774372790271536, 0.8494512032072281, -0.23549059135564554, 0.4657603971191082]
v: 0.125
- phase: A
O: [0.7403304241587872, 0.5263542470691401, -0.1337642666383515, -0.39619337529526266]
v: 0.125
- phase: A
O: [0.3862324493674717, -0.7790044035666689, 0.2926767318730596, 0.39789064439798927]
v: 0.125
- phase: A
O: [0.12717246708576896, 0.2790094483039878, 0.8667779585068207, -0.39328979394229346]
v: 0.125
homogenization: RGC
homogenization:
RGC:
N_constituents: 8
mechanical:
type: RGC
D_alpha: [4.0e-06, 4.0e-06, 2.0e-06]
a_g: [0.0, 0.0, 0.0]
c_alpha: 2.0
cluster_size: [2, 2, 2]
output: [M, Delta_V, avg_a_dot, max_a_dot]
xi_alpha: 10.0
phase:
A:
lattice: cF
mechanical:
output: [F, F_e, F_p, L_p]
elastic: {C_11: 106.75e9, C_12: 60.41e9, C_44: 28.34e9, type: hooke}
plastic:
N_sl: [12]
a_sl: 2.25
atol_xi: 1.0
dot_gamma_0_sl: 0.001
h_0_sl_sl: 75e6
h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4]
n_sl: 20
output: [xi_sl]
type: phenopowerlaw
xi_0_sl: [31e6]
xi_inf_sl: [63e6]
B:
lattice: cI
mechanical:
output: [F, P, F_e, F_p, L_p, O]
elastic: {C_11: 106.75e9, C_12: 60.41e9, C_44: 28.34e9, type: hooke}
plastic:
N_sl: [12]
a_sl: 2.25
atol_xi: 1.0
dot_gamma_0_sl: 0.001
h_0_sl_sl: 75e6
h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4]
n_sl: 20
output: [xi_sl]
type: phenopowerlaw
xi_0_sl: [31e6]
xi_inf_sl: [63e6]
C:
lattice: cI
mechanical:
output: [F]
elastic: {C_11: 106.75e9, C_12: 60.41e9, C_44: 28.34e9, type: hooke}

View File

@ -0,0 +1,30 @@
<?xml version="1.0"?>
<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
<RectilinearGrid WholeExtent="0 2 0 4 0 3">
<FieldData>
<Array type="String" Name="comments" NumberOfTuples="0" format="binary">
AAAAAACAAAAAAAAA
</Array>
</FieldData>
<Piece Extent="0 2 0 4 0 3">
<PointData>
</PointData>
<CellData>
<DataArray type="Int64" Name="material" format="binary" RangeMin="0" RangeMax="3">
AQAAAACAAADAAAAAIgAAAA==eF5jZIAARjSamQAfXRwdoMszoYnD+DBAyBx0GgYADegAHw==
</DataArray>
</CellData>
<Coordinates>
<DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="1">
AQAAAACAAAAYAAAAEQAAAA==eF5jYEAGD+wh9Ad7AA0uAk8=
</DataArray>
<DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="1">
AQAAAACAAAAoAAAAFwAAAA==eF5jYEAGF+wh9AMo/QJKf7AHADzEBIU=
</DataArray>
<DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="1">
AQAAAACAAAAgAAAAGAAAAA==eF5jYICAUDC4ag+hn9pDRD/YAwBmSwdk
</DataArray>
</Coordinates>
</Piece>
</RectilinearGrid>
</VTKFile>

View File

@ -0,0 +1,17 @@
---
solver:
mechanical: spectral_basic
loadstep:
- boundary_conditions:
mechanical:
dot_F: [x, 0, 0,
0, -1.0e-3, 0,
0, 0, x]
P: [0, x, x,
x, x, x,
x, x, 0]
discretization:
t: 5
N: 10
f_out: 2

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

View File

@ -0,0 +1 @@
3b83384def67552ab7dd211efc0d54fd

View File

@ -0,0 +1 @@
c32c86ed50dbb39a93ca2a2ebe47d9cb

View File

@ -0,0 +1 @@
ead4f6fcaff174fddc041d701e54ac60

View File

@ -0,0 +1 @@
bde8b728110c2c05a6a4740f7c5f9c06

View File

@ -0,0 +1 @@
e09bfa9248283fc390003ad28d15d36e

View File

@ -0,0 +1 @@
3f21254164f96de8ee4a28249ae72cc6

View File

@ -136,9 +136,9 @@ class TestConfigMaterial:
def test_load_DREAM3D_reference(self,tmp_path,ref_path,update): def test_load_DREAM3D_reference(self,tmp_path,ref_path,update):
cur = ConfigMaterial.load_DREAM3D(ref_path/'measured.dream3d') cur = ConfigMaterial.load_DREAM3D(ref_path/'measured.dream3d')
ref = ConfigMaterial.load(ref_path/'measured.material_yaml') ref = ConfigMaterial.load(ref_path/'measured.material.yaml')
if update: if update:
cur.save(ref_path/'measured.material_yaml') cur.save(ref_path/'measured.material.yaml')
for i,m in enumerate(ref['material']): for i,m in enumerate(ref['material']):
assert Rotation(m['constituents'][0]['O']).isclose(Rotation(cur['material'][i]['constituents'][0]['O'])) assert Rotation(m['constituents'][0]['O']).isclose(Rotation(cur['material'][i]['constituents'][0]['O']))
assert cur.is_valid and cur['phase'] == ref['phase'] and cur['homogenization'] == ref['homogenization'] assert cur.is_valid and cur['phase'] == ref['phase'] and cur['homogenization'] == ref['homogenization']

View File

@ -14,8 +14,7 @@ from damask import grid_filters
def grid_equal(a,b): def grid_equal(a,b):
return np.all(a.material == b.material) and \ return np.all(a.material == b.material) and \
np.all(a.cells == b.cells) and \ np.all(a.cells == b.cells) and \
np.allclose(a.size, b.size) and \ np.allclose(a.size, b.size)
str(a.diff(b)) == str(b.diff(a))
@pytest.fixture @pytest.fixture
def default(): def default():
@ -42,13 +41,9 @@ class TestGrid:
def _patch_datetime_now(self, patch_datetime_now): def _patch_datetime_now(self, patch_datetime_now):
print('patched datetime.datetime.now') print('patched datetime.datetime.now')
def test_diff_equal(self,default):
assert str(default.diff(default)) == ''
def test_equal(self,default):
def test_diff_not_equal(self,default): assert default == default
new = Grid(default.material[1:,1:,1:]+1,default.size*.9,np.ones(3)-default.origin,comments=['modified'])
assert str(default.diff(new)) != ''
def test_repr(self,default): def test_repr(self,default):
print(default) print(default)

View File

@ -1,12 +1,14 @@
import bz2
import pickle
import time import time
import shutil import shutil
import os import os
import sys import sys
import hashlib
from datetime import datetime from datetime import datetime
import pytest import pytest
import numpy as np import numpy as np
import h5py
from damask import Result from damask import Result
from damask import Rotation from damask import Rotation
@ -21,8 +23,7 @@ def default(tmp_path,ref_path):
fname = '12grains6x7x8_tensionY.hdf5' fname = '12grains6x7x8_tensionY.hdf5'
shutil.copy(ref_path/fname,tmp_path) shutil.copy(ref_path/fname,tmp_path)
f = Result(tmp_path/fname) f = Result(tmp_path/fname)
f.view('times',20.0) return f.view('times',20.0)
return f
@pytest.fixture @pytest.fixture
def single_phase(tmp_path,ref_path): def single_phase(tmp_path,ref_path):
@ -36,6 +37,17 @@ def ref_path(ref_path_base):
"""Directory containing reference results.""" """Directory containing reference results."""
return ref_path_base/'Result' return ref_path_base/'Result'
def dict_equal(d1, d2):
for k in d1:
if (k not in d2):
return False
else:
if type(d1[k]) is dict:
return dict_equal(d1[k],d2[k])
else:
if not np.allclose(d1[k],d2[k]):
return False
return True
class TestResult: class TestResult:
@ -44,51 +56,39 @@ class TestResult:
def test_view_all(self,default): def test_view_all(self,default):
default.view('increments',True) a = default.view('increments',True).get('F')
a = default.get_dataset_location('F')
default.view('increments','*')
b = default.get_dataset_location('F')
default.view('increments',default.increments_in_range(0,np.iinfo(int).max))
c = default.get_dataset_location('F')
default.view('times',True) assert dict_equal(a,default.view('increments','*').get('F'))
d = default.get_dataset_location('F') assert dict_equal(a,default.view('increments',default.increments_in_range(0,np.iinfo(int).max)).get('F'))
default.view('times','*')
e = default.get_dataset_location('F') assert dict_equal(a,default.view('times',True).get('F'))
default.view('times',default.times_in_range(0.0,np.inf)) assert dict_equal(a,default.view('times','*').get('F'))
f = default.get_dataset_location('F') assert dict_equal(a,default.view('times',default.times_in_range(0.0,np.inf)).get('F'))
assert a == b == c == d == e ==f
@pytest.mark.parametrize('what',['increments','times','phases']) # ToDo: discuss homogenizations @pytest.mark.parametrize('what',['increments','times','phases']) # ToDo: discuss homogenizations
def test_view_none(self,default,what): def test_view_none(self,default,what):
default.view(what,False) a = default.view(what,False).get('F')
a = default.get_dataset_location('F') b = default.view(what,[]).get('F')
default.view(what,[])
b = default.get_dataset_location('F')
assert a == b == [] assert a == b == {}
@pytest.mark.parametrize('what',['increments','times','phases']) # ToDo: discuss homogenizations @pytest.mark.parametrize('what',['increments','times','phases']) # ToDo: discuss homogenizations
def test_view_more(self,default,what): def test_view_more(self,default,what):
default.view(what,False) empty = default.view(what,False)
default.view_more(what,'*')
a = default.get_dataset_location('F')
default.view(what,True) a = empty.view_more(what,'*').get('F')
b = default.get_dataset_location('F') b = empty.view_more(what,True).get('F')
assert a == b assert dict_equal(a,b)
@pytest.mark.parametrize('what',['increments','times','phases']) # ToDo: discuss homogenizations @pytest.mark.parametrize('what',['increments','times','phases']) # ToDo: discuss homogenizations
def test_view_less(self,default,what): def test_view_less(self,default,what):
default.view(what,True) full = default.view(what,True)
default.view_less(what,'*')
a = default.get_dataset_location('F')
default.view(what,False) a = full.view_less(what,'*').get('F')
b = default.get_dataset_location('F') b = full.view_less(what,True).get('F')
assert a == b == [] assert a == b == {}
def test_view_invalid(self,default): def test_view_invalid(self,default):
with pytest.raises(AttributeError): with pytest.raises(AttributeError):
@ -96,10 +96,8 @@ class TestResult:
def test_add_absolute(self,default): def test_add_absolute(self,default):
default.add_absolute('F_e') default.add_absolute('F_e')
loc = {'F_e': default.get_dataset_location('F_e'), in_memory = np.abs(default.place('F_e'))
'|F_e|': default.get_dataset_location('|F_e|')} in_file = default.place('|F_e|')
in_memory = np.abs(default.read_dataset(loc['F_e'],0))
in_file = default.read_dataset(loc['|F_e|'],0)
assert np.allclose(in_memory,in_file) assert np.allclose(in_memory,in_file)
@pytest.mark.parametrize('mode',['direct','function']) @pytest.mark.parametrize('mode',['direct','function'])
@ -115,77 +113,59 @@ class TestResult:
default.enable_user_function(f.my_func) default.enable_user_function(f.my_func)
default.add_calculation('x','my_func(#F#)','-','my notes') default.add_calculation('x','my_func(#F#)','-','my notes')
loc = {'F': default.get_dataset_location('F'), in_memory = 2.0*np.abs(default.place('F'))-1.0
'x': default.get_dataset_location('x')} in_file = default.place('x')
in_memory = 2.0*np.abs(default.read_dataset(loc['F'],0))-1.0
in_file = default.read_dataset(loc['x'],0)
assert np.allclose(in_memory,in_file) assert np.allclose(in_memory,in_file)
def test_add_stress_Cauchy(self,default): def test_add_stress_Cauchy(self,default):
default.add_stress_Cauchy('P','F') default.add_stress_Cauchy('P','F')
loc = {'F': default.get_dataset_location('F'), in_memory = mechanics.stress_Cauchy(default.place('P'), default.place('F'))
'P': default.get_dataset_location('P'), in_file = default.place('sigma')
'sigma':default.get_dataset_location('sigma')}
in_memory = mechanics.stress_Cauchy(default.read_dataset(loc['P'],0),
default.read_dataset(loc['F'],0))
in_file = default.read_dataset(loc['sigma'],0)
assert np.allclose(in_memory,in_file) assert np.allclose(in_memory,in_file)
def test_add_determinant(self,default): def test_add_determinant(self,default):
default.add_determinant('P') default.add_determinant('P')
loc = {'P': default.get_dataset_location('P'), in_memory = np.linalg.det(default.place('P'))
'det(P)':default.get_dataset_location('det(P)')} in_file = default.place('det(P)')
in_memory = np.linalg.det(default.read_dataset(loc['P'],0)).reshape(-1,1)
in_file = default.read_dataset(loc['det(P)'],0)
assert np.allclose(in_memory,in_file) assert np.allclose(in_memory,in_file)
def test_add_deviator(self,default): def test_add_deviator(self,default):
default.add_deviator('P') default.add_deviator('P')
loc = {'P' :default.get_dataset_location('P'), in_memory = tensor.deviatoric(default.place('P'))
's_P':default.get_dataset_location('s_P')} in_file = default.place('s_P')
in_memory = tensor.deviatoric(default.read_dataset(loc['P'],0))
in_file = default.read_dataset(loc['s_P'],0)
assert np.allclose(in_memory,in_file) assert np.allclose(in_memory,in_file)
@pytest.mark.parametrize('eigenvalue,function',[('max',np.amax),('min',np.amin)]) @pytest.mark.parametrize('eigenvalue,function',[('max',np.amax),('min',np.amin)])
def test_add_eigenvalue(self,default,eigenvalue,function): def test_add_eigenvalue(self,default,eigenvalue,function):
default.add_stress_Cauchy('P','F') default.add_stress_Cauchy('P','F')
default.add_eigenvalue('sigma',eigenvalue) default.add_eigenvalue('sigma',eigenvalue)
loc = {'sigma' :default.get_dataset_location('sigma'), in_memory = function(tensor.eigenvalues(default.place('sigma')),axis=1)
'lambda':default.get_dataset_location(f'lambda_{eigenvalue}(sigma)')} in_file = default.place(f'lambda_{eigenvalue}(sigma)')
in_memory = function(tensor.eigenvalues(default.read_dataset(loc['sigma'],0)),axis=1,keepdims=True)
in_file = default.read_dataset(loc['lambda'],0)
assert np.allclose(in_memory,in_file) assert np.allclose(in_memory,in_file)
@pytest.mark.parametrize('eigenvalue,idx',[('max',2),('mid',1),('min',0)]) @pytest.mark.parametrize('eigenvalue,idx',[('max',2),('mid',1),('min',0)])
def test_add_eigenvector(self,default,eigenvalue,idx): def test_add_eigenvector(self,default,eigenvalue,idx):
default.add_stress_Cauchy('P','F') default.add_stress_Cauchy('P','F')
default.add_eigenvector('sigma',eigenvalue) default.add_eigenvector('sigma',eigenvalue)
loc = {'sigma' :default.get_dataset_location('sigma'), in_memory = tensor.eigenvectors(default.place('sigma'))[:,idx]
'v(sigma)':default.get_dataset_location(f'v_{eigenvalue}(sigma)')} in_file = default.place(f'v_{eigenvalue}(sigma)')
in_memory = tensor.eigenvectors(default.read_dataset(loc['sigma'],0))[:,idx]
in_file = default.read_dataset(loc['v(sigma)'],0)
assert np.allclose(in_memory,in_file) assert np.allclose(in_memory,in_file)
@pytest.mark.parametrize('d',[[1,0,0],[0,1,0],[0,0,1]]) @pytest.mark.parametrize('d',[[1,0,0],[0,1,0],[0,0,1]])
def test_add_IPF_color(self,default,d): def test_add_IPF_color(self,default,d):
default.add_IPF_color(d,'O') default.add_IPF_color(d,'O')
loc = {'O': default.get_dataset_location('O'), qu = default.place('O')
'color': default.get_dataset_location('IPFcolor_[{} {} {}]'.format(*d))} crystal_structure = qu.dtype.metadata['lattice']
qu = default.read_dataset(loc['O']).view(np.double).squeeze()
crystal_structure = default._get_attribute(default.get_dataset_location('O')[0],'lattice')
c = Orientation(rotation=qu,lattice=crystal_structure) c = Orientation(rotation=qu,lattice=crystal_structure)
in_memory = np.uint8(c.IPF_color(np.array(d))*255) in_memory = np.uint8(c.IPF_color(np.array(d))*255)
in_file = default.read_dataset(loc['color']) in_file = default.place('IPFcolor_({} {} {})'.format(*d))
assert np.allclose(in_memory,in_file) assert np.allclose(in_memory,in_file)
def test_add_maximum_shear(self,default): def test_add_maximum_shear(self,default):
default.add_stress_Cauchy('P','F') default.add_stress_Cauchy('P','F')
default.add_maximum_shear('sigma') default.add_maximum_shear('sigma')
loc = {'sigma' :default.get_dataset_location('sigma'), in_memory = mechanics.maximum_shear(default.place('sigma'))
'max_shear(sigma)':default.get_dataset_location('max_shear(sigma)')} in_file = default.place('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) assert np.allclose(in_memory,in_file)
def test_add_Mises_strain(self,default): def test_add_Mises_strain(self,default):
@ -194,26 +174,22 @@ class TestResult:
default.add_strain('F',t,m) default.add_strain('F',t,m)
label = f'epsilon_{t}^{m}(F)' label = f'epsilon_{t}^{m}(F)'
default.add_equivalent_Mises(label) default.add_equivalent_Mises(label)
loc = {label :default.get_dataset_location(label), in_memory = mechanics.equivalent_strain_Mises(default.place(label))
label+'_vM':default.get_dataset_location(label+'_vM')} in_file = default.place(label+'_vM')
in_memory = mechanics.equivalent_strain_Mises(default.read_dataset(loc[label],0)).reshape(-1,1)
in_file = default.read_dataset(loc[label+'_vM'],0)
assert np.allclose(in_memory,in_file) assert np.allclose(in_memory,in_file)
def test_add_Mises_stress(self,default): def test_add_Mises_stress(self,default):
default.add_stress_Cauchy('P','F') default.add_stress_Cauchy('P','F')
default.add_equivalent_Mises('sigma') default.add_equivalent_Mises('sigma')
loc = {'sigma' :default.get_dataset_location('sigma'), in_memory = mechanics.equivalent_stress_Mises(default.place('sigma'))
'sigma_vM':default.get_dataset_location('sigma_vM')} in_file = default.place('sigma_vM')
in_memory = mechanics.equivalent_stress_Mises(default.read_dataset(loc['sigma'],0)).reshape(-1,1)
in_file = default.read_dataset(loc['sigma_vM'],0)
assert np.allclose(in_memory,in_file) assert np.allclose(in_memory,in_file)
def test_add_Mises_invalid(self,default): def test_add_Mises_invalid(self,default):
default.add_stress_Cauchy('P','F') default.add_stress_Cauchy('P','F')
default.add_calculation('sigma_y','#sigma#',unit='y') default.add_calculation('sigma_y','#sigma#',unit='y')
default.add_equivalent_Mises('sigma_y') default.add_equivalent_Mises('sigma_y')
assert default.get_dataset_location('sigma_y_vM') == [] assert default.get('sigma_y_vM') == {}
def test_add_Mises_stress_strain(self,default): def test_add_Mises_stress_strain(self,default):
default.add_stress_Cauchy('P','F') default.add_stress_Cauchy('P','F')
@ -221,26 +197,20 @@ class TestResult:
default.add_calculation('sigma_x','#sigma#',unit='x') default.add_calculation('sigma_x','#sigma#',unit='x')
default.add_equivalent_Mises('sigma_y',kind='strain') default.add_equivalent_Mises('sigma_y',kind='strain')
default.add_equivalent_Mises('sigma_x',kind='stress') default.add_equivalent_Mises('sigma_x',kind='stress')
loc = {'y' :default.get_dataset_location('sigma_y_vM'), assert not np.allclose(default.place('sigma_y_vM'),default.place('sigma_x_vM'))
'x' :default.get_dataset_location('sigma_x_vM')}
assert not np.allclose(default.read_dataset(loc['y'],0),default.read_dataset(loc['x'],0))
def test_add_norm(self,default): @pytest.mark.parametrize('ord',[1,2])
default.add_norm('F',1) @pytest.mark.parametrize('dataset,axis',[('F',(1,2)),('xi_sl',(1,))])
loc = {'F': default.get_dataset_location('F'), def test_add_norm(self,default,ord,dataset,axis):
'|F|_1':default.get_dataset_location('|F|_1')} default.add_norm(dataset,ord)
in_memory = np.linalg.norm(default.read_dataset(loc['F'],0),ord=1,axis=(1,2),keepdims=True) in_memory = np.linalg.norm(default.place(dataset),ord=ord,axis=axis,keepdims=True)
in_file = default.read_dataset(loc['|F|_1'],0) in_file = default.place(f'|{dataset}|_{ord}')
assert np.allclose(in_memory,in_file) assert np.allclose(in_memory,in_file)
def test_add_stress_second_Piola_Kirchhoff(self,default): def test_add_stress_second_Piola_Kirchhoff(self,default):
default.add_stress_second_Piola_Kirchhoff('P','F') default.add_stress_second_Piola_Kirchhoff('P','F')
loc = {'F':default.get_dataset_location('F'), in_memory = mechanics.stress_second_Piola_Kirchhoff(default.place('P'),default.place('F'))
'P':default.get_dataset_location('P'), in_file = default.place('S')
'S':default.get_dataset_location('S')}
in_memory = mechanics.stress_second_Piola_Kirchhoff(default.read_dataset(loc['P'],0),
default.read_dataset(loc['F'],0))
in_file = default.read_dataset(loc['S'],0)
assert np.allclose(in_memory,in_file) assert np.allclose(in_memory,in_file)
@pytest.mark.skip(reason='requires rework of lattice.f90') @pytest.mark.skip(reason='requires rework of lattice.f90')
@ -248,30 +218,24 @@ class TestResult:
def test_add_pole(self,default,polar): def test_add_pole(self,default,polar):
pole = np.array([1.,0.,0.]) pole = np.array([1.,0.,0.])
default.add_pole('O',pole,polar) default.add_pole('O',pole,polar)
loc = {'O': default.get_dataset_location('O'), rot = Rotation(default.place('O'))
'pole': default.get_dataset_location('p^{}_[1 0 0)'.format(u'' if polar else 'xy'))}
rot = Rotation(default.read_dataset(loc['O']).view(np.double))
rotated_pole = rot * np.broadcast_to(pole,rot.shape+(3,)) rotated_pole = rot * np.broadcast_to(pole,rot.shape+(3,))
xy = rotated_pole[:,0:2]/(1.+abs(pole[2])) xy = rotated_pole[:,0:2]/(1.+abs(pole[2]))
in_memory = xy if not polar else \ in_memory = xy if not polar else \
np.block([np.sqrt(xy[:,0:1]*xy[:,0:1]+xy[:,1:2]*xy[:,1:2]),np.arctan2(xy[:,1:2],xy[:,0:1])]) np.block([np.sqrt(xy[:,0:1]*xy[:,0:1]+xy[:,1:2]*xy[:,1:2]),np.arctan2(xy[:,1:2],xy[:,0:1])])
in_file = default.read_dataset(loc['pole']) in_file = default.place('p^{}_[1 0 0)'.format(u'' if polar else 'xy'))
assert np.allclose(in_memory,in_file) assert np.allclose(in_memory,in_file)
def test_add_rotation(self,default): def test_add_rotation(self,default):
default.add_rotation('F') default.add_rotation('F')
loc = {'F': default.get_dataset_location('F'), in_memory = mechanics.rotation(default.place('F')).as_matrix()
'R(F)': default.get_dataset_location('R(F)')} in_file = default.place('R(F)')
in_memory = mechanics.rotation(default.read_dataset(loc['F'],0)).as_matrix()
in_file = default.read_dataset(loc['R(F)'],0)
assert np.allclose(in_memory,in_file) assert np.allclose(in_memory,in_file)
def test_add_spherical(self,default): def test_add_spherical(self,default):
default.add_spherical('P') default.add_spherical('P')
loc = {'P': default.get_dataset_location('P'), in_memory = tensor.spherical(default.place('P'),False)
'p_P': default.get_dataset_location('p_P')} in_file = default.place('p_P')
in_memory = tensor.spherical(default.read_dataset(loc['P'],0),False).reshape(-1,1)
in_file = default.read_dataset(loc['p_P'],0)
assert np.allclose(in_memory,in_file) assert np.allclose(in_memory,in_file)
def test_add_strain(self,default): def test_add_strain(self,default):
@ -279,26 +243,20 @@ class TestResult:
m = np.random.random()*2.0 - 1.0 m = np.random.random()*2.0 - 1.0
default.add_strain('F',t,m) default.add_strain('F',t,m)
label = f'epsilon_{t}^{m}(F)' label = f'epsilon_{t}^{m}(F)'
loc = {'F': default.get_dataset_location('F'), in_memory = mechanics.strain(default.place('F'),t,m)
label: default.get_dataset_location(label)} in_file = default.place(label)
in_memory = mechanics.strain(default.read_dataset(loc['F'],0),t,m)
in_file = default.read_dataset(loc[label],0)
assert np.allclose(in_memory,in_file) assert np.allclose(in_memory,in_file)
def test_add_stretch_right(self,default): def test_add_stretch_right(self,default):
default.add_stretch_tensor('F','U') default.add_stretch_tensor('F','U')
loc = {'F': default.get_dataset_location('F'), in_memory = mechanics.stretch_right(default.place('F'))
'U(F)': default.get_dataset_location('U(F)')} in_file = default.place('U(F)')
in_memory = mechanics.stretch_right(default.read_dataset(loc['F'],0))
in_file = default.read_dataset(loc['U(F)'],0)
assert np.allclose(in_memory,in_file) assert np.allclose(in_memory,in_file)
def test_add_stretch_left(self,default): def test_add_stretch_left(self,default):
default.add_stretch_tensor('F','V') default.add_stretch_tensor('F','V')
loc = {'F': default.get_dataset_location('F'), in_memory = mechanics.stretch_left(default.place('F'))
'V(F)': default.get_dataset_location('V(F)')} in_file = default.place('V(F)')
in_memory = mechanics.stretch_left(default.read_dataset(loc['F'],0))
in_file = default.read_dataset(loc['V(F)'],0)
assert np.allclose(in_memory,in_file) assert np.allclose(in_memory,in_file)
def test_add_invalid(self,default): def test_add_invalid(self,default):
@ -307,48 +265,40 @@ class TestResult:
@pytest.mark.parametrize('overwrite',['off','on']) @pytest.mark.parametrize('overwrite',['off','on'])
def test_add_overwrite(self,default,overwrite): def test_add_overwrite(self,default,overwrite):
default.view('times',default.times_in_range(0,np.inf)[-1]) last = default.view('times',default.times_in_range(0,np.inf)[-1])
default.add_stress_Cauchy() last.add_stress_Cauchy()
loc = default.get_dataset_location('sigma')
with h5py.File(default.fname,'r') as f: created_first = last.place('sigma').dtype.metadata['created']
# h5py3 compatibility
try:
created_first = f[loc[0]].attrs['created'].decode()
except AttributeError:
created_first = f[loc[0]].attrs['created']
created_first = datetime.strptime(created_first,'%Y-%m-%d %H:%M:%S%z') created_first = datetime.strptime(created_first,'%Y-%m-%d %H:%M:%S%z')
if overwrite == 'on': if overwrite == 'on':
default.allow_modification() last = last.allow_modification()
else: else:
default.disallow_modification() last = last.disallow_modification()
time.sleep(2.) time.sleep(2.)
try: try:
default.add_calculation('sigma','#sigma#*0.0+311.','not the Cauchy stress') last.add_calculation('sigma','#sigma#*0.0+311.','not the Cauchy stress')
except ValueError: except ValueError:
pass pass
with h5py.File(default.fname,'r') as f:
# h5py3 compatibility created_second = last.place('sigma').dtype.metadata['created']
try:
created_second = f[loc[0]].attrs['created'].decode()
except AttributeError:
created_second = f[loc[0]].attrs['created']
created_second = datetime.strptime(created_second,'%Y-%m-%d %H:%M:%S%z') created_second = datetime.strptime(created_second,'%Y-%m-%d %H:%M:%S%z')
if overwrite == 'on': if overwrite == 'on':
assert created_first < created_second and np.allclose(default.read_dataset(loc),311.) assert created_first < created_second and np.allclose(last.place('sigma'),311.)
else: else:
assert created_first == created_second and not np.allclose(default.read_dataset(loc),311.) assert created_first == created_second and not np.allclose(last.place('sigma'),311.)
@pytest.mark.parametrize('allowed',['off','on']) @pytest.mark.parametrize('allowed',['off','on'])
def test_rename(self,default,allowed): def test_rename(self,default,allowed):
if allowed == 'on': if allowed == 'on':
F = default.read_dataset(default.get_dataset_location('F')) F = default.place('F')
default.allow_modification() default = default.allow_modification()
default.rename('F','new_name') default.rename('F','new_name')
assert np.all(F == default.read_dataset(default.get_dataset_location('new_name'))) assert np.all(F == default.place('new_name'))
default.disallow_modification() default = default.disallow_modification()
with pytest.raises(PermissionError): with pytest.raises(PermissionError):
default.rename('P','another_new_name') default.rename('P','another_new_name')
@ -363,10 +313,30 @@ class TestResult:
b = default.coordinates0_node.reshape(tuple(default.cells+1)+(3,),order='F') b = default.coordinates0_node.reshape(tuple(default.cells+1)+(3,),order='F')
assert np.allclose(a,b) assert np.allclose(a,b)
@pytest.mark.parametrize('output',['F',[],['F','P']]) # need to wait for writing in parallel, output order might change if select more then one
def test_vtk(self,tmp_path,default,output): @pytest.mark.parametrize('output',['F','*',['P']],ids=range(3))
@pytest.mark.parametrize('fname',['12grains6x7x8_tensionY.hdf5'],ids=range(1))
@pytest.mark.parametrize('inc',[4,0],ids=range(2))
def test_vtk(self,request,tmp_path,ref_path,update,output,fname,inc):
result = Result(ref_path/fname).view('increments',inc)
os.chdir(tmp_path) os.chdir(tmp_path)
default.save_VTK(output) result.save_VTK(output)
fname = fname.split('.')[0]+f'_inc{(inc if type(inc) == int else inc[0]):0>2}.vtr'
last = ''
for i in range(10):
if os.path.isfile(tmp_path/fname):
with open(fname) as f:
cur = hashlib.md5(f.read().encode()).hexdigest()
if cur == last:
break
else:
last = cur
time.sleep(.5)
if update:
with open((ref_path/'save_VTK'/request.node.name).with_suffix('.md5'),'w') as f:
f.write(cur)
with open((ref_path/'save_VTK'/request.node.name).with_suffix('.md5')) as f:
assert cur == f.read()
@pytest.mark.parametrize('mode',['point','cell']) @pytest.mark.parametrize('mode',['point','cell'])
def test_vtk_mode(self,tmp_path,single_phase,mode): def test_vtk_mode(self,tmp_path,single_phase,mode):
@ -387,3 +357,52 @@ class TestResult:
def test_XDMF_invalid(self,default): def test_XDMF_invalid(self,default):
with pytest.raises(TypeError): with pytest.raises(TypeError):
default.save_XDMF() default.save_XDMF()
@pytest.mark.parametrize('view,output,flatten,prune',
[({},['F','P','F','L_p','F_e','F_p'],True,True),
({'increments':3},'F',True,True),
({'increments':[1,8,3,4,5,6,7]},['F','P'],True,True),
({'phases':['A','B']},['F','P'],True,True),
({'phases':['A','C'],'homogenizations':False},['F','P','O'],True,True),
({'phases':False,'homogenizations':False},['F','P','O'],True,True),
({'phases':False},['Delta_V'],True,True),
({},['u_p','u_n'],False,False)],
ids=list(range(8)))
def test_get(self,update,request,ref_path,view,output,flatten,prune):
result = Result(ref_path/'4grains2x4x3_compressionY.hdf5')
for key,value in view.items():
result = result.view(key,value)
fname = request.node.name
cur = result.get(output,flatten,prune)
if update:
with bz2.BZ2File((ref_path/'get'/fname).with_suffix('.pbz2'),'w') as f:
pickle.dump(cur,f)
with bz2.BZ2File((ref_path/'get'/fname).with_suffix('.pbz2')) as f:
assert dict_equal(cur,pickle.load(f))
@pytest.mark.parametrize('view,output,flatten,constituents,prune',
[({},['F','P','F','L_p','F_e','F_p'],True,True,None),
({'increments':3},'F',True,True,[0,1,2,3,4,5,6,7]),
({'increments':[1,8,3,4,5,6,7]},['F','P'],True,True,1),
({'phases':['A','B']},['F','P'],True,True,[1,2]),
({'phases':['A','C'],'homogenizations':False},['F','P','O'],True,True,[0,7]),
({'phases':False,'homogenizations':False},['F','P','O'],True,True,[1,2,3,4]),
({'phases':False},['Delta_V'],True,True,[1,2,4]),
({},['u_p','u_n'],False,False,None)],
ids=list(range(8)))
def test_place(self,update,request,ref_path,view,output,flatten,prune,constituents):
result = Result(ref_path/'4grains2x4x3_compressionY.hdf5')
for key,value in view.items():
result = result.view(key,value)
fname = request.node.name
cur = result.place(output,flatten,prune,constituents)
if update:
with bz2.BZ2File((ref_path/'place'/fname).with_suffix('.pbz2'),'w') as f:
pickle.dump(cur,f)
with bz2.BZ2File((ref_path/'place'/fname).with_suffix('.pbz2')) as f:
assert dict_equal(cur,pickle.load(f))

View File

@ -4,6 +4,7 @@ import time
import pytest import pytest
import numpy as np import numpy as np
import numpy.ma as ma
from damask import VTK from damask import VTK
from damask import grid_filters from damask import grid_filters
@ -134,6 +135,15 @@ class TestVTK:
with pytest.raises(TypeError): with pytest.raises(TypeError):
default.add('invalid_type','valid') default.add('invalid_type','valid')
def test_add_masked(self,default):
data = np.random.rand(5*6*7,3)
masked = ma.MaskedArray(data,mask=data<.4,fill_value=42.)
default.add(masked,'D')
result_masked = str(default)
default.add(np.where(masked.mask,masked.fill_value,masked),'D')
assert result_masked == str(default)
def test_comments(self,tmp_path,default): def test_comments(self,tmp_path,default):
default.add_comments(['this is a comment']) default.add_comments(['this is a comment'])
default.save(tmp_path/'with_comments',parallel=False) default.save(tmp_path/'with_comments',parallel=False)

View File

@ -136,3 +136,25 @@ class TestUtil:
else: else:
with pytest.raises(ValueError): with pytest.raises(ValueError):
util.DREAM3D_cell_data_group(tmp_path/'cell_data_group.dream3d') util.DREAM3D_cell_data_group(tmp_path/'cell_data_group.dream3d')
@pytest.mark.parametrize('full,reduced',[({}, {}),
({'A':{}}, {}),
({'A':{'B':{}}}, {}),
({'A':{'B':'C'}},)*2,
({'A':{'B':{},'C':'D'}}, {'A':{'C':'D'}})])
def test_prune(self,full,reduced):
assert util.dict_prune(full) == reduced
@pytest.mark.parametrize('full,reduced',[({}, {}),
({'A':{}}, {}),
({'A':'F'}, 'F'),
({'A':{'B':{}}}, {}),
({'A':{'B':'C'}}, 'C'),
({'A':1,'B':2},)*2,
({'A':{'B':'C','D':'E'}}, {'B':'C','D':'E'}),
({'B':'C','D':'E'},)*2,
({'A':{'B':{},'C':'D'}}, {'B':{},'C':'D'})])
def test_flatten(self,full,reduced):
assert util.dict_flatten(full) == reduced

View File

@ -604,7 +604,7 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg)
msg = 'read only the first document' msg = 'read only the first document'
case default case default
msg = 'unknown warning number' msg = 'unknown warning number'
end select end select
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(IO_STDERR,'(/,a)') ' ┌'//IO_DIVIDER//'┐' write(IO_STDERR,'(/,a)') ' ┌'//IO_DIVIDER//'┐'
@ -658,7 +658,7 @@ subroutine selfTest
if(any([1,1,1] /= IO_stringPos('a'))) error stop 'IO_stringPos' if(any([1,1,1] /= IO_stringPos('a'))) error stop 'IO_stringPos'
if(any([2,2,3,5,5] /= IO_stringPos(' aa b'))) error stop 'IO_stringPos' if(any([2,2,3,5,5] /= IO_stringPos(' aa b'))) error stop 'IO_stringPos'
str=' 1.0 xxx' str = ' 1.0 xxx'
chunkPos = IO_stringPos(str) chunkPos = IO_stringPos(str)
if(dNeq(1.0_pReal,IO_floatValue(str,chunkPos,1))) error stop 'IO_floatValue' if(dNeq(1.0_pReal,IO_floatValue(str,chunkPos,1))) error stop 'IO_floatValue'

View File

@ -266,7 +266,7 @@ subroutine selfTest
if(any(dNeq(l2%get_as1dFloat(1),[2.0_pReal,3.0_pReal]))) error stop 'byIndex_as1dFloat' if(any(dNeq(l2%get_as1dFloat(1),[2.0_pReal,3.0_pReal]))) error stop 'byIndex_as1dFloat'
call l2%append(l3) call l2%append(l3)
x = l2%as2dFloat() x = l2%as2dFloat()
if(x(2,1)/= 4.0_pReal) error stop 'byKey_as2dFloat' if(dNeq(x(2,1),4.0_pReal)) error stop 'byKey_as2dFloat'
if(any(dNeq(pack(l2%as2dFloat(),.true.),& if(any(dNeq(pack(l2%as2dFloat(),.true.),&
[2.0_pReal,4.0_pReal,3.0_pReal,5.0_pReal]))) error stop 'byKey_as2dFloat' [2.0_pReal,4.0_pReal,3.0_pReal,5.0_pReal]))) error stop 'byKey_as2dFloat'
n => l2 n => l2

View File

@ -47,5 +47,8 @@
#include "homogenization_mechanical_isostrain.f90" #include "homogenization_mechanical_isostrain.f90"
#include "homogenization_mechanical_RGC.f90" #include "homogenization_mechanical_RGC.f90"
#include "homogenization_thermal.f90" #include "homogenization_thermal.f90"
#include "homogenization_thermal_pass.f90"
#include "homogenization_thermal_isotemperature.f90"
#include "homogenization_damage.f90" #include "homogenization_damage.f90"
#include "homogenization_damage_pass.f90"
#include "CPFEM.f90" #include "CPFEM.f90"

View File

@ -196,7 +196,7 @@ function grid_damage_spectral_solution(timeinc) result(solution)
ce = 0 ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1 ce = ce + 1
call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),ce) call homogenization_set_phi(phi_current(i,j,k),ce)
enddo; enddo; enddo enddo; enddo; enddo
call VecMin(solution_vec,devNull,phi_min,ierr); CHKERRQ(ierr) call VecMin(solution_vec,devNull,phi_min,ierr); CHKERRQ(ierr)
@ -233,7 +233,7 @@ subroutine grid_damage_spectral_forward(cutBack)
call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr)
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1 ce = ce + 1
call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),ce) call homogenization_set_phi(phi_current(i,j,k),ce)
enddo; enddo; enddo enddo; enddo; enddo
else else
phi_lastInc = phi_current phi_lastInc = phi_current
@ -259,7 +259,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
PetscObject :: dummy PetscObject :: dummy
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer :: i, j, k, ce integer :: i, j, k, ce
real(pReal) :: phiDot, dPhiDot_dPhi, mobility real(pReal) :: phiDot, mobility
phi_current = x_scal phi_current = x_scal
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -281,7 +281,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
ce = 0 ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1 ce = ce + 1
call damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi_current(i,j,k),ce) call damage_nonlocal_getSourceAndItsTangent(phiDot, phi_current(i,j,k),ce)
mobility = damage_nonlocal_getMobility(ce) mobility = damage_nonlocal_getMobility(ce)
scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + phiDot) & scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + phiDot) &
+ mobility*(phi_lastInc(i,j,k) - phi_current(i,j,k)) & + mobility*(phi_lastInc(i,j,k) - phi_current(i,j,k)) &

View File

@ -282,9 +282,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
ce = ce + 1 ce = ce + 1
call thermal_conduction_getSource(Tdot,1,ce) call thermal_conduction_getSource(Tdot,1,ce)
scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + Tdot) & scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + Tdot) &
+ thermal_conduction_getMassDensity (ce)* & + homogenization_thermal_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) &
thermal_conduction_getSpecificHeat(ce)*(T_lastInc(i,j,k) - &
T_current(i,j,k))&
+ mu_ref*T_current(i,j,k) + mu_ref*T_current(i,j,k)
enddo; enddo; enddo enddo; enddo; enddo
@ -314,7 +312,7 @@ subroutine updateReference
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1 ce = ce + 1
K_ref = K_ref + thermal_conduction_getConductivity(ce) K_ref = K_ref + thermal_conduction_getConductivity(ce)
mu_ref = mu_ref + thermal_conduction_getMassDensity(ce)* thermal_conduction_getSpecificHeat(ce) mu_ref = mu_ref + homogenization_thermal_mu_T(ce)
enddo; enddo; enddo enddo; enddo; enddo
K_ref = K_ref*wgt K_ref = K_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)

View File

@ -39,8 +39,6 @@ module homogenization
thermal_type !< thermal transport model thermal_type !< thermal transport model
integer(kind(DAMAGE_none_ID)), dimension(:), allocatable :: & integer(kind(DAMAGE_none_ID)), dimension(:), allocatable :: &
damage_type !< nonlocal damage model damage_type !< nonlocal damage model
integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:), allocatable :: &
homogenization_type !< type of each homogenization
type, private :: tNumerics_damage type, private :: tNumerics_damage
real(pReal) :: & real(pReal) :: &
@ -117,6 +115,16 @@ module homogenization
integer, intent(in) :: ho integer, intent(in) :: ho
end subroutine mechanical_results end subroutine mechanical_results
module subroutine damage_results(ho,group)
integer, intent(in) :: ho
character(len=*), intent(in) :: group
end subroutine damage_results
module subroutine thermal_results(ho,group)
integer, intent(in) :: ho
character(len=*), intent(in) :: group
end subroutine thermal_results
module function mechanical_updateState(subdt,subF,ce) result(doneAndHappy) module function mechanical_updateState(subdt,subF,ce) result(doneAndHappy)
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
subdt !< current time step subdt !< current time step
@ -133,26 +141,16 @@ module homogenization
real(pReal), dimension(3,3) :: K real(pReal), dimension(3,3) :: K
end function thermal_conduction_getConductivity end function thermal_conduction_getConductivity
module function thermal_conduction_getSpecificHeat(ce) result(c_P) module function homogenization_thermal_mu_T(ce) result(mu_T)
integer, intent(in) :: ce integer, intent(in) :: ce
real(pReal) :: c_P real(pReal) :: mu_T
end function thermal_conduction_getSpecificHeat end function homogenization_thermal_mu_T
module function thermal_conduction_getMassDensity(ce) result(rho)
integer, intent(in) :: ce
real(pReal) :: rho
end function thermal_conduction_getMassDensity
module subroutine homogenization_thermal_setField(T,dot_T, ce) module subroutine homogenization_thermal_setField(T,dot_T, ce)
integer, intent(in) :: ce integer, intent(in) :: ce
real(pReal), intent(in) :: T, dot_T real(pReal), intent(in) :: T, dot_T
end subroutine homogenization_thermal_setField end subroutine homogenization_thermal_setField
module subroutine thermal_conduction_results(ho,group)
integer, intent(in) :: ho
character(len=*), intent(in) :: group
end subroutine thermal_conduction_results
module function homogenization_thermal_T(ce) result(T) module function homogenization_thermal_T(ce) result(T)
integer, intent(in) :: ce integer, intent(in) :: ce
real(pReal) :: T real(pReal) :: T
@ -170,37 +168,31 @@ module homogenization
real(pReal) :: M real(pReal) :: M
end function damage_nonlocal_getMobility end function damage_nonlocal_getMobility
module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ce) module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, phi, ce)
integer, intent(in) :: ce integer, intent(in) :: ce
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi phi
real(pReal) :: & real(pReal), intent(out) :: &
phiDot, dPhiDot_dPhi phiDot
end subroutine damage_nonlocal_getSourceAndItsTangent end subroutine damage_nonlocal_getSourceAndItsTangent
module subroutine damage_nonlocal_putNonLocalDamage(phi,ce) module subroutine homogenization_set_phi(phi,ce)
integer, intent(in) :: ce integer, intent(in) :: ce
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi phi
end subroutine damage_nonlocal_putNonLocalDamage end subroutine homogenization_set_phi
module subroutine damage_nonlocal_results(ho,group)
integer, intent(in) :: ho
character(len=*), intent(in) :: group
end subroutine damage_nonlocal_results
end interface end interface
public :: & public :: &
homogenization_init, & homogenization_init, &
materialpoint_stressAndItsTangent, & materialpoint_stressAndItsTangent, &
thermal_conduction_getSpecificHeat, & homogenization_thermal_mu_T, &
thermal_conduction_getConductivity, & thermal_conduction_getConductivity, &
thermal_conduction_getMassDensity, &
thermal_conduction_getSource, & thermal_conduction_getSource, &
damage_nonlocal_getMobility, & damage_nonlocal_getMobility, &
damage_nonlocal_getSourceAndItsTangent, & damage_nonlocal_getSourceAndItsTangent, &
damage_nonlocal_putNonLocalDamage, & homogenization_set_phi, &
homogenization_thermal_setfield, & homogenization_thermal_setfield, &
homogenization_thermal_T, & homogenization_thermal_T, &
homogenization_forward, & homogenization_forward, &
@ -211,7 +203,6 @@ module homogenization
DAMAGE_NONLOCAL_ID DAMAGE_NONLOCAL_ID
public :: & public :: &
damage_nonlocal_init, &
damage_nonlocal_getDiffusion damage_nonlocal_getDiffusion
contains contains
@ -242,8 +233,6 @@ subroutine homogenization_init()
call mechanical_init(num_homog) call mechanical_init(num_homog)
call thermal_init() call thermal_init()
call damage_init() call damage_init()
call damage_nonlocal_init()
end subroutine homogenization_init end subroutine homogenization_init
@ -259,25 +248,25 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
NiterationMPstate, & NiterationMPstate, &
ip, & !< integration point number ip, & !< integration point number
el, & !< element number el, & !< element number
myNgrains, co, ce, ho, me, ph myNgrains, co, ce, ho, en, ph
logical :: & logical :: &
converged converged
logical, dimension(2) :: & logical, dimension(2) :: &
doneAndHappy doneAndHappy
!$OMP PARALLEL !$OMP PARALLEL
!$OMP DO PRIVATE(ce,me,ho,myNgrains,NiterationMPstate,converged,doneAndHappy) !$OMP DO PRIVATE(ce,en,ho,myNgrains,NiterationMPstate,converged,doneAndHappy)
do el = FEsolving_execElem(1),FEsolving_execElem(2) do el = FEsolving_execElem(1),FEsolving_execElem(2)
ho = material_homogenizationAt(el) ho = material_homogenizationAt(el)
myNgrains = homogenization_Nconstituents(ho) myNgrains = homogenization_Nconstituents(ho)
do ip = FEsolving_execIP(1),FEsolving_execIP(2) do ip = FEsolving_execIP(1),FEsolving_execIP(2)
ce = (el-1)*discretization_nIPs + ip ce = (el-1)*discretization_nIPs + ip
me = material_homogenizationMemberAt2(ce) en = material_homogenizationEntry(ce)
call phase_restore(ce,.false.) ! wrong name (is more a forward function) call phase_restore(ce,.false.) ! wrong name (is more a forward function)
if(homogState(ho)%sizeState > 0) homogState(ho)%state(:,me) = homogState(ho)%state0(:,me) if(homogState(ho)%sizeState > 0) homogState(ho)%state(:,en) = homogState(ho)%state0(:,en)
if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%state(:,me) = damageState_h(ho)%state0(:,me) if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%state(:,en) = damageState_h(ho)%state0(:,en)
call damage_partition(ce) call damage_partition(ce)
doneAndHappy = [.false.,.true.] doneAndHappy = [.false.,.true.]
@ -371,14 +360,14 @@ subroutine homogenization_results
case(DAMAGE_NONLOCAL_ID) case(DAMAGE_NONLOCAL_ID)
group = trim(group_base)//'/damage' group = trim(group_base)//'/damage'
call results_closeGroup(results_addGroup(group)) call results_closeGroup(results_addGroup(group))
call damage_nonlocal_results(ho,group) call damage_results(ho,group)
end select end select
select case(thermal_type(ho)) select case(thermal_type(ho))
case(THERMAL_CONDUCTION_ID) case(THERMAL_CONDUCTION_ID)
group = trim(group_base)//'/thermal' group = trim(group_base)//'/thermal'
call results_closeGroup(results_addGroup(group)) call results_closeGroup(results_addGroup(group))
call thermal_conduction_results(ho,group) call thermal_results(ho,group)
end select end select
enddo enddo
@ -458,41 +447,6 @@ subroutine homogenization_restartRead(fileHandle)
end subroutine homogenization_restartRead end subroutine homogenization_restartRead
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine damage_nonlocal_init
integer :: Ninstances,Nmaterialpoints,h
class(tNode), pointer :: &
num_generic, &
material_homogenization
print'(/,a)', ' <<<+- damage_nonlocal init -+>>>'; flush(6)
!------------------------------------------------------------------------------------
! read numerics parameter
num_generic => config_numerics%get('generic',defaultVal= emptyDict)
num_damage%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal)
Ninstances = count(damage_type == DAMAGE_nonlocal_ID)
material_homogenization => config_material%get('homogenization')
do h = 1, material_homogenization%length
if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle
Nmaterialpoints = count(material_homogenizationAt == h)
damageState_h(h)%sizeState = 1
allocate(damageState_h(h)%state0 (1,Nmaterialpoints), source=1.0_pReal)
allocate(damageState_h(h)%state (1,Nmaterialpoints), source=1.0_pReal)
enddo
end subroutine damage_nonlocal_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns homogenized non local damage diffusion tensor in reference configuration !> @brief returns homogenized non local damage diffusion tensor in reference configuration
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -505,12 +459,12 @@ function damage_nonlocal_getDiffusion(ce)
ho, & ho, &
co co
ho = material_homogenizationAt2(ce) ho = material_homogenizationID(ce)
damage_nonlocal_getDiffusion = 0.0_pReal damage_nonlocal_getDiffusion = 0.0_pReal
do co = 1, homogenization_Nconstituents(ho) do co = 1, homogenization_Nconstituents(ho)
damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + & damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + &
crystallite_push33ToRef(co,ce,lattice_D(1:3,1:3,material_phaseAt2(co,ce))) crystallite_push33ToRef(co,ce,lattice_D(1:3,1:3,material_phaseID(co,ce)))
enddo enddo
damage_nonlocal_getDiffusion = & damage_nonlocal_getDiffusion = &
@ -521,14 +475,12 @@ end function damage_nonlocal_getDiffusion
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief parses the homogenization part from the material configuration !> @brief parses the homogenization part from the material configuration
! ToDo: This should be done in homogenization
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine material_parseHomogenization subroutine material_parseHomogenization
class(tNode), pointer :: & class(tNode), pointer :: &
material_homogenization, & material_homogenization, &
homog, & homog, &
homogMech, &
homogThermal, & homogThermal, &
homogDamage homogDamage
@ -536,23 +488,11 @@ subroutine material_parseHomogenization
material_homogenization => config_material%get('homogenization') material_homogenization => config_material%get('homogenization')
allocate(homogenization_type(size(material_name_homogenization)), source=HOMOGENIZATION_undefined_ID)
allocate(thermal_type(size(material_name_homogenization)), source=THERMAL_isothermal_ID) allocate(thermal_type(size(material_name_homogenization)), source=THERMAL_isothermal_ID)
allocate(damage_type (size(material_name_homogenization)), source=DAMAGE_none_ID) allocate(damage_type (size(material_name_homogenization)), source=DAMAGE_none_ID)
do h=1, size(material_name_homogenization) do h=1, size(material_name_homogenization)
homog => material_homogenization%get(h) homog => material_homogenization%get(h)
homogMech => homog%get('mechanical')
select case (homogMech%get_asString('type'))
case('pass')
homogenization_type(h) = HOMOGENIZATION_NONE_ID
case('isostrain')
homogenization_type(h) = HOMOGENIZATION_ISOSTRAIN_ID
case('RGC')
homogenization_type(h) = HOMOGENIZATION_RGC_ID
case default
call IO_error(500,ext_msg=homogMech%get_asString('type'))
end select
if (homog%contains('thermal')) then if (homog%contains('thermal')) then
homogThermal => homog%get('thermal') homogThermal => homog%get('thermal')

View File

@ -1,10 +1,17 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, KU Leuven !> @author Martin Diehl, KU Leuven
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule(homogenization) homogenization_damage submodule(homogenization) damage
use lattice use lattice
interface
module subroutine pass_init
end subroutine pass_init
end interface
type :: tDataContainer type :: tDataContainer
real(pReal), dimension(:), allocatable :: phi real(pReal), dimension(:), allocatable :: phi
end type tDataContainer end type tDataContainer
@ -30,19 +37,22 @@ module subroutine damage_init()
class(tNode), pointer :: & class(tNode), pointer :: &
configHomogenizations, & configHomogenizations, &
configHomogenization, & configHomogenization, &
configHomogenizationDamage configHomogenizationDamage, &
num_generic, &
material_homogenization
integer :: ho integer :: ho
integer :: Ninstances,Nmaterialpoints,h
print'(/,a)', ' <<<+- homogenization:damage init -+>>>' print'(/,a)', ' <<<+- homogenization:damage init -+>>>'
print'(/,a)', ' <<<+- homogenization:damage:isodamage init -+>>>' print'(/,a)', ' <<<+- homogenization:damage:pass init -+>>>'
configHomogenizations => config_material%get('homogenization') configHomogenizations => config_material%get('homogenization')
allocate(param(configHomogenizations%length)) allocate(param(configHomogenizations%length))
allocate(current(configHomogenizations%length)) allocate(current(configHomogenizations%length))
do ho = 1, configHomogenizations%length do ho = 1, configHomogenizations%length
allocate(current(ho)%phi(count(material_homogenizationAt2==ho)), source=1.0_pReal) allocate(current(ho)%phi(count(material_homogenizationID==ho)), source=1.0_pReal)
configHomogenization => configHomogenizations%get(ho) configHomogenization => configHomogenizations%get(ho)
associate(prm => param(ho)) associate(prm => param(ho))
if (configHomogenization%contains('damage')) then if (configHomogenization%contains('damage')) then
@ -58,6 +68,24 @@ module subroutine damage_init()
end associate end associate
enddo enddo
!------------------------------------------------------------------------------------
! read numerics parameter
num_generic => config_numerics%get('generic',defaultVal= emptyDict)
num_damage%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal)
Ninstances = count(damage_type == DAMAGE_nonlocal_ID)
material_homogenization => config_material%get('homogenization')
do h = 1, material_homogenization%length
if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle
Nmaterialpoints = count(material_homogenizationAt == h)
damageState_h(h)%sizeState = 1
allocate(damageState_h(h)%state0 (1,Nmaterialpoints), source=1.0_pReal)
allocate(damageState_h(h)%state (1,Nmaterialpoints), source=1.0_pReal)
enddo
end subroutine damage_init end subroutine damage_init
@ -69,14 +97,10 @@ module subroutine damage_partition(ce)
real(pReal) :: phi real(pReal) :: phi
integer, intent(in) :: ce integer, intent(in) :: ce
integer :: co
if(damageState_h(material_homogenizationID(ce))%sizeState < 1) return
if(damageState_h(material_homogenizationAt2(ce))%sizeState < 1) return phi = damagestate_h(material_homogenizationID(ce))%state(1,material_homogenizationEntry(ce))
phi = damagestate_h(material_homogenizationAt2(ce))%state(1,material_homogenizationMemberAt2(ce)) call phase_damage_set_phi(phi,1,ce)
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce))
call phase_damage_set_phi(phi,co,ce)
enddo
end subroutine damage_partition end subroutine damage_partition
@ -88,17 +112,9 @@ end subroutine damage_partition
module function damage_nonlocal_getMobility(ce) result(M) module function damage_nonlocal_getMobility(ce) result(M)
integer, intent(in) :: ce integer, intent(in) :: ce
integer :: &
co
real(pReal) :: M real(pReal) :: M
M = 0.0_pReal M = lattice_M(material_phaseID(1,ce))
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce))
M = M + lattice_M(material_phaseAt2(co,ce))
enddo
M = M/real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal)
end function damage_nonlocal_getMobility end function damage_nonlocal_getMobility
@ -106,20 +122,15 @@ end function damage_nonlocal_getMobility
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates homogenized damage driving forces !> @brief calculates homogenized damage driving forces
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ce) module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, phi, ce)
integer, intent(in) :: ce integer, intent(in) :: ce
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi phi
real(pReal) :: & real(pReal), intent(out) :: &
phiDot, dPhiDot_dPhi phiDot
phiDot = 0.0_pReal phiDot = phase_damage_phi_dot(phi, 1, ce)
dPhiDot_dPhi = 0.0_pReal
call phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce)
phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal)
dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal)
end subroutine damage_nonlocal_getSourceAndItsTangent end subroutine damage_nonlocal_getSourceAndItsTangent
@ -127,26 +138,27 @@ end subroutine damage_nonlocal_getSourceAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief updated nonlocal damage field with solution from damage phase field PDE !> @brief updated nonlocal damage field with solution from damage phase field PDE
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine damage_nonlocal_putNonLocalDamage(phi,ce) module subroutine homogenization_set_phi(phi,ce)
integer, intent(in) :: ce integer, intent(in) :: ce
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi phi
integer :: & integer :: &
ho, & ho, &
me en
ho = material_homogenizationAt2(ce) ho = material_homogenizationID(ce)
me = material_homogenizationMemberAt2(ce) en = material_homogenizationEntry(ce)
damagestate_h(ho)%state(1,me) = phi damagestate_h(ho)%state(1,en) = phi
current(ho)%phi(en) = phi
end subroutine damage_nonlocal_putNonLocalDamage end subroutine homogenization_set_phi
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine damage_nonlocal_results(ho,group) module subroutine damage_results(ho,group)
integer, intent(in) :: ho integer, intent(in) :: ho
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
@ -163,6 +175,6 @@ module subroutine damage_nonlocal_results(ho,group)
enddo outputsLoop enddo outputsLoop
end associate end associate
end subroutine damage_nonlocal_results end subroutine damage_results
end submodule homogenization_damage end submodule damage

View File

@ -0,0 +1,14 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, KU Leuven
!> @brief Dummy homogenization scheme for 1 constituent per material point
!--------------------------------------------------------------------------------------------------
submodule(homogenization:damage) damage_pass
contains
module subroutine pass_init
end subroutine pass_init
end submodule damage_pass

View File

@ -7,51 +7,51 @@ submodule(homogenization) mechanical
interface interface
module subroutine mechanical_pass_init module subroutine pass_init
end subroutine mechanical_pass_init end subroutine pass_init
module subroutine mechanical_isostrain_init module subroutine isostrain_init
end subroutine mechanical_isostrain_init end subroutine isostrain_init
module subroutine mechanical_RGC_init(num_homogMech) module subroutine RGC_init(num_homogMech)
class(tNode), pointer, intent(in) :: & class(tNode), pointer, intent(in) :: &
num_homogMech !< pointer to mechanical homogenization numerics data num_homogMech !< pointer to mechanical homogenization numerics data
end subroutine mechanical_RGC_init end subroutine RGC_init
module subroutine mechanical_isostrain_partitionDeformation(F,avgF) module subroutine isostrain_partitionDeformation(F,avgF)
real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient
real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point
end subroutine mechanical_isostrain_partitionDeformation end subroutine isostrain_partitionDeformation
module subroutine mechanical_RGC_partitionDeformation(F,avgF,ce) module subroutine RGC_partitionDeformation(F,avgF,ce)
real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient
real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point
integer, intent(in) :: & integer, intent(in) :: &
ce ce
end subroutine mechanical_RGC_partitionDeformation end subroutine RGC_partitionDeformation
module subroutine mechanical_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho) module subroutine isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho)
real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point
real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses
real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
integer, intent(in) :: ho integer, intent(in) :: ho
end subroutine mechanical_isostrain_averageStressAndItsTangent end subroutine isostrain_averageStressAndItsTangent
module subroutine mechanical_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho) module subroutine RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho)
real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point
real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses
real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
integer, intent(in) :: ho integer, intent(in) :: ho
end subroutine mechanical_RGC_averageStressAndItsTangent end subroutine RGC_averageStressAndItsTangent
module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
logical, dimension(2) :: doneAndHappy logical, dimension(2) :: doneAndHappy
real(pReal), dimension(:,:,:), intent(in) :: & real(pReal), dimension(:,:,:), intent(in) :: &
P,& !< partitioned stresses P,& !< partitioned stresses
@ -61,16 +61,19 @@ submodule(homogenization) mechanical
real(pReal), intent(in) :: dt !< time increment real(pReal), intent(in) :: dt !< time increment
integer, intent(in) :: & integer, intent(in) :: &
ce !< cell ce !< cell
end function mechanical_RGC_updateState end function RGC_updateState
module subroutine mechanical_RGC_results(ho,group) module subroutine RGC_results(ho,group)
integer, intent(in) :: ho !< homogenization type integer, intent(in) :: ho !< homogenization type
character(len=*), intent(in) :: group !< group name in HDF5 file character(len=*), intent(in) :: group !< group name in HDF5 file
end subroutine mechanical_RGC_results end subroutine RGC_results
end interface end interface
integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:), allocatable :: &
homogenization_type !< type of each homogenization
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -86,15 +89,17 @@ module subroutine mechanical_init(num_homog)
print'(/,a)', ' <<<+- homogenization:mechanical init -+>>>' print'(/,a)', ' <<<+- homogenization:mechanical init -+>>>'
call material_parseHomogenization2()
allocate(homogenization_dPdF(3,3,3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal) allocate(homogenization_dPdF(3,3,3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal)
homogenization_F0 = spread(math_I3,3,discretization_nIPs*discretization_Nelems) ! initialize to identity homogenization_F0 = spread(math_I3,3,discretization_nIPs*discretization_Nelems) ! initialize to identity
homogenization_F = homogenization_F0 ! initialize to identity homogenization_F = homogenization_F0 ! initialize to identity
allocate(homogenization_P(3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal) allocate(homogenization_P(3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal)
num_homogMech => num_homog%get('mech',defaultVal=emptyDict) num_homogMech => num_homog%get('mech',defaultVal=emptyDict)
if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call mechanical_pass_init if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call pass_init
if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call mechanical_isostrain_init if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call isostrain_init
if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call mechanical_RGC_init(num_homogMech) if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call RGC_init(num_homogMech)
end subroutine mechanical_init end subroutine mechanical_init
@ -110,24 +115,24 @@ module subroutine mechanical_partition(subF,ce)
ce ce
integer :: co integer :: co
real(pReal), dimension (3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) :: Fs real(pReal), dimension (3,3,homogenization_Nconstituents(material_homogenizationID(ce))) :: Fs
chosenHomogenization: select case(homogenization_type(material_homogenizationAt2(ce))) chosenHomogenization: select case(homogenization_type(material_homogenizationID(ce)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization case (HOMOGENIZATION_NONE_ID) chosenHomogenization
Fs(1:3,1:3,1) = subF Fs(1:3,1:3,1) = subF
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
call mechanical_isostrain_partitionDeformation(Fs,subF) call isostrain_partitionDeformation(Fs,subF)
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
call mechanical_RGC_partitionDeformation(Fs,subF,ce) call RGC_partitionDeformation(Fs,subF,ce)
end select chosenHomogenization end select chosenHomogenization
do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
call phase_mechanical_setF(Fs(1:3,1:3,co),co,ce) call phase_set_F(Fs(1:3,1:3,co),co,ce)
enddo enddo
@ -143,37 +148,37 @@ module subroutine mechanical_homogenize(dt,ce)
integer, intent(in) :: ce integer, intent(in) :: ce
integer :: co integer :: co
real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationID(ce)))
real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationID(ce)))
chosenHomogenization: select case(homogenization_type(material_homogenizationAt2(ce))) chosenHomogenization: select case(homogenization_type(material_homogenizationID(ce)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization case (HOMOGENIZATION_NONE_ID) chosenHomogenization
homogenization_P(1:3,1:3,ce) = phase_mechanical_getP(1,ce) homogenization_P(1:3,1:3,ce) = phase_P(1,ce)
homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(dt,1,ce) homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(dt,1,ce)
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ce) dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ce)
Ps(:,:,co) = phase_mechanical_getP(co,ce) Ps(:,:,co) = phase_P(co,ce)
enddo enddo
call mechanical_isostrain_averageStressAndItsTangent(& call isostrain_averageStressAndItsTangent(&
homogenization_P(1:3,1:3,ce), & homogenization_P(1:3,1:3,ce), &
homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& homogenization_dPdF(1:3,1:3,1:3,1:3,ce),&
Ps,dPdFs, & Ps,dPdFs, &
material_homogenizationAt2(ce)) material_homogenizationID(ce))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ce) dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ce)
Ps(:,:,co) = phase_mechanical_getP(co,ce) Ps(:,:,co) = phase_P(co,ce)
enddo enddo
call mechanical_RGC_averageStressAndItsTangent(& call RGC_averageStressAndItsTangent(&
homogenization_P(1:3,1:3,ce), & homogenization_P(1:3,1:3,ce), &
homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& homogenization_dPdF(1:3,1:3,1:3,1:3,ce),&
Ps,dPdFs, & Ps,dPdFs, &
material_homogenizationAt2(ce)) material_homogenizationID(ce))
end select chosenHomogenization end select chosenHomogenization
@ -195,18 +200,18 @@ module function mechanical_updateState(subdt,subF,ce) result(doneAndHappy)
logical, dimension(2) :: doneAndHappy logical, dimension(2) :: doneAndHappy
integer :: co integer :: co
real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationID(ce)))
real(pReal) :: Fs(3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) real(pReal) :: Fs(3,3,homogenization_Nconstituents(material_homogenizationID(ce)))
real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationID(ce)))
if (homogenization_type(material_homogenizationAt2(ce)) == HOMOGENIZATION_RGC_ID) then if (homogenization_type(material_homogenizationID(ce)) == HOMOGENIZATION_RGC_ID) then
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(subdt,co,ce) dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(subdt,co,ce)
Fs(:,:,co) = phase_mechanical_getF(co,ce) Fs(:,:,co) = phase_F(co,ce)
Ps(:,:,co) = phase_mechanical_getP(co,ce) Ps(:,:,co) = phase_P(co,ce)
enddo enddo
doneAndHappy = mechanical_RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ce) doneAndHappy = RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ce)
else else
doneAndHappy = .true. doneAndHappy = .true.
endif endif
@ -230,7 +235,7 @@ module subroutine mechanical_results(group_base,ho)
select case(homogenization_type(ho)) select case(homogenization_type(ho))
case(HOMOGENIZATION_rgc_ID) case(HOMOGENIZATION_rgc_ID)
call mechanical_RGC_results(ho,group) call RGC_results(ho,group)
end select end select
@ -244,4 +249,38 @@ module subroutine mechanical_results(group_base,ho)
end subroutine mechanical_results end subroutine mechanical_results
!--------------------------------------------------------------------------------------------------
!> @brief parses the homogenization part from the material configuration
!--------------------------------------------------------------------------------------------------
subroutine material_parseHomogenization2()
class(tNode), pointer :: &
material_homogenization, &
homog, &
homogMech
integer :: h
material_homogenization => config_material%get('homogenization')
allocate(homogenization_type(size(material_name_homogenization)), source=HOMOGENIZATION_undefined_ID)
do h=1, size(material_name_homogenization)
homog => material_homogenization%get(h)
homogMech => homog%get('mechanical')
select case (homogMech%get_asString('type'))
case('pass')
homogenization_type(h) = HOMOGENIZATION_NONE_ID
case('isostrain')
homogenization_type(h) = HOMOGENIZATION_ISOSTRAIN_ID
case('RGC')
homogenization_type(h) = HOMOGENIZATION_RGC_ID
case default
call IO_error(500,ext_msg=homogMech%get_asString('type'))
end select
enddo
end subroutine material_parseHomogenization2
end submodule mechanical end submodule mechanical

View File

@ -71,7 +71,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields, reads information from material configuration file !> @brief allocates all necessary fields, reads information from material configuration file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mechanical_RGC_init(num_homogMech) module subroutine RGC_init(num_homogMech)
class(tNode), pointer, intent(in) :: & class(tNode), pointer, intent(in) :: &
num_homogMech !< pointer to mechanical homogenization numerics data num_homogMech !< pointer to mechanical homogenization numerics data
@ -152,7 +152,7 @@ module subroutine mechanical_RGC_init(num_homogMech)
prm%N_constituents = homogMech%get_as1dInt('cluster_size',requiredSize=3) prm%N_constituents = homogMech%get_as1dInt('cluster_size',requiredSize=3)
if (homogenization_Nconstituents(ho) /= product(prm%N_constituents)) & if (homogenization_Nconstituents(ho) /= product(prm%N_constituents)) &
call IO_error(211,ext_msg='N_constituents (mechanical_RGC)') call IO_error(211,ext_msg='N_constituents (RGC)')
prm%xi_alpha = homogMech%get_asFloat('xi_alpha') prm%xi_alpha = homogMech%get_asFloat('xi_alpha')
prm%c_alpha = homogMech%get_asFloat('c_alpha') prm%c_alpha = homogMech%get_asFloat('c_alpha')
@ -187,13 +187,13 @@ module subroutine mechanical_RGC_init(num_homogMech)
enddo enddo
end subroutine mechanical_RGC_init end subroutine RGC_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief partitions the deformation gradient onto the constituents !> @brief partitions the deformation gradient onto the constituents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mechanical_RGC_partitionDeformation(F,avgF,ce) module subroutine RGC_partitionDeformation(F,avgF,ce)
real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned F per grain real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned F per grain
@ -204,12 +204,12 @@ module subroutine mechanical_RGC_partitionDeformation(F,avgF,ce)
real(pReal), dimension(3) :: aVect,nVect real(pReal), dimension(3) :: aVect,nVect
integer, dimension(4) :: intFace integer, dimension(4) :: intFace
integer, dimension(3) :: iGrain3 integer, dimension(3) :: iGrain3
integer :: iGrain,iFace,i,j,ho,me integer :: iGrain,iFace,i,j,ho,en
associate(prm => param(material_homogenizationAt2(ce))) associate(prm => param(material_homogenizationID(ce)))
ho = material_homogenizationAt2(ce) ho = material_homogenizationID(ce)
me = material_homogenizationMemberAt2(ce) en = material_homogenizationEntry(ce)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! compute the deformation gradient of individual grains due to relaxations ! compute the deformation gradient of individual grains due to relaxations
F = 0.0_pReal F = 0.0_pReal
@ -217,8 +217,8 @@ module subroutine mechanical_RGC_partitionDeformation(F,avgF,ce)
iGrain3 = grain1to3(iGrain,prm%N_constituents) iGrain3 = grain1to3(iGrain,prm%N_constituents)
do iFace = 1,6 do iFace = 1,6
intFace = getInterface(iFace,iGrain3) ! identifying 6 interfaces of each grain intFace = getInterface(iFace,iGrain3) ! identifying 6 interfaces of each grain
aVect = relaxationVector(intFace,ho,me) ! get the relaxation vectors for each interface from global relaxation vector array aVect = relaxationVector(intFace,ho,en) ! get the relaxation vectors for each interface from global relaxation vector array
nVect = interfaceNormal(intFace,ho,me) nVect = interfaceNormal(intFace,ho,en)
forall (i=1:3,j=1:3) & forall (i=1:3,j=1:3) &
F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! calculating deformation relaxations due to interface relaxation F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! calculating deformation relaxations due to interface relaxation
enddo enddo
@ -227,14 +227,14 @@ module subroutine mechanical_RGC_partitionDeformation(F,avgF,ce)
end associate end associate
end subroutine mechanical_RGC_partitionDeformation end subroutine RGC_partitionDeformation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief update the internal state of the homogenization scheme and tell whether "done" and !> @brief update the internal state of the homogenization scheme and tell whether "done" and
! "happy" with result ! "happy" with result
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
logical, dimension(2) :: doneAndHappy logical, dimension(2) :: doneAndHappy
real(pReal), dimension(:,:,:), intent(in) :: & real(pReal), dimension(:,:,:), intent(in) :: &
P,& !< partitioned stresses P,& !< partitioned stresses
@ -247,7 +247,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
integer, dimension(4) :: intFaceN,intFaceP,faceID integer, dimension(4) :: intFaceN,intFaceP,faceID
integer, dimension(3) :: nGDim,iGr3N,iGr3P integer, dimension(3) :: nGDim,iGr3N,iGr3P
integer :: ho,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ipert,iGrain,nGrain, me integer :: ho,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ipert,nGrain, en
real(pReal), dimension(3,3,size(P,3)) :: R,pF,pR,D,pD real(pReal), dimension(3,3,size(P,3)) :: R,pF,pR,D,pD
real(pReal), dimension(3,size(P,3)) :: NN,devNull real(pReal), dimension(3,size(P,3)) :: NN,devNull
real(pReal), dimension(3) :: normP,normN,mornP,mornN real(pReal), dimension(3) :: normP,normN,mornP,mornN
@ -261,9 +261,9 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
return return
endif zeroTimeStep endif zeroTimeStep
ho = material_homogenizationAt2(ce) ho = material_homogenizationID(ce)
en = material_homogenizationEntry(ce)
me = material_homogenizationMemberAt2(ce)
associate(stt => state(ho), st0 => state0(ho), dst => dependentState(ho), prm => param(ho)) associate(stt => state(ho), st0 => state0(ho), dst => dependentState(ho), prm => param(ho))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -278,16 +278,16 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
! allocate the size of the global relaxation arrays/jacobian matrices depending on the size of the cluster ! allocate the size of the global relaxation arrays/jacobian matrices depending on the size of the cluster
allocate(resid(3*nIntFaceTot), source=0.0_pReal) allocate(resid(3*nIntFaceTot), source=0.0_pReal)
allocate(tract(nIntFaceTot,3), source=0.0_pReal) allocate(tract(nIntFaceTot,3), source=0.0_pReal)
relax = stt%relaxationVector(:,me) relax = stt%relaxationVector(:,en)
drelax = stt%relaxationVector(:,me) - st0%relaxationVector(:,me) drelax = stt%relaxationVector(:,en) - st0%relaxationVector(:,en)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! computing interface mismatch and stress penalty tensor for all interfaces of all grains ! computing interface mismatch and stress penalty tensor for all interfaces of all grains
call stressPenalty(R,NN,avgF,F,ho,me) call stressPenalty(R,NN,avgF,F,ho,en)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculating volume discrepancy and stress penalty related to overall volume discrepancy ! calculating volume discrepancy and stress penalty related to overall volume discrepancy
call volumePenalty(D,dst%volumeDiscrepancy(me),avgF,F,nGrain) call volumePenalty(D,dst%volumeDiscrepancy(en),avgF,F,nGrain)
!------------------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------------------
! computing the residual stress from the balance of traction at all (interior) interfaces ! computing the residual stress from the balance of traction at all (interior) interfaces
@ -299,7 +299,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate system (3-dimensional index) iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate system (3-dimensional index)
iGrN = grain3to1(iGr3N,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) iGrN = grain3to1(iGr3N,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceN = getInterface(2*faceID(1),iGr3N) intFaceN = getInterface(2*faceID(1),iGr3N)
normN = interfaceNormal(intFaceN,ho,me) normN = interfaceNormal(intFaceN,ho,en)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! identify the right/up/front grain (+|P) ! identify the right/up/front grain (+|P)
@ -307,7 +307,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate system (3-dimensional index) iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate system (3-dimensional index)
iGrP = grain3to1(iGr3P,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) iGrP = grain3to1(iGr3P,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceP = getInterface(2*faceID(1)-1,iGr3P) intFaceP = getInterface(2*faceID(1)-1,iGr3P)
normP = interfaceNormal(intFaceP,ho,me) normP = interfaceNormal(intFaceP,ho,en)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! compute the residual of traction at the interface (in local system, 4-dimensional index) ! compute the residual of traction at the interface (in local system, 4-dimensional index)
@ -335,9 +335,9 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
if (residMax < num%rtol*stresMax .or. residMax < num%atol) then if (residMax < num%rtol*stresMax .or. residMax < num%atol) then
doneAndHappy = .true. doneAndHappy = .true.
dst%mismatch(1:3,me) = sum(NN,2)/real(nGrain,pReal) dst%mismatch(1:3,en) = sum(NN,2)/real(nGrain,pReal)
dst%relaxationRate_avg(me) = sum(abs(drelax))/dt/real(3*nIntFaceTot,pReal) dst%relaxationRate_avg(en) = sum(abs(drelax))/dt/real(3*nIntFaceTot,pReal)
dst%relaxationRate_max(me) = maxval(abs(drelax))/dt dst%relaxationRate_max(en) = maxval(abs(drelax))/dt
return return
@ -363,10 +363,10 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate sytem iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate sytem
iGrN = grain3to1(iGr3N,param(ho)%N_constituents) ! translate into global grain ID iGrN = grain3to1(iGr3N,param(ho)%N_constituents) ! translate into global grain ID
intFaceN = getInterface(2*faceID(1),iGr3N) ! identifying the connecting interface in local coordinate system intFaceN = getInterface(2*faceID(1),iGr3N) ! identifying the connecting interface in local coordinate system
normN = interfaceNormal(intFaceN,ho,me) normN = interfaceNormal(intFaceN,ho,en)
do iFace = 1,6 do iFace = 1,6
intFaceN = getInterface(iFace,iGr3N) ! identifying all interfaces that influence relaxation of the above interface intFaceN = getInterface(iFace,iGr3N) ! identifying all interfaces that influence relaxation of the above interface
mornN = interfaceNormal(intFaceN,ho,me) mornN = interfaceNormal(intFaceN,ho,en)
iMun = interface4to1(intFaceN,param(ho)%N_constituents) ! translate the interfaces ID into local 4-dimensional index iMun = interface4to1(intFaceN,param(ho)%N_constituents) ! translate the interfaces ID into local 4-dimensional index
if (iMun > 0) then ! get the corresponding tangent if (iMun > 0) then ! get the corresponding tangent
do i=1,3; do j=1,3; do k=1,3; do l=1,3 do i=1,3; do j=1,3; do k=1,3; do l=1,3
@ -384,10 +384,10 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate sytem iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate sytem
iGrP = grain3to1(iGr3P,param(ho)%N_constituents) ! translate into global grain ID iGrP = grain3to1(iGr3P,param(ho)%N_constituents) ! translate into global grain ID
intFaceP = getInterface(2*faceID(1)-1,iGr3P) ! identifying the connecting interface in local coordinate system intFaceP = getInterface(2*faceID(1)-1,iGr3P) ! identifying the connecting interface in local coordinate system
normP = interfaceNormal(intFaceP,ho,me) normP = interfaceNormal(intFaceP,ho,en)
do iFace = 1,6 do iFace = 1,6
intFaceP = getInterface(iFace,iGr3P) ! identifying all interfaces that influence relaxation of the above interface intFaceP = getInterface(iFace,iGr3P) ! identifying all interfaces that influence relaxation of the above interface
mornP = interfaceNormal(intFaceP,ho,me) mornP = interfaceNormal(intFaceP,ho,en)
iMun = interface4to1(intFaceP,param(ho)%N_constituents) ! translate the interfaces ID into local 4-dimensional index iMun = interface4to1(intFaceP,param(ho)%N_constituents) ! translate the interfaces ID into local 4-dimensional index
if (iMun > 0) then ! get the corresponding tangent if (iMun > 0) then ! get the corresponding tangent
do i=1,3; do j=1,3; do k=1,3; do l=1,3 do i=1,3; do j=1,3; do k=1,3; do l=1,3
@ -408,9 +408,9 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
do ipert = 1,3*nIntFaceTot do ipert = 1,3*nIntFaceTot
p_relax = relax p_relax = relax
p_relax(ipert) = relax(ipert) + num%pPert ! perturb the relaxation vector p_relax(ipert) = relax(ipert) + num%pPert ! perturb the relaxation vector
stt%relaxationVector(:,me) = p_relax stt%relaxationVector(:,en) = p_relax
call grainDeformation(pF,avgF,ho,me) ! rain deformation from perturbed state call grainDeformation(pF,avgF,ho,en) ! rain deformation from perturbed state
call stressPenalty(pR,DevNull, avgF,pF,ho,me) ! stress penalty due to interface mismatch from perturbed state call stressPenalty(pR,DevNull, avgF,pF,ho,en) ! stress penalty due to interface mismatch from perturbed state
call volumePenalty(pD,devNull(1,1), avgF,pF,nGrain) ! stress penalty due to volume discrepancy from perturbed state call volumePenalty(pD,devNull(1,1), avgF,pF,nGrain) ! stress penalty due to volume discrepancy from perturbed state
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -424,7 +424,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
iGr3N = faceID(2:4) ! identify the grain ID in local coordinate system (3-dimensional index) iGr3N = faceID(2:4) ! identify the grain ID in local coordinate system (3-dimensional index)
iGrN = grain3to1(iGr3N,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) iGrN = grain3to1(iGr3N,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceN = getInterface(2*faceID(1),iGr3N) ! identify the interface ID of the grain intFaceN = getInterface(2*faceID(1),iGr3N) ! identify the interface ID of the grain
normN = interfaceNormal(intFaceN,ho,me) normN = interfaceNormal(intFaceN,ho,en)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! identify the right/up/front grain (+|P) ! identify the right/up/front grain (+|P)
@ -432,7 +432,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identify the grain ID in local coordinate system (3-dimensional index) iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identify the grain ID in local coordinate system (3-dimensional index)
iGrP = grain3to1(iGr3P,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) iGrP = grain3to1(iGr3P,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceP = getInterface(2*faceID(1)-1,iGr3P) ! identify the interface ID of the grain intFaceP = getInterface(2*faceID(1)-1,iGr3P) ! identify the interface ID of the grain
normP = interfaceNormal(intFaceP,ho,me) normP = interfaceNormal(intFaceP,ho,en)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! compute the residual stress (contribution of mismatch and volume penalties) from perturbed state ! compute the residual stress (contribution of mismatch and volume penalties) from perturbed state
@ -472,7 +472,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
do i = 1,3*nIntFaceTot;do j = 1,3*nIntFaceTot do i = 1,3*nIntFaceTot;do j = 1,3*nIntFaceTot
drelax(i) = drelax(i) - jnverse(i,j)*resid(j) ! Calculate the correction for the state variable drelax(i) = drelax(i) - jnverse(i,j)*resid(j) ! Calculate the correction for the state variable
enddo; enddo enddo; enddo
stt%relaxationVector(:,me) = relax + drelax ! Updateing the state variable for the next iteration stt%relaxationVector(:,en) = relax + drelax ! Updateing the state variable for the next iteration
if (any(abs(drelax) > num%maxdRelax)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large if (any(abs(drelax) > num%maxdRelax)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large
doneAndHappy = [.true.,.false.] doneAndHappy = [.true.,.false.]
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
@ -488,14 +488,14 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
!------------------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------------------
!> @brief calculate stress-like penalty due to deformation mismatch !> @brief calculate stress-like penalty due to deformation mismatch
!------------------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------------------
subroutine stressPenalty(rPen,nMis,avgF,fDef,ho,me) subroutine stressPenalty(rPen,nMis,avgF,fDef,ho,en)
real(pReal), dimension (:,:,:), intent(out) :: rPen !< stress-like penalty real(pReal), dimension (:,:,:), intent(out) :: rPen !< stress-like penalty
real(pReal), dimension (:,:), intent(out) :: nMis !< total amount of mismatch real(pReal), dimension (:,:), intent(out) :: nMis !< total amount of mismatch
real(pReal), dimension (:,:,:), intent(in) :: fDef !< deformation gradients real(pReal), dimension (:,:,:), intent(in) :: fDef !< deformation gradients
real(pReal), dimension (3,3), intent(in) :: avgF !< initial effective stretch tensor real(pReal), dimension (3,3), intent(in) :: avgF !< initial effective stretch tensor
integer, intent(in) :: ho, me integer, intent(in) :: ho, en
integer, dimension (4) :: intFace integer, dimension (4) :: intFace
integer, dimension (3) :: iGrain3,iGNghb3,nGDim integer, dimension (3) :: iGrain3,iGNghb3,nGDim
@ -515,7 +515,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
! get the correction factor the modulus of penalty stress representing the evolution of area of ! get the correction factor the modulus of penalty stress representing the evolution of area of
! the interfaces due to deformations ! the interfaces due to deformations
surfCorr = surfaceCorrection(avgF,ho,me) surfCorr = surfaceCorrection(avgF,ho,en)
associate(prm => param(ho)) associate(prm => param(ho))
@ -527,7 +527,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
interfaceLoop: do iFace = 1,6 interfaceLoop: do iFace = 1,6
intFace = getInterface(iFace,iGrain3) ! get the 4-dimensional index of the interface in local numbering system of the grain intFace = getInterface(iFace,iGrain3) ! get the 4-dimensional index of the interface in local numbering system of the grain
nVect = interfaceNormal(intFace,ho,me) nVect = interfaceNormal(intFace,ho,en)
iGNghb3 = iGrain3 ! identify the neighboring grain across the interface iGNghb3 = iGrain3 ! identify the neighboring grain across the interface
iGNghb3(abs(intFace(1))) = iGNghb3(abs(intFace(1))) & iGNghb3(abs(intFace(1))) = iGNghb3(abs(intFace(1))) &
+ int(real(intFace(1),pReal)/real(abs(intFace(1)),pReal)) + int(real(intFace(1),pReal)/real(abs(intFace(1)),pReal))
@ -611,14 +611,14 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
!> @brief compute the correction factor accouted for surface evolution (area change) due to !> @brief compute the correction factor accouted for surface evolution (area change) due to
! deformation ! deformation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function surfaceCorrection(avgF,ho,me) function surfaceCorrection(avgF,ho,en)
real(pReal), dimension(3) :: surfaceCorrection real(pReal), dimension(3) :: surfaceCorrection
real(pReal), dimension(3,3), intent(in) :: avgF !< average F real(pReal), dimension(3,3), intent(in) :: avgF !< average F
integer, intent(in) :: & integer, intent(in) :: &
ho, & ho, &
me en
real(pReal), dimension(3,3) :: invC real(pReal), dimension(3,3) :: invC
real(pReal), dimension(3) :: nVect real(pReal), dimension(3) :: nVect
real(pReal) :: detF real(pReal) :: detF
@ -629,7 +629,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
surfaceCorrection = 0.0_pReal surfaceCorrection = 0.0_pReal
do iBase = 1,3 do iBase = 1,3
nVect = interfaceNormal([iBase,1,1,1],ho,me) nVect = interfaceNormal([iBase,1,1,1],ho,en)
do i = 1,3; do j = 1,3 do i = 1,3; do j = 1,3
surfaceCorrection(iBase) = surfaceCorrection(iBase) + invC(i,j)*nVect(i)*nVect(j) ! compute the component of (the inverse of) the stretch in the direction of the normal surfaceCorrection(iBase) = surfaceCorrection(iBase) + invC(i,j)*nVect(i)*nVect(j) ! compute the component of (the inverse of) the stretch in the direction of the normal
enddo; enddo enddo; enddo
@ -651,7 +651,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
real(pReal), dimension(6,6) :: C real(pReal), dimension(6,6) :: C
C = phase_homogenizedC(material_phaseAt2(grainID,ce),material_phaseMemberAt2(grainID,ce)) C = phase_homogenizedC(material_phaseID(grainID,ce),material_phaseEntry(grainID,ce))
equivalentMu = lattice_equivalent_mu(C,'voigt') equivalentMu = lattice_equivalent_mu(C,'voigt')
end function equivalentMu end function equivalentMu
@ -661,14 +661,14 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
!> @brief calculating the grain deformation gradient (the same with !> @brief calculating the grain deformation gradient (the same with
! homogenization_RGC_partitionDeformation, but used only for perturbation scheme) ! homogenization_RGC_partitionDeformation, but used only for perturbation scheme)
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
subroutine grainDeformation(F, avgF, ho, me) subroutine grainDeformation(F, avgF, ho, en)
real(pReal), dimension(:,:,:), intent(out) :: F !< partitioned F per grain real(pReal), dimension(:,:,:), intent(out) :: F !< partitioned F per grain
real(pReal), dimension(:,:), intent(in) :: avgF !< averaged F real(pReal), dimension(:,:), intent(in) :: avgF !< averaged F
integer, intent(in) :: & integer, intent(in) :: &
ho, & ho, &
me en
real(pReal), dimension(3) :: aVect,nVect real(pReal), dimension(3) :: aVect,nVect
integer, dimension(4) :: intFace integer, dimension(4) :: intFace
@ -685,8 +685,8 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
iGrain3 = grain1to3(iGrain,prm%N_constituents) iGrain3 = grain1to3(iGrain,prm%N_constituents)
do iFace = 1,6 do iFace = 1,6
intFace = getInterface(iFace,iGrain3) intFace = getInterface(iFace,iGrain3)
aVect = relaxationVector(intFace,ho,me) aVect = relaxationVector(intFace,ho,en)
nVect = interfaceNormal(intFace,ho,me) nVect = interfaceNormal(intFace,ho,en)
forall (i=1:3,j=1:3) & forall (i=1:3,j=1:3) &
F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations
enddo enddo
@ -697,13 +697,13 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
end subroutine grainDeformation end subroutine grainDeformation
end function mechanical_RGC_updateState end function RGC_updateState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief derive average stress and stiffness from constituent quantities !> @brief derive average stress and stiffness from constituent quantities
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mechanical_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho) module subroutine RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho)
real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point
@ -715,13 +715,13 @@ module subroutine mechanical_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dP
avgP = sum(P,3) /real(product(param(ho)%N_constituents),pReal) avgP = sum(P,3) /real(product(param(ho)%N_constituents),pReal)
dAvgPdAvgF = sum(dPdF,5)/real(product(param(ho)%N_constituents),pReal) dAvgPdAvgF = sum(dPdF,5)/real(product(param(ho)%N_constituents),pReal)
end subroutine mechanical_RGC_averageStressAndItsTangent end subroutine RGC_averageStressAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mechanical_RGC_results(ho,group) module subroutine RGC_results(ho,group)
integer, intent(in) :: ho integer, intent(in) :: ho
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
@ -747,17 +747,17 @@ module subroutine mechanical_RGC_results(ho,group)
enddo outputsLoop enddo outputsLoop
end associate end associate
end subroutine mechanical_RGC_results end subroutine RGC_results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief collect relaxation vectors of an interface !> @brief collect relaxation vectors of an interface
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function relaxationVector(intFace,ho,me) pure function relaxationVector(intFace,ho,en)
real(pReal), dimension (3) :: relaxationVector real(pReal), dimension (3) :: relaxationVector
integer, intent(in) :: ho,me integer, intent(in) :: ho,en
integer, dimension(4), intent(in) :: intFace !< set of interface ID in 4D array (normal and position) integer, dimension(4), intent(in) :: intFace !< set of interface ID in 4D array (normal and position)
integer :: iNum integer :: iNum
@ -770,7 +770,7 @@ pure function relaxationVector(intFace,ho,me)
iNum = interface4to1(intFace,prm%N_constituents) ! identify the position of the interface in global state array iNum = interface4to1(intFace,prm%N_constituents) ! identify the position of the interface in global state array
if (iNum > 0) then if (iNum > 0) then
relaxationVector = stt%relaxationVector((3*iNum-2):(3*iNum),me) relaxationVector = stt%relaxationVector((3*iNum-2):(3*iNum),en)
else else
relaxationVector = 0.0_pReal relaxationVector = 0.0_pReal
endif endif
@ -783,14 +783,14 @@ end function relaxationVector
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief identify the normal of an interface !> @brief identify the normal of an interface
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function interfaceNormal(intFace,ho,me) pure function interfaceNormal(intFace,ho,en)
real(pReal), dimension(3) :: interfaceNormal real(pReal), dimension(3) :: interfaceNormal
integer, dimension(4), intent(in) :: intFace !< interface ID in 4D array (normal and position) integer, dimension(4), intent(in) :: intFace !< interface ID in 4D array (normal and position)
integer, intent(in) :: & integer, intent(in) :: &
ho, & ho, &
me en
integer :: nPos integer :: nPos
associate (dst => dependentState(ho)) associate (dst => dependentState(ho))
@ -801,7 +801,7 @@ pure function interfaceNormal(intFace,ho,me)
nPos = abs(intFace(1)) ! identify the position of the interface in global state array nPos = abs(intFace(1)) ! identify the position of the interface in global state array
interfaceNormal(nPos) = real(intFace(1)/abs(intFace(1)),pReal) ! get the normal vector w.r.t. cluster axis interfaceNormal(nPos) = real(intFace(1)/abs(intFace(1)),pReal) ! get the normal vector w.r.t. cluster axis
interfaceNormal = matmul(dst%orientation(1:3,1:3,me),interfaceNormal) ! map the normal vector into sample coordinate system (basis) interfaceNormal = matmul(dst%orientation(1:3,1:3,en),interfaceNormal) ! map the normal vector into sample coordinate system (basis)
end associate end associate

View File

@ -26,7 +26,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields, reads information from material configuration file !> @brief allocates all neccessary fields, reads information from material configuration file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mechanical_isostrain_init module subroutine isostrain_init
integer :: & integer :: &
h, & h, &
@ -56,7 +56,7 @@ module subroutine mechanical_isostrain_init
case ('avg') case ('avg')
prm%mapping = average_ID prm%mapping = average_ID
case default case default
call IO_error(211,ext_msg='sum'//' (mechanical_isostrain)') call IO_error(211,ext_msg='sum'//' (isostrain)')
end select end select
Nmaterialpoints = count(material_homogenizationAt == h) Nmaterialpoints = count(material_homogenizationAt == h)
@ -68,13 +68,13 @@ module subroutine mechanical_isostrain_init
enddo enddo
end subroutine mechanical_isostrain_init end subroutine isostrain_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief partitions the deformation gradient onto the constituents !> @brief partitions the deformation gradient onto the constituents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mechanical_isostrain_partitionDeformation(F,avgF) module subroutine isostrain_partitionDeformation(F,avgF)
real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient
@ -82,13 +82,13 @@ module subroutine mechanical_isostrain_partitionDeformation(F,avgF)
F = spread(avgF,3,size(F,3)) F = spread(avgF,3,size(F,3))
end subroutine mechanical_isostrain_partitionDeformation end subroutine isostrain_partitionDeformation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief derive average stress and stiffness from constituent quantities !> @brief derive average stress and stiffness from constituent quantities
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mechanical_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho) module subroutine isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho)
real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point
@ -110,6 +110,6 @@ module subroutine mechanical_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvg
end associate end associate
end subroutine mechanical_isostrain_averageStressAndItsTangent end subroutine isostrain_averageStressAndItsTangent
end submodule isostrain end submodule isostrain

View File

@ -4,14 +4,14 @@
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief dummy homogenization homogenization scheme for 1 constituent per material point !> @brief dummy homogenization homogenization scheme for 1 constituent per material point
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule(homogenization:mechanical) none submodule(homogenization:mechanical) mechanical_pass
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields, reads information from material configuration file !> @brief allocates all necessary fields, reads information from material configuration file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mechanical_pass_init module subroutine pass_init
integer :: & integer :: &
Ninstances, & Ninstances, &
@ -27,7 +27,7 @@ module subroutine mechanical_pass_init
if(homogenization_type(h) /= HOMOGENIZATION_NONE_ID) cycle if(homogenization_type(h) /= HOMOGENIZATION_NONE_ID) cycle
if(homogenization_Nconstituents(h) /= 1) & if(homogenization_Nconstituents(h) /= 1) &
call IO_error(211,ext_msg='N_constituents (mechanical_pass)') call IO_error(211,ext_msg='N_constituents (pass)')
Nmaterialpoints = count(material_homogenizationAt == h) Nmaterialpoints = count(material_homogenizationAt == h)
homogState(h)%sizeState = 0 homogState(h)%sizeState = 0
@ -36,6 +36,6 @@ module subroutine mechanical_pass_init
enddo enddo
end subroutine mechanical_pass_init end subroutine pass_init
end submodule none end submodule mechanical_pass

View File

@ -5,6 +5,16 @@ submodule(homogenization) thermal
use lattice use lattice
interface
module subroutine pass_init
end subroutine pass_init
module subroutine isotemperature_init
end subroutine isotemperature_init
end interface
type :: tDataContainer type :: tDataContainer
real(pReal), dimension(:), allocatable :: T, dot_T real(pReal), dimension(:), allocatable :: T, dot_T
end type tDataContainer end type tDataContainer
@ -35,7 +45,7 @@ module subroutine thermal_init()
print'(/,a)', ' <<<+- homogenization:thermal init -+>>>' print'(/,a)', ' <<<+- homogenization:thermal init -+>>>'
print'(/,a)', ' <<<+- homogenization:thermal:isotemperature init -+>>>' print'(/,a)', ' <<<+- homogenization:thermal:pass init -+>>>'
@ -44,8 +54,8 @@ module subroutine thermal_init()
allocate(current(configHomogenizations%length)) allocate(current(configHomogenizations%length))
do ho = 1, configHomogenizations%length do ho = 1, configHomogenizations%length
allocate(current(ho)%T(count(material_homogenizationAt2==ho)), source=300.0_pReal) allocate(current(ho)%T(count(material_homogenizationID==ho)), source=300.0_pReal)
allocate(current(ho)%dot_T(count(material_homogenizationAt2==ho)), source=0.0_pReal) allocate(current(ho)%dot_T(count(material_homogenizationID==ho)), source=0.0_pReal)
configHomogenization => configHomogenizations%get(ho) configHomogenization => configHomogenizations%get(ho)
associate(prm => param(ho)) associate(prm => param(ho))
if (configHomogenization%contains('thermal')) then if (configHomogenization%contains('thermal')) then
@ -75,9 +85,9 @@ module subroutine thermal_partition(ce)
integer :: co integer :: co
T = current(material_homogenizationAt2(ce))%T(material_homogenizationMemberAt2(ce)) T = current(material_homogenizationID(ce))%T(material_homogenizationEntry(ce))
dot_T = current(material_homogenizationAt2(ce))%dot_T(material_homogenizationMemberAt2(ce)) dot_T = current(material_homogenizationID(ce))%dot_T(material_homogenizationEntry(ce))
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
call phase_thermal_setField(T,dot_T,co,ce) call phase_thermal_setField(T,dot_T,co,ce)
enddo enddo
@ -109,19 +119,29 @@ module function thermal_conduction_getConductivity(ce) result(K)
K = 0.0_pReal K = 0.0_pReal
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
K = K + crystallite_push33ToRef(co,ce,lattice_K(:,:,material_phaseAt2(co,ce))) K = K + crystallite_push33ToRef(co,ce,lattice_K(:,:,material_phaseID(co,ce)))
enddo enddo
K = K / real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) K = K / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
end function thermal_conduction_getConductivity end function thermal_conduction_getConductivity
module function homogenization_thermal_mu_T(ce) result(mu_T)
integer, intent(in) :: ce
real(pReal) :: mu_T
mu_T = c_P(ce) * rho(ce)
end function homogenization_thermal_mu_T
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns homogenized specific heat capacity !> @brief returns homogenized specific heat capacity
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function thermal_conduction_getSpecificHeat(ce) result(c_P) function c_P(ce)
integer, intent(in) :: ce integer, intent(in) :: ce
real(pReal) :: c_P real(pReal) :: c_P
@ -129,21 +149,20 @@ module function thermal_conduction_getSpecificHeat(ce) result(c_P)
integer :: co integer :: co
c_P = 0.0_pReal c_P = lattice_c_p(material_phaseID(1,ce))
do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) c_P = c_P + lattice_c_p(material_phaseID(co,ce))
c_P = c_P + lattice_c_p(material_phaseAt2(co,ce))
enddo enddo
c_P = c_P / real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) c_P = c_P / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
end function thermal_conduction_getSpecificHeat end function c_P
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns homogenized mass density !> @brief returns homogenized mass density
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function thermal_conduction_getMassDensity(ce) result(rho) function rho(ce)
integer, intent(in) :: ce integer, intent(in) :: ce
real(pReal) :: rho real(pReal) :: rho
@ -151,15 +170,14 @@ module function thermal_conduction_getMassDensity(ce) result(rho)
integer :: co integer :: co
rho = 0.0_pReal rho = lattice_rho(material_phaseID(1,ce))
do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) rho = rho + lattice_rho(material_phaseID(co,ce))
rho = rho + lattice_rho(material_phaseAt2(co,ce))
enddo enddo
rho = rho / real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) rho = rho / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
end function thermal_conduction_getMassDensity end function rho
@ -172,8 +190,8 @@ module subroutine homogenization_thermal_setField(T,dot_T, ce)
real(pReal), intent(in) :: T, dot_T real(pReal), intent(in) :: T, dot_T
current(material_homogenizationAt2(ce))%T(material_homogenizationMemberAt2(ce)) = T current(material_homogenizationID(ce))%T(material_homogenizationEntry(ce)) = T
current(material_homogenizationAt2(ce))%dot_T(material_homogenizationMemberAt2(ce)) = dot_T current(material_homogenizationID(ce))%dot_T(material_homogenizationEntry(ce)) = dot_T
end subroutine homogenization_thermal_setField end subroutine homogenization_thermal_setField
@ -183,7 +201,7 @@ end subroutine homogenization_thermal_setField
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine thermal_conduction_results(ho,group) module subroutine thermal_results(ho,group)
integer, intent(in) :: ho integer, intent(in) :: ho
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
@ -199,7 +217,7 @@ module subroutine thermal_conduction_results(ho,group)
enddo outputsLoop enddo outputsLoop
end associate end associate
end subroutine thermal_conduction_results end subroutine thermal_results
module function homogenization_thermal_T(ce) result(T) module function homogenization_thermal_T(ce) result(T)
@ -207,7 +225,7 @@ module function homogenization_thermal_T(ce) result(T)
integer, intent(in) :: ce integer, intent(in) :: ce
real(pReal) :: T real(pReal) :: T
T = current(material_homogenizationAt2(ce))%T(material_homogenizationMemberAt2(ce)) T = current(material_homogenizationID(ce))%T(material_homogenizationEntry(ce))
end function homogenization_thermal_T end function homogenization_thermal_T

View File

@ -0,0 +1,14 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, KU Leuven
!> @brief Dummy homogenization scheme for 1 constituent per material point
!--------------------------------------------------------------------------------------------------
submodule(homogenization:thermal) isotemperature
contains
module subroutine isotemperature_init
end subroutine isotemperature_init
end submodule isotemperature

View File

@ -0,0 +1,14 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, KU Leuven
!> @brief Dummy homogenization scheme for 1 constituent per material point
!--------------------------------------------------------------------------------------------------
submodule(homogenization:thermal) thermal_pass
contains
module subroutine pass_init
end subroutine pass_init
end submodule thermal_pass

View File

@ -28,14 +28,14 @@ module material
integer, dimension(:), allocatable, public, protected :: & ! (elem) integer, dimension(:), allocatable, public, protected :: & ! (elem)
material_homogenizationAt, & !< homogenization ID of each element material_homogenizationAt, & !< homogenization ID of each element
material_homogenizationAt2, & !< per cell material_homogenizationID, & !< per cell
material_homogenizationMemberAt2 !< cell material_homogenizationEntry !< cell
integer, dimension(:,:), allocatable :: & ! (ip,elem) integer, dimension(:,:), allocatable :: & ! (ip,elem)
material_homogenizationMemberAt !< position of the element within its homogenization instance material_homogenizationMemberAt !< position of the element within its homogenization instance
integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem) integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem)
material_phaseAt, & !< phase ID of each element material_phaseAt, & !< phase ID of each element
material_phaseAt2, & !< per constituent,cell material_phaseID, & !< per constituent,cell
material_phaseMemberAt2 !< per constituent, cell material_phaseEntry !< per constituent, cell
integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,IP,elem) integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,IP,elem)
material_phaseMemberAt !< position of the element within its phase instance material_phaseMemberAt !< position of the element within its phase instance
@ -117,10 +117,10 @@ subroutine material_parseMaterial
allocate(material_phaseMemberAt(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems),source=0) allocate(material_phaseMemberAt(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems),source=0)
allocate(material_homogenizationAt2(discretization_nIPs*discretization_Nelems),source=0) allocate(material_homogenizationID(discretization_nIPs*discretization_Nelems),source=0)
allocate(material_homogenizationMemberAt2(discretization_nIPs*discretization_Nelems),source=0) allocate(material_homogenizationEntry(discretization_nIPs*discretization_Nelems),source=0)
allocate(material_phaseAt2(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0) allocate(material_phaseID(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0)
allocate(material_phaseMemberAt2(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0) allocate(material_phaseEntry(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0)
do el = 1, discretization_Nelems do el = 1, discretization_Nelems
material => materials%get(discretization_materialAt(el)) material => materials%get(discretization_materialAt(el))
@ -131,8 +131,8 @@ subroutine material_parseMaterial
ce = (el-1)*discretization_nIPs + ip ce = (el-1)*discretization_nIPs + ip
counterHomogenization(material_homogenizationAt(el)) = counterHomogenization(material_homogenizationAt(el)) + 1 counterHomogenization(material_homogenizationAt(el)) = counterHomogenization(material_homogenizationAt(el)) + 1
material_homogenizationMemberAt(ip,el) = counterHomogenization(material_homogenizationAt(el)) material_homogenizationMemberAt(ip,el) = counterHomogenization(material_homogenizationAt(el))
material_homogenizationAt2(ce) = material_homogenizationAt(el) material_homogenizationID(ce) = material_homogenizationAt(el)
material_homogenizationMemberAt2(ce) = material_homogenizationMemberAt(ip,el) material_homogenizationEntry(ce) = material_homogenizationMemberAt(ip,el)
enddo enddo
frac = 0.0_pReal frac = 0.0_pReal
@ -146,8 +146,8 @@ subroutine material_parseMaterial
counterPhase(material_phaseAt(co,el)) = counterPhase(material_phaseAt(co,el)) + 1 counterPhase(material_phaseAt(co,el)) = counterPhase(material_phaseAt(co,el)) + 1
material_phaseMemberAt(co,ip,el) = counterPhase(material_phaseAt(co,el)) material_phaseMemberAt(co,ip,el) = counterPhase(material_phaseAt(co,el))
material_phaseAt2(co,ce) = material_phaseAt(co,el) material_phaseID(co,ce) = material_phaseAt(co,el)
material_phaseMemberAt2(co,ce) = material_phaseMemberAt(co,ip,el) material_phaseEntry(co,ce) = material_phaseMemberAt(co,ip,el)
enddo enddo
enddo enddo

View File

@ -85,8 +85,8 @@ program DAMASK_mesh
stagItMax = num_mesh%get_asInt('maxStaggeredIter',defaultVal=10) stagItMax = num_mesh%get_asInt('maxStaggeredIter',defaultVal=10)
maxCutBack = num_mesh%get_asInt('maxCutBack',defaultVal=3) maxCutBack = num_mesh%get_asInt('maxCutBack',defaultVal=3)
if (stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter') if (stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter')
if (maxCutBack < 0) call IO_error(301,ext_msg='maxCutBack') if (maxCutBack < 0) call IO_error(301,ext_msg='maxCutBack')
! reading basic information from load case file and allocate data structure containing load cases ! reading basic information from load case file and allocate data structure containing load cases
call DMGetDimension(geomMesh,dimPlex,ierr) !< dimension of mesh (2D or 3D) call DMGetDimension(geomMesh,dimPlex,ierr) !< dimension of mesh (2D or 3D)

View File

@ -100,14 +100,12 @@ module phase
integer, intent(in) :: ph integer, intent(in) :: ph
end subroutine damage_results end subroutine damage_results
module subroutine mechanical_windForward(ph,me)
integer, intent(in) :: ph, me
end subroutine mechanical_windForward
module subroutine mechanical_forward() module subroutine mechanical_forward()
end subroutine mechanical_forward end subroutine mechanical_forward
module subroutine damage_forward()
end subroutine damage_forward
module subroutine thermal_forward() module subroutine thermal_forward()
end subroutine thermal_forward end subroutine thermal_forward
@ -147,20 +145,20 @@ module phase
real(pReal), dimension(3,3) :: L_p real(pReal), dimension(3,3) :: L_p
end function mechanical_L_p end function mechanical_L_p
module function phase_mechanical_getF(co,ce) result(F) module function phase_F(co,ce) result(F)
integer, intent(in) :: co, ce integer, intent(in) :: co, ce
real(pReal), dimension(3,3) :: F real(pReal), dimension(3,3) :: F
end function phase_mechanical_getF end function phase_F
module function mechanical_F_e(ph,me) result(F_e) module function mechanical_F_e(ph,me) result(F_e)
integer, intent(in) :: ph,me integer, intent(in) :: ph,me
real(pReal), dimension(3,3) :: F_e real(pReal), dimension(3,3) :: F_e
end function mechanical_F_e end function mechanical_F_e
module function phase_mechanical_getP(co,ce) result(P) module function phase_P(co,ce) result(P)
integer, intent(in) :: co, ce integer, intent(in) :: co, ce
real(pReal), dimension(3,3) :: P real(pReal), dimension(3,3) :: P
end function phase_mechanical_getP end function phase_P
module function phase_damage_get_phi(co,ip,el) result(phi) module function phase_damage_get_phi(co,ip,el) result(phi)
integer, intent(in) :: co, ip, el integer, intent(in) :: co, ip, el
@ -183,10 +181,10 @@ module phase
end function damage_phi end function damage_phi
module subroutine phase_mechanical_setF(F,co,ce) module subroutine phase_set_F(F,co,ce)
real(pReal), dimension(3,3), intent(in) :: F real(pReal), dimension(3,3), intent(in) :: F
integer, intent(in) :: co, ce integer, intent(in) :: co, ce
end subroutine phase_mechanical_setF end subroutine phase_set_F
module subroutine phase_thermal_setField(T,dot_T, co,ce) module subroutine phase_thermal_setField(T,dot_T, co,ce)
real(pReal), intent(in) :: T, dot_T real(pReal), intent(in) :: T, dot_T
@ -229,14 +227,13 @@ module phase
end function phase_homogenizedC end function phase_homogenizedC
module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce) module function phase_damage_phi_dot(phi,co,ce) result(phi_dot)
integer, intent(in) :: ce integer, intent(in) :: ce,co
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi !< damage parameter phi !< damage parameter
real(pReal), intent(inout) :: & real(pReal) :: &
phiDot, & phi_dot
dPhiDot_dPhi end function phase_damage_phi_dot
end subroutine phase_damage_getRateAndItsTangents
module subroutine phase_thermal_getRate(TDot, ph,me) module subroutine phase_thermal_getRate(TDot, ph,me)
integer, intent(in) :: ph, me integer, intent(in) :: ph, me
@ -261,7 +258,7 @@ module phase
end subroutine plastic_dependentState end subroutine plastic_dependentState
module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me)
integer, intent(in) :: ph, me integer, intent(in) :: ph, me
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
S S
@ -269,9 +266,9 @@ module phase
Ld !< damage velocity gradient Ld !< damage velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: & real(pReal), intent(out), dimension(3,3,3,3) :: &
dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor)
end subroutine kinematics_cleavage_opening_LiAndItsTangent end subroutine damage_anisobrittle_LiAndItsTangent
module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) module subroutine damage_isoductile_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me)
integer, intent(in) :: ph, me integer, intent(in) :: ph, me
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
S S
@ -279,7 +276,7 @@ module phase
Ld !< damage velocity gradient Ld !< damage velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: & real(pReal), intent(out), dimension(3,3,3,3) :: &
dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor)
end subroutine kinematics_slipplane_opening_LiAndItsTangent end subroutine damage_isoductile_LiAndItsTangent
end interface end interface
@ -303,7 +300,7 @@ module phase
public :: & public :: &
phase_init, & phase_init, &
phase_homogenizedC, & phase_homogenizedC, &
phase_damage_getRateAndItsTangents, & phase_damage_phi_dot, &
phase_thermal_getRate, & phase_thermal_getRate, &
phase_results, & phase_results, &
phase_allocateState, & phase_allocateState, &
@ -323,21 +320,19 @@ module phase
phase_thermal_setField, & phase_thermal_setField, &
phase_damage_set_phi, & phase_damage_set_phi, &
phase_damage_get_phi, & phase_damage_get_phi, &
phase_mechanical_getP, & phase_P, &
phase_mechanical_setF, & phase_set_F, &
phase_mechanical_getF, & phase_F
phase_windForward
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Initialze constitutive models for individual physics !> @brief Initialize constitutive models for individual physics
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine phase_init subroutine phase_init
integer :: & integer :: &
ph, & !< counter in phase loop ph
so !< counter in source loop
class (tNode), pointer :: & class (tNode), pointer :: &
debug_constitutive, & debug_constitutive, &
materials, & materials, &
@ -382,12 +377,12 @@ end subroutine phase_init
!> @brief Allocate the components of the state structure for a given phase !> @brief Allocate the components of the state structure for a given phase
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine phase_allocateState(state, & subroutine phase_allocateState(state, &
Nconstituents,sizeState,sizeDotState,sizeDeltaState) NEntries,sizeState,sizeDotState,sizeDeltaState)
class(tState), intent(out) :: & class(tState), intent(out) :: &
state state
integer, intent(in) :: & integer, intent(in) :: &
Nconstituents, & NEntries, &
sizeState, & sizeState, &
sizeDotState, & sizeDotState, &
sizeDeltaState sizeDeltaState
@ -398,13 +393,13 @@ subroutine phase_allocateState(state, &
state%sizeDeltaState = sizeDeltaState state%sizeDeltaState = sizeDeltaState
state%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition state%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition
allocate(state%atol (sizeState), source=0.0_pReal) allocate(state%atol (sizeState), source=0.0_pReal)
allocate(state%state0 (sizeState,Nconstituents), source=0.0_pReal) allocate(state%state0 (sizeState,NEntries), source=0.0_pReal)
allocate(state%state (sizeState,Nconstituents), source=0.0_pReal) allocate(state%state (sizeState,NEntries), source=0.0_pReal)
allocate(state%dotState (sizeDotState,Nconstituents), source=0.0_pReal) allocate(state%dotState (sizeDotState,NEntries), source=0.0_pReal)
allocate(state%deltaState (sizeDeltaState,Nconstituents), source=0.0_pReal) allocate(state%deltaState (sizeDeltaState,NEntries), source=0.0_pReal)
end subroutine phase_allocateState end subroutine phase_allocateState
@ -422,10 +417,10 @@ subroutine phase_restore(ce,includeL)
co co
do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
if (damageState(material_phaseAt2(co,ce))%sizeState > 0) & if (damageState(material_phaseID(co,ce))%sizeState > 0) &
damageState(material_phaseAt2(co,ce))%state( :,material_phasememberAt2(co,ce)) = & damageState(material_phaseID(co,ce))%state( :,material_phaseEntry(co,ce)) = &
damageState(material_phaseAt2(co,ce))%state0(:,material_phasememberAt2(co,ce)) damageState(material_phaseID(co,ce))%state0(:,material_phaseEntry(co,ce))
enddo enddo
call mechanical_restore(ce,includeL) call mechanical_restore(ce,includeL)
@ -435,21 +430,13 @@ end subroutine phase_restore
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Forward data after successful increment. !> @brief Forward data after successful increment.
! ToDo: Any guessing for the current states possible?
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine phase_forward() subroutine phase_forward()
integer :: ph
call mechanical_forward() call mechanical_forward()
call damage_forward()
call thermal_forward() call thermal_forward()
do ph = 1, size(damageState)
if (damageState(ph)%sizeState > 0) &
damageState(ph)%state0 = damageState(ph)%state
enddo
end subroutine phase_forward end subroutine phase_forward
@ -484,11 +471,9 @@ subroutine crystallite_init()
integer :: & integer :: &
ph, & ph, &
me, &
co, & !< counter in integration point component loop co, & !< counter in integration point component loop
ip, & !< counter in integration point loop ip, & !< counter in integration point loop
el, & !< counter in element loop el, & !< counter in element loop
so, &
cMax, & !< maximum number of integration point components cMax, & !< maximum number of integration point components
iMax, & !< maximum number of integration points iMax, & !< maximum number of integration points
eMax !< maximum number of elements eMax !< maximum number of elements
@ -550,12 +535,10 @@ subroutine crystallite_init()
flush(IO_STDOUT) flush(IO_STDOUT)
!$OMP PARALLEL DO PRIVATE(ph,me) !$OMP PARALLEL DO
do el = 1, size(material_phaseMemberAt,3) do el = 1, size(material_phaseMemberAt,3)
do ip = 1, size(material_phaseMemberAt,2) do ip = 1, size(material_phaseMemberAt,2)
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
call crystallite_orientations(co,ip,el) call crystallite_orientations(co,ip,el)
call plastic_dependentState(co,ip,el) ! update dependent state variables to be consistent with basic states call plastic_dependentState(co,ip,el) ! update dependent state variables to be consistent with basic states
enddo enddo
@ -567,34 +550,6 @@ subroutine crystallite_init()
end subroutine crystallite_init end subroutine crystallite_init
!--------------------------------------------------------------------------------------------------
!> @brief Wind homog inc forward.
!--------------------------------------------------------------------------------------------------
subroutine phase_windForward(ip,el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
integer :: &
co, & !< constituent number
so, ph, me
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
call mechanical_windForward(ph,me)
if(damageState(ph)%sizeState > 0) damageState(ph)%state0(:,me) = damageState(ph)%state(:,me)
enddo
end subroutine phase_windForward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates orientations !> @brief calculates orientations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -629,11 +584,11 @@ function crystallite_push33ToRef(co,ce, tensor33)
real(pReal), dimension(3,3) :: crystallite_push33ToRef real(pReal), dimension(3,3) :: crystallite_push33ToRef
real(pReal), dimension(3,3) :: T real(pReal), dimension(3,3) :: T
integer :: ph, me integer :: ph, en
ph = material_phaseAt2(co,ce) ph = material_phaseID(co,ce)
me = material_phaseMemberAt2(co,ce) en = material_phaseEntry(co,ce)
T = matmul(material_orientation0(co,ph,me)%asMatrix(),transpose(math_inv33(phase_mechanical_getF(co,ce)))) ! ToDo: initial orientation correct? T = matmul(material_orientation0(co,ph,en)%asMatrix(),transpose(math_inv33(phase_F(co,ce)))) ! ToDo: initial orientation correct?
crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T)) crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T))
@ -658,8 +613,7 @@ end function converged
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Write current restart information (Field and constitutive data) to file. !> @brief Write restart data to file.
! ToDo: Merge data into one file for MPI
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine phase_restartWrite(fileHandle) subroutine phase_restartWrite(fileHandle)
@ -687,8 +641,7 @@ end subroutine phase_restartWrite
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Read data for restart !> @brief Read restart data from file.
! ToDo: Merge data into one file for MPI
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine phase_restartRead(fileHandle) subroutine phase_restartRead(fileHandle)

View File

@ -1,7 +1,7 @@
!---------------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------------
!> @brief internal microstructure state for all damage sources and kinematics constitutive models !> @brief internal microstructure state for all damage sources and kinematics constitutive models
!---------------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------------
submodule(phase) damagee submodule(phase) damage
enum, bind(c); enumerator :: & enum, bind(c); enumerator :: &
DAMAGE_UNDEFINED_ID, & DAMAGE_UNDEFINED_ID, &
DAMAGE_ISOBRITTLE_ID, & DAMAGE_ISOBRITTLE_ID, &
@ -65,43 +65,6 @@ submodule(phase) damagee
integer, intent(in) :: ph,me integer, intent(in) :: ph,me
end subroutine isoductile_dotState end subroutine isoductile_dotState
module subroutine anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me)
integer, intent(in) :: ph,me
real(pReal), intent(in) :: &
phi !< damage parameter
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
end subroutine anisobrittle_getRateAndItsTangent
module subroutine anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me)
integer, intent(in) :: ph,me
real(pReal), intent(in) :: &
phi !< damage parameter
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
end subroutine anisoductile_getRateAndItsTangent
module subroutine isobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me)
integer, intent(in) :: ph,me
real(pReal), intent(in) :: &
phi !< damage parameter
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
end subroutine isobrittle_getRateAndItsTangent
module subroutine isoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me)
integer, intent(in) :: ph,me
real(pReal), intent(in) :: &
phi !< damage parameter
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
end subroutine isoductile_getRateAndItsTangent
module subroutine anisobrittle_results(phase,group) module subroutine anisobrittle_results(phase,group)
integer, intent(in) :: phase integer, intent(in) :: phase
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
@ -151,7 +114,7 @@ module subroutine damage_init
do ph = 1,phases%length do ph = 1,phases%length
Nmembers = count(material_phaseAt2 == ph) Nmembers = count(material_phaseID == ph)
allocate(current(ph)%phi(Nmembers),source=1.0_pReal) allocate(current(ph)%phi(Nmembers),source=1.0_pReal)
allocate(current(ph)%d_phi_d_dot_phi(Nmembers),source=0.0_pReal) allocate(current(ph)%d_phi_d_dot_phi(Nmembers),source=0.0_pReal)
@ -179,53 +142,30 @@ end subroutine damage_init
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief returns local part of nonlocal damage driving force !< @brief returns local part of nonlocal damage driving force
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce) module function phase_damage_phi_dot(phi,co,ce) result(phi_dot)
integer, intent(in) :: ce integer, intent(in) :: ce,co
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi !< damage parameter phi !< damage parameter
real(pReal), intent(inout) :: &
phiDot, &
dPhiDot_dPhi
real(pReal) :: & real(pReal) :: &
localphiDot, & phi_dot
dLocalphiDot_dPhi
integer :: & integer :: &
ph, & ph, &
co, & en
me
phiDot = 0.0_pReal ph = material_phaseID(co,ce)
dPhiDot_dPhi = 0.0_pReal en = material_phaseEntry(co,ce)
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) select case(phase_source(ph))
ph = material_phaseAt2(co,ce) case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ISODUCTILE_ID,DAMAGE_ANISOBRITTLE_ID,DAMAGE_ANISODUCTILE_ID)
me = material_phasememberAt2(co,ce) phi_dot = 1.0_pReal &
- phi*damageState(ph)%state(1,en)
case default
phi_dot = 0.0_pReal
end select
select case(phase_source(ph)) end function phase_damage_phi_dot
case (DAMAGE_ISOBRITTLE_ID)
call isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, ph, me)
case (DAMAGE_ISODUCTILE_ID)
call isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, ph, me)
case (DAMAGE_ANISOBRITTLE_ID)
call anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me)
case (DAMAGE_ANISODUCTILE_ID)
call anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me)
case default
localphiDot = 0.0_pReal
dLocalphiDot_dPhi = 0.0_pReal
end select
phiDot = phiDot + localphiDot
dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi
enddo
end subroutine phase_damage_getRateAndItsTangents
@ -477,7 +417,7 @@ module subroutine phase_damage_set_phi(phi,co,ce)
integer, intent(in) :: ce, co integer, intent(in) :: ce, co
current(material_phaseAt2(co,ce))%phi(material_phaseMemberAt2(co,ce)) = phi current(material_phaseID(co,ce))%phi(material_phaseEntry(co,ce)) = phi
end subroutine phase_damage_set_phi end subroutine phase_damage_set_phi
@ -503,4 +443,21 @@ module function damage_phi(ph,me) result(phi)
end function damage_phi end function damage_phi
end submodule damagee
!--------------------------------------------------------------------------------------------------
!> @brief Forward data after successful increment.
!--------------------------------------------------------------------------------------------------
module subroutine damage_forward()
integer :: ph
do ph = 1, size(damageState)
if (damageState(ph)%sizeState > 0) &
damageState(ph)%state0 = damageState(ph)%state
enddo
end subroutine damage_forward
end submodule damage

View File

@ -4,7 +4,7 @@
!> @brief material subroutine incorporating anisotropic brittle damage source mechanism !> @brief material subroutine incorporating anisotropic brittle damage source mechanism
!> @details to be done !> @details to be done
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule (phase:damagee) anisobrittle submodule (phase:damage) anisobrittle
type :: tParameters !< container type for internal constitutive parameters type :: tParameters !< container type for internal constitutive parameters
real(pReal) :: & real(pReal) :: &
@ -120,9 +120,6 @@ module subroutine anisobrittle_dotState(S, ph,me)
S S
integer :: & integer :: &
sourceOffset, &
damageOffset, &
homog, &
i i
real(pReal) :: & real(pReal) :: &
traction_d, traction_t, traction_n, traction_crit traction_d, traction_t, traction_n, traction_crit
@ -148,29 +145,6 @@ module subroutine anisobrittle_dotState(S, ph,me)
end subroutine anisobrittle_dotState end subroutine anisobrittle_dotState
!--------------------------------------------------------------------------------------------------
!> @brief returns local part of nonlocal damage driving force
!--------------------------------------------------------------------------------------------------
module subroutine anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me)
integer, intent(in) :: &
ph, &
me
real(pReal), intent(in) :: &
phi
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
dLocalphiDot_dPhi = -damageState(ph)%state(1,me)
localphiDot = 1.0_pReal &
+ dLocalphiDot_dPhi*phi
end subroutine anisobrittle_getRateAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -197,7 +171,7 @@ end subroutine anisobrittle_results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient !> @brief contains the constitutive equation for calculating the velocity gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me)
integer, intent(in) :: & integer, intent(in) :: &
ph,me ph,me
@ -253,6 +227,6 @@ module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S,
enddo enddo
end associate end associate
end subroutine kinematics_cleavage_opening_LiAndItsTangent end subroutine damage_anisobrittle_LiAndItsTangent
end submodule anisobrittle end submodule anisobrittle

View File

@ -4,7 +4,7 @@
!> @brief material subroutine incorporating anisotropic ductile damage source mechanism !> @brief material subroutine incorporating anisotropic ductile damage source mechanism
!> @details to be done !> @details to be done
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule(phase:damagee) anisoductile submodule(phase:damage) anisoductile
type :: tParameters !< container type for internal constitutive parameters type :: tParameters !< container type for internal constitutive parameters
real(pReal) :: & real(pReal) :: &
@ -35,7 +35,7 @@ module function anisoductile_init() result(mySources)
pl, & pl, &
sources, & sources, &
src src
integer :: Ninstances,Nmembers,p integer :: Ninstances,Nmembers,ph
integer, dimension(:), allocatable :: N_sl integer, dimension(:), allocatable :: N_sl
character(len=pStringLen) :: extmsg = '' character(len=pStringLen) :: extmsg = ''
@ -50,15 +50,15 @@ module function anisoductile_init() result(mySources)
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(phases%length)) allocate(param(phases%length))
do p = 1, phases%length do ph = 1, phases%length
if(mySources(p)) then if(mySources(ph)) then
phase => phases%get(p) phase => phases%get(ph)
mech => phase%get('mechanical') mech => phase%get('mechanical')
pl => mech%get('plastic') pl => mech%get('plastic')
sources => phase%get('damage') sources => phase%get('damage')
associate(prm => param(p)) associate(prm => param(ph))
src => sources%get(1) src => sources%get(1)
N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
@ -78,10 +78,10 @@ module function anisoductile_init() result(mySources)
if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit' if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit'
Nmembers=count(material_phaseAt2==p) Nmembers=count(material_phaseID==ph)
call phase_allocateState(damageState(p),Nmembers,1,1,0) call phase_allocateState(damageState(ph),Nmembers,1,1,0)
damageState(p)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal) damageState(ph)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol'
end associate end associate
@ -113,29 +113,6 @@ module subroutine anisoductile_dotState(ph,me)
end subroutine anisoductile_dotState end subroutine anisoductile_dotState
!--------------------------------------------------------------------------------------------------
!> @brief returns local part of nonlocal damage driving force
!--------------------------------------------------------------------------------------------------
module subroutine anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me)
integer, intent(in) :: &
ph, &
me
real(pReal), intent(in) :: &
phi
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
dLocalphiDot_dPhi = -damageState(ph)%state(1,me)
localphiDot = 1.0_pReal &
+ dLocalphiDot_dPhi*phi
end subroutine anisoductile_getRateAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

View File

@ -4,7 +4,7 @@
!> @brief material subroutine incoprorating isotropic brittle damage source mechanism !> @brief material subroutine incoprorating isotropic brittle damage source mechanism
!> @details to be done !> @details to be done
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule(phase:damagee) isobrittle submodule(phase:damage) isobrittle
type :: tParameters !< container type for internal constitutive parameters type :: tParameters !< container type for internal constitutive parameters
real(pReal) :: & real(pReal) :: &
@ -31,7 +31,7 @@ module function isobrittle_init() result(mySources)
phase, & phase, &
sources, & sources, &
src src
integer :: Nmembers,p integer :: Nmembers,ph
character(len=pStringLen) :: extmsg = '' character(len=pStringLen) :: extmsg = ''
@ -45,12 +45,12 @@ module function isobrittle_init() result(mySources)
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(phases%length)) allocate(param(phases%length))
do p = 1, phases%length do ph = 1, phases%length
if(mySources(p)) then if(mySources(ph)) then
phase => phases%get(p) phase => phases%get(ph)
sources => phase%get('damage') sources => phase%get('damage')
associate(prm => param(p)) associate(prm => param(ph))
src => sources%get(1) src => sources%get(1)
prm%W_crit = src%get_asFloat('W_crit') prm%W_crit = src%get_asFloat('W_crit')
@ -64,10 +64,10 @@ module function isobrittle_init() result(mySources)
! sanity checks ! sanity checks
if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit' if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit'
Nmembers = count(material_phaseAt2==p) Nmembers = count(material_phaseID==ph)
call phase_allocateState(damageState(p),Nmembers,1,1,1) call phase_allocateState(damageState(ph),Nmembers,1,1,1)
damageState(p)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal) damageState(ph)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol'
end associate end associate
@ -113,29 +113,6 @@ module subroutine isobrittle_deltaState(C, Fe, ph,me)
end subroutine isobrittle_deltaState end subroutine isobrittle_deltaState
!--------------------------------------------------------------------------------------------------
!> @brief returns local part of nonlocal damage driving force
!--------------------------------------------------------------------------------------------------
module subroutine isobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me)
integer, intent(in) :: &
ph, me
real(pReal), intent(in) :: &
phi
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
associate(prm => param(ph))
localphiDot = 1.0_pReal &
- phi*damageState(ph)%state(1,me)
dLocalphiDot_dPhi = - damageState(ph)%state(1,me)
end associate
end subroutine isobrittle_getRateAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

View File

@ -4,7 +4,7 @@
!> @brief material subroutine incorporating isotropic ductile damage source mechanism !> @brief material subroutine incorporating isotropic ductile damage source mechanism
!> @details to be done !> @details to be done
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule(phase:damagee) isoductile submodule(phase:damage) isoductile
type:: tParameters !< container type for internal constitutive parameters type:: tParameters !< container type for internal constitutive parameters
real(pReal) :: & real(pReal) :: &
@ -33,7 +33,7 @@ module function isoductile_init() result(mySources)
phase, & phase, &
sources, & sources, &
src src
integer :: Ninstances,Nmembers,p integer :: Ninstances,Nmembers,ph
character(len=pStringLen) :: extmsg = '' character(len=pStringLen) :: extmsg = ''
@ -47,12 +47,12 @@ module function isoductile_init() result(mySources)
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(phases%length)) allocate(param(phases%length))
do p = 1, phases%length do ph = 1, phases%length
if(mySources(p)) then if(mySources(ph)) then
phase => phases%get(p) phase => phases%get(ph)
sources => phase%get('damage') sources => phase%get('damage')
associate(prm => param(p)) associate(prm => param(ph))
src => sources%get(1) src => sources%get(1)
prm%q = src%get_asFloat('q') prm%q = src%get_asFloat('q')
@ -68,10 +68,10 @@ module function isoductile_init() result(mySources)
if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit' if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit'
Nmembers=count(material_phaseAt2==p) Nmembers=count(material_phaseID==ph)
call phase_allocateState(damageState(p),Nmembers,1,1,0) call phase_allocateState(damageState(ph),Nmembers,1,1,0)
damageState(p)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal) damageState(ph)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol'
end associate end associate
@ -103,29 +103,6 @@ module subroutine isoductile_dotState(ph, me)
end subroutine isoductile_dotState end subroutine isoductile_dotState
!--------------------------------------------------------------------------------------------------
!> @brief returns local part of nonlocal damage driving force
!--------------------------------------------------------------------------------------------------
module subroutine isoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me)
integer, intent(in) :: &
ph, &
me
real(pReal), intent(in) :: &
phi
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
dLocalphiDot_dPhi = -damageState(ph)%state(1,me)
localphiDot = 1.0_pReal &
+ dLocalphiDot_dPhi*phi
end subroutine isoductile_getRateAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

View File

@ -1060,26 +1060,6 @@ subroutine crystallite_results(group,ph)
end subroutine crystallite_results end subroutine crystallite_results
!--------------------------------------------------------------------------------------------------
!> @brief Wind homog inc forward.
!--------------------------------------------------------------------------------------------------
module subroutine mechanical_windForward(ph,me)
integer, intent(in) :: ph, me
phase_mechanical_Fp0(ph)%data(1:3,1:3,me) = phase_mechanical_Fp(ph)%data(1:3,1:3,me)
phase_mechanical_Fi0(ph)%data(1:3,1:3,me) = phase_mechanical_Fi(ph)%data(1:3,1:3,me)
phase_mechanical_F0(ph)%data(1:3,1:3,me) = phase_mechanical_F(ph)%data(1:3,1:3,me)
phase_mechanical_Li0(ph)%data(1:3,1:3,me) = phase_mechanical_Li(ph)%data(1:3,1:3,me)
phase_mechanical_Lp0(ph)%data(1:3,1:3,me) = phase_mechanical_Lp(ph)%data(1:3,1:3,me)
phase_mechanical_S0(ph)%data(1:3,1:3,me) = phase_mechanical_S(ph)%data(1:3,1:3,me)
plasticState(ph)%State0(:,me) = plasticState(ph)%state(:,me)
end subroutine mechanical_windForward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Forward data after successful increment. !> @brief Forward data after successful increment.
! ToDo: Any guessing for the current states possible? ! ToDo: Any guessing for the current states possible?
@ -1213,8 +1193,6 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
if (todo) then if (todo) then
subF = subF0 & subF = subF0 &
+ subStep * (phase_mechanical_F(ph)%data(1:3,1:3,me) - phase_mechanical_F0(ph)%data(1:3,1:3,me)) + subStep * (phase_mechanical_F(ph)%data(1:3,1:3,me) - phase_mechanical_F0(ph)%data(1:3,1:3,me))
phase_mechanical_Fe(ph)%data(1:3,1:3,me) = matmul(subF,math_inv33(matmul(phase_mechanical_Fi(ph)%data(1:3,1:3,me), &
phase_mechanical_Fp(ph)%data(1:3,1:3,me))))
converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el) converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el)
converged_ = converged_ .and. .not. integrateDamageState(subStep * dt,co,ip,el) converged_ = converged_ .and. .not. integrateDamageState(subStep * dt,co,ip,el)
endif endif
@ -1237,9 +1215,9 @@ module subroutine mechanical_restore(ce,includeL)
co, ph, me co, ph, me
do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
ph = material_phaseAt2(co,ce) ph = material_phaseID(co,ce)
me = material_phaseMemberAt2(co,ce) me = material_phaseEntry(co,ce)
if (includeL) then if (includeL) then
phase_mechanical_Lp(ph)%data(1:3,1:3,me) = phase_mechanical_Lp0(ph)%data(1:3,1:3,me) phase_mechanical_Lp(ph)%data(1:3,1:3,me) = phase_mechanical_Lp0(ph)%data(1:3,1:3,me)
phase_mechanical_Li(ph)%data(1:3,1:3,me) = phase_mechanical_Li0(ph)%data(1:3,1:3,me) phase_mechanical_Li(ph)%data(1:3,1:3,me) = phase_mechanical_Li0(ph)%data(1:3,1:3,me)
@ -1287,8 +1265,8 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF)
logical :: error logical :: error
ph = material_phaseAt2(co,ce) ph = material_phaseID(co,ce)
me = material_phaseMemberAt2(co,ce) me = material_phaseEntry(co,ce)
call phase_hooke_SandItsTangents(devNull,dSdFe,dSdFi, & call phase_hooke_SandItsTangents(devNull,dSdFe,dSdFi, &
phase_mechanical_Fe(ph)%data(1:3,1:3,me), & phase_mechanical_Fe(ph)%data(1:3,1:3,me), &
@ -1443,20 +1421,6 @@ module function mechanical_L_p(ph,me) result(L_p)
end function mechanical_L_p end function mechanical_L_p
!----------------------------------------------------------------------------------------------
!< @brief Get deformation gradient (for use by homogenization)
!----------------------------------------------------------------------------------------------
module function phase_mechanical_getF(co,ce) result(F)
integer, intent(in) :: co, ce
real(pReal), dimension(3,3) :: F
F = phase_mechanical_F(material_phaseAt2(co,ce))%data(1:3,1:3,material_phaseMemberAt2(co,ce))
end function phase_mechanical_getF
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief Get elastic deformation gradient (for use by non-mech physics) !< @brief Get elastic deformation gradient (for use by non-mech physics)
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
@ -1471,31 +1435,46 @@ module function mechanical_F_e(ph,me) result(F_e)
end function mechanical_F_e end function mechanical_F_e
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief Get second Piola-Kichhoff stress (for use by homogenization) !< @brief Get second Piola-Kichhoff stress (for use by homogenization)
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module function phase_mechanical_getP(co,ce) result(P) module function phase_P(co,ce) result(P)
integer, intent(in) :: co, ce integer, intent(in) :: co, ce
real(pReal), dimension(3,3) :: P real(pReal), dimension(3,3) :: P
P = phase_mechanical_P(material_phaseAt2(co,ce))%data(1:3,1:3,material_phaseMemberAt2(co,ce)) P = phase_mechanical_P(material_phaseID(co,ce))%data(1:3,1:3,material_phaseEntry(co,ce))
end function phase_mechanical_getP end function phase_P
! setter for homogenization !----------------------------------------------------------------------------------------------
module subroutine phase_mechanical_setF(F,co,ce) !< @brief Get deformation gradient (for use by homogenization)
!----------------------------------------------------------------------------------------------
module function phase_F(co,ce) result(F)
integer, intent(in) :: co, ce
real(pReal), dimension(3,3) :: F
F = phase_mechanical_F(material_phaseID(co,ce))%data(1:3,1:3,material_phaseEntry(co,ce))
end function phase_F
!----------------------------------------------------------------------------------------------
!< @brief Set deformation gradient (for use by homogenization)
!----------------------------------------------------------------------------------------------
module subroutine phase_set_F(F,co,ce)
real(pReal), dimension(3,3), intent(in) :: F real(pReal), dimension(3,3), intent(in) :: F
integer, intent(in) :: co, ce integer, intent(in) :: co, ce
phase_mechanical_F(material_phaseAt2(co,ce))%data(1:3,1:3,material_phaseMemberAt2(co,ce)) = F phase_mechanical_F(material_phaseID(co,ce))%data(1:3,1:3,material_phaseEntry(co,ce)) = F
end subroutine phase_mechanical_setF end subroutine phase_set_F
end submodule mechanical end submodule mechanical

View File

@ -9,13 +9,13 @@ submodule(phase:mechanical) eigen
model_damage model_damage
interface interface
module function kinematics_cleavage_opening_init() result(myKinematics) module function damage_anisobrittle_init() result(myKinematics)
logical, dimension(:), allocatable :: myKinematics logical, dimension(:), allocatable :: myKinematics
end function kinematics_cleavage_opening_init end function damage_anisobrittle_init
module function kinematics_slipplane_opening_init() result(myKinematics) module function damage_isoductile_init() result(myKinematics)
logical, dimension(:), allocatable :: myKinematics logical, dimension(:), allocatable :: myKinematics
end function kinematics_slipplane_opening_init end function damage_isoductile_init
module function thermalexpansion_init(kinematics_length) result(myKinematics) module function thermalexpansion_init(kinematics_length) result(myKinematics)
integer, intent(in) :: kinematics_length integer, intent(in) :: kinematics_length
@ -46,7 +46,6 @@ module subroutine eigendeformation_init(phases)
class(tNode), pointer :: & class(tNode), pointer :: &
phase, & phase, &
kinematics, & kinematics, &
damage, &
mechanics mechanics
print'(/,a)', ' <<<+- phase:mechanical:eigen init -+>>>' print'(/,a)', ' <<<+- phase:mechanical:eigen init -+>>>'
@ -70,8 +69,8 @@ module subroutine eigendeformation_init(phases)
allocate(model_damage(phases%length), source = KINEMATICS_UNDEFINED_ID) allocate(model_damage(phases%length), source = KINEMATICS_UNDEFINED_ID)
where(kinematics_cleavage_opening_init()) model_damage = KINEMATICS_cleavage_opening_ID where(damage_anisobrittle_init()) model_damage = KINEMATICS_cleavage_opening_ID
where(kinematics_slipplane_opening_init()) model_damage = KINEMATICS_slipplane_opening_ID where(damage_isoductile_init()) model_damage = KINEMATICS_slipplane_opening_ID
end subroutine eigendeformation_init end subroutine eigendeformation_init
@ -198,12 +197,12 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
select case (model_damage(ph)) select case (model_damage(ph))
case (KINEMATICS_cleavage_opening_ID) case (KINEMATICS_cleavage_opening_ID)
call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me) call damage_anisobrittle_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me)
Li = Li + my_Li Li = Li + my_Li
dLi_dS = dLi_dS + my_dLi_dS dLi_dS = dLi_dS + my_dLi_dS
active = .true. active = .true.
case (KINEMATICS_slipplane_opening_ID) case (KINEMATICS_slipplane_opening_ID)
call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me) call damage_isoductile_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me)
Li = Li + my_Li Li = Li + my_Li
dLi_dS = dLi_dS + my_dLi_dS dLi_dS = dLi_dS + my_dLi_dS
active = .true. active = .true.

View File

@ -13,7 +13,7 @@ contains
!> @brief module initialization !> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function kinematics_cleavage_opening_init() result(myKinematics) module function damage_anisobrittle_init() result(myKinematics)
logical, dimension(:), allocatable :: myKinematics logical, dimension(:), allocatable :: myKinematics
@ -24,7 +24,7 @@ module function kinematics_cleavage_opening_init() result(myKinematics)
print'(/,a)', ' <<<+- phase:mechanical:eigen:cleavageopening init -+>>>' print'(/,a)', ' <<<+- phase:mechanical:eigen:cleavageopening init -+>>>'
print'(a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT) print'(a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT)
end function kinematics_cleavage_opening_init end function damage_anisobrittle_init
end submodule cleavageopening end submodule cleavageopening

View File

@ -6,7 +6,7 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule(phase:eigen) slipplaneopening submodule(phase:eigen) slipplaneopening
integer, dimension(:), allocatable :: kinematics_slipplane_opening_instance integer, dimension(:), allocatable :: damage_isoductile_instance
type :: tParameters !< container type for internal constitutive parameters type :: tParameters !< container type for internal constitutive parameters
integer :: & integer :: &
@ -32,7 +32,7 @@ contains
!> @brief module initialization !> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function kinematics_slipplane_opening_init() result(myKinematics) module function damage_isoductile_init() result(myKinematics)
logical, dimension(:), allocatable :: myKinematics logical, dimension(:), allocatable :: myKinematics
@ -107,13 +107,13 @@ module function kinematics_slipplane_opening_init() result(myKinematics)
enddo enddo
end function kinematics_slipplane_opening_init end function damage_isoductile_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient !> @brief contains the constitutive equation for calculating the velocity gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) module subroutine damage_isoductile_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me)
integer, intent(in) :: & integer, intent(in) :: &
ph, me ph, me
@ -179,6 +179,6 @@ module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S
end associate end associate
end subroutine kinematics_slipplane_opening_LiAndItsTangent end subroutine damage_isoductile_LiAndItsTangent
end submodule slipplaneopening end submodule slipplaneopening

View File

@ -220,7 +220,7 @@ module function plastic_dislotungsten_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
Nmembers = count(material_phaseAt2 == ph) Nmembers = count(material_phaseID == ph)
sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl
sizeState = sizeDotState sizeState = sizeDotState

View File

@ -406,7 +406,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
Nmembers = count(material_phaseAt2 == ph) Nmembers = count(material_phaseID == ph)
sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl & sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl &
+ size(['f_tw']) * prm%sum_N_tw & + size(['f_tw']) * prm%sum_N_tw &
+ size(['f_tr']) * prm%sum_N_tr + size(['f_tr']) * prm%sum_N_tr

View File

@ -119,7 +119,7 @@ module function plastic_isotropic_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
Nmembers = count(material_phaseAt2 == ph) Nmembers = count(material_phaseID == ph)
sizeDotState = size(['xi ','gamma']) sizeDotState = size(['xi ','gamma'])
sizeState = sizeDotState sizeState = sizeDotState

View File

@ -165,7 +165,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
Nmembers = count(material_phaseAt2 == ph) Nmembers = count(material_phaseID == ph)
sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%sum_N_sl !ToDo: adjust names like in material.yaml sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%sum_N_sl !ToDo: adjust names like in material.yaml
sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names like in material.yaml sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names like in material.yaml
sizeState = sizeDotState + sizeDeltaState sizeState = sizeDotState + sizeDeltaState

View File

@ -30,7 +30,7 @@ module function plastic_none_init() result(myPlasticity)
phases => config_material%get('phase') phases => config_material%get('phase')
do ph = 1, phases%length do ph = 1, phases%length
if(.not. myPlasticity(ph)) cycle if(.not. myPlasticity(ph)) cycle
call phase_allocateState(plasticState(ph),count(material_phaseAt2 == ph),0,0,0) call phase_allocateState(plasticState(ph),count(material_phaseID == ph),0,0,0)
enddo enddo
end function plastic_none_init end function plastic_none_init

View File

@ -398,7 +398,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
Nmembers = count(material_phaseAt2 == ph) Nmembers = count(material_phaseID == ph)
sizeDotState = size([ 'rhoSglEdgePosMobile ','rhoSglEdgeNegMobile ', & sizeDotState = size([ 'rhoSglEdgePosMobile ','rhoSglEdgeNegMobile ', &
'rhoSglScrewPosMobile ','rhoSglScrewNegMobile ', & 'rhoSglScrewPosMobile ','rhoSglScrewNegMobile ', &
'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', & 'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', &
@ -527,7 +527,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
if(.not. myPlasticity(ph)) cycle if(.not. myPlasticity(ph)) cycle
phase => phases%get(ph) phase => phases%get(ph)
Nmembers = count(material_phaseAt2 == ph) Nmembers = count(material_phaseID == ph)
l = 0 l = 0
do t = 1,4 do t = 1,4
do s = 1,param(ph)%sum_N_sl do s = 1,param(ph)%sum_N_sl
@ -1824,9 +1824,9 @@ subroutine storeGeometry(ph)
V = reshape(IPvolume,[product(shape(IPvolume))]) V = reshape(IPvolume,[product(shape(IPvolume))])
do ce = 1, size(material_homogenizationMemberAt2,1) do ce = 1, size(material_homogenizationEntry,1)
do co = 1, homogenization_maxNconstituents do co = 1, homogenization_maxNconstituents
if (material_phaseAt2(co,ce) == ph) geom(ph)%V_0(material_phaseMemberAt2(co,ce)) = V(ce) if (material_phaseID(co,ce) == ph) geom(ph)%V_0(material_phaseEntry(co,ce)) = V(ce)
enddo enddo
enddo enddo

View File

@ -223,7 +223,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
Nmembers = count(material_phaseAt2 == ph) Nmembers = count(material_phaseID == ph)
sizeDotState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl & sizeDotState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl &
+ size(['xi_tw ','gamma_tw']) * prm%sum_N_tw + size(['xi_tw ','gamma_tw']) * prm%sum_N_tw
sizeState = sizeDotState sizeState = sizeDotState

View File

@ -89,7 +89,7 @@ module subroutine thermal_init(phases)
allocate(thermal_Nsources(phases%length),source = 0) allocate(thermal_Nsources(phases%length),source = 0)
do ph = 1, phases%length do ph = 1, phases%length
Nmembers = count(material_phaseAt2 == ph) Nmembers = count(material_phaseID == ph)
allocate(current(ph)%T(Nmembers),source=300.0_pReal) allocate(current(ph)%T(Nmembers),source=300.0_pReal)
allocate(current(ph)%dot_T(Nmembers),source=0.0_pReal) allocate(current(ph)%dot_T(Nmembers),source=0.0_pReal)
phase => phases%get(ph) phase => phases%get(ph)
@ -268,8 +268,8 @@ module subroutine phase_thermal_setField(T,dot_T, co,ce)
integer, intent(in) :: ce, co integer, intent(in) :: ce, co
current(material_phaseAt2(co,ce))%T(material_phaseMemberAt2(co,ce)) = T current(material_phaseID(co,ce))%T(material_phaseEntry(co,ce)) = T
current(material_phaseAt2(co,ce))%dot_T(material_phaseMemberAt2(co,ce)) = dot_T current(material_phaseID(co,ce))%dot_T(material_phaseEntry(co,ce)) = dot_T
end subroutine phase_thermal_setField end subroutine phase_thermal_setField

View File

@ -54,7 +54,7 @@ module function dissipation_init(source_length) result(mySources)
src => sources%get(so) src => sources%get(so)
prm%kappa = src%get_asFloat('kappa') prm%kappa = src%get_asFloat('kappa')
Nmembers = count(material_phaseAt2 == ph) Nmembers = count(material_phaseID == ph)
call phase_allocateState(thermalState(ph)%p(so),Nmembers,0,0,0) call phase_allocateState(thermalState(ph)%p(so),Nmembers,0,0,0)
end associate end associate

View File

@ -67,7 +67,7 @@ module function externalheat_init(source_length) result(mySources)
prm%f_T = src%get_as1dFloat('f_T',requiredSize = size(prm%t_n)) prm%f_T = src%get_as1dFloat('f_T',requiredSize = size(prm%t_n))
Nmembers = count(material_phaseAt2 == ph) Nmembers = count(material_phaseID == ph)
call phase_allocateState(thermalState(ph)%p(so),Nmembers,1,1,0) call phase_allocateState(thermalState(ph)%p(so),Nmembers,1,1,0)
end associate end associate
endif endif