Merge branch 'development' into fix-docstrings-for-sphinx

This commit is contained in:
Vitesh Shah 2020-09-26 17:56:42 +02:00
commit d0028b2eea
58 changed files with 1278 additions and 1183 deletions

View File

@ -344,6 +344,7 @@ Phenopowerlaw_singleSlip:
Pytest_grid:
stage: grid
script:
- module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel
- cd pytest
- pytest
except:

@ -1 +1 @@
Subproject commit 7d2aaef2341b83db53b58c04c4bdd88a64d47765
Subproject commit 68837540cab7435d8e2a06ae4c74e069e9386f35

View File

@ -1 +1 @@
v3.0.0-alpha-292-g14bfaa60c
v3.0.0-alpha-335-gdb8f6400f

View File

@ -52,16 +52,11 @@ parser.add_option('-q', '--quaternion',
type = 'string',
metavar='string',
help = 'name of the dataset containing pointwise/average orientation as quaternion [%default]')
parser.add_option('--homogenization',
dest = 'homogenization',
type = 'int', metavar = 'int',
help = 'homogenization index to be used [%default]')
parser.set_defaults(pointwise = 'CellData',
quaternion = 'Quats',
phase = 'Phases',
microstructure = 'FeatureIds',
homogenization = 1,
)
(options, filenames) = parser.parse_args()
@ -150,8 +145,7 @@ for name in filenames:
header = [scriptID + ' ' + ' '.join(sys.argv[1:])]\
+ config_header
geom = damask.Geom(microstructure,size,origin,
homogenization=options.homogenization,comments=header)
geom = damask.Geom(microstructure,size,origin,comments=header)
damask.util.croak(geom)
geom.save_ASCII(os.path.splitext(name)[0]+'.geom',compress=False)

View File

@ -52,10 +52,6 @@ parser.add_option('-p', '--periods',
dest = 'periods',
type = 'int', metavar = 'int',
help = 'number of repetitions of unit cell [%default]')
parser.add_option('--homogenization',
dest = 'homogenization',
type = 'int', metavar = 'int',
help = 'homogenization index to be used [%default]')
parser.add_option('--m',
dest = 'microstructure',
type = 'int', nargs = 2, metavar = 'int int',
@ -66,7 +62,6 @@ parser.set_defaults(type = minimal_surfaces[0],
periods = 1,
grid = (16,16,16),
size = (1.0,1.0,1.0),
homogenization = 1,
microstructure = (1,2),
)
@ -85,7 +80,6 @@ microstructure = np.where(options.threshold < surface[options.type](x,y,z),
options.microstructure[1],options.microstructure[0])
geom=damask.Geom(microstructure,options.size,
homogenization=options.homogenization,
comments=[scriptID + ' ' + ' '.join(sys.argv[1:])])
damask.util.croak(geom)

View File

@ -57,10 +57,6 @@ parser.add_option('-w', '--omega',
dest='omega',
type='float', metavar = 'float',
help='rotation angle around normal of osteon [%default]')
parser.add_option( '--homogenization',
dest='homogenization',
type='int', metavar = 'int',
help='homogenization index to be used [%default]')
parser.set_defaults(canal = 25e-6,
osteon = 100e-6,
@ -70,7 +66,7 @@ parser.set_defaults(canal = 25e-6,
amplitude = 60,
size = (300e-6,300e-6),
grid = (512,512),
homogenization = 1)
)
(options,filename) = parser.parse_args()
@ -139,7 +135,7 @@ header = [scriptID + ' ' + ' '.join(sys.argv[1:])]\
+ config_header
geom = damask.Geom(microstructure.reshape(grid),
size,-size/2,
homogenization=options.homogenization,comments=header)
comments=header)
damask.util.croak(geom)
geom.save_ASCII(sys.stdout if name is None else name,compress=False)

View File

@ -44,14 +44,9 @@ parser.add_option('--axes',
dest = 'axes',
type = 'string', nargs = 3, metavar = ' '.join(['string']*3),
help = 'orientation coordinate frame in terms of position coordinate frame [+x +y +z]')
parser.add_option('--homogenization',
dest = 'homogenization',
type = 'int', metavar = 'int',
help = 'homogenization index to be used [%default]')
parser.set_defaults(homogenization = 1,
pos = 'pos',
parser.set_defaults(pos = 'pos',
)
(options,filenames) = parser.parse_args()
@ -102,7 +97,7 @@ for name in filenames:
header = [scriptID + ' ' + ' '.join(sys.argv[1:])]\
+ config_header
geom = damask.Geom(microstructure,size,origin,
homogenization=options.homogenization,comments=header)
comments=header)
damask.util.croak(geom)
geom.save_ASCII(sys.stdout if name is None else os.path.splitext(name)[0]+'.geom',compress=False)

View File

@ -142,10 +142,6 @@ group.add_option('--without-config',
dest = 'config',
action = 'store_false',
help = 'omit material configuration header')
group.add_option('--homogenization',
dest = 'homogenization',
type = 'int', metavar = 'int',
help = 'homogenization index to be used [%default]')
group.add_option('--phase',
dest = 'phase',
type = 'int', metavar = 'int',
@ -157,7 +153,6 @@ parser.set_defaults(pos = 'pos',
weight = 'weight',
microstructure = 'microstructure',
eulers = 'euler',
homogenization = 1,
phase = 1,
cpus = 2,
laguerre = False,
@ -225,7 +220,7 @@ for name in filenames:
header = [scriptID + ' ' + ' '.join(sys.argv[1:])]\
+ config_header
geom = damask.Geom(indices.reshape(grid),size,origin,
homogenization=options.homogenization,comments=header)
comments=header)
damask.util.croak(geom)
geom.save_ASCII(sys.stdout if name is None else os.path.splitext(name)[0]+'.geom',compress=False)

View File

@ -41,7 +41,7 @@ parser.add_option('-N', '--iterations',
help = 'curvature flow iterations [%default]')
parser.add_option('-i', '--immutable',
action = 'extend', dest = 'immutable', metavar = '<int LIST>',
help = 'list of immutable microstructure indices')
help = 'list of immutable material indices')
parser.add_option('--ndimage',
dest = 'ndimage', action='store_true',
help = 'use ndimage.gaussian_filter in lieu of explicit FFT')
@ -64,15 +64,15 @@ for name in filenames:
geom = damask.Geom.load_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
grid_original = geom.get_grid()
grid_original = geom.grid
damask.util.croak(geom)
microstructure = np.tile(geom.microstructure,np.where(grid_original == 1, 2,1)) # make one copy along dimensions with grid == 1
grid = np.array(microstructure.shape)
material = np.tile(geom.material,np.where(grid_original == 1, 2,1)) # make one copy along dimensions with grid == 1
grid = np.array(material.shape)
# --- initialize support data ---------------------------------------------------------------------
# store a copy the initial microstructure to find locations of immutable indices
microstructure_original = np.copy(microstructure)
# store a copy of the initial material indices to find locations of immutable indices
material_original = np.copy(material)
if not options.ndimage:
X,Y,Z = np.mgrid[0:grid[0],0:grid[1],0:grid[2]]
@ -88,14 +88,14 @@ for name in filenames:
for smoothIter in range(options.N):
interfaceEnergy = np.zeros(microstructure.shape,dtype=np.float32)
interfaceEnergy = np.zeros(material.shape,dtype=np.float32)
for i in (-1,0,1):
for j in (-1,0,1):
for k in (-1,0,1):
# assign interfacial energy to all voxels that have a differing neighbor (in Moore neighborhood)
interfaceEnergy = np.maximum(interfaceEnergy,
getInterfaceEnergy(microstructure,np.roll(np.roll(np.roll(
microstructure,i,axis=0), j,axis=1), k,axis=2)))
getInterfaceEnergy(material,np.roll(np.roll(np.roll(
material,i,axis=0), j,axis=1), k,axis=2)))
# periodically extend interfacial energy array by half a grid size in positive and negative directions
periodic_interfaceEnergy = np.tile(interfaceEnergy,(3,3,3))[grid[0]//2:-grid[0]//2,
@ -129,13 +129,13 @@ for name in filenames:
iterations = int(round(options.d*2.))-1),# fat boundary
periodic_bulkEnergy[grid[0]//2:-grid[0]//2, # retain filled energy on fat boundary...
grid[1]//2:-grid[1]//2,
grid[2]//2:-grid[2]//2], # ...and zero everywhere else
grid[2]//2:-grid[2]//2], # ...and zero everywhere else
0.)).astype(np.complex64) *
gauss).astype(np.float32)
periodic_diffusedEnergy = np.tile(diffusedEnergy,(3,3,3))[grid[0]//2:-grid[0]//2,
grid[1]//2:-grid[1]//2,
grid[2]//2:-grid[2]//2] # periodically extend the smoothed bulk energy
grid[2]//2:-grid[2]//2] # periodically extend the smoothed bulk energy
# transform voxels close to interface region
@ -143,33 +143,35 @@ for name in filenames:
return_distances = False,
return_indices = True) # want index of closest bulk grain
periodic_microstructure = np.tile(microstructure,(3,3,3))[grid[0]//2:-grid[0]//2,
grid[1]//2:-grid[1]//2,
grid[2]//2:-grid[2]//2] # periodically extend the microstructure
periodic_material = np.tile(material,(3,3,3))[grid[0]//2:-grid[0]//2,
grid[1]//2:-grid[1]//2,
grid[2]//2:-grid[2]//2] # periodically extend the geometry
microstructure = periodic_microstructure[index[0],
index[1],
index[2]].reshape(2*grid)[grid[0]//2:-grid[0]//2,
grid[1]//2:-grid[1]//2,
grid[2]//2:-grid[2]//2] # extent grains into interface region
material = periodic_material[index[0],
index[1],
index[2]].reshape(2*grid)[grid[0]//2:-grid[0]//2,
grid[1]//2:-grid[1]//2,
grid[2]//2:-grid[2]//2] # extent grains into interface region
# replace immutable microstructures with closest mutable ones
index = ndimage.morphology.distance_transform_edt(np.in1d(microstructure,options.immutable).reshape(grid),
# replace immutable materials with closest mutable ones
index = ndimage.morphology.distance_transform_edt(np.in1d(material,options.immutable).reshape(grid),
return_distances = False,
return_indices = True)
microstructure = microstructure[index[0],
index[1],
index[2]]
material = material[index[0],
index[1],
index[2]]
immutable = np.zeros(microstructure.shape, dtype=np.bool)
# find locations where immutable microstructures have been in original structure
immutable = np.zeros(material.shape, dtype=np.bool)
# find locations where immutable materials have been in original structure
for micro in options.immutable:
immutable += microstructure_original == micro
immutable += material_original == micro
# undo any changes involving immutable microstructures
microstructure = np.where(immutable, microstructure_original,microstructure)
# undo any changes involving immutable materials
material = np.where(immutable, material_original,material)
geom = geom.duplicate(microstructure[0:grid_original[0],0:grid_original[1],0:grid_original[2]])
geom.add_comments(scriptID + ' ' + ' '.join(sys.argv[1:]))
geom.save_ASCII(sys.stdout if name is None else name,compress=False)
damask.Geom(material = material[0:grid_original[0],0:grid_original[1],0:grid_original[2]],
size = geom.size,
origin = geom.origin,
comments = geom.comments + [scriptID + ' ' + ' '.join(sys.argv[1:])],
)\
.save_ASCII(sys.stdout if name is None else name,compress=False)

View File

@ -100,7 +100,7 @@ def mesh(r,d):
#-------------------------------------------------------------------------------------------------
def material():
def materials():
return [\
"*new_mater standard",
"*mater_option general:state:solid",
@ -130,10 +130,10 @@ def geometry():
#-------------------------------------------------------------------------------------------------
def initial_conditions(microstructures):
def initial_conditions(material):
elements = []
element = 0
for id in microstructures:
for id in material:
element += 1
if len(elements) < id:
for i in range(id-len(elements)):
@ -153,7 +153,7 @@ def initial_conditions(microstructures):
for grain,elementList in enumerate(elements):
cmds.append([\
"*new_icond",
"*icond_name microstructure_%i"%(grain+1),
"*icond_name material_%i"%(grain+1),
"*icond_type state_variable",
"*icond_param_value state_var_id 2",
"*icond_dof_value var %i"%(grain+1),
@ -197,14 +197,14 @@ for name in filenames:
damask.util.report(scriptName,name)
geom = damask.Geom.load_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
microstructure = geom.get_microstructure().flatten(order='F')
material = geom.material.flatten(order='F')
cmds = [\
init(),
mesh(geom.grid,geom.size),
material(),
materials(),
geometry(),
initial_conditions(microstructure),
initial_conditions(material),
'*identify_sets',
'*show_model',
'*redraw',

View File

@ -30,7 +30,7 @@ class myThread (threading.Thread):
def run(self):
global bestSeedsUpdate
global bestSeedsVFile
global nMicrostructures
global nMaterials
global delta
global points
global target
@ -70,7 +70,7 @@ class myThread (threading.Thread):
selectedMs = []
direction = []
for i in range(NmoveGrains):
selectedMs.append(random.randrange(1,nMicrostructures))
selectedMs.append(random.randrange(1,nMaterials))
direction.append((np.random.random()-0.5)*delta)
@ -101,11 +101,11 @@ class myThread (threading.Thread):
#--- evaluate current seeds file ------------------------------------------------------------------
perturbedGeom = damask.Geom.load_ASCII(perturbedGeomVFile)
myNmicrostructures = len(np.unique(perturbedGeom.microstructure))
currentData=np.bincount(perturbedGeom.microstructure.ravel())[1:]/points
myNmaterials = len(np.unique(perturbedGeom.material))
currentData = np.bincount(perturbedGeom.material.ravel())[1:]/points
currentError=[]
currentHist=[]
for i in range(nMicrostructures): # calculate the deviation in all bins per histogram
for i in range(nMaterials): # calculate the deviation in all bins per histogram
currentHist.append(np.histogram(currentData,bins=target[i]['bins'])[0])
currentError.append(np.sqrt(np.square(np.array(target[i]['histogram']-currentHist[i])).sum()))
@ -117,12 +117,12 @@ class myThread (threading.Thread):
bestMatch = match
#--- count bin classes with no mismatch ----------------------------------------------------------------------
myMatch=0
for i in range(nMicrostructures):
for i in range(nMaterials):
if currentError[i] > 0.0: break
myMatch = i+1
if myNmicrostructures == nMicrostructures:
for i in range(min(nMicrostructures,myMatch+options.bins)):
if myNmaterials == nMaterials:
for i in range(min(nMaterials,myMatch+options.bins)):
if currentError[i] > target[i]['error']: # worse fitting, next try
randReset = True
break
@ -141,25 +141,25 @@ class myThread (threading.Thread):
for line in perturbedSeedsVFile:
currentSeedsFile.write(line)
bestSeedsVFile.write(line)
for j in range(nMicrostructures): # save new errors for all bins
for j in range(nMaterials): # save new errors for all bins
target[j]['error'] = currentError[j]
if myMatch > match: # one or more new bins have no deviation
damask.util.croak( 'Stage {:d} cleared'.format(myMatch))
match=myMatch
sys.stdout.flush()
break
if i == min(nMicrostructures,myMatch+options.bins)-1: # same quality as before: take it to keep on moving
if i == min(nMaterials,myMatch+options.bins)-1: # same quality as before: take it to keep on moving
bestSeedsUpdate = time.time()
perturbedSeedsVFile.seek(0)
bestSeedsVFile.close()
bestSeedsVFile = StringIO()
bestSeedsVFile.writelines(perturbedSeedsVFile.readlines())
for j in range(nMicrostructures):
for j in range(nMaterials):
target[j]['error'] = currentError[j]
randReset = True
else: #--- not all grains are tessellated
damask.util.croak('Thread {:d}: Microstructure mismatch ({:d} microstructures mapped)'\
.format(self.threadID,myNmicrostructures))
damask.util.croak('Thread {:d}: Material mismatch ({:d} material indices mapped)'\
.format(self.threadID,myNmaterials))
randReset = True
@ -217,10 +217,10 @@ points = np.array(options.grid).prod().astype('float')
# ----------- calculate target distribution and bin edges
targetGeom = damask.Geom.load_ASCII(os.path.splitext(os.path.basename(options.target))[0]+'.geom')
nMicrostructures = len(np.unique(targetGeom.microstructure))
targetVolFrac = np.bincount(targetGeom.microstructure.flatten())/targetGeom.grid.prod().astype(np.float)
nMaterials = len(np.unique(targetGeom.material))
targetVolFrac = np.bincount(targetGeom.material.flatten())/targetGeom.grid.prod().astype(np.float)
target = []
for i in range(1,nMicrostructures+1):
for i in range(1,nMaterials+1):
targetHist,targetBins = np.histogram(targetVolFrac,bins=i) #bin boundaries
target.append({'histogram':targetHist,'bins':targetBins})
@ -233,7 +233,7 @@ else:
bestSeedsVFile.write(damask.util.execute('seeds_fromRandom'+\
' -g '+' '.join(list(map(str, options.grid)))+\
' -r {:d}'.format(options.randomSeed)+\
' -N '+str(nMicrostructures))[0])
' -N '+str(nMaterials))[0])
bestSeedsUpdate = time.time()
# ----------- tessellate initial seed file to get and evaluate geom file
@ -244,11 +244,11 @@ initialGeomVFile.write(damask.util.execute('geom_fromVoronoiTessellation '+
initialGeomVFile.seek(0)
initialGeom = damask.Geom.load_ASCII(initialGeomVFile)
if len(np.unique(targetGeom.microstructure)) != nMicrostructures:
damask.util.croak('error. Microstructure count mismatch')
if len(np.unique(targetGeom.material)) != nMaterials:
damask.util.croak('error. Material count mismatch')
initialData = np.bincount(initialGeom.microstructure.flatten())/points
for i in range(nMicrostructures):
initialData = np.bincount(initialGeom.material.flatten())/points
for i in range(nMaterials):
initialHist = np.histogram(initialData,bins=target[i]['bins'])[0]
target[i]['error']=np.sqrt(np.square(np.array(target[i]['histogram']-initialHist)).sum())
@ -257,13 +257,13 @@ if target[0]['error'] > 0.0:
target[0]['error'] *=((target[0]['bins'][0]-np.min(initialData))**2.0+
(target[0]['bins'][1]-np.max(initialData))**2.0)**0.5
match=0
for i in range(nMicrostructures):
for i in range(nMaterials):
if target[i]['error'] > 0.0: break
match = i+1
if options.maxseeds < 1:
maxSeeds = len(np.unique(initialGeom.microstructure))
maxSeeds = len(np.unique(initialGeom.material))
else:
maxSeeds = options.maxseeds

View File

@ -17,7 +17,7 @@ scriptID = ' '.join([scriptName,damask.version])
#--------------------------------------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
Create seed file taking microstructure indices from given geom file.
Create seed file taking material indices from given geom file.
Indices can be black-listed or white-listed.
""", version = scriptID)
@ -47,11 +47,11 @@ for name in filenames:
damask.util.report(scriptName,name)
geom = damask.Geom.load_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
microstructure = geom.get_microstructure().reshape((-1,1),order='F')
material = geom.material.reshape((-1,1),order='F')
mask = np.logical_and(np.in1d(microstructure,options.whitelist,invert=False) if options.whitelist else \
mask = np.logical_and(np.in1d(material,options.whitelist,invert=False) if options.whitelist else \
np.full(geom.grid.prod(),True,dtype=bool),
np.in1d(microstructure,options.blacklist,invert=True) if options.blacklist else \
np.in1d(material,options.blacklist,invert=True) if options.blacklist else \
np.full(geom.grid.prod(),True,dtype=bool))
seeds = damask.grid_filters.cell_coord0(geom.grid,geom.size).reshape(-1,3,order='F')
@ -61,8 +61,8 @@ for name in filenames:
'grid\ta {}\tb {}\tc {}'.format(*geom.grid),
'size\tx {}\ty {}\tz {}'.format(*geom.size),
'origin\tx {}\ty {}\tz {}'.format(*geom.origin),
'homogenization\t{}'.format(geom.homogenization)]
]
damask.Table(seeds[mask],{'pos':(3,)},comments)\
.add('microstructure',microstructure[mask].astype(int))\
.add('material',material[mask].astype(int))\
.save(sys.stdout if name is None else os.path.splitext(name)[0]+'.seeds',legacy=True)

View File

@ -76,7 +76,7 @@ for name in filenames:
g[2] = k + offset[2]
g %= geom.grid
seeds[n,0:3] = (g+0.5)/geom.grid # normalize coordinates to box
seeds[n, 3] = geom.microstructure[g[0],g[1],g[2]]
seeds[n, 3] = geom.material[g[0],g[1],g[2]]
if options.x: g[0] += 1
if options.y: g[1] += 1
n += 1
@ -88,9 +88,9 @@ for name in filenames:
'grid\ta {}\tb {}\tc {}'.format(*geom.grid),
'size\tx {}\ty {}\tz {}'.format(*geom.size),
'origin\tx {}\ty {}\tz {}'.format(*geom.origin),
'homogenization\t{}'.format(geom.homogenization)]
]
table = damask.Table(seeds,{'pos':(3,),'microstructure':(1,)},comments)
table.set('microstructure',table.get('microstructure').astype(np.int))\
table = damask.Table(seeds,{'pos':(3,),'material':(1,)},comments)
table.set('material',table.get('material').astype(np.int))\
.save(sys.stdout if name is None else \
os.path.splitext(name)[0]+f'_poked_{options.N}.seeds',legacy=True)

View File

@ -20,6 +20,12 @@ from ._result import Result # noqa
from ._geom import Geom # noqa
from ._material import Material # noqa
from . import solver # noqa
from . import util # noqa
from . import seeds # noqa
from . import grid_filters # noqa
from . import mechanics # noqa
# deprecated
Environment = _

View File

@ -15,39 +15,54 @@ from . import grid_filters
class Geom:
"""Geometry definition for grid solvers."""
def __init__(self,microstructure,size,origin=[0.0,0.0,0.0],homogenization=1,comments=[]):
def __init__(self,material,size,origin=[0.0,0.0,0.0],comments=[]):
"""
New geometry definition from array of microstructures and size.
New geometry definition from array of material, size, and origin.
Parameters
----------
microstructure : numpy.ndarray
Microstructure array (3D)
material : numpy.ndarray
Material index array (3D).
size : list or numpy.ndarray
Physical size of the microstructure in meter.
Physical size of the geometry in meter.
origin : list or numpy.ndarray, optional
Physical origin of the microstructure in meter.
homogenization : int, optional
Homogenization index.
Physical origin of the geometry in meter.
comments : list of str, optional
Comment lines.
"""
self.set_microstructure(microstructure)
self.set_size(size)
self.set_origin(origin)
self.set_homogenization(homogenization)
self.set_comments(comments)
if len(material.shape) != 3:
raise ValueError(f'Invalid material shape {material.shape}.')
elif material.dtype not in np.sctypes['float'] + np.sctypes['int']:
raise TypeError(f'Invalid material data type {material.dtype}.')
else:
self.material = np.copy(material)
if self.material.dtype in np.sctypes['float'] and \
np.all(self.material == self.material.astype(int).astype(float)):
self.material = self.material.astype(int)
if len(size) != 3 or any(np.array(size) <= 0):
raise ValueError(f'Invalid size {size}.')
else:
self.size = np.array(size)
if len(origin) != 3:
raise ValueError(f'Invalid origin {origin}.')
else:
self.origin = np.array(origin)
self.comments = [str(c) for c in comments] if isinstance(comments,list) else [str(comments)]
def __repr__(self):
"""Basic information on geometry definition."""
return util.srepr([
f'grid a b c: {util.srepr(self.get_grid ()," x ")}',
f'size x y z: {util.srepr(self.get_size ()," x ")}',
f'origin x y z: {util.srepr(self.get_origin()," ")}',
f'# materialpoints: {self.N_microstructure}',
f'max materialpoint: {np.nanmax(self.microstructure)}',
f'grid a b c: {util.srepr(self.grid, " x ")}',
f'size x y z: {util.srepr(self.size, " x ")}',
f'origin x y z: {util.srepr(self.origin," ")}',
f'# materials: {self.N_materials}',
f'max material: {np.nanmax(self.material)}',
])
@ -61,43 +76,6 @@ class Geom:
return self.__copy__()
def duplicate(self,microstructure=None,size=None,origin=None,comments=None,autosize=False):
"""
Create a duplicate having updated microstructure, size, and origin.
Parameters
----------
microstructure : numpy.ndarray, optional
Microstructure array (3D).
size : list or numpy.ndarray, optional
Physical size of the microstructure in meter.
origin : list or numpy.ndarray, optional
Physical origin of the microstructure in meter.
comments : list of str, optional
Comment lines.
autosize : bool, optional
Ignore size parameter and rescale according to change of grid points.
"""
if size is not None and autosize:
raise ValueError('Auto-sizing conflicts with explicit size parameter.')
grid_old = self.get_grid()
dup = self.copy()
dup.set_microstructure(microstructure)
dup.set_origin(origin)
if comments is not None:
dup.set_comments(comments)
if size is not None:
dup.set_size(size)
elif autosize:
dup.set_size(dup.get_grid()/grid_old*self.get_size())
return dup
def diff(self,other):
"""
Report property differences of self relative to other.
@ -109,176 +87,37 @@ class Geom:
"""
message = []
if np.any(other.get_grid() != self.get_grid()):
message.append(util.delete(f'grid a b c: {util.srepr(other.get_grid()," x ")}'))
message.append(util.emph( f'grid a b c: {util.srepr( self.get_grid()," x ")}'))
if np.any(other.grid != self.grid):
message.append(util.delete(f'grid a b c: {util.srepr(other.grid," x ")}'))
message.append(util.emph( f'grid a b c: {util.srepr( self.grid," x ")}'))
if np.any(other.get_size() != self.get_size()):
message.append(util.delete(f'size x y z: {util.srepr(other.get_size()," x ")}'))
message.append(util.emph( f'size x y z: {util.srepr( self.get_size()," x ")}'))
if not np.allclose(other.size,self.size):
message.append(util.delete(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 np.any(other.get_origin() != self.get_origin()):
message.append(util.delete(f'origin x y z: {util.srepr(other.get_origin()," ")}'))
message.append(util.emph( f'origin x y z: {util.srepr( self.get_origin()," ")}'))
if not np.allclose(other.origin,self.origin):
message.append(util.delete(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_microstructure != self.N_microstructure:
message.append(util.delete(f'# materialpoints: {other.N_microstructure}'))
message.append(util.emph( f'# materialpoints: { self.N_microstructure}'))
if other.N_materials != self.N_materials:
message.append(util.delete(f'# materials: {other.N_materials}'))
message.append(util.emph( f'# materials: { self.N_materials}'))
if np.nanmax(other.microstructure) != np.nanmax(self.microstructure):
message.append(util.delete(f'max materialpoint: {np.nanmax(other.microstructure)}'))
message.append(util.emph( f'max materialpoint: {np.nanmax( self.microstructure)}'))
if np.nanmax(other.material) != np.nanmax(self.material):
message.append(util.delete(f'max material: {np.nanmax(other.material)}'))
message.append(util.emph( f'max material: {np.nanmax( self.material)}'))
return util.return_message(message)
def set_comments(self,comments):
"""
Replace all existing comments.
Parameters
----------
comments : list of str
All comments.
"""
self.comments = []
self.add_comments(comments)
def add_comments(self,comments):
"""
Append comments to existing comments.
Parameters
----------
comments : list of str
New comments.
"""
self.comments += [str(c) for c in comments] if isinstance(comments,list) else [str(comments)]
def set_microstructure(self,microstructure):
"""
Replace the existing microstructure representation.
The complete microstructure is replaced (indcluding grid definition),
unless a masked array is provided in which case the grid dimensions
need to match and masked entries are not replaced.
Parameters
----------
microstructure : numpy.ndarray or numpy.ma.core.MaskedArray of shape (:,:,:)
Microstructure indices.
"""
if microstructure is not None:
if isinstance(microstructure,np.ma.core.MaskedArray):
self.microstructure = np.where(microstructure.mask,
self.microstructure,microstructure.data)
else:
self.microstructure = np.copy(microstructure)
if self.microstructure.dtype in np.sctypes['float'] and \
np.all(self.microstructure == self.microstructure.astype(int).astype(float)):
self.microstructure = self.microstructure.astype(int)
if len(self.microstructure.shape) != 3:
raise ValueError(f'Invalid microstructure shape {microstructure.shape}')
elif self.microstructure.dtype not in np.sctypes['float'] + np.sctypes['int']:
raise TypeError(f'Invalid microstructure data type {microstructure.dtype}')
def set_size(self,size):
"""
Replace the existing size information.
Parameters
----------
size : list or numpy.ndarray
Physical size of the microstructure in meter.
"""
if size is not None:
if len(size) != 3 or any(np.array(size) <= 0):
raise ValueError(f'Invalid size {size}')
else:
self.size = np.array(size)
def set_origin(self,origin):
"""
Replace the existing origin information.
Parameters
----------
origin : list or numpy.ndarray
Physical origin of the microstructure in meter.
"""
if origin is not None:
if len(origin) != 3:
raise ValueError(f'Invalid origin {origin}')
else:
self.origin = np.array(origin)
def set_homogenization(self,homogenization):
"""
Replace the existing homogenization index.
Parameters
----------
homogenization : int
Homogenization index.
"""
if homogenization is not None:
if not isinstance(homogenization,int) or homogenization < 1:
raise TypeError(f'Invalid homogenization {homogenization}.')
else:
self.homogenization = homogenization
@property
def grid(self):
return self.get_grid()
return np.asarray(self.material.shape)
@property
def N_microstructure(self):
return np.unique(self.microstructure).size
def get_microstructure(self):
"""Return the microstructure representation."""
return np.copy(self.microstructure)
def get_size(self):
"""Return the physical size in meter."""
return np.copy(self.size)
def get_origin(self):
"""Return the origin in meter."""
return np.copy(self.origin)
def get_grid(self):
"""Return the grid discretization."""
return np.asarray(self.microstructure.shape)
def get_homogenization(self):
"""Return the homogenization index."""
return self.homogenization
def get_comments(self):
"""Return the comments."""
return self.comments[:]
def N_materials(self):
return np.unique(self.material).size
@staticmethod
@ -318,12 +157,10 @@ class Geom:
size = np.array([float(dict(zip(items[1::2],items[2::2]))[i]) for i in ['x','y','z']])
elif key == 'origin':
origin = np.array([float(dict(zip(items[1::2],items[2::2]))[i]) for i in ['x','y','z']])
elif key == 'homogenization':
homogenization = int(items[1])
else:
comments.append(line.strip())
microstructure = np.empty(grid.prod()) # initialize as flat array
material = np.empty(grid.prod()) # initialize as flat array
i = 0
for line in content[header_length:]:
items = line.split('#')[0].split()
@ -335,16 +172,16 @@ class Geom:
abs(int(items[2])-int(items[0]))+1,dtype=float)
else: items = list(map(float,items))
else: items = list(map(float,items))
microstructure[i:i+len(items)] = items
material[i:i+len(items)] = items
i += len(items)
if i != grid.prod():
raise TypeError(f'Invalid file: expected {grid.prod()} entries, found {i}')
if not np.any(np.mod(microstructure,1) != 0.0): # no float present
microstructure = microstructure.astype('int')
if not np.any(np.mod(material,1) != 0.0): # no float present
material = material.astype('int')
return Geom(microstructure.reshape(grid,order='F'),size,origin,homogenization,comments)
return Geom(material.reshape(grid,order='F'),size,origin,comments)
@staticmethod
@ -363,9 +200,11 @@ class Geom:
comments = v.get_comments()
grid = np.array(v.vtk_data.GetDimensions())-1
bbox = np.array(v.vtk_data.GetBounds()).reshape(3,2).T
size = bbox[1] - bbox[0]
return Geom(v.get('materialpoint').reshape(grid,order='F'),size,bbox[0],comments=comments)
return Geom(material = v.get('material').reshape(grid,order='F'),
size = bbox[1] - bbox[0],
origin = bbox[0],
comments=comments)
@staticmethod
@ -373,7 +212,7 @@ class Geom:
return np.argmin(np.sum((np.broadcast_to(point,(len(seeds),3))-seeds)**2,axis=1) - weights)
@staticmethod
def from_Laguerre_tessellation(grid,size,seeds,weights,periodic=True):
def from_Laguerre_tessellation(grid,size,seeds,weights,material=None,periodic=True):
"""
Generate geometry from Laguerre tessellation.
@ -382,11 +221,14 @@ class Geom:
grid : int numpy.ndarray of shape (3)
Number of grid points in x,y,z direction.
size : list or numpy.ndarray of shape (3)
Physical size of the microstructure in meter.
Physical size of the geometry in meter.
seeds : numpy.ndarray of shape (:,3)
Position of the seed points in meter. All points need to lay within the box.
weights : numpy.ndarray of shape (seeds.shape[0])
Weights of the seeds. Setting all weights to 1.0 gives a standard Voronoi tessellation.
material : numpy.ndarray of shape (seeds.shape[0]), optional
Material ID of the seeds. Defaults to None, in which case materials are
consecutively numbered.
periodic : Boolean, optional
Perform a periodic tessellation. Defaults to True.
@ -406,21 +248,22 @@ class Geom:
result = pool.map_async(partial(Geom._find_closest_seed,seeds_p,weights_p), [coord for coord in coords])
pool.close()
pool.join()
microstructure = np.array(result.get())
material_ = np.array(result.get())
if periodic:
microstructure = microstructure.reshape(grid*3)
microstructure = microstructure[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0]
material_ = material_.reshape(grid*3)
material_ = material_[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0]
else:
microstructure = microstructure.reshape(grid)
material_ = material_.reshape(grid)
return Geom(microstructure+1,size,homogenization=1,
comments=util.execution_stamp('Geom','from_Laguerre_tessellation'),
return Geom(material = material_+1 if material is None else material[material_],
size = size,
comments = util.execution_stamp('Geom','from_Laguerre_tessellation'),
)
@staticmethod
def from_Voronoi_tessellation(grid,size,seeds,periodic=True):
def from_Voronoi_tessellation(grid,size,seeds,material=None,periodic=True):
"""
Generate geometry from Voronoi tessellation.
@ -429,19 +272,23 @@ class Geom:
grid : int numpy.ndarray of shape (3)
Number of grid points in x,y,z direction.
size : list or numpy.ndarray of shape (3)
Physical size of the microstructure in meter.
Physical size of the geometry in meter.
seeds : numpy.ndarray of shape (:,3)
Position of the seed points in meter. All points need to lay within the box.
material : numpy.ndarray of shape (seeds.shape[0]), optional
Material ID of the seeds. Defaults to None, in which case materials are
consecutively numbered.
periodic : Boolean, optional
Perform a periodic tessellation. Defaults to True.
"""
coords = grid_filters.cell_coord0(grid,size).reshape(-1,3)
KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds)
devNull,microstructure = KDTree.query(coords)
devNull,material_ = KDTree.query(coords)
return Geom(microstructure.reshape(grid)+1,size,homogenization=1,
comments=util.execution_stamp('Geom','from_Voronoi_tessellation'),
return Geom(material = (material_+1 if material is None else material[material_]).reshape(grid),
size = size,
comments = util.execution_stamp('Geom','from_Voronoi_tessellation'),
)
@ -457,24 +304,25 @@ class Geom:
Compress geometry with 'x of y' and 'a to b'.
"""
header = [f'{len(self.comments)+4} header'] + self.comments
header.append('grid a {} b {} c {}'.format(*self.get_grid()))
header.append('size x {} y {} z {}'.format(*self.get_size()))
header.append('origin x {} y {} z {}'.format(*self.get_origin()))
header.append(f'homogenization {self.get_homogenization()}')
header = [f'{len(self.comments)+4} header'] + self.comments \
+ ['grid a {} b {} c {}'.format(*self.grid),
'size x {} y {} z {}'.format(*self.size),
'origin x {} y {} z {}'.format(*self.origin),
'homogenization 1',
]
grid = self.get_grid()
grid = self.grid
if compress is None:
plain = grid.prod()/self.N_microstructure < 250
plain = grid.prod()/self.N_materials < 250
else:
plain = not compress
if plain:
format_string = '%g' if self.microstructure.dtype in np.sctypes['float'] else \
'%{}i'.format(1+int(np.floor(np.log10(np.nanmax(self.microstructure)))))
format_string = '%g' if self.material.dtype in np.sctypes['float'] else \
'%{}i'.format(1+int(np.floor(np.log10(np.nanmax(self.material)))))
np.savetxt(fname,
self.microstructure.reshape([grid[0],np.prod(grid[1:])],order='F').T,
self.material.reshape([grid[0],np.prod(grid[1:])],order='F').T,
header='\n'.join(header), fmt=format_string, comments='')
else:
try:
@ -485,7 +333,7 @@ class Geom:
compressType = None
former = start = -1
reps = 0
for current in self.microstructure.flatten('F'):
for current in self.material.flatten('F'):
if abs(current - former) == 1 and (start - current) == reps*(former - current):
compressType = 'to'
reps += 1
@ -530,7 +378,7 @@ class Geom:
"""
v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin)
v.add(self.microstructure.flatten(order='F'),'materialpoint')
v.add(self.material.flatten(order='F'),'material')
v.add_comments(self.comments)
v.save(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr',parallel=False,compress=compress)
@ -562,11 +410,11 @@ class Geom:
0 gives octahedron (ǀxǀ^(2^0) + ǀyǀ^(2^0) + ǀzǀ^(2^0) < 1)
1 gives a sphere (ǀxǀ^(2^1) + ǀyǀ^(2^1) + ǀzǀ^(2^1) < 1)
fill : int, optional
Fill value for primitive. Defaults to microstructure.max() + 1.
Fill value for primitive. Defaults to material.max() + 1.
R : damask.Rotation, optional
Rotation of primitive. Defaults to no rotation.
inverse : Boolean, optional
Retain original microstructure within primitive and fill outside.
Retain original materials within primitive and fill outside.
Defaults to False.
periodic : Boolean, optional
Repeat primitive over boundaries. Defaults to True.
@ -588,22 +436,23 @@ class Geom:
if periodic: # translate back to center
mask = np.roll(mask,((c-np.ones(3)*.5)*self.grid).astype(int),(0,1,2))
fill_ = np.full_like(self.microstructure,np.nanmax(self.microstructure)+1 if fill is None else fill)
ms = np.ma.MaskedArray(fill_,np.logical_not(mask) if inverse else mask)
fill_ = np.full_like(self.material,np.nanmax(self.material)+1 if fill is None else fill)
return self.duplicate(ms,
comments=self.get_comments()+[util.execution_stamp('Geom','add_primitive')],
)
return Geom(material = np.where(np.logical_not(mask) if inverse else mask, self.material,fill_),
size = self.size,
origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','add_primitive')],
)
def mirror(self,directions,reflect=False):
"""
Mirror microstructure along given directions.
Mirror geometry along given directions.
Parameters
----------
directions : iterable containing str
Direction(s) along which the microstructure is mirrored.
Direction(s) along which the geometry is mirrored.
Valid entries are 'x', 'y', 'z'.
reflect : bool, optional
Reflect (include) outermost layers. Defaults to False.
@ -614,28 +463,30 @@ class Geom:
raise ValueError(f'Invalid direction {set(directions).difference(valid)} specified.')
limits = [None,None] if reflect else [-2,0]
ms = self.get_microstructure()
mat = self.material.copy()
if 'z' in directions:
ms = np.concatenate([ms,ms[:,:,limits[0]:limits[1]:-1]],2)
if 'y' in directions:
ms = np.concatenate([ms,ms[:,limits[0]:limits[1]:-1,:]],1)
if 'x' in directions:
ms = np.concatenate([ms,ms[limits[0]:limits[1]:-1,:,:]],0)
mat = np.concatenate([mat,mat[limits[0]:limits[1]:-1,:,:]],0)
if 'y' in directions:
mat = np.concatenate([mat,mat[:,limits[0]:limits[1]:-1,:]],1)
if 'z' in directions:
mat = np.concatenate([mat,mat[:,:,limits[0]:limits[1]:-1]],2)
return self.duplicate(ms,
comments=self.get_comments()+[util.execution_stamp('Geom','mirror')],
autosize=True)
return Geom(material = mat,
size = self.size/self.grid*np.asarray(mat.shape),
origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','mirror')],
)
def flip(self,directions):
"""
Flip microstructure along given directions.
Flip geometry along given directions.
Parameters
----------
directions : iterable containing str
Direction(s) along which the microstructure is flipped.
Direction(s) along which the geometry is flipped.
Valid entries are 'x', 'y', 'z'.
"""
@ -643,16 +494,18 @@ class Geom:
if not set(directions).issubset(valid):
raise ValueError(f'Invalid direction {set(directions).difference(valid)} specified.')
ms = np.flip(self.microstructure, (valid.index(d) for d in directions if d in valid))
mat = np.flip(self.material, (valid.index(d) for d in directions if d in valid))
return self.duplicate(ms,
comments=self.get_comments()+[util.execution_stamp('Geom','flip')],
)
return Geom(material = mat,
size = self.size,
origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','flip')],
)
def scale(self,grid,periodic=True):
"""
Scale microstructure to new grid.
Scale geometry to new grid.
Parameters
----------
@ -662,21 +515,23 @@ class Geom:
Assume geometry to be periodic. Defaults to True.
"""
return self.duplicate(ndimage.interpolation.zoom(
self.microstructure,
grid/self.get_grid(),
output=self.microstructure.dtype,
order=0,
mode=('wrap' if periodic else 'nearest'),
prefilter=False
),
comments=self.get_comments()+[util.execution_stamp('Geom','scale')],
)
return Geom(material = ndimage.interpolation.zoom(
self.material,
grid/self.grid,
output=self.material.dtype,
order=0,
mode=('wrap' if periodic else 'nearest'),
prefilter=False
),
size = self.size,
origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','scale')],
)
def clean(self,stencil=3,selection=None,periodic=True):
"""
Smooth microstructure by selecting most frequent index within given stencil at each location.
Smooth geometry by selecting most frequent material index within given stencil at each location.
Parameters
----------
@ -696,83 +551,87 @@ class Geom:
else:
return me
return self.duplicate(ndimage.filters.generic_filter(
self.microstructure,
mostFrequent,
size=(stencil if selection is None else stencil//2*2+1,)*3,
mode=('wrap' if periodic else 'nearest'),
extra_keywords=dict(selection=selection),
).astype(self.microstructure.dtype),
comments=self.get_comments()+[util.execution_stamp('Geom','clean')],
)
return Geom(material = ndimage.filters.generic_filter(
self.material,
mostFrequent,
size=(stencil if selection is None else stencil//2*2+1,)*3,
mode=('wrap' if periodic else 'nearest'),
extra_keywords=dict(selection=selection),
).astype(self.material.dtype),
size = self.size,
origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','clean')],
)
def renumber(self):
"""Renumber sorted microstructure indices to 1,...,N."""
renumbered = np.empty(self.get_grid(),dtype=self.microstructure.dtype)
for i, oldID in enumerate(np.unique(self.microstructure)):
renumbered = np.where(self.microstructure == oldID, i+1, renumbered)
"""Renumber sorted material indices to 1,...,N."""
renumbered = np.empty(self.grid,dtype=self.material.dtype)
for i, oldID in enumerate(np.unique(self.material)):
renumbered = np.where(self.material == oldID, i+1, renumbered)
return self.duplicate(renumbered,
comments=self.get_comments()+[util.execution_stamp('Geom','renumber')],
)
return Geom(material = renumbered,
size = self.size,
origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','renumber')],
)
def rotate(self,R,fill=None):
"""
Rotate microstructure (pad if required).
Rotate geometry (pad if required).
Parameters
----------
R : damask.Rotation
Rotation to apply to the microstructure.
Rotation to apply to the geometry.
fill : int or float, optional
Microstructure index to fill the corners. Defaults to microstructure.max() + 1.
Material index to fill the corners. Defaults to material.max() + 1.
"""
if fill is None: fill = np.nanmax(self.microstructure) + 1
dtype = float if np.isnan(fill) or int(fill) != fill or self.microstructure.dtype==np.float else int
if fill is None: fill = np.nanmax(self.material) + 1
dtype = float if np.isnan(fill) or int(fill) != fill or self.material.dtype==np.float else int
Eulers = R.as_Eulers(degrees=True)
microstructure_in = self.get_microstructure()
material_in = self.material.copy()
# These rotations are always applied in the reference coordinate system, i.e. (z,x,z) not (z,x',z'')
# see https://www.cs.utexas.edu/~theshark/courses/cs354/lectures/cs354-14.pdf
for angle,axes in zip(Eulers[::-1], [(0,1),(1,2),(0,1)]):
microstructure_out = ndimage.rotate(microstructure_in,angle,axes,order=0,
material_out = ndimage.rotate(material_in,angle,axes,order=0,
prefilter=False,output=dtype,cval=fill)
if np.prod(microstructure_in.shape) == np.prod(microstructure_out.shape):
if np.prod(material_in.shape) == np.prod(material_out.shape):
# avoid scipy interpolation errors for rotations close to multiples of 90°
microstructure_in = np.rot90(microstructure_in,k=np.rint(angle/90.).astype(int),axes=axes)
material_in = np.rot90(material_in,k=np.rint(angle/90.).astype(int),axes=axes)
else:
microstructure_in = microstructure_out
material_in = material_out
origin = self.origin-(np.asarray(microstructure_in.shape)-self.grid)*.5 * self.size/self.grid
origin = self.origin-(np.asarray(material_in.shape)-self.grid)*.5 * self.size/self.grid
return self.duplicate(microstructure_in,
origin=origin,
comments=self.get_comments()+[util.execution_stamp('Geom','rotate')],
autosize=True,
)
return Geom(material = material_in,
size = self.size/self.grid*np.asarray(material_in.shape),
origin = origin,
comments = self.comments+[util.execution_stamp('Geom','rotate')],
)
def canvas(self,grid=None,offset=None,fill=None):
"""
Crop or enlarge/pad microstructure.
Crop or enlarge/pad geometry.
Parameters
----------
grid : numpy.ndarray of shape (3)
Number of grid points in x,y,z direction.
offset : numpy.ndarray of shape (3)
Offset (measured in grid points) from old to new microstructure[0,0,0].
Offset (measured in grid points) from old to new geometry [0,0,0].
fill : int or float, optional
Microstructure index to fill the background. Defaults to microstructure.max() + 1.
Material index to fill the background. Defaults to material.max() + 1.
"""
if offset is None: offset = 0
if fill is None: fill = np.nanmax(self.microstructure) + 1
dtype = float if int(fill) != fill or self.microstructure.dtype in np.sctypes['float'] else int
if fill is None: fill = np.nanmax(self.material) + 1
dtype = float if int(fill) != fill or self.material.dtype in np.sctypes['float'] else int
canvas = np.full(self.grid if grid is None else grid,fill,dtype)
@ -781,39 +640,41 @@ class Geom:
ll = np.clip(-offset, 0,np.minimum( grid,self.grid-offset))
ur = np.clip(-offset+self.grid,0,np.minimum( grid,self.grid-offset))
canvas[ll[0]:ur[0],ll[1]:ur[1],ll[2]:ur[2]] = self.microstructure[LL[0]:UR[0],LL[1]:UR[1],LL[2]:UR[2]]
canvas[ll[0]:ur[0],ll[1]:ur[1],ll[2]:ur[2]] = self.material[LL[0]:UR[0],LL[1]:UR[1],LL[2]:UR[2]]
return self.duplicate(canvas,
origin=self.origin+offset*self.size/self.grid,
comments=self.get_comments()+[util.execution_stamp('Geom','canvas')],
autosize=True,
)
return Geom(material = canvas,
size = self.size/self.grid*np.asarray(canvas.shape),
origin = self.origin+offset*self.size/self.grid,
comments = self.comments+[util.execution_stamp('Geom','canvas')],
)
def substitute(self,from_microstructure,to_microstructure):
def substitute(self,from_material,to_material):
"""
Substitute microstructure indices.
Substitute material indices.
Parameters
----------
from_microstructure : iterable of ints
Microstructure indices to be substituted.
to_microstructure : iterable of ints
New microstructure indices.
from_material : iterable of ints
Material indices to be substituted.
to_material : iterable of ints
New material indices.
"""
substituted = self.get_microstructure()
for from_ms,to_ms in zip(from_microstructure,to_microstructure):
substituted[self.microstructure==from_ms] = to_ms
substituted = self.material.copy()
for from_ms,to_ms in zip(from_material,to_material):
substituted[self.material==from_ms] = to_ms
return self.duplicate(substituted,
comments=self.get_comments()+[util.execution_stamp('Geom','substitute')],
)
return Geom(material = substituted,
size = self.size,
origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','substitute')],
)
def vicinity_offset(self,vicinity=1,offset=None,trigger=[],periodic=True):
"""
Offset microstructure index of points in the vicinity of xxx.
Offset material index of points in the vicinity of xxx.
Different from themselves (or listed as triggers) within a given (cubic) vicinity,
i.e. within the region close to a grain/phase boundary.
@ -822,14 +683,14 @@ class Geom:
Parameters
----------
vicinity : int, optional
Voxel distance checked for presence of other microstructure.
Voxel distance checked for presence of other materials.
Defaults to 1.
offset : int, optional
Offset (positive or negative) to tag microstructure indices,
defaults to microstructure.max() + 1.
Offset (positive or negative) to tag material indices,
defaults to material.max() + 1.
trigger : list of ints, optional
List of microstructure indices triggering a change.
Defaults to [], meaning that different neigboors trigger a change.
List of material indices that trigger a change.
Defaults to [], meaning that any different neighbor triggers a change.
periodic : Boolean, optional
Assume geometry to be periodic. Defaults to True.
@ -845,14 +706,15 @@ class Geom:
trigger = list(trigger)
return np.any(np.in1d(stencil,np.array(trigger)))
offset_ = np.nanmax(self.microstructure) if offset is None else offset
mask = ndimage.filters.generic_filter(self.microstructure,
offset_ = np.nanmax(self.material) if offset is None else offset
mask = ndimage.filters.generic_filter(self.material,
tainted_neighborhood,
size=1+2*vicinity,
mode='wrap' if periodic else 'nearest',
extra_keywords={'trigger':trigger})
microstructure = np.ma.MaskedArray(self.microstructure + offset_, np.logical_not(mask))
return self.duplicate(microstructure,
comments=self.get_comments()+[util.execution_stamp('Geom','vicinity_offset')],
)
return Geom(material = np.where(mask, self.material + offset_,self.material),
size = self.size,
origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','vicinity_offset')],
)

113
python/damask/seeds.py Normal file
View File

@ -0,0 +1,113 @@
from scipy import spatial as _spatial
import numpy as _np
from . import util
from . import grid_filters
def from_random(size,N_seeds,grid=None,seed=None):
"""
Random seeding in space.
Parameters
----------
size : numpy.ndarray of shape (3)
Physical size of the seeding domain.
N_seeds : int
Number of seeds.
grid : numpy.ndarray of shape (3), optional.
If given, ensures that all seeds initiate one grain if using a
standard Voronoi tessellation.
seed : {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional
A seed to initialize the BitGenerator. Defaults to None.
If None, then fresh, unpredictable entropy will be pulled from the OS.
"""
rng = _np.random.default_rng(seed)
if grid is None:
coords = rng.random((N_seeds,3)) * size
else:
grid_coords = grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F')
coords = grid_coords[rng.choice(_np.prod(grid),N_seeds, replace=False)] \
+ _np.broadcast_to(size/grid,(N_seeds,3))*(rng.random((N_seeds,3))*.5-.25) # wobble without leaving grid
return coords
def from_Poisson_disc(size,N_seeds,N_candidates,distance,periodic=True,seed=None):
"""
Seeding in space according to a Poisson disc distribution.
Parameters
----------
size : numpy.ndarray of shape (3)
Physical size of the seeding domain.
N_seeds : int
Number of seeds.
N_candidates : int
Number of candidates to consider for finding best candidate.
distance : float
Minimum acceptable distance to other seeds.
periodic : boolean, optional
Calculate minimum distance for periodically repeated grid.
seed : {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional
A seed to initialize the BitGenerator. Defaults to None.
If None, then fresh, unpredictable entropy will be pulled from the OS.
"""
rng = _np.random.default_rng(seed)
coords = _np.empty((N_seeds,3))
coords[0] = rng.random(3) * size
i = 1
progress = util._ProgressBar(N_seeds+1,'',50)
while i < N_seeds:
candidates = rng.random((N_candidates,3))*_np.broadcast_to(size,(N_candidates,3))
tree = _spatial.cKDTree(coords[:i],boxsize=size) if periodic else \
_spatial.cKDTree(coords[:i])
distances, dev_null = tree.query(candidates)
best = distances.argmax()
if distances[best] > distance: # require minimum separation
coords[i] = candidates[best] # maximum separation to existing point cloud
i += 1
progress.update(i)
return coords
def from_geom(geom,selection=None,invert=False,average=False,periodic=True):
"""
Create seed from existing geometry description.
Parameters
----------
geom : damask.Geom
Geometry, from which the material IDs are used as seeds.
selection : iterable of integers, optional
Material IDs to consider.
invert : boolean, false
Do not consider the material IDs given in selection. Defaults to False.
average : boolean, optional
Seed corresponds to center of gravity of material ID cloud.
periodic : boolean, optional
Center of gravity with periodic boundaries.
"""
material = geom.material.reshape((-1,1),order='F')
mask = _np.full(geom.grid.prod(),True,dtype=bool) if selection is None else \
_np.isin(material,selection,invert=invert)
coords = grid_filters.cell_coord0(geom.grid,geom.size).reshape(-1,3,order='F')
if not average:
return (coords[mask],material[mask])
else:
materials = _np.unique(material[mask])
coords_ = _np.zeros((materials.size,3),dtype=float)
for i,mat in enumerate(materials):
pc = (2*_np.pi*coords[material[:,0]==mat,:]-geom.origin)/geom.size
coords_[i] = geom.origin + geom.size / 2 / _np.pi * (_np.pi +
_np.arctan2(-_np.average(_np.sin(pc),axis=0),
-_np.average(_np.cos(pc),axis=0))) \
if periodic else \
_np.average(coords[material[:,0]==mat,:],axis=0)
return (coords_,materials)

View File

@ -173,7 +173,7 @@ def scale_to_coprime(v):
m = (np.array(v) * reduce(lcm, map(lambda x: int(get_square_denominator(x)),v)) ** 0.5).astype(np.int)
m = m//reduce(np.gcd,m)
with np.errstate(divide='ignore'):
with np.errstate(invalid='ignore'):
if not np.allclose(np.ma.masked_invalid(v/m),v[np.argmax(abs(v))]/m[np.argmax(abs(v))]):
raise ValueError(f'Invalid result {m} for input {v}. Insufficient precision?')

View File

@ -1,11 +1,16 @@
<?xml version="1.0"?>
<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
<RectilinearGrid WholeExtent="0 8 0 5 0 4">
<FieldData>
<Array type="String" Name="comments" NumberOfTuples="1" format="binary">
AQAAAACAAABJAAAATgAAAA==eF4FwUEKgCAUBNCO4rIWX8ZJsbxA5/iUFqQVBJ2/9zZt+p52yXeza816mW+0sBCtz6HCGGSPE1wJjMX0BCGYhTQuJLrkKfDA0P0d3xK6
</Array>
</FieldData>
<Piece Extent="0 8 0 5 0 4">
<PointData>
</PointData>
<CellData>
<DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
<DataArray type="Int64" Name="material" format="binary" RangeMin="1" RangeMax="41">
AQAAAACAAAAABQAAZwAAAA==eF7t0rcOgmAAhVEgNmyo2AuoWN//BR04EwsJcfzvcvabL47qxcFOJg177HPAIUdMOeaEU844Z8YFl1wx55obbrnjngceeeKZFxYseeWNd1Z88MkX3/zwy+Z/wf8YOqzX1uEPlgwHCA==
</DataArray>
</CellData>

View File

@ -1,11 +1,16 @@
<?xml version="1.0"?>
<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
<RectilinearGrid WholeExtent="0 8 0 5 0 4">
<FieldData>
<Array type="String" Name="comments" NumberOfTuples="1" format="binary">
AQAAAACAAABJAAAATwAAAA==eF5LScxNLM7Wc0/Nz9VLzklNzFMoM9Yz0DPQTcwpyEjUNTI31U03tzAwTDM1Mk9T0DAyMDLQNbDUNTJSMDS1MjK0MgFyTQwMNBkAHc8SuA==
</Array>
</FieldData>
<Piece Extent="0 8 0 5 0 4">
<PointData>
</PointData>
<CellData>
<DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
<DataArray type="Int64" Name="material" format="binary" RangeMin="1" RangeMax="41">
AQAAAACAAAAABQAAagAAAA==eF7t0rkOglAARFExLrgCKuKuqLj8/w9acCoSY7B+05x+cqNOvSj4l92GPfY54JAxRxxzwilnnDNhyowLLrlizjULbrjljnseeOSJZ15Y8sob76z44JMvvtn8L9jObz2GDuv96vADk5QHBg==
</DataArray>
</CellData>

View File

@ -1,11 +1,16 @@
<?xml version="1.0"?>
<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
<RectilinearGrid WholeExtent="0 8 0 5 0 4">
<FieldData>
<Array type="String" Name="comments" NumberOfTuples="1" format="binary">
AQAAAACAAABJAAAATwAAAA==eF4FwdEJgDAMBUBH6ad+JLzElmoXcI6grYKtCoLze7dZs/fkJd+N15rtct/IYJDV5zDSGGiPE6QEjcX1CgVhJlUnIakkLwQPDN0PHdcSuQ==
</Array>
</FieldData>
<Piece Extent="0 8 0 5 0 4">
<PointData>
</PointData>
<CellData>
<DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
<DataArray type="Int64" Name="material" format="binary" RangeMin="1" RangeMax="41">
AQAAAACAAAAABQAAZAAAAA==eF7t0scRglAAQEEBAyZUMCuomPtv0ANbgMNw/O+yDbyo1xQFWxkzYZ8DDjliyjEnnHLGOTMuuOSKOQuuueGWO+554JEnnlmy4oVX3ljzzgeffPHND7+Mg50aPmz698MfmvQHCg==
</DataArray>
</CellData>

View File

@ -1,11 +1,16 @@
<?xml version="1.0"?>
<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
<RectilinearGrid WholeExtent="0 8 0 5 0 4">
<FieldData>
<Array type="String" Name="comments" NumberOfTuples="1" format="binary">
AQAAAACAAABJAAAATwAAAA==eF4FwVEKgCAQBcCO4md97PJcE9MLdI6ltCCtIOj8zuza9Lt4zU/jrWa9ze8YDNL6nkoSPB1hgS1eQjGjQECIJGKsT2KTi4QZmIYOHg4SwA==
</Array>
</FieldData>
<Piece Extent="0 8 0 5 0 4">
<PointData>
</PointData>
<CellData>
<DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
<DataArray type="Int64" Name="material" format="binary" RangeMin="1" RangeMax="41">
AQAAAACAAAAABQAAYwAAAA==eF7t0scBgkAAAEHBgBEwgDmBsf8GfTANCN/bzzSwUa8pCrYyZp8DDjliwjEnnHLGORdMmTHnkiuuuWHBklvuuOeBR5545oVX3nhnxZoPPvnimx9+GQc7GT5sqvjvhz+ZtAcJ
</DataArray>
</CellData>

View File

@ -1,11 +1,16 @@
<?xml version="1.0"?>
<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
<RectilinearGrid WholeExtent="0 8 0 5 0 4">
<FieldData>
<Array type="String" Name="comments" NumberOfTuples="1" format="binary">
AQAAAACAAABJAAAATwAAAA==eF4FwdEJgDAMBUBH6ad+JLzElmoXcI6grYKtCoLze7dZs/fkJd+N15rtct/IYJDV5zDSGGiPE6QEjcX1CgVhJlUnIakkLwQPDN0PHdcSuQ==
</Array>
</FieldData>
<Piece Extent="0 8 0 5 0 4">
<PointData>
</PointData>
<CellData>
<DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="2">
<DataArray type="Int64" Name="material" format="binary" RangeMin="1" RangeMax="2">
AQAAAACAAAAABQAAGwAAAA==eF5jZIAAxlF6lB4AmmmUpogeDUfKaAD7jwDw
</DataArray>
</CellData>

View File

@ -1,11 +1,16 @@
<?xml version="1.0"?>
<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
<RectilinearGrid WholeExtent="0 8 0 5 0 4">
<FieldData>
<Array type="String" Name="comments" NumberOfTuples="1" format="binary">
AQAAAACAAABJAAAATwAAAA==eF4FwVEKgCAQBcCO4md97PJcE9MLdI6ltCCtIOj8zuza9Lt4zU/jrWa9ze8YDNL6nkoSPB1hgS1eQjGjQECIJGKsT2KTi4QZmIYOHg4SwA==
</Array>
</FieldData>
<Piece Extent="0 8 0 5 0 4">
<PointData>
</PointData>
<CellData>
<DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="2">
<DataArray type="Int64" Name="material" format="binary" RangeMin="1" RangeMax="2">
AQAAAACAAAAABQAAGQAAAA==eF5jZIAAxlF6lB4AmmmUHqUHkAYA/M8A8Q==
</DataArray>
</CellData>

View File

@ -1,11 +1,16 @@
<?xml version="1.0"?>
<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
<RectilinearGrid WholeExtent="0 8 0 5 0 4">
<FieldData>
<Array type="String" Name="comments" NumberOfTuples="1" format="binary">
AQAAAACAAABJAAAATgAAAA==eF4FwUEKgCAUBNCO4rIWX8ZJsbxA5/iUFqQVBJ2/9zZt+p52yXeza816mW+0sBCtz6HCGGSPE1wJjMX0BCGYhTQuJLrkKfDA0P0d3xK6
</Array>
</FieldData>
<Piece Extent="0 8 0 5 0 4">
<PointData>
</PointData>
<CellData>
<DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
<DataArray type="Int64" Name="material" format="binary" RangeMin="1" RangeMax="41">
AQAAAACAAAAABQAAZwAAAA==eF7t0rcOgmAAhVEgNmyo2AuoWN//BR04EwsJcfzvcvabL47qxcFOJg177HPAIUdMOeaEU844Z8YFl1wx55obbrnjngceeeKZFxYseeWNd1Z88MkX3/zwy+Z/wf8YOqzX1uEPlgwHCA==
</DataArray>
</CellData>

View File

@ -1,11 +1,16 @@
<?xml version="1.0"?>
<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
<RectilinearGrid WholeExtent="0 8 0 5 0 4">
<FieldData>
<Array type="String" Name="comments" NumberOfTuples="1" format="binary">
AQAAAACAAABJAAAATwAAAA==eF5LScxNLM7Wc0/Nz9VLzklNzFMoM9Yz0DPQTcwpyEjUNTI31U03tzAwTDM1Mk9T0DAyMDLQNbDUNTJSMDS1MjK0MgFyTQwMNBkAHc8SuA==
</Array>
</FieldData>
<Piece Extent="0 8 0 5 0 4">
<PointData>
</PointData>
<CellData>
<DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
<DataArray type="Int64" Name="material" format="binary" RangeMin="1" RangeMax="41">
AQAAAACAAAAABQAAagAAAA==eF7t0rkOglAARFExLrgCKuKuqLj8/w9acCoSY7B+05x+cqNOvSj4l92GPfY54JAxRxxzwilnnDNhyowLLrlizjULbrjljnseeOSJZ15Y8sob76z44JMvvtn8L9jObz2GDuv96vADk5QHBg==
</DataArray>
</CellData>

View File

@ -1,11 +1,16 @@
<?xml version="1.0"?>
<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
<RectilinearGrid WholeExtent="0 8 0 5 0 4">
<FieldData>
<Array type="String" Name="comments" NumberOfTuples="1" format="binary">
AQAAAACAAABJAAAATwAAAA==eF4FwdEJgDAMBUBH6ad+JLzElmoXcI6grYKtCoLze7dZs/fkJd+N15rtct/IYJDV5zDSGGiPE6QEjcX1CgVhJlUnIakkLwQPDN0PHdcSuQ==
</Array>
</FieldData>
<Piece Extent="0 8 0 5 0 4">
<PointData>
</PointData>
<CellData>
<DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
<DataArray type="Int64" Name="material" format="binary" RangeMin="1" RangeMax="41">
AQAAAACAAAAABQAAZAAAAA==eF7t0scRglAAQEEBAyZUMCuomPtv0ANbgMNw/O+yDbyo1xQFWxkzYZ8DDjliyjEnnHLGOTMuuOSKOQuuueGWO+554JEnnlmy4oVX3ljzzgeffPHND7+Mg50aPmz698MfmvQHCg==
</DataArray>
</CellData>

View File

@ -1,11 +1,16 @@
<?xml version="1.0"?>
<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
<RectilinearGrid WholeExtent="0 8 0 5 0 4">
<FieldData>
<Array type="String" Name="comments" NumberOfTuples="1" format="binary">
AQAAAACAAABJAAAATwAAAA==eF5LScxNLM7Wc0/Nz9VLzklNzFMoM9Yz0DPQTcwpyEjUNTI31U03tzAwTDM1Mk9T0DAyMDLQNbDUNTJSMDS1MjK0MgFyTQwMNBkAHc8SuA==
</Array>
</FieldData>
<Piece Extent="0 8 0 5 0 4">
<PointData>
</PointData>
<CellData>
<DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
<DataArray type="Int64" Name="material" format="binary" RangeMin="1" RangeMax="41">
AQAAAACAAAAABQAAYwAAAA==eF7t0scBgkAAAEHBgBEwgDmBsf8GfTANCN/bzzSwUa8pCrYyZp8DDjliwjEnnHLGORdMmTHnkiuuuWHBklvuuOeBR5545oVX3nhnxZoPPvnimx9+GQc7GT5sqvjvhz+ZtAcJ
</DataArray>
</CellData>

View File

@ -1,11 +1,16 @@
<?xml version="1.0"?>
<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
<RectilinearGrid WholeExtent="0 8 0 5 0 4">
<FieldData>
<Array type="String" Name="comments" NumberOfTuples="1" format="binary">
AQAAAACAAABJAAAATwAAAA==eF4FwdEJgDAMBUBH6ad+JLzElmoXcI6grYKtCoLze7dZs/fkJd+N15rtct/IYJDV5zDSGGiPE6QEjcX1CgVhJlUnIakkLwQPDN0PHdcSuQ==
</Array>
</FieldData>
<Piece Extent="0 8 0 5 0 4">
<PointData>
</PointData>
<CellData>
<DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="2">
<DataArray type="Int64" Name="material" format="binary" RangeMin="1" RangeMax="2">
AQAAAACAAAAABQAAIgAAAA==eF5jZIAAxlGaLJoJjSakntr6hzqN7v9RepSmJw0AC04A9Q==
</DataArray>
</CellData>

View File

@ -1,11 +1,16 @@
<?xml version="1.0"?>
<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
<RectilinearGrid WholeExtent="0 8 0 5 0 4">
<FieldData>
<Array type="String" Name="comments" NumberOfTuples="1" format="binary">
AQAAAACAAABJAAAATwAAAA==eF4FwVEKgCAQBcCO4md97PJcE9MLdI6ltCCtIOj8zuza9Lt4zU/jrWa9ze8YDNL6nkoSPB1hgS1eQjGjQECIJGKsT2KTi4QZmIYOHg4SwA==
</Array>
</FieldData>
<Piece Extent="0 8 0 5 0 4">
<PointData>
</PointData>
<CellData>
<DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="2">
<DataArray type="Int64" Name="material" format="binary" RangeMin="1" RangeMax="2">
AQAAAACAAAAABQAALwAAAA==eF5jZIAAxlGaLJoJjSakHpc+cvUTUkdrmlL3j9KU0dROF5TqH2iaVPcDAALOANU=
</DataArray>
</CellData>

View File

@ -1,11 +1,16 @@
<?xml version="1.0"?>
<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
<RectilinearGrid WholeExtent="0 8 0 5 0 4">
<FieldData>
<Array type="String" Name="comments" NumberOfTuples="1" format="binary">
AQAAAACAAABJAAAATgAAAA==eF4FwUEKgCAUBNCO4rIWX8ZJsbxA5/iUFqQVBJ2/9zZt+p52yXeza816mW+0sBCtz6HCGGSPE1wJjMX0BCGYhTQuJLrkKfDA0P0d3xK6
</Array>
</FieldData>
<Piece Extent="0 8 0 5 0 4">
<PointData>
</PointData>
<CellData>
<DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
<DataArray type="Int64" Name="material" format="binary" RangeMin="1" RangeMax="41">
AQAAAACAAAAABQAAcQAAAA==eF7t0rkOglAUBFAxKu6igvsKrv//gxYcm9fQGEPBNKe6yc1kolaZqPEndthljzH7HHDIEceccMoZE8654JIpM6645oZb7rjngUeeeOaFV+YseOOdDz754pthf+3Aqr7rdv9vw3+/NjssU7XDD0/8BuQ=
</DataArray>
</CellData>

View File

@ -1,11 +1,16 @@
<?xml version="1.0"?>
<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
<RectilinearGrid WholeExtent="0 8 0 5 0 4">
<FieldData>
<Array type="String" Name="comments" NumberOfTuples="1" format="binary">
AQAAAACAAABJAAAATwAAAA==eF5LScxNLM7Wc0/Nz9VLzklNzFMoM9Yz0DPQTcwpyEjUNTI31U03tzAwTDM1Mk9T0DAyMDLQNbDUNTJSMDS1MjK0MgFyTQwMNBkAHc8SuA==
</Array>
</FieldData>
<Piece Extent="0 8 0 5 0 4">
<PointData>
</PointData>
<CellData>
<DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
<DataArray type="Int64" Name="material" format="binary" RangeMin="1" RangeMax="41">
AQAAAACAAAAABQAAYQAAAA==eF7t0scVglAAAEHgqZgBA2ZExdR/gx6YCpDj38s0sEnUlgR7ccAhR0w55oRTzjjngktmzFlwxTU33LLkjnseeOSJZ15Y8cqaN975YMMnX3zzwy/j4F+GD9u6fvgD+gwHCA==
</DataArray>
</CellData>

View File

@ -1,11 +1,16 @@
<?xml version="1.0"?>
<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
<RectilinearGrid WholeExtent="0 8 0 5 0 4">
<FieldData>
<Array type="String" Name="comments" NumberOfTuples="1" format="binary">
AQAAAACAAABJAAAATgAAAA==eF4FwUEKgCAUBNCO4rIWX8ZJsbxA5/iUFqQVBJ2/9zZt+p52yXeza816mW+0sBCtz6HCGGSPE1wJjMX0BCGYhTQuJLrkKfDA0P0d3xK6
</Array>
</FieldData>
<Piece Extent="0 8 0 5 0 4">
<PointData>
</PointData>
<CellData>
<DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
<DataArray type="Int64" Name="material" format="binary" RangeMin="1" RangeMax="41">
AQAAAACAAAAABQAAZAAAAA==eF7t0scRglAAQEEBAyZUMCuomPtv0ANbgMNw/O+yDbyo1xQFWxkzYZ8DDjliyjEnnHLGOTMuuOSKOQuuueGWO+554JEnnlmy4oVX3ljzzgeffPHND7+Mg50aPmz698MfmvQHCg==
</DataArray>
</CellData>

View File

@ -1,11 +1,16 @@
<?xml version="1.0"?>
<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
<RectilinearGrid WholeExtent="0 8 0 5 0 4">
<FieldData>
<Array type="String" Name="comments" NumberOfTuples="1" format="binary">
AQAAAACAAABJAAAATwAAAA==eF5LScxNLM7Wc0/Nz9VLzklNzFMoM9Yz0DPQTcwpyEjUNTI31U03tzAwTDM1Mk9T0DAyMDLQNbDUNTJSMDS1MjK0MgFyTQwMNBkAHc8SuA==
</Array>
</FieldData>
<Piece Extent="0 8 0 5 0 4">
<PointData>
</PointData>
<CellData>
<DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="2" RangeMax="41">
<DataArray type="Int64" Name="material" format="binary" RangeMin="2" RangeMax="41">
AQAAAACAAAAABQAAZAAAAA==eF7t0rcSglAARFEHE0bAgBkE8///oAWnF8b2bXP6nRv1mkXBv+xzwCFHHDPmhFPOOOeCSyZMmXHFNTfcMueOex545IlnXliw5JUVa95454NPvvjmh79+DXYzdNisbYdfSqMHMg==
</DataArray>
</CellData>

View File

@ -1,11 +1,16 @@
<?xml version="1.0"?>
<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
<RectilinearGrid WholeExtent="0 8 0 5 0 4">
<FieldData>
<Array type="String" Name="comments" NumberOfTuples="1" format="binary">
AQAAAACAAABJAAAATwAAAA==eF4FwdEJgDAMBUBH6ad+JLzElmoXcI6grYKtCoLze7dZs/fkJd+N15rtct/IYJDV5zDSGGiPE6QEjcX1CgVhJlUnIakkLwQPDN0PHdcSuQ==
</Array>
</FieldData>
<Piece Extent="0 8 0 5 0 4">
<PointData>
</PointData>
<CellData>
<DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="2">
<DataArray type="Int64" Name="material" format="binary" RangeMin="1" RangeMax="2">
AQAAAACAAAAABQAAIAAAAA==eF5jZIAAxlF6lB4AmokAPdj1DzRNyP2jNH4aAMufANU=
</DataArray>
</CellData>

View File

@ -1,11 +1,16 @@
<?xml version="1.0"?>
<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
<RectilinearGrid WholeExtent="0 8 0 5 0 4">
<FieldData>
<Array type="String" Name="comments" NumberOfTuples="1" format="binary">
AQAAAACAAABJAAAATwAAAA==eF4FwVEKgCAQBcCO4md97PJcE9MLdI6ltCCtIOj8zuza9Lt4zU/jrWa9ze8YDNL6nkoSPB1hgS1eQjGjQECIJGKsT2KTi4QZmIYOHg4SwA==
</Array>
</FieldData>
<Piece Extent="0 8 0 5 0 4">
<PointData>
</PointData>
<CellData>
<DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="2">
<DataArray type="Int64" Name="material" format="binary" RangeMin="1" RangeMax="2">
AQAAAACAAAAABQAAMAAAAA==eF5jYoAAJhw0IwEalz566aeUptT+oa6fUppS+4e6fkppSu0f6voppSm1HwBAngDh
</DataArray>
</CellData>

View File

@ -11,9 +11,9 @@ from damask import util
def geom_equal(a,b):
return np.all(a.get_microstructure() == b.get_microstructure()) and \
np.all(a.get_grid() == b.get_grid()) and \
np.allclose(a.get_size(), b.get_size()) and \
return np.all(a.material == b.material) and \
np.all(a.grid == b.grid) and \
np.allclose(a.size, b.size) and \
str(a.diff(b)) == str(b.diff(a))
@pytest.fixture
@ -33,40 +33,21 @@ def reference_dir(reference_dir_base):
class TestGeom:
@pytest.mark.parametrize('flavor',['plain','explicit'])
def test_duplicate(self,default,flavor):
if flavor == 'plain':
modified = default.duplicate()
elif flavor == 'explicit':
modified = default.duplicate(
default.get_microstructure(),
default.get_size(),
default.get_origin()
)
print(modified)
assert geom_equal(default,modified)
def test_diff_equal(self,default):
assert str(default.diff(default)) == ''
def test_diff_not_equal(self,default):
new = Geom(default.microstructure[1:,1:,1:]+1,default.size*.9,np.ones(3)-default.origin,comments=['modified'])
new = Geom(default.material[1:,1:,1:]+1,default.size*.9,np.ones(3)-default.origin,comments=['modified'])
assert str(default.diff(new)) != ''
@pytest.mark.parametrize('masked',[True,False])
def test_set_microstructure(self,default,masked):
old = default.get_microstructure()
new = np.random.randint(200,size=default.grid)
default.set_microstructure(np.ma.MaskedArray(new,np.full_like(new,masked)))
assert np.all(default.microstructure==(old if masked else new))
def test_write_read_str(self,default,tmpdir):
default.save_ASCII(str(tmpdir/'default.geom'))
new = Geom.load_ASCII(str(tmpdir/'default.geom'))
assert geom_equal(default,new)
def test_write_read_file(self,default,tmpdir):
with open(tmpdir/'default.geom','w') as f:
default.save_ASCII(f,compress=True)
@ -74,6 +55,7 @@ class TestGeom:
new = Geom.load_ASCII(f)
assert geom_equal(default,new)
def test_read_write_vtr(self,default,tmpdir):
default.save(tmpdir/'default')
for _ in range(10):
@ -83,6 +65,7 @@ class TestGeom:
new = Geom.load(tmpdir/'default.vtr')
assert geom_equal(new,default)
def test_invalid_geom(self,tmpdir):
with open('invalid_file','w') as f:
f.write('this is not a valid header')
@ -90,6 +73,7 @@ class TestGeom:
with pytest.raises(TypeError):
Geom.load_ASCII(f)
def test_invalid_vtr(self,tmpdir):
v = VTK.from_rectilinearGrid(np.random.randint(5,10,3)*2,np.random.random(3) + 1.0)
v.save(tmpdir/'no_materialpoint.vtr')
@ -99,6 +83,13 @@ class TestGeom:
with pytest.raises(ValueError):
Geom.load(tmpdir/'no_materialpoint.vtr')
def test_invalid_material(self):
with pytest.raises(TypeError):
Geom(np.zeros((3,3,3),dtype='complex'),np.ones(3))
def test_cast_to_int(self):
g = Geom(np.zeros((3,3,3)),np.ones(3))
assert g.material.dtype in np.sctypes['int']
@pytest.mark.parametrize('compress',[True,False])
def test_compress(self,default,tmpdir,compress):
@ -106,31 +97,32 @@ class TestGeom:
new = Geom.load_ASCII(tmpdir/'default.geom')
assert geom_equal(new,default)
def test_invalid_combination(self,default):
with pytest.raises(ValueError):
default.duplicate(default.microstructure[1:,1:,1:],size=np.ones(3), autosize=True)
def test_invalid_size(self,default):
with pytest.raises(ValueError):
default.duplicate(default.microstructure[1:,1:,1:],size=np.ones(2))
Geom(default.material[1:,1:,1:],
size=np.ones(2))
def test_invalid_origin(self,default):
with pytest.raises(ValueError):
default.duplicate(default.microstructure[1:,1:,1:],origin=np.ones(4))
Geom(default.material[1:,1:,1:],
size=np.ones(3),
origin=np.ones(4))
def test_invalid_microstructure_size(self,default):
microstructure = np.ones((3,3))
def test_invalid_materials_shape(self,default):
material = np.ones((3,3))
with pytest.raises(ValueError):
default.duplicate(microstructure)
Geom(material,
size=np.ones(3))
def test_invalid_microstructure_type(self,default):
microstructure = np.random.randint(1,300,(3,4,5))==1
with pytest.raises(TypeError):
default.duplicate(microstructure)
def test_invalid_homogenization(self,default):
def test_invalid_materials_type(self,default):
material = np.random.randint(1,300,(3,4,5))==1
with pytest.raises(TypeError):
default.set_homogenization(homogenization=0)
Geom(material)
@pytest.mark.parametrize('directions,reflect',[
(['x'], False),
@ -147,6 +139,7 @@ class TestGeom:
assert geom_equal(Geom.load_ASCII(reference),
modified)
@pytest.mark.parametrize('directions',[(1,2,'y'),('a','b','x'),[1]])
def test_mirror_invalid(self,default,directions):
with pytest.raises(ValueError):
@ -168,13 +161,16 @@ class TestGeom:
assert geom_equal(Geom.load_ASCII(reference),
modified)
def test_flip_invariant(self,default):
assert geom_equal(default,default.flip([]))
@pytest.mark.parametrize('direction',[['x'],['x','y']])
def test_flip_double(self,default,direction):
assert geom_equal(default,default.flip(direction).flip(direction))
@pytest.mark.parametrize('directions',[(1,2,'y'),('a','b','x'),[1]])
def test_flip_invalid(self,default,directions):
with pytest.raises(ValueError):
@ -196,6 +192,7 @@ class TestGeom:
current
)
@pytest.mark.parametrize('grid',[
(10,11,10),
[10,13,10],
@ -213,22 +210,29 @@ class TestGeom:
assert geom_equal(Geom.load_ASCII(reference),
modified)
def test_renumber(self,default):
microstructure = default.get_microstructure()
for m in np.unique(microstructure):
microstructure[microstructure==m] = microstructure.max() + np.random.randint(1,30)
modified = default.duplicate(microstructure)
material = default.material.copy()
for m in np.unique(material):
material[material==m] = material.max() + np.random.randint(1,30)
modified = Geom(material,
default.size,
default.origin)
assert not geom_equal(modified,default)
assert geom_equal(default,
modified.renumber())
def test_substitute(self,default):
offset = np.random.randint(1,500)
modified = default.duplicate(default.get_microstructure() + offset)
modified = Geom(default.material + offset,
default.size,
default.origin)
assert not geom_equal(modified,default)
assert geom_equal(default,
modified.substitute(np.arange(default.microstructure.max())+1+offset,
np.arange(default.microstructure.max())+1))
modified.substitute(np.arange(default.material.max())+1+offset,
np.arange(default.material.max())+1))
@pytest.mark.parametrize('axis_angle',[np.array([1,0,0,86.7]), np.array([0,1,0,90.4]), np.array([0,0,1,90]),
np.array([1,0,0,175]),np.array([0,-1,0,178]),np.array([0,0,1,180])])
@ -238,6 +242,7 @@ class TestGeom:
modified.rotate(Rotation.from_axis_angle(axis_angle,degrees=True))
assert geom_equal(default,modified)
@pytest.mark.parametrize('Eulers',[[32.0,68.0,21.0],
[0.0,32.0,240.0]])
def test_rotate(self,default,update,reference_dir,Eulers):
@ -248,11 +253,13 @@ class TestGeom:
assert geom_equal(Geom.load_ASCII(reference),
modified)
def test_canvas(self,default):
grid = default.grid
grid_add = np.random.randint(0,30,(3))
modified = default.canvas(grid + grid_add)
assert np.all(modified.microstructure[:grid[0],:grid[1],:grid[2]] == default.microstructure)
assert np.all(modified.material[:grid[0],:grid[1],:grid[2]] == default.material)
@pytest.mark.parametrize('center1,center2',[(np.random.random(3)*.5,np.random.random()*8),
(np.random.randint(4,8,(3)),np.random.randint(9,12,(3)))])
@ -265,13 +272,14 @@ class TestGeom:
np.random.rand()*4,
np.random.randint(20)])
def test_add_primitive_shift(self,center1,center2,diameter,exponent):
"""Same volume fraction for periodic microstructures and different center."""
"""Same volume fraction for periodic geometries and different center."""
o = np.random.random(3)-.5
g = np.random.randint(8,32,(3))
s = np.random.random(3)+.5
G_1 = Geom(np.ones(g,'i'),s,o).add_primitive(diameter,center1,exponent)
G_2 = Geom(np.ones(g,'i'),s,o).add_primitive(diameter,center2,exponent)
assert np.count_nonzero(G_1.microstructure!=2) == np.count_nonzero(G_2.microstructure!=2)
assert np.count_nonzero(G_1.material!=2) == np.count_nonzero(G_2.material!=2)
@pytest.mark.parametrize('center',[np.random.randint(4,10,(3)),
np.random.randint(2,10),
@ -288,6 +296,7 @@ class TestGeom:
G_2 = Geom(np.ones(g,'i'),[1.,1.,1.]).add_primitive(.3,center,1,fill,Rotation.from_Eulers(eu),inverse,periodic=periodic)
assert geom_equal(G_1,G_2)
@pytest.mark.parametrize('trigger',[[1],[]])
def test_vicinity_offset(self,trigger):
offset = np.random.randint(2,4)
@ -306,13 +315,15 @@ class TestGeom:
geom = Geom(m,np.random.rand(3)).vicinity_offset(vicinity,offset,trigger=trigger)
assert np.all(m2==geom.microstructure)
assert np.all(m2==geom.material)
@pytest.mark.parametrize('periodic',[True,False])
def test_vicinity_offset_invariant(self,default,periodic):
old = default.get_microstructure()
default.vicinity_offset(trigger=[old.max()+1,old.min()-1])
assert np.all(old==default.microstructure)
offset = default.vicinity_offset(trigger=[default.material.max()+1,
default.material.min()-1])
assert np.all(offset.material==default.material)
@pytest.mark.parametrize('periodic',[True,False])
def test_tessellation_approaches(self,periodic):
@ -320,10 +331,11 @@ class TestGeom:
size = np.random.random(3) + 1.0
N_seeds= np.random.randint(10,30)
seeds = np.random.rand(N_seeds,3) * np.broadcast_to(size,(N_seeds,3))
Voronoi = Geom.from_Voronoi_tessellation( grid,size,seeds, periodic)
Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(N_seeds),periodic)
Voronoi = Geom.from_Voronoi_tessellation( grid,size,seeds, np.arange(N_seeds)+5,periodic)
Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(N_seeds),np.arange(N_seeds)+5,periodic)
assert geom_equal(Laguerre,Voronoi)
def test_Laguerre_weights(self):
grid = np.random.randint(10,20,3)
size = np.random.random(3) + 1.0
@ -332,18 +344,19 @@ class TestGeom:
weights= np.full((N_seeds),-np.inf)
ms = np.random.randint(1, N_seeds+1)
weights[ms-1] = np.random.random()
Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,weights,np.random.random()>0.5)
assert np.all(Laguerre.microstructure == ms)
Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,weights,periodic=np.random.random()>0.5)
assert np.all(Laguerre.material == ms)
@pytest.mark.parametrize('approach',['Laguerre','Voronoi'])
def test_tessellate_bicrystal(self,approach):
grid = np.random.randint(5,10,3)*2
size = grid.astype(np.float)
seeds = np.vstack((size*np.array([0.5,0.25,0.5]),size*np.array([0.5,0.75,0.5])))
microstructure = np.ones(grid)
microstructure[:,grid[1]//2:,:] = 2
material = np.ones(grid)
material[:,grid[1]//2:,:] = 2
if approach == 'Laguerre':
geom = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(2),np.random.random()>0.5)
geom = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(2),periodic=np.random.random()>0.5)
elif approach == 'Voronoi':
geom = Geom.from_Voronoi_tessellation(grid,size,seeds, np.random.random()>0.5)
assert np.all(geom.microstructure == microstructure)
geom = Geom.from_Voronoi_tessellation(grid,size,seeds, periodic=np.random.random()>0.5)
assert np.all(geom.material == material)

View File

@ -0,0 +1,36 @@
import pytest
import numpy as np
from scipy.spatial import cKDTree
from damask import seeds
from damask import Geom
class TestSeeds:
@pytest.mark.parametrize('grid',[None,np.ones(3,dtype='i')*10])
def test_from_random(self,grid):
N_seeds = np.random.randint(30,300)
size = np.ones(3) + np.random.random(3)
coords = seeds.from_random(size,N_seeds,grid)
assert (0<=coords).all() and (coords<size).all()
@pytest.mark.parametrize('periodic',[True,False])
def test_from_Poisson_disc(self,periodic):
N_seeds = np.random.randint(30,300)
N_candidates = N_seeds//15
distance = np.random.random()
size = np.ones(3)*distance*N_seeds
coords = seeds.from_Poisson_disc(size,N_seeds,N_candidates,distance,periodic=periodic)
min_dists, _ = cKDTree(coords,boxsize=size).query(coords, 2) if periodic else \
cKDTree(coords).query(coords, 2)
assert (0<= coords).all() and (coords<size).all() and np.min(min_dists[:,1])>=distance
def test_from_geom(self):
grid = np.random.randint(10,20,3)
N_seeds = np.random.randint(30,300)
size = np.ones(3) + np.random.random(3)
coords = seeds.from_random(size,N_seeds,grid)
geom_1 = Geom.from_Voronoi_tessellation(grid,size,coords)
coords,material = seeds.from_geom(geom_1)
geom_2 = Geom.from_Voronoi_tessellation(grid,size,coords,material)
assert (geom_2.material==geom_1.material).all()

View File

@ -17,18 +17,18 @@ submodule(constitutive:constitutive_plastic) plastic_disloTungsten
D_0 = 1.0_pReal, & !< prefactor for self-diffusion coefficient
Q_cl = 1.0_pReal !< activation energy for dislocation climb
real(pReal), allocatable, dimension(:) :: &
b_sl, & !< magnitude of burgers vector [m]
b_sl, & !< magnitude of Burgers vector [m]
D_a, &
i_sl, & !< Adj. parameter for distance between 2 forest dislocations
atomicVolume, & !< factor to calculate atomic volume
tau_0, & !< Peierls stress
f_at, & !< factor to calculate atomic volume
tau_Peierls, & !< Peierls stress
!* mobility law parameters
delta_F, & !< activation energy for glide [J]
v0, & !< dislocation velocity prefactor [m/s]
Q_s, & !< activation energy for glide [J]
v_0, & !< dislocation velocity prefactor [m/s]
p, & !< p-exponent in glide velocity
q, & !< q-exponent in glide velocity
B, & !< friction coefficient
kink_height, & !< height of the kink pair
h, & !< height of the kink pair
w, & !< width of the kink pair
omega !< attempt frequency for kink pair nucleation
real(pReal), allocatable, dimension(:,:) :: &
@ -142,7 +142,7 @@ module function plastic_disloTungsten_init() result(myPlasticity)
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
if(trim(phase%get_asString('lattice')) == 'bcc') then
a = pl%get_asFloats('nonSchmid_coefficients',defaultVal = emptyRealArray)
a = pl%get_asFloats('a_nonSchmid',defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,a,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,a,-1)
else
@ -158,17 +158,17 @@ module function plastic_disloTungsten_init() result(myPlasticity)
rho_mob_0 = pl%get_asFloats('rho_mob_0', requiredSize=size(N_sl))
rho_dip_0 = pl%get_asFloats('rho_dip_0', requiredSize=size(N_sl))
prm%v0 = pl%get_asFloats('v_0', requiredSize=size(N_sl))
prm%v_0 = pl%get_asFloats('v_0', requiredSize=size(N_sl))
prm%b_sl = pl%get_asFloats('b_sl', requiredSize=size(N_sl))
prm%delta_F = pl%get_asFloats('Q_s', requiredSize=size(N_sl))
prm%Q_s = pl%get_asFloats('Q_s', requiredSize=size(N_sl))
prm%i_sl = pl%get_asFloats('i_sl', requiredSize=size(N_sl))
prm%tau_0 = pl%get_asFloats('tau_peierls', requiredSize=size(N_sl))
prm%tau_Peierls = pl%get_asFloats('tau_Peierls', requiredSize=size(N_sl))
prm%p = pl%get_asFloats('p_sl', requiredSize=size(N_sl), &
defaultVal=[(1.0_pReal,i=1,size(N_sl))])
prm%q = pl%get_asFloats('q_sl', requiredSize=size(N_sl), &
defaultVal=[(1.0_pReal,i=1,size(N_sl))])
prm%kink_height = pl%get_asFloats('h', requiredSize=size(N_sl))
prm%h = pl%get_asFloats('h', requiredSize=size(N_sl))
prm%w = pl%get_asFloats('w', requiredSize=size(N_sl))
prm%omega = pl%get_asFloats('omega', requiredSize=size(N_sl))
prm%B = pl%get_asFloats('B', requiredSize=size(N_sl))
@ -176,7 +176,7 @@ module function plastic_disloTungsten_init() result(myPlasticity)
prm%D = pl%get_asFloat('D')
prm%D_0 = pl%get_asFloat('D_0')
prm%Q_cl = pl%get_asFloat('Q_cl')
prm%atomicVolume = pl%get_asFloat('f_at') * prm%b_sl**3.0_pReal
prm%f_at = pl%get_asFloat('f_at') * prm%b_sl**3.0_pReal
prm%D_a = pl%get_asFloat('D_a') * prm%b_sl
prm%dipoleformation = pl%get_asBool('dipole_formation_factor', defaultVal = .true.)
@ -186,16 +186,16 @@ module function plastic_disloTungsten_init() result(myPlasticity)
rho_dip_0 = math_expand(rho_dip_0, N_sl)
prm%q = math_expand(prm%q, N_sl)
prm%p = math_expand(prm%p, N_sl)
prm%delta_F = math_expand(prm%delta_F, N_sl)
prm%Q_s = math_expand(prm%Q_s, N_sl)
prm%b_sl = math_expand(prm%b_sl, N_sl)
prm%kink_height = math_expand(prm%kink_height, N_sl)
prm%h = math_expand(prm%h, N_sl)
prm%w = math_expand(prm%w, N_sl)
prm%omega = math_expand(prm%omega, N_sl)
prm%tau_0 = math_expand(prm%tau_0, N_sl)
prm%v0 = math_expand(prm%v0, N_sl)
prm%tau_Peierls = math_expand(prm%tau_Peierls, N_sl)
prm%v_0 = math_expand(prm%v_0, N_sl)
prm%B = math_expand(prm%B, N_sl)
prm%i_sl = math_expand(prm%i_sl, N_sl)
prm%atomicVolume = math_expand(prm%atomicVolume, N_sl)
prm%f_at = math_expand(prm%f_at, N_sl)
prm%D_a = math_expand(prm%D_a, N_sl)
! sanity checks
@ -203,17 +203,17 @@ module function plastic_disloTungsten_init() result(myPlasticity)
if ( prm%Q_cl <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl'
if (any(rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_mob_0'
if (any(rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_dip_0'
if (any(prm%v0 < 0.0_pReal)) extmsg = trim(extmsg)//' v_0'
if (any(prm%v_0 < 0.0_pReal)) extmsg = trim(extmsg)//' v_0'
if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl'
if (any(prm%delta_F <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_s'
if (any(prm%tau_0 < 0.0_pReal)) extmsg = trim(extmsg)//' tau_peierls'
if (any(prm%Q_s <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_s'
if (any(prm%tau_Peierls < 0.0_pReal)) extmsg = trim(extmsg)//' tau_Peierls'
if (any(prm%D_a <= 0.0_pReal)) extmsg = trim(extmsg)//' D_a or b_sl'
if (any(prm%atomicVolume <= 0.0_pReal)) extmsg = trim(extmsg)//' f_at or b_sl'
if (any(prm%f_at <= 0.0_pReal)) extmsg = trim(extmsg)//' f_at or b_sl'
else slipActive
rho_mob_0= emptyRealArray; rho_dip_0 = emptyRealArray
allocate(prm%b_sl,prm%D_a,prm%i_sl,prm%atomicVolume,prm%tau_0, &
prm%delta_F,prm%v0,prm%p,prm%q,prm%B,prm%kink_height,prm%w,prm%omega, &
allocate(prm%b_sl,prm%D_a,prm%i_sl,prm%f_at,prm%tau_Peierls, &
prm%Q_s,prm%v_0,prm%p,prm%q,prm%B,prm%h,prm%w,prm%omega, &
source = emptyRealArray)
allocate(prm%forestProjection(0,0))
allocate(prm%h_sl_sl (0,0))
@ -354,7 +354,7 @@ module subroutine plastic_disloTungsten_dotState(Mp,T,instance,of)
dot_rho_dip_formation = merge(2.0_pReal*dip_distance* stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of))/prm%b_sl, & ! ToDo: ignore region of spontaneous annihilation
0.0_pReal, &
prm%dipoleformation)
v_cl = (3.0_pReal*prm%mu*VacancyDiffusion*prm%atomicVolume/(2.0_pReal*pi*kB*T)) &
v_cl = (3.0_pReal*prm%mu*VacancyDiffusion*prm%f_at/(2.0_pReal*pi*kB*T)) &
* (1.0_pReal/(dip_distance+prm%D_a))
dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,of))/(dip_distance-prm%D_a) ! ToDo: Discuss with Franz: Stress dependency?
end where
@ -477,12 +477,12 @@ pure subroutine kinetics(Mp,T,instance,of, &
if (present(tau_pos_out)) tau_pos_out = tau_pos
if (present(tau_neg_out)) tau_neg_out = tau_neg
associate(BoltzmannRatio => prm%delta_F/(kB*T), &
dot_gamma_0 => stt%rho_mob(:,of)*prm%b_sl*prm%v0, &
associate(BoltzmannRatio => prm%Q_s/(kB*T), &
dot_gamma_0 => stt%rho_mob(:,of)*prm%b_sl*prm%v_0, &
effectiveLength => dst%Lambda_sl(:,of) - prm%w)
significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check)
StressRatio = (abs(tau_pos)-dst%threshold_stress(:,of))/prm%tau_0
StressRatio = (abs(tau_pos)-dst%threshold_stress(:,of))/prm%tau_Peierls
StressRatio_p = StressRatio** prm%p
StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal)
needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)
@ -490,7 +490,7 @@ pure subroutine kinetics(Mp,T,instance,of, &
t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength)
t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos)
vel = prm%kink_height/(t_n + t_k)
vel = prm%h/(t_n + t_k)
dot_gamma_pos = dot_gamma_0 * sign(vel,tau_pos) * 0.5_pReal
else where significantPositiveTau
@ -500,10 +500,10 @@ pure subroutine kinetics(Mp,T,instance,of, &
if (present(ddot_gamma_dtau_pos)) then
significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check)
dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) &
* (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_0
* (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls
dtk = -1.0_pReal * t_k / tau_pos
dvel = -1.0_pReal * prm%kink_height * (dtk + dtn) / (t_n + t_k)**2.0_pReal
dvel = -1.0_pReal * prm%h * (dtk + dtn) / (t_n + t_k)**2.0_pReal
ddot_gamma_dtau_pos = dot_gamma_0 * dvel* 0.5_pReal
else where significantPositiveTau2
@ -512,7 +512,7 @@ pure subroutine kinetics(Mp,T,instance,of, &
endif
significantNegativeTau: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check)
StressRatio = (abs(tau_neg)-dst%threshold_stress(:,of))/prm%tau_0
StressRatio = (abs(tau_neg)-dst%threshold_stress(:,of))/prm%tau_Peierls
StressRatio_p = StressRatio** prm%p
StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal)
needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)
@ -520,7 +520,7 @@ pure subroutine kinetics(Mp,T,instance,of, &
t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength)
t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos)
vel = prm%kink_height/(t_n + t_k)
vel = prm%h/(t_n + t_k)
dot_gamma_neg = dot_gamma_0 * sign(vel,tau_neg) * 0.5_pReal
else where significantNegativeTau
@ -530,10 +530,10 @@ pure subroutine kinetics(Mp,T,instance,of, &
if (present(ddot_gamma_dtau_neg)) then
significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check)
dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) &
* (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_0
* (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls
dtk = -1.0_pReal * t_k / tau_neg
dvel = -1.0_pReal * prm%kink_height * (dtk + dtn) / (t_n + t_k)**2.0_pReal
dvel = -1.0_pReal * prm%h * (dtk + dtn) / (t_n + t_k)**2.0_pReal
ddot_gamma_dtau_neg = dot_gamma_0 * dvel * 0.5_pReal
else where significantNegativeTau2

View File

@ -16,38 +16,38 @@ submodule(constitutive:constitutive_plastic) plastic_dislotwin
real(pReal) :: &
mu = 1.0_pReal, & !< equivalent shear modulus
nu = 1.0_pReal, & !< equivalent shear Poisson's ratio
D0 = 1.0_pReal, & !< prefactor for self-diffusion coefficient
Qsd = 1.0_pReal, & !< activation energy for dislocation climb
D_0 = 1.0_pReal, & !< prefactor for self-diffusion coefficient
Q_cl = 1.0_pReal, & !< activation energy for dislocation climb
omega = 1.0_pReal, & !< frequency factor for dislocation climb
D = 1.0_pReal, & !< grain size
p_sb = 1.0_pReal, & !< p-exponent in shear band velocity
q_sb = 1.0_pReal, & !< q-exponent in shear band velocity
CEdgeDipMinDistance = 1.0_pReal, & !< adjustment parameter to calculate minimum dipole distance
D_a = 1.0_pReal, & !< adjustment parameter to calculate minimum dipole distance
i_tw = 1.0_pReal, & !< adjustment parameter to calculate MFP for twinning
tau_0 = 1.0_pReal, & !< strength due to elements in solid solution
L_tw = 1.0_pReal, & !< Length of twin nuclei in Burgers vectors
L_tr = 1.0_pReal, & !< Length of trans nuclei in Burgers vectors
xc_twin = 1.0_pReal, & !< critical distance for formation of twin nucleus
xc_trans = 1.0_pReal, & !< critical distance for formation of trans nucleus
x_c_tw = 1.0_pReal, & !< critical distance for formation of twin nucleus
x_c_tr = 1.0_pReal, & !< critical distance for formation of trans nucleus
V_cs = 1.0_pReal, & !< cross slip volume
sbResistance = 1.0_pReal, & !< value for shearband resistance
sbVelocity = 1.0_pReal, & !< value for shearband velocity_0
xi_sb = 1.0_pReal, & !< value for shearband resistance
v_sb = 1.0_pReal, & !< value for shearband velocity_0
E_sb = 1.0_pReal, & !< activation energy for shear bands
SFE_0K = 1.0_pReal, & !< stacking fault energy at zero K
dSFE_dT = 1.0_pReal, & !< temperature dependence of stacking fault energy
gamma_fcc_hex = 1.0_pReal, & !< Free energy difference between austensite and martensite
Gamma_sf_0K = 1.0_pReal, & !< stacking fault energy at zero K
dGamma_sf_dT = 1.0_pReal, & !< temperature dependence of stacking fault energy
delta_G = 1.0_pReal, & !< Free energy difference between austensite and martensite
i_tr = 1.0_pReal, & !< adjustment parameter to calculate MFP for transformation
h = 1.0_pReal !< Stack height of hex nucleus
real(pReal), allocatable, dimension(:) :: &
b_sl, & !< absolute length of burgers vector [m] for each slip system
b_tw, & !< absolute length of burgers vector [m] for each twin system
b_tr, & !< absolute length of burgers vector [m] for each transformation system
Delta_F,& !< activation energy for glide [J] for each slip system
v0, & !< dislocation velocity prefactor [m/s] for each slip system
b_sl, & !< absolute length of Burgers vector [m] for each slip system
b_tw, & !< absolute length of Burgers vector [m] for each twin system
b_tr, & !< absolute length of Burgers vector [m] for each transformation system
Q_s,& !< activation energy for glide [J] for each slip system
v_0, & !< dislocation velocity prefactor [m/s] for each slip system
dot_N_0_tw, & !< twin nucleation rate [1/m³s] for each twin system
dot_N_0_tr, & !< trans nucleation rate [1/m³s] for each trans system
t_tw, & !< twin thickness [m] for each twin system
CLambdaSlip, & !< Adj. parameter for distance between 2 forest dislocations for each slip system
i_sl, & !< Adj. parameter for distance between 2 forest dislocations for each slip system
t_tr, & !< martensite lamellar thickness [m] for each trans system and instance
p, & !< p-exponent in glide velocity
q, & !< q-exponent in glide velocity
@ -209,23 +209,23 @@ module function plastic_dislotwin_init() result(myPlasticity)
rho_mob_0 = pl%get_asFloats('rho_mob_0', requiredSize=size(N_sl))
rho_dip_0 = pl%get_asFloats('rho_dip_0', requiredSize=size(N_sl))
prm%v0 = pl%get_asFloats('v_0', requiredSize=size(N_sl))
prm%b_sl = pl%get_asFloats('b_sl', requiredSize=size(N_sl))
prm%Delta_F = pl%get_asFloats('Q_s', requiredSize=size(N_sl))
prm%CLambdaSlip = pl%get_asFloats('i_sl', requiredSize=size(N_sl))
prm%p = pl%get_asFloats('p_sl', requiredSize=size(N_sl))
prm%q = pl%get_asFloats('q_sl', requiredSize=size(N_sl))
prm%B = pl%get_asFloats('B', requiredSize=size(N_sl), &
prm%v_0 = pl%get_asFloats('v_0', requiredSize=size(N_sl))
prm%b_sl = pl%get_asFloats('b_sl', requiredSize=size(N_sl))
prm%Q_s = pl%get_asFloats('Q_s', requiredSize=size(N_sl))
prm%i_sl = pl%get_asFloats('i_sl', requiredSize=size(N_sl))
prm%p = pl%get_asFloats('p_sl', requiredSize=size(N_sl))
prm%q = pl%get_asFloats('q_sl', requiredSize=size(N_sl))
prm%B = pl%get_asFloats('B', requiredSize=size(N_sl), &
defaultVal=[(0.0_pReal, i=1,size(N_sl))])
prm%tau_0 = pl%get_asFloat('tau_0')
prm%CEdgeDipMinDistance = pl%get_asFloat('D_a')
prm%D0 = pl%get_asFloat('D_0')
prm%Qsd = pl%get_asFloat('Q_cl')
prm%D_a = pl%get_asFloat('D_a')
prm%D_0 = pl%get_asFloat('D_0')
prm%Q_cl = pl%get_asFloat('Q_cl')
prm%ExtendedDislocations = pl%get_asBool('extend_dislocations',defaultVal = .false.)
if (prm%ExtendedDislocations) then
prm%SFE_0K = pl%get_asFloat('Gamma_sf_0K')
prm%dSFE_dT = pl%get_asFloat('dGamma_sf_dT')
prm%Gamma_sf_0K = pl%get_asFloat('Gamma_sf_0K')
prm%dGamma_sf_dT = pl%get_asFloat('dGamma_sf_dT')
endif
prm%dipoleformation = .not. pl%get_asBool('no_dipole_formation',defaultVal = .false.)
@ -238,29 +238,29 @@ module function plastic_dislotwin_init() result(myPlasticity)
! expand: family => system
rho_mob_0 = math_expand(rho_mob_0, N_sl)
rho_dip_0 = math_expand(rho_dip_0, N_sl)
prm%v0 = math_expand(prm%v0, N_sl)
prm%v_0 = math_expand(prm%v_0, N_sl)
prm%b_sl = math_expand(prm%b_sl, N_sl)
prm%Delta_F = math_expand(prm%Delta_F, N_sl)
prm%CLambdaSlip = math_expand(prm%CLambdaSlip, N_sl)
prm%Q_s = math_expand(prm%Q_s, N_sl)
prm%i_sl = math_expand(prm%i_sl, N_sl)
prm%p = math_expand(prm%p, N_sl)
prm%q = math_expand(prm%q, N_sl)
prm%B = math_expand(prm%B, N_sl)
! sanity checks
if ( prm%D0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0'
if ( prm%Qsd <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl'
if (any(rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_mob_0'
if (any(rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_dip_0'
if (any(prm%v0 < 0.0_pReal)) extmsg = trim(extmsg)//' v_0'
if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl'
if (any(prm%Delta_F <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_s'
if (any(prm%CLambdaSlip <= 0.0_pReal)) extmsg = trim(extmsg)//' i_sl'
if (any(prm%B < 0.0_pReal)) extmsg = trim(extmsg)//' B'
if (any(prm%p<=0.0_pReal .or. prm%p>1.0_pReal)) extmsg = trim(extmsg)//' p_sl'
if (any(prm%q< 1.0_pReal .or. prm%q>2.0_pReal)) extmsg = trim(extmsg)//' q_sl'
if ( prm%D_0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0'
if ( prm%Q_cl <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl'
if (any(rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_mob_0'
if (any(rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_dip_0'
if (any(prm%v_0 < 0.0_pReal)) extmsg = trim(extmsg)//' v_0'
if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl'
if (any(prm%Q_s <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_s'
if (any(prm%i_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' i_sl'
if (any(prm%B < 0.0_pReal)) extmsg = trim(extmsg)//' B'
if (any(prm%p<=0.0_pReal .or. prm%p>1.0_pReal)) extmsg = trim(extmsg)//' p_sl'
if (any(prm%q< 1.0_pReal .or. prm%q>2.0_pReal)) extmsg = trim(extmsg)//' q_sl'
else slipActive
rho_mob_0 = emptyRealArray; rho_dip_0 = emptyRealArray
allocate(prm%b_sl,prm%Delta_F,prm%v0,prm%CLambdaSlip,prm%p,prm%q,prm%B,source=emptyRealArray)
allocate(prm%b_sl,prm%Q_s,prm%v_0,prm%i_sl,prm%p,prm%q,prm%B,source=emptyRealArray)
allocate(prm%forestProjection(0,0),prm%h_sl_sl(0,0))
endif slipActive
@ -279,7 +279,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%t_tw = pl%get_asFloats('t_tw', requiredSize=size(N_tw))
prm%r = pl%get_asFloats('p_tw', requiredSize=size(N_tw))
prm%xc_twin = pl%get_asFloat('x_c_tw')
prm%x_c_tw = pl%get_asFloat('x_c_tw')
prm%L_tw = pl%get_asFloat('L_tw')
prm%i_tw = pl%get_asFloat('i_tw')
@ -300,7 +300,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%r = math_expand(prm%r,N_tw)
! sanity checks
if ( prm%xc_twin < 0.0_pReal) extmsg = trim(extmsg)//' x_c_twin'
if ( prm%x_c_tw < 0.0_pReal) extmsg = trim(extmsg)//' x_c_twin'
if ( prm%L_tw < 0.0_pReal) extmsg = trim(extmsg)//' L_tw'
if ( prm%i_tw < 0.0_pReal) extmsg = trim(extmsg)//' i_tw'
if (any(prm%b_tw < 0.0_pReal)) extmsg = trim(extmsg)//' b_tw'
@ -324,8 +324,8 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%h = pl%get_asFloat('h', defaultVal=0.0_pReal) ! ToDo: How to handle that???
prm%i_tr = pl%get_asFloat('i_tr', defaultVal=0.0_pReal) ! ToDo: How to handle that???
prm%gamma_fcc_hex = pl%get_asFloat('delta_G')
prm%xc_trans = pl%get_asFloat('x_c_tr', defaultVal=0.0_pReal) ! ToDo: How to handle that???
prm%delta_G = pl%get_asFloat('delta_G')
prm%x_c_tr = pl%get_asFloat('x_c_tr', defaultVal=0.0_pReal) ! ToDo: How to handle that???
prm%L_tr = pl%get_asFloat('L_tr')
prm%h_tr_tr = lattice_interaction_TransByTrans(N_tr,pl%get_asFloats('h_tr_tr'), &
@ -351,7 +351,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%s = math_expand(prm%s,N_tr)
! sanity checks
if ( prm%xc_trans < 0.0_pReal) extmsg = trim(extmsg)//' x_c_trans'
if ( prm%x_c_tr < 0.0_pReal) extmsg = trim(extmsg)//' x_c_trans'
if ( prm%L_tr < 0.0_pReal) extmsg = trim(extmsg)//' L_tr'
if ( prm%i_tr < 0.0_pReal) extmsg = trim(extmsg)//' i_tr'
if (any(prm%t_tr < 0.0_pReal)) extmsg = trim(extmsg)//' t_tr'
@ -366,15 +366,15 @@ module function plastic_dislotwin_init() result(myPlasticity)
!--------------------------------------------------------------------------------------------------
! shearband related parameters
prm%sbVelocity = pl%get_asFloat('v_sb',defaultVal=0.0_pReal)
if (prm%sbVelocity > 0.0_pReal) then
prm%sbResistance = pl%get_asFloat('xi_sb')
prm%v_sb = pl%get_asFloat('v_sb',defaultVal=0.0_pReal)
if (prm%v_sb > 0.0_pReal) then
prm%xi_sb = pl%get_asFloat('xi_sb')
prm%E_sb = pl%get_asFloat('Q_sb')
prm%p_sb = pl%get_asFloat('p_sb')
prm%q_sb = pl%get_asFloat('q_sb')
! sanity checks
if (prm%sbResistance < 0.0_pReal) extmsg = trim(extmsg)//' xi_sb'
if (prm%xi_sb < 0.0_pReal) extmsg = trim(extmsg)//' xi_sb'
if (prm%E_sb < 0.0_pReal) extmsg = trim(extmsg)//' Q_sb'
if (prm%p_sb <= 0.0_pReal) extmsg = trim(extmsg)//' p_sb'
if (prm%q_sb <= 0.0_pReal) extmsg = trim(extmsg)//' q_sb'
@ -386,8 +386,8 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%D = pl%get_asFloat('D')
twinOrSlipActive: if (prm%sum_N_tw + prm%sum_N_tr > 0) then
prm%SFE_0K = pl%get_asFloat('Gamma_sf_0K')
prm%dSFE_dT = pl%get_asFloat('dGamma_sf_dT')
prm%Gamma_sf_0K = pl%get_asFloat('Gamma_sf_0K')
prm%dGamma_sf_dT = pl%get_asFloat('dGamma_sf_dT')
prm%V_cs = pl%get_asFloat('V_cs')
endif twinOrSlipActive
@ -602,7 +602,7 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
Lp = Lp * f_unrotated
dLp_dMp = dLp_dMp * f_unrotated
shearBandingContribution: if(dNeq0(prm%sbVelocity)) then
shearBandingContribution: if(dNeq0(prm%v_sb)) then
BoltzmannRatio = prm%E_sb/(kB*T)
call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design?
@ -613,10 +613,10 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
tau = math_tensordot(Mp,P_sb)
significantShearBandStress: if (abs(tau) > tol_math_check) then
StressRatio_p = (abs(tau)/prm%sbResistance)**prm%p_sb
dot_gamma_sb = sign(prm%sbVelocity*exp(-BoltzmannRatio*(1-StressRatio_p)**prm%q_sb), tau)
ddot_gamma_dtau = abs(dot_gamma_sb)*BoltzmannRatio* prm%p_sb*prm%q_sb/ prm%sbResistance &
* (abs(tau)/prm%sbResistance)**(prm%p_sb-1.0_pReal) &
StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb
dot_gamma_sb = sign(prm%v_sb*exp(-BoltzmannRatio*(1-StressRatio_p)**prm%q_sb), tau)
ddot_gamma_dtau = abs(dot_gamma_sb)*BoltzmannRatio* prm%p_sb*prm%q_sb/ prm%xi_sb &
* (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) &
* (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal)
Lp = Lp + dot_gamma_sb * P_sb
@ -654,7 +654,7 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of)
Gamma, & !< stacking fault energy
tau, &
sigma_cl, & !< climb stress
b_d !< ratio of burgers vector to stacking fault width
b_d !< ratio of Burgers vector to stacking fault width
real(pReal), dimension(param(instance)%sum_N_sl) :: &
dot_rho_dip_formation, &
dot_rho_dip_climb, &
@ -675,7 +675,7 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of)
call kinetics_slip(Mp,T,instance,of,dot_gamma_sl)
dot%gamma_sl(:,of) = abs(dot_gamma_sl)
rho_dip_distance_min = prm%CEdgeDipMinDistance*prm%b_sl
rho_dip_distance_min = prm%D_a*prm%b_sl
slipState: do i = 1, prm%sum_N_sl
tau = math_tensordot(Mp,prm%P_sl(1:3,1:3,i))
@ -701,12 +701,12 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of)
!@details: Refer: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981
sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i)))
if (prm%ExtendedDislocations) then
Gamma = prm%SFE_0K + prm%dSFE_dT * T
Gamma = prm%Gamma_sf_0K + prm%dGamma_sf_dT * T
b_d = 24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu)* Gamma/(prm%mu*prm%b_sl(i))
else
b_d = 1.0_pReal
endif
v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Qsd/(kB*T)) &
v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(kB*T)) &
* (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal)
dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,of) &
@ -768,7 +768,7 @@ module subroutine plastic_dislotwin_dependentState(T,instance,of)
sumf_twin = sum(stt%f_tw(1:prm%sum_N_tw,of))
sumf_trans = sum(stt%f_tr(1:prm%sum_N_tr,of))
Gamma = prm%SFE_0K + prm%dSFE_dT * T
Gamma = prm%Gamma_sf_0K + prm%dGamma_sf_dT * T
!* rescaled volume fraction for topology
f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,of)/prm%t_tw ! this is per system ...
@ -776,7 +776,7 @@ module subroutine plastic_dislotwin_dependentState(T,instance,of)
! ToDo ...Physically correct, but naming could be adjusted
inv_lambda_sl_sl = sqrt(matmul(prm%forestProjection, &
stt%rho_mob(:,of)+stt%rho_dip(:,of)))/prm%CLambdaSlip
stt%rho_mob(:,of)+stt%rho_dip(:,of)))/prm%i_sl
if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) &
inv_lambda_sl_tw = matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_twin)
@ -805,21 +805,21 @@ module subroutine plastic_dislotwin_dependentState(T,instance,of)
!* threshold stress for growing twin/martensite
if(prm%sum_N_tw == prm%sum_N_sl) &
dst%tau_hat_tw(:,of) = Gamma/(3.0_pReal*prm%b_tw) &
+ 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_sl) ! slip burgers here correct?
+ 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_sl) ! slip Burgers here correct?
if(prm%sum_N_tr == prm%sum_N_sl) &
dst%tau_hat_tr(:,of) = Gamma/(3.0_pReal*prm%b_tr) &
+ 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip burgers here correct?
+ prm%h*prm%gamma_fcc_hex/ (3.0_pReal*prm%b_tr)
+ 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip Burgers here correct?
+ prm%h*prm%delta_G/ (3.0_pReal*prm%b_tr)
dst%V_tw(:,of) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,of)**2.0_pReal
dst%V_tr(:,of) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,of)**2.0_pReal
x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the burgers vector for slip and is the same for twin and trans
dst%tau_r_tw(:,of) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%xc_twin)+cos(pi/3.0_pReal)/x0)
x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip and is the same for twin and trans
dst%tau_r_tw(:,of) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0)
x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the burgers vector for slip
dst%tau_r_tr(:,of) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%xc_trans)+cos(pi/3.0_pReal)/x0)
x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip
dst%tau_r_tr(:,of) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0)
end associate
@ -927,8 +927,8 @@ pure subroutine kinetics_slip(Mp,T,instance,of, &
significantStress: where(tau_eff > tol_math_check)
stressRatio = tau_eff/prm%tau_0
StressRatio_p = stressRatio** prm%p
BoltzmannRatio = prm%Delta_F/(kB*T)
v_wait_inverse = prm%v0**(-1.0_pReal) * exp(BoltzmannRatio*(1.0_pReal-StressRatio_p)** prm%q)
BoltzmannRatio = prm%Q_s/(kB*T)
v_wait_inverse = prm%v_0**(-1.0_pReal) * exp(BoltzmannRatio*(1.0_pReal-StressRatio_p)** prm%q)
v_run_inverse = prm%B/(tau_eff*prm%b_sl)
dot_gamma_sl = sign(stt%rho_mob(:,of)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau)

View File

@ -14,7 +14,7 @@ submodule(constitutive:constitutive_plastic) plastic_isotropic
M, & !< Taylor factor
dot_gamma_0, & !< reference strain rate
n, & !< stress exponent
h0, &
h_0, &
h_ln, &
xi_inf, & !< maximum critical stress
a, &
@ -109,7 +109,7 @@ module function plastic_isotropic_init() result(myPlasticity)
prm%xi_inf = pl%get_asFloat('xi_inf')
prm%dot_gamma_0 = pl%get_asFloat('dot_gamma_0')
prm%n = pl%get_asFloat('n')
prm%h0 = pl%get_asFloat('h_0')
prm%h_0 = pl%get_asFloat('h_0')
prm%M = pl%get_asFloat('M')
prm%h_ln = pl%get_asFloat('h_ln', defaultVal=0.0_pReal)
prm%c_1 = pl%get_asFloat('c_1', defaultVal=0.0_pReal)
@ -310,7 +310,7 @@ module subroutine plastic_isotropic_dotState(Mp,instance,of)
/ prm%c_4 * (dot_gamma / prm%dot_gamma_0)**(1.0_pReal / prm%n)
endif
dot%xi(of) = dot_gamma &
* ( prm%h0 + prm%h_ln * log(dot_gamma) ) &
* ( prm%h_0 + prm%h_ln * log(dot_gamma) ) &
* abs( 1.0_pReal - stt%xi(of)/xi_inf_star )**prm%a &
* sign(1.0_pReal, 1.0_pReal - stt%xi(of)/xi_inf_star)
else

View File

@ -9,15 +9,15 @@ submodule(constitutive:constitutive_plastic) plastic_kinehardening
type :: tParameters
real(pReal) :: &
gdot0 = 1.0_pReal, & !< reference shear strain rate for slip
n = 1.0_pReal !< stress exponent for slip
n = 1.0_pReal, & !< stress exponent for slip
dot_gamma_0 = 1.0_pReal !< reference shear strain rate for slip
real(pReal), allocatable, dimension(:) :: &
theta0, & !< initial hardening rate of forward stress for each slip
theta1, & !< asymptotic hardening rate of forward stress for each slip
theta0_b, & !< initial hardening rate of back stress for each slip
theta1_b, & !< asymptotic hardening rate of back stress for each slip
tau1, &
tau1_b
h_0_f, & !< initial hardening rate of forward stress for each slip
h_inf_f, & !< asymptotic hardening rate of forward stress for each slip
h_0_b, & !< initial hardening rate of back stress for each slip
h_inf_b, & !< asymptotic hardening rate of back stress for each slip
xi_inf_f, &
xi_inf_b
real(pReal), allocatable, dimension(:,:) :: &
interaction_slipslip !< slip resistance from slip activity
real(pReal), allocatable, dimension(:,:,:) :: &
@ -125,7 +125,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
if(trim(phase%get_asString('lattice')) == 'bcc') then
a = pl%get_asFloats('nonSchmid_coefficients',defaultVal = emptyRealArray)
a = pl%get_asFloats('a_nonSchmid',defaultVal = emptyRealArray)
if(size(a) > 0) prm%nonSchmidActive = .true.
prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,a,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,a,-1)
@ -137,38 +137,38 @@ module function plastic_kinehardening_init() result(myPlasticity)
pl%get_asFloats('h_sl_sl'), &
phase%get_asString('lattice'))
xi_0 = pl%get_asFloats('xi_0', requiredSize=size(N_sl))
prm%tau1 = pl%get_asFloats('xi_inf_f', requiredSize=size(N_sl))
prm%tau1_b = pl%get_asFloats('xi_inf_b', requiredSize=size(N_sl))
prm%theta0 = pl%get_asFloats('h_0_f', requiredSize=size(N_sl))
prm%theta1 = pl%get_asFloats('h_inf_f', requiredSize=size(N_sl))
prm%theta0_b = pl%get_asFloats('h_0_b', requiredSize=size(N_sl))
prm%theta1_b = pl%get_asFloats('h_inf_b', requiredSize=size(N_sl))
xi_0 = pl%get_asFloats('xi_0', requiredSize=size(N_sl))
prm%xi_inf_f = pl%get_asFloats('xi_inf_f', requiredSize=size(N_sl))
prm%xi_inf_b = pl%get_asFloats('xi_inf_b', requiredSize=size(N_sl))
prm%h_0_f = pl%get_asFloats('h_0_f', requiredSize=size(N_sl))
prm%h_inf_f = pl%get_asFloats('h_inf_f', requiredSize=size(N_sl))
prm%h_0_b = pl%get_asFloats('h_0_b', requiredSize=size(N_sl))
prm%h_inf_b = pl%get_asFloats('h_inf_b', requiredSize=size(N_sl))
prm%gdot0 = pl%get_asFloat('dot_gamma_0')
prm%n = pl%get_asFloat('n')
prm%dot_gamma_0 = pl%get_asFloat('dot_gamma_0')
prm%n = pl%get_asFloat('n')
! expand: family => system
xi_0 = math_expand(xi_0, N_sl)
prm%tau1 = math_expand(prm%tau1, N_sl)
prm%tau1_b = math_expand(prm%tau1_b, N_sl)
prm%theta0 = math_expand(prm%theta0, N_sl)
prm%theta1 = math_expand(prm%theta1, N_sl)
prm%theta0_b = math_expand(prm%theta0_b,N_sl)
prm%theta1_b = math_expand(prm%theta1_b,N_sl)
xi_0 = math_expand(xi_0, N_sl)
prm%xi_inf_f = math_expand(prm%xi_inf_f, N_sl)
prm%xi_inf_b = math_expand(prm%xi_inf_b, N_sl)
prm%h_0_f = math_expand(prm%h_0_f, N_sl)
prm%h_inf_f = math_expand(prm%h_inf_f, N_sl)
prm%h_0_b = math_expand(prm%h_0_b, N_sl)
prm%h_inf_b = math_expand(prm%h_inf_b, N_sl)
!--------------------------------------------------------------------------------------------------
! sanity checks
if ( prm%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' dot_gamma_0'
if ( prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' n'
if (any(xi_0 <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_0'
if (any(prm%tau1 <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_inf_f'
if (any(prm%tau1_b <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_inf_b'
if ( prm%dot_gamma_0 <= 0.0_pReal) extmsg = trim(extmsg)//' dot_gamma_0'
if ( prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' n'
if (any(xi_0 <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_0'
if (any(prm%xi_inf_f <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_inf_f'
if (any(prm%xi_inf_b <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_inf_b'
!ToDo: Any sensible checks for theta?
else slipActive
xi_0 = emptyRealArray
allocate(prm%tau1,prm%tau1_b,prm%theta0,prm%theta1,prm%theta0_b,prm%theta1_b,source=emptyRealArray)
allocate(prm%xi_inf_f,prm%xi_inf_b,prm%h_0_f,prm%h_inf_f,prm%h_0_b,prm%h_inf_b,source=emptyRealArray)
allocate(prm%interaction_SlipSlip(0,0))
endif slipActive
@ -303,16 +303,16 @@ module subroutine plastic_kinehardening_dotState(Mp,instance,of)
dot%crss(:,of) = matmul(prm%interaction_SlipSlip,dot%accshear(:,of)) &
* ( prm%theta1 &
+ (prm%theta0 - prm%theta1 + prm%theta0*prm%theta1*sumGamma/prm%tau1) &
* exp(-sumGamma*prm%theta0/prm%tau1) &
* ( prm%h_inf_f &
+ (prm%h_0_f - prm%h_inf_f + prm%h_0_f*prm%h_inf_f*sumGamma/prm%xi_inf_f) &
* exp(-sumGamma*prm%h_0_f/prm%xi_inf_f) &
)
dot%crss_back(:,of) = stt%sense(:,of)*dot%accshear(:,of) * &
( prm%theta1_b + &
(prm%theta0_b - prm%theta1_b &
+ prm%theta0_b*prm%theta1_b/(prm%tau1_b+stt%chi0(:,of))*(stt%accshear(:,of)-stt%gamma0(:,of))&
) *exp(-(stt%accshear(:,of)-stt%gamma0(:,of)) *prm%theta0_b/(prm%tau1_b+stt%chi0(:,of))) &
( prm%h_inf_b + &
(prm%h_0_b - prm%h_inf_b &
+ prm%h_0_b*prm%h_inf_b/(prm%xi_inf_b+stt%chi0(:,of))*(stt%accshear(:,of)-stt%gamma0(:,of))&
) *exp(-(stt%accshear(:,of)-stt%gamma0(:,of)) *prm%h_0_b/(prm%xi_inf_b+stt%chi0(:,of))) &
)
end associate
@ -442,14 +442,14 @@ pure subroutine kinetics(Mp,instance,of, &
enddo
where(dNeq0(tau_pos))
gdot_pos = prm%gdot0 * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active
gdot_pos = prm%dot_gamma_0 * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active
* sign(abs(tau_pos/stt%crss(:,of))**prm%n, tau_pos)
else where
gdot_pos = 0.0_pReal
end where
where(dNeq0(tau_neg))
gdot_neg = prm%gdot0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2
gdot_neg = prm%dot_gamma_0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2
* sign(abs(tau_neg/stt%crss(:,of))**prm%n, tau_neg)
else where
gdot_neg = 0.0_pReal

View File

@ -46,58 +46,58 @@ submodule(constitutive:constitutive_plastic) plastic_nonlocal
type :: tInitialParameters !< container type for internal constitutive parameters
real(pReal) :: &
rhoSglScatter, & !< standard deviation of scatter in initial dislocation density
rhoSglRandom, &
rhoSglRandomBinning
sigma_rho_u, & !< standard deviation of scatter in initial dislocation density
random_rho_u, &
random_rho_u_binning
real(pReal), dimension(:), allocatable :: &
rhoSglEdgePos0, & !< initial edge_pos dislocation density
rhoSglEdgeNeg0, & !< initial edge_neg dislocation density
rhoSglScrewPos0, & !< initial screw_pos dislocation density
rhoSglScrewNeg0, & !< initial screw_neg dislocation density
rhoDipEdge0, & !< initial edge dipole dislocation density
rhoDipScrew0 !< initial screw dipole dislocation density
rho_u_ed_pos_0, & !< initial edge_pos dislocation density
rho_u_ed_neg_0, & !< initial edge_neg dislocation density
rho_u_sc_pos_0, & !< initial screw_pos dislocation density
rho_u_sc_neg_0, & !< initial screw_neg dislocation density
rho_d_ed_0, & !< initial edge dipole dislocation density
rho_d_sc_0 !< initial screw dipole dislocation density
integer, dimension(:), allocatable :: &
N_sl
end type tInitialParameters
type :: tParameters !< container type for internal constitutive parameters
real(pReal) :: &
atomicVolume, & !< atomic volume
Dsd0, & !< prefactor for self-diffusion coefficient
selfDiffusionEnergy, & !< activation enthalpy for diffusion
V_at, & !< atomic volume
D_0, & !< prefactor for self-diffusion coefficient
Q_cl, & !< activation enthalpy for diffusion
atol_rho, & !< absolute tolerance for dislocation density in state integration
significantRho, & !< density considered significant
significantN, & !< number of dislocations considered significant
doublekinkwidth, & !< width of a doubkle kink in multiples of the burgers vector length b
solidSolutionEnergy, & !< activation energy for solid solution in J
solidSolutionSize, & !< solid solution obstacle size in multiples of the burgers vector length
solidSolutionConcentration, & !< concentration of solid solution in atomic parts
rho_significant, & !< density considered significant
rho_min, & !< number of dislocations considered significant
w, & !< width of a doubkle kink in multiples of the Burgers vector length b
Q_sol, & !< activation energy for solid solution in J
f_sol, & !< solid solution obstacle size in multiples of the Burgers vector length
c_sol, & !< concentration of solid solution in atomic parts
p, & !< parameter for kinetic law (Kocks,Argon,Ashby)
q, & !< parameter for kinetic law (Kocks,Argon,Ashby)
viscosity, & !< viscosity for dislocation glide in Pa s
fattack, & !< attack frequency in Hz
surfaceTransmissivity, & !< transmissivity at free surface
grainboundaryTransmissivity, & !< transmissivity at grain boundary (identified by different texture)
CFLfactor, & !< safety factor for CFL flux condition
fEdgeMultiplication, & !< factor that determines how much edge dislocations contribute to multiplication (0...1)
linetensionEffect, &
edgeJogFactor, &
eta, & !< viscosity for dislocation glide in Pa s
nu_a, & !< attack frequency in Hz
chi_surface, & !< transmissivity at free surface
chi_GB, & !< transmissivity at grain boundary (identified by different texture)
f_c, & !< safety factor for CFL flux condition
f_ed_mult, & !< factor that determines how much edge dislocations contribute to multiplication (0...1)
f_F, &
f_ed, &
mu, &
nu
real(pReal), dimension(:), allocatable :: &
minDipoleHeight_edge, & !< minimum stable edge dipole height
minDipoleHeight_screw, & !< minimum stable screw dipole height
peierlsstress_edge, &
peierlsstress_screw, &
lambda0, & !< mean free path prefactor for each
burgers !< absolute length of burgers vector [m]
d_ed, & !< minimum stable edge dipole height
d_sc, & !< minimum stable screw dipole height
tau_Peierls_ed, &
tau_Peierls_sc, &
i_sl, & !< mean free path prefactor for each
b_sl !< absolute length of Burgers vector [m]
real(pReal), dimension(:,:), allocatable :: &
slip_normal, &
slip_direction, &
slip_transverse, &
minDipoleHeight, & ! edge and screw
peierlsstress, & ! edge and screw
interactionSlipSlip ,& !< coefficients for slip-slip interaction
h_sl_sl ,& !< coefficients for slip-slip interaction
forestProjection_Edge, & !< matrix of forest projections of edge dislocations
forestProjection_Screw !< matrix of forest projections of screw dislocations
real(pReal), dimension(:,:,:), allocatable :: &
@ -244,7 +244,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
if(trim(phase%get_asString('lattice')) == 'bcc') then
a = pl%get_asFloats('nonSchmid_coefficients',defaultVal = emptyRealArray)
a = pl%get_asFloats('a_nonSchmid',defaultVal = emptyRealArray)
if(size(a) > 0) prm%nonSchmidActive = .true.
prm%nonSchmid_pos = lattice_nonSchmidMatrix(ini%N_sl,a,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(ini%N_sl,a,-1)
@ -253,9 +253,9 @@ module function plastic_nonlocal_init() result(myPlasticity)
prm%nonSchmid_neg = prm%Schmid
endif
prm%interactionSlipSlip = lattice_interaction_SlipBySlip(ini%N_sl, &
pl%get_asFloats('h_sl_sl'), &
phase%get_asString('lattice'))
prm%h_sl_sl = lattice_interaction_SlipBySlip(ini%N_sl, &
pl%get_asFloats('h_sl_sl'), &
phase%get_asString('lattice'))
prm%forestProjection_edge = lattice_forestProjection_edge (ini%N_sl,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
@ -279,113 +279,113 @@ module function plastic_nonlocal_init() result(myPlasticity)
enddo
enddo
ini%rhoSglEdgePos0 = pl%get_asFloats('rho_u_ed_pos_0', requiredSize=size(ini%N_sl))
ini%rhoSglEdgeNeg0 = pl%get_asFloats('rho_u_ed_neg_0', requiredSize=size(ini%N_sl))
ini%rhoSglScrewPos0 = pl%get_asFloats('rho_u_sc_pos_0', requiredSize=size(ini%N_sl))
ini%rhoSglScrewNeg0 = pl%get_asFloats('rho_u_sc_neg_0', requiredSize=size(ini%N_sl))
ini%rhoDipEdge0 = pl%get_asFloats('rho_d_ed_0', requiredSize=size(ini%N_sl))
ini%rhoDipScrew0 = pl%get_asFloats('rho_d_sc_0', requiredSize=size(ini%N_sl))
ini%rho_u_ed_pos_0 = pl%get_asFloats('rho_u_ed_pos_0', requiredSize=size(ini%N_sl))
ini%rho_u_ed_neg_0 = pl%get_asFloats('rho_u_ed_neg_0', requiredSize=size(ini%N_sl))
ini%rho_u_sc_pos_0 = pl%get_asFloats('rho_u_sc_pos_0', requiredSize=size(ini%N_sl))
ini%rho_u_sc_neg_0 = pl%get_asFloats('rho_u_sc_neg_0', requiredSize=size(ini%N_sl))
ini%rho_d_ed_0 = pl%get_asFloats('rho_d_ed_0', requiredSize=size(ini%N_sl))
ini%rho_d_sc_0 = pl%get_asFloats('rho_d_sc_0', requiredSize=size(ini%N_sl))
prm%lambda0 = pl%get_asFloats('i_sl', requiredSize=size(ini%N_sl))
prm%burgers = pl%get_asFloats('b_sl', requiredSize=size(ini%N_sl))
prm%i_sl = pl%get_asFloats('i_sl', requiredSize=size(ini%N_sl))
prm%b_sl = pl%get_asFloats('b_sl', requiredSize=size(ini%N_sl))
prm%lambda0 = math_expand(prm%lambda0,ini%N_sl)
prm%burgers = math_expand(prm%burgers,ini%N_sl)
prm%i_sl = math_expand(prm%i_sl,ini%N_sl)
prm%b_sl = math_expand(prm%b_sl,ini%N_sl)
prm%minDipoleHeight_edge = pl%get_asFloats('d_ed', requiredSize=size(ini%N_sl))
prm%minDipoleHeight_screw = pl%get_asFloats('d_sc', requiredSize=size(ini%N_sl))
prm%minDipoleHeight_edge = math_expand(prm%minDipoleHeight_edge, ini%N_sl)
prm%minDipoleHeight_screw = math_expand(prm%minDipoleHeight_screw,ini%N_sl)
prm%d_ed = pl%get_asFloats('d_ed', requiredSize=size(ini%N_sl))
prm%d_sc = pl%get_asFloats('d_sc', requiredSize=size(ini%N_sl))
prm%d_ed = math_expand(prm%d_ed,ini%N_sl)
prm%d_sc = math_expand(prm%d_sc,ini%N_sl)
allocate(prm%minDipoleHeight(prm%sum_N_sl,2))
prm%minDipoleHeight(:,1) = prm%minDipoleHeight_edge
prm%minDipoleHeight(:,2) = prm%minDipoleHeight_screw
prm%minDipoleHeight(:,1) = prm%d_ed
prm%minDipoleHeight(:,2) = prm%d_sc
prm%peierlsstress_edge = pl%get_asFloats('tau_peierls_ed', requiredSize=size(ini%N_sl))
prm%peierlsstress_screw = pl%get_asFloats('tau_peierls_sc', requiredSize=size(ini%N_sl))
prm%peierlsstress_edge = math_expand(prm%peierlsstress_edge, ini%N_sl)
prm%peierlsstress_screw = math_expand(prm%peierlsstress_screw,ini%N_sl)
prm%tau_Peierls_ed = pl%get_asFloats('tau_Peierls_ed', requiredSize=size(ini%N_sl))
prm%tau_Peierls_sc = pl%get_asFloats('tau_Peierls_sc', requiredSize=size(ini%N_sl))
prm%tau_Peierls_ed = math_expand(prm%tau_Peierls_ed,ini%N_sl)
prm%tau_Peierls_sc = math_expand(prm%tau_Peierls_sc,ini%N_sl)
allocate(prm%peierlsstress(prm%sum_N_sl,2))
prm%peierlsstress(:,1) = prm%peierlsstress_edge
prm%peierlsstress(:,2) = prm%peierlsstress_screw
prm%peierlsstress(:,1) = prm%tau_Peierls_ed
prm%peierlsstress(:,2) = prm%tau_Peierls_sc
prm%significantRho = pl%get_asFloat('rho_significant')
prm%significantN = pl%get_asFloat('rho_num_significant', 0.0_pReal)
prm%CFLfactor = pl%get_asFloat('f_c',defaultVal=2.0_pReal)
prm%rho_significant = pl%get_asFloat('rho_significant')
prm%rho_min = pl%get_asFloat('rho_min', 0.0_pReal)
prm%f_c = pl%get_asFloat('f_c',defaultVal=2.0_pReal)
prm%atomicVolume = pl%get_asFloat('V_at')
prm%Dsd0 = pl%get_asFloat('D_0') !,'dsd0'
prm%selfDiffusionEnergy = pl%get_asFloat('Q_cl') !,'qsd'
prm%linetensionEffect = pl%get_asFloat('f_F')
prm%edgeJogFactor = pl%get_asFloat('f_ed') !,'edgejogs'
prm%doublekinkwidth = pl%get_asFloat('w')
prm%solidSolutionEnergy = pl%get_asFloat('Q_sol')
prm%solidSolutionSize = pl%get_asFloat('f_sol')
prm%solidSolutionConcentration = pl%get_asFloat('c_sol')
prm%V_at = pl%get_asFloat('V_at')
prm%D_0 = pl%get_asFloat('D_0')
prm%Q_cl = pl%get_asFloat('Q_cl')
prm%f_F = pl%get_asFloat('f_F')
prm%f_ed = pl%get_asFloat('f_ed') !,'edgejogs'
prm%w = pl%get_asFloat('w')
prm%Q_sol = pl%get_asFloat('Q_sol')
prm%f_sol = pl%get_asFloat('f_sol')
prm%c_sol = pl%get_asFloat('c_sol')
prm%p = pl%get_asFloat('p_sl')
prm%q = pl%get_asFloat('q_sl')
prm%viscosity = pl%get_asFloat('eta')
prm%fattack = pl%get_asFloat('nu_a')
prm%p = pl%get_asFloat('p_sl')
prm%q = pl%get_asFloat('q_sl')
prm%eta = pl%get_asFloat('eta')
prm%nu_a = pl%get_asFloat('nu_a')
! ToDo: discuss logic
ini%rhoSglScatter = pl%get_asFloat('sigma_rho_u')
ini%rhoSglRandom = pl%get_asFloat('random_rho_u',defaultVal= 0.0_pReal)
ini%sigma_rho_u = pl%get_asFloat('sigma_rho_u')
ini%random_rho_u = pl%get_asFloat('random_rho_u',defaultVal= 0.0_pReal)
if (pl%contains('random_rho_u')) &
ini%rhoSglRandomBinning = pl%get_asFloat('random_rho_u_binning',defaultVal=0.0_pReal) !ToDo: useful default?
ini%random_rho_u_binning = pl%get_asFloat('random_rho_u_binning',defaultVal=0.0_pReal) !ToDo: useful default?
! if (rhoSglRandom(instance) < 0.0_pReal) &
! if (rhoSglRandomBinning(instance) <= 0.0_pReal) &
prm%surfaceTransmissivity = pl%get_asFloat('chi_surface',defaultVal=1.0_pReal)
prm%grainboundaryTransmissivity = pl%get_asFloat('chi_GB', defaultVal=-1.0_pReal)
prm%fEdgeMultiplication = pl%get_asFloat('f_ed_mult')
prm%chi_surface = pl%get_asFloat('chi_surface',defaultVal=1.0_pReal)
prm%chi_GB = pl%get_asFloat('chi_GB', defaultVal=-1.0_pReal)
prm%f_ed_mult = pl%get_asFloat('f_ed_mult')
prm%shortRangeStressCorrection = pl%get_asBool('short_range_stress_correction', defaultVal = .false.)
!--------------------------------------------------------------------------------------------------
! sanity checks
if (any(prm%burgers < 0.0_pReal)) extmsg = trim(extmsg)//' b_sl'
if (any(prm%lambda0 <= 0.0_pReal)) extmsg = trim(extmsg)//' i_sl'
if (any(prm%b_sl < 0.0_pReal)) extmsg = trim(extmsg)//' b_sl'
if (any(prm%i_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' i_sl'
if (any(ini%rhoSglEdgePos0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_ed_pos_0'
if (any(ini%rhoSglEdgeNeg0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_ed_neg_0'
if (any(ini%rhoSglScrewPos0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_sc_pos_0'
if (any(ini%rhoSglScrewNeg0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_sc_neg_0'
if (any(ini%rhoDipEdge0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_d_ed_0'
if (any(ini%rhoDipScrew0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_d_sc_0'
if (any(ini%rho_u_ed_pos_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_ed_pos_0'
if (any(ini%rho_u_ed_neg_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_ed_neg_0'
if (any(ini%rho_u_sc_pos_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_sc_pos_0'
if (any(ini%rho_u_sc_neg_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_sc_neg_0'
if (any(ini%rho_d_ed_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_d_ed_0'
if (any(ini%rho_d_sc_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_d_sc_0'
if (any(prm%peierlsstress < 0.0_pReal)) extmsg = trim(extmsg)//' tau_peierls'
if (any(prm%minDipoleHeight < 0.0_pReal)) extmsg = trim(extmsg)//' d_ed or d_sc'
if (prm%viscosity <= 0.0_pReal) extmsg = trim(extmsg)//' eta'
if (prm%selfDiffusionEnergy <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl'
if (prm%fattack <= 0.0_pReal) extmsg = trim(extmsg)//' nu_a'
if (prm%doublekinkwidth <= 0.0_pReal) extmsg = trim(extmsg)//' w'
if (prm%Dsd0 < 0.0_pReal) extmsg = trim(extmsg)//' D_0'
if (prm%atomicVolume <= 0.0_pReal) extmsg = trim(extmsg)//' V_at' ! ToDo: in disloTungsten, the atomic volume is given as a factor
if (prm%eta <= 0.0_pReal) extmsg = trim(extmsg)//' eta'
if (prm%Q_cl <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl'
if (prm%nu_a <= 0.0_pReal) extmsg = trim(extmsg)//' nu_a'
if (prm%w <= 0.0_pReal) extmsg = trim(extmsg)//' w'
if (prm%D_0 < 0.0_pReal) extmsg = trim(extmsg)//' D_0'
if (prm%V_at <= 0.0_pReal) extmsg = trim(extmsg)//' V_at' ! ToDo: in disloTungsten, the atomic volume is given as a factor
if (prm%significantN < 0.0_pReal) extmsg = trim(extmsg)//' rho_num_significant'
if (prm%significantrho < 0.0_pReal) extmsg = trim(extmsg)//' rho_significant'
if (prm%atol_rho < 0.0_pReal) extmsg = trim(extmsg)//' atol_rho'
if (prm%CFLfactor < 0.0_pReal) extmsg = trim(extmsg)//' f_c'
if (prm%rho_min < 0.0_pReal) extmsg = trim(extmsg)//' rho_min'
if (prm%rho_significant < 0.0_pReal) extmsg = trim(extmsg)//' rho_significant'
if (prm%atol_rho < 0.0_pReal) extmsg = trim(extmsg)//' atol_rho'
if (prm%f_c < 0.0_pReal) extmsg = trim(extmsg)//' f_c'
if (prm%p <= 0.0_pReal .or. prm%p > 1.0_pReal) extmsg = trim(extmsg)//' p_sl'
if (prm%q < 1.0_pReal .or. prm%q > 2.0_pReal) extmsg = trim(extmsg)//' q_sl'
if (prm%p <= 0.0_pReal .or. prm%p > 1.0_pReal) extmsg = trim(extmsg)//' p_sl'
if (prm%q < 1.0_pReal .or. prm%q > 2.0_pReal) extmsg = trim(extmsg)//' q_sl'
if (prm%linetensionEffect < 0.0_pReal .or. prm%linetensionEffect > 1.0_pReal) &
if (prm%f_F < 0.0_pReal .or. prm%f_F > 1.0_pReal) &
extmsg = trim(extmsg)//' f_F'
if (prm%edgeJogFactor < 0.0_pReal .or. prm%edgeJogFactor > 1.0_pReal) &
if (prm%f_ed < 0.0_pReal .or. prm%f_ed > 1.0_pReal) &
extmsg = trim(extmsg)//' f_ed'
if (prm%solidSolutionEnergy <= 0.0_pReal) extmsg = trim(extmsg)//' Q_sol'
if (prm%solidSolutionSize <= 0.0_pReal) extmsg = trim(extmsg)//' f_sol'
if (prm%solidSolutionConcentration <= 0.0_pReal) extmsg = trim(extmsg)//' c_sol'
if (prm%Q_sol <= 0.0_pReal) extmsg = trim(extmsg)//' Q_sol'
if (prm%f_sol <= 0.0_pReal) extmsg = trim(extmsg)//' f_sol'
if (prm%c_sol <= 0.0_pReal) extmsg = trim(extmsg)//' c_sol'
if (prm%grainboundaryTransmissivity > 1.0_pReal) extmsg = trim(extmsg)//' chi_GB'
if (prm%surfaceTransmissivity < 0.0_pReal .or. prm%surfaceTransmissivity > 1.0_pReal) &
extmsg = trim(extmsg)//' chi_surface'
if (prm%chi_GB > 1.0_pReal) extmsg = trim(extmsg)//' chi_GB'
if (prm%chi_surface < 0.0_pReal .or. prm%chi_surface > 1.0_pReal) &
extmsg = trim(extmsg)//' chi_surface'
if (prm%fEdgeMultiplication < 0.0_pReal .or. prm%fEdgeMultiplication > 1.0_pReal) &
extmsg = trim(extmsg)//' f_ed_mult'
if (prm%f_ed_mult < 0.0_pReal .or. prm%f_ed_mult > 1.0_pReal) &
extmsg = trim(extmsg)//' f_ed_mult'
endif slipActive
@ -412,7 +412,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
call IO_error(212,ext_msg='IPneighborhood does not exist')
plasticState(p)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention
plasticState(p)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention
st0%rho => plasticState(p)%state0 (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
stt%rho => plasticState(p)%state (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
@ -620,16 +620,16 @@ module subroutine plastic_nonlocal_dependentState(F, Fp, instance, of, ip, el)
! coefficients are corrected for the line tension effect
! (see Kubin,Devincre,Hoc; 2008; Modeling dislocation storage rates and mean free paths in face-centered cubic crystals)
if (any(lattice_structure(material_phaseAt(1,el)) == [LATTICE_bcc_ID,LATTICE_fcc_ID])) then
myInteractionMatrix = prm%interactionSlipSlip &
* spread(( 1.0_pReal - prm%linetensionEffect &
+ prm%linetensionEffect &
* log(0.35_pReal * prm%burgers * sqrt(max(stt%rho_forest(:,of),prm%significantRho))) &
/ log(0.35_pReal * prm%burgers * 1e6_pReal))** 2.0_pReal,2,prm%sum_N_sl)
myInteractionMatrix = prm%h_sl_sl &
* spread(( 1.0_pReal - prm%f_F &
+ prm%f_F &
* log(0.35_pReal * prm%b_sl * sqrt(max(stt%rho_forest(:,of),prm%rho_significant))) &
/ log(0.35_pReal * prm%b_sl * 1e6_pReal))** 2.0_pReal,2,prm%sum_N_sl)
else
myInteractionMatrix = prm%interactionSlipSlip
myInteractionMatrix = prm%h_sl_sl
endif
dst%tau_pass(:,of) = prm%mu * prm%burgers &
dst%tau_pass(:,of) = prm%mu * prm%b_sl &
* sqrt(matmul(myInteractionMatrix,sum(abs(rho),2)))
!*** calculate the dislocation stress of the neighboring excess dislocation densities
@ -731,7 +731,7 @@ module subroutine plastic_nonlocal_dependentState(F, Fp, instance, of, ip, el)
where(rhoTotal > 0.0_pReal) rhoExcessGradient_over_rho = rhoExcessGradient / rhoTotal
! ... gives the local stress correction when multiplied with a factor
dst%tau_back(s,of) = - prm%mu * prm%burgers(s) / (2.0_pReal * PI) &
dst%tau_back(s,of) = - prm%mu * prm%b_sl(s) / (2.0_pReal * PI) &
* ( rhoExcessGradient_over_rho(1) / (1.0_pReal - prm%nu) &
+ rhoExcessGradient_over_rho(2))
enddo
@ -841,7 +841,7 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) &
rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t))
gdotTotal = sum(rhoSgl(:,1:4) * v, 2) * prm%burgers
gdotTotal = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
@ -850,10 +850,10 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
forall (i=1:3,j=1:3,k=1:3,l=1:3) &
dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) &
+ prm%Schmid(i,j,s) * prm%Schmid(k,l,s) &
* sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * prm%burgers(s) &
* sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * prm%b_sl(s) &
+ prm%Schmid(i,j,s) &
* ( prm%nonSchmid_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) &
- prm%nonSchmid_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%burgers(s)
- prm%nonSchmid_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%b_sl(s)
enddo
end associate
@ -929,8 +929,8 @@ module subroutine plastic_nonlocal_deltaState(Mp,instance,of,ip,el)
if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal
enddo
dUpper(:,1) = prm%mu * prm%burgers/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau))
dUpper(:,2) = prm%mu * prm%burgers/(4.0_pReal * PI * abs(tau))
dUpper(:,1) = prm%mu * prm%b_sl/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau))
dUpper(:,2) = prm%mu * prm%b_sl/(4.0_pReal * PI * abs(tau))
where(dNeq0(sqrt(sum(abs(rho(:,edg)),2)))) &
dUpper(:,1) = min(1.0_pReal/sqrt(sum(abs(rho(:,edg)),2)),dUpper(:,1))
@ -1041,7 +1041,7 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, &
my_rhoSgl0 = rho0(:,sgl)
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),of)
gdot = rhoSgl(:,1:4) * v * spread(prm%burgers,2,4)
gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
#ifdef DEBUG
if (debugConstitutive%basic &
@ -1060,8 +1060,8 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, &
enddo
dLower = prm%minDipoleHeight
dUpper(:,1) = prm%mu * prm%burgers/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau))
dUpper(:,2) = prm%mu * prm%burgers/(4.0_pReal * PI * abs(tau))
dUpper(:,1) = prm%mu * prm%b_sl/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau))
dUpper(:,2) = prm%mu * prm%b_sl/(4.0_pReal * PI * abs(tau))
where(dNeq0(sqrt(sum(abs(rho(:,edg)),2)))) &
dUpper(:,1) = min(1.0_pReal/sqrt(sum(abs(rho(:,edg)),2)),dUpper(:,1))
@ -1075,18 +1075,18 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, &
rhoDotMultiplication = 0.0_pReal
isBCC: if (lattice_structure(ph) == LATTICE_bcc_ID) then
forall (s = 1:ns, sum(abs(v(s,1:4))) > 0.0_pReal)
rhoDotMultiplication(s,1:2) = sum(abs(gdot(s,3:4))) / prm%burgers(s) & ! assuming double-cross-slip of screws to be decisive for multiplication
* sqrt(stt%rho_forest(s,of)) / prm%lambda0(s) ! & ! mean free path
rhoDotMultiplication(s,1:2) = sum(abs(gdot(s,3:4))) / prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication
* sqrt(stt%rho_forest(s,of)) / prm%i_sl(s) ! & ! mean free path
! * 2.0_pReal * sum(abs(v(s,3:4))) / sum(abs(v(s,1:4))) ! ratio of screw to overall velocity determines edge generation
rhoDotMultiplication(s,3:4) = sum(abs(gdot(s,3:4))) /prm%burgers(s) & ! assuming double-cross-slip of screws to be decisive for multiplication
* sqrt(stt%rho_forest(s,of)) / prm%lambda0(s) ! & ! mean free path
rhoDotMultiplication(s,3:4) = sum(abs(gdot(s,3:4))) /prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication
* sqrt(stt%rho_forest(s,of)) / prm%i_sl(s) ! & ! mean free path
! * 2.0_pReal * sum(abs(v(s,1:2))) / sum(abs(v(s,1:4))) ! ratio of edge to overall velocity determines screw generation
endforall
else isBCC
rhoDotMultiplication(:,1:4) = spread( &
(sum(abs(gdot(:,1:2)),2) * prm%fEdgeMultiplication + sum(abs(gdot(:,3:4)),2)) &
* sqrt(stt%rho_forest(:,of)) / prm%lambda0 / prm%burgers, 2, 4)
(sum(abs(gdot(:,1:2)),2) * prm%f_ed_mult + sum(abs(gdot(:,3:4)),2)) &
* sqrt(stt%rho_forest(:,of)) / prm%i_sl / prm%b_sl, 2, 4)
endif isBCC
forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,instance),of)
@ -1097,20 +1097,20 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, &
!*** formation by glide
do c = 1,2
rhoDotSingle2DipoleGlide(:,2*c-1) = -2.0_pReal * dUpper(:,c) / prm%burgers &
rhoDotSingle2DipoleGlide(:,2*c-1) = -2.0_pReal * dUpper(:,c) / prm%b_sl &
* ( rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) & ! negative mobile --> positive mobile
+ rhoSgl(:,2*c) * abs(gdot(:,2*c-1)) & ! positive mobile --> negative mobile
+ abs(rhoSgl(:,2*c+4)) * abs(gdot(:,2*c-1))) ! positive mobile --> negative immobile
rhoDotSingle2DipoleGlide(:,2*c) = -2.0_pReal * dUpper(:,c) / prm%burgers &
rhoDotSingle2DipoleGlide(:,2*c) = -2.0_pReal * dUpper(:,c) / prm%b_sl &
* ( rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) & ! negative mobile --> positive mobile
+ rhoSgl(:,2*c) * abs(gdot(:,2*c-1)) & ! positive mobile --> negative mobile
+ abs(rhoSgl(:,2*c+3)) * abs(gdot(:,2*c))) ! negative mobile --> positive immobile
rhoDotSingle2DipoleGlide(:,2*c+3) = -2.0_pReal * dUpper(:,c) / prm%burgers &
rhoDotSingle2DipoleGlide(:,2*c+3) = -2.0_pReal * dUpper(:,c) / prm%b_sl &
* rhoSgl(:,2*c+3) * abs(gdot(:,2*c)) ! negative mobile --> positive immobile
rhoDotSingle2DipoleGlide(:,2*c+4) = -2.0_pReal * dUpper(:,c) / prm%burgers &
rhoDotSingle2DipoleGlide(:,2*c+4) = -2.0_pReal * dUpper(:,c) / prm%b_sl &
* rhoSgl(:,2*c+4) * abs(gdot(:,2*c-1)) ! positive mobile --> negative immobile
rhoDotSingle2DipoleGlide(:,c+8) = abs(rhoDotSingle2DipoleGlide(:,2*c+3)) &
@ -1123,27 +1123,27 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, &
!*** athermal annihilation
rhoDotAthermalAnnihilation = 0.0_pReal
forall (c=1:2) &
rhoDotAthermalAnnihilation(:,c+8) = -2.0_pReal * dLower(:,c) / prm%burgers &
* ( 2.0_pReal * (rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) + rhoSgl(:,2*c) * abs(gdot(:,2*c-1))) & ! was single hitting single
rhoDotAthermalAnnihilation(:,c+8) = -2.0_pReal * dLower(:,c) / prm%b_sl &
* ( 2.0_pReal * (rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) + rhoSgl(:,2*c) * abs(gdot(:,2*c-1))) & ! was single hitting single
+ 2.0_pReal * (abs(rhoSgl(:,2*c+3)) * abs(gdot(:,2*c)) + abs(rhoSgl(:,2*c+4)) * abs(gdot(:,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single
+ rhoDip(:,c) * (abs(gdot(:,2*c-1)) + abs(gdot(:,2*c)))) ! single knocks dipole constituent
+ rhoDip(:,c) * (abs(gdot(:,2*c-1)) + abs(gdot(:,2*c)))) ! single knocks dipole constituent
! annihilated screw dipoles leave edge jogs behind on the colinear system
if (lattice_structure(ph) == LATTICE_fcc_ID) &
forall (s = 1:ns, prm%colinearSystem(s) > 0) &
rhoDotAthermalAnnihilation(prm%colinearSystem(s),1:2) = - rhoDotAthermalAnnihilation(s,10) &
* 0.25_pReal * sqrt(stt%rho_forest(s,of)) * (dUpper(s,2) + dLower(s,2)) * prm%edgeJogFactor
* 0.25_pReal * sqrt(stt%rho_forest(s,of)) * (dUpper(s,2) + dLower(s,2)) * prm%f_ed
!*** thermally activated annihilation of edge dipoles by climb
rhoDotThermalAnnihilation = 0.0_pReal
selfDiffusion = prm%Dsd0 * exp(-prm%selfDiffusionEnergy / (kB * Temperature))
vClimb = prm%atomicVolume * selfDiffusion * prm%mu &
selfDiffusion = prm%D_0 * exp(-prm%Q_cl / (kB * Temperature))
vClimb = prm%V_at * selfDiffusion * prm%mu &
/ ( kB * Temperature * PI * (1.0_pReal-prm%nu) * (dUpper(:,1) + dLower(:,1)))
forall (s = 1:ns, dUpper(s,1) > dLower(s,1)) &
rhoDotThermalAnnihilation(s,9) = max(- 4.0_pReal * rhoDip(s,1) * vClimb(s) / (dUpper(s,1) - dLower(s,1)), &
- rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) &
- rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have
- rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have
rhoDot = rhoDotFlux(F,Fp,timestep, instance,of,ip,el) &
+ rhoDotMultiplication &
@ -1252,7 +1252,7 @@ function rhoDotFlux(F,Fp,timestep, instance,of,ip,el)
my_rhoSgl0 = rho0(:,sgl)
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),of) !ToDo: MD: I think we should use state0 here
gdot = rhoSgl(:,1:4) * v * spread(prm%burgers,2,4)
gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,instance),of)
@ -1264,14 +1264,14 @@ function rhoDotFlux(F,Fp,timestep, instance,of,ip,el)
!*** check CFL (Courant-Friedrichs-Lewy) condition for flux
if (any( abs(gdot) > 0.0_pReal & ! any active slip system ...
.and. prm%CFLfactor * abs(v0) * timestep &
.and. prm%f_c * abs(v0) * timestep &
> IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here)
#ifdef DEBUG
if (debugConstitutive%extensive) then
print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip
print'(a,e10.3,a,e10.3)', '<< CONST >> velocity is at ', &
maxval(abs(v0), abs(gdot) > 0.0_pReal &
.and. prm%CFLfactor * abs(v0) * timestep &
.and. prm%f_c * abs(v0) * timestep &
> IPvolume(ip,el) / maxval(IParea(:,ip,el))), &
' at a timestep of ',timestep
print*, '<< CONST >> enforcing cutback !!!'
@ -1334,8 +1334,8 @@ function rhoDotFlux(F,Fp,timestep, instance,of,ip,el)
neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,neighbor_instance),no),0.0_pReal)
endforall
where (neighbor_rhoSgl0 * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%significantN &
.or. neighbor_rhoSgl0 < prm%significantRho) &
where (neighbor_rhoSgl0 * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%rho_min &
.or. neighbor_rhoSgl0 < prm%rho_significant) &
neighbor_rhoSgl0 = 0.0_pReal
normal_neighbor2me_defConf = math_det33(Favg) * matmul(math_inv33(transpose(Favg)), &
IPareaNormal(1:3,neighbor_n,neighbor_ip,neighbor_el)) ! normal of the interface in (average) deformed configuration (pointing neighbor => me)
@ -1460,7 +1460,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e)
if (neighbor_e <= 0 .or. neighbor_i <= 0) then
!* FREE SURFACE
!* Set surface transmissivity to the value specified in the material.config
forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%surfaceTransmissivity)
forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_surface)
elseif (neighbor_phase /= ph) then
!* PHASE BOUNDARY
!* If we encounter a different nonlocal phase at the neighbor,
@ -1469,13 +1469,13 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e)
!* we do not consider this to be a phase boundary, so completely compatible.
if (.not. phase_localPlasticity(neighbor_phase) .and. .not. phase_localPlasticity(ph)) &
forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = 0.0_pReal
elseif (prm%grainboundaryTransmissivity >= 0.0_pReal) then
elseif (prm%chi_GB >= 0.0_pReal) then
!* GRAIN BOUNDARY !
!* fixed transmissivity for adjacent ips with different texture (only if explicitly given in material.config)
if (any(dNeq(material_orientation0(1,i,e)%asQuaternion(), &
material_orientation0(1,neighbor_i,neighbor_e)%asQuaternion())) .and. &
(.not. phase_localPlasticity(neighbor_phase))) &
forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%grainboundaryTransmissivity)
forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB)
else
!* GRAIN BOUNDARY ?
!* Compatibility defined by relative orientation of slip systems:
@ -1631,7 +1631,7 @@ subroutine stateInit(ini,phase,NipcMyPhase,instance)
associate(stt => state(instance))
if (ini%rhoSglRandom > 0.0_pReal) then ! randomly distribute dislocation segments on random slip system and of random type in the volume
if (ini%random_rho_u > 0.0_pReal) then ! randomly distribute dislocation segments on random slip system and of random type in the volume
do e = 1,discretization_nElem
do i = 1,discretization_nIP
if (material_phaseAt(1,e) == phase) volume(material_phasememberAt(1,i,e)) = IPvolume(i,e)
@ -1639,11 +1639,11 @@ subroutine stateInit(ini,phase,NipcMyPhase,instance)
enddo
totalVolume = sum(volume)
minimumIPVolume = minval(volume)
densityBinning = ini%rhoSglRandomBinning / minimumIpVolume ** (2.0_pReal / 3.0_pReal)
densityBinning = ini%random_rho_u_binning / minimumIpVolume ** (2.0_pReal / 3.0_pReal)
! fill random material points with dislocation segments until the desired overall density is reached
meanDensity = 0.0_pReal
do while(meanDensity < ini%rhoSglRandom)
do while(meanDensity < ini%random_rho_u)
call random_number(rnd)
phasemember = nint(rnd(1)*real(NipcMyPhase,pReal) + 0.5_pReal)
s = nint(rnd(2)*real(sum(ini%N_sl),pReal)*4.0_pReal + 0.5_pReal)
@ -1656,15 +1656,15 @@ subroutine stateInit(ini,phase,NipcMyPhase,instance)
from = 1 + sum(ini%N_sl(1:f-1))
upto = sum(ini%N_sl(1:f))
do s = from,upto
noise = [math_sampleGaussVar(0.0_pReal, ini%rhoSglScatter), &
math_sampleGaussVar(0.0_pReal, ini%rhoSglScatter)]
stt%rho_sgl_mob_edg_pos(s,e) = ini%rhoSglEdgePos0(f) + noise(1)
stt%rho_sgl_mob_edg_neg(s,e) = ini%rhoSglEdgeNeg0(f) + noise(1)
stt%rho_sgl_mob_scr_pos(s,e) = ini%rhoSglScrewPos0(f) + noise(2)
stt%rho_sgl_mob_scr_neg(s,e) = ini%rhoSglScrewNeg0(f) + noise(2)
noise = [math_sampleGaussVar(0.0_pReal, ini%sigma_rho_u), &
math_sampleGaussVar(0.0_pReal, ini%sigma_rho_u)]
stt%rho_sgl_mob_edg_pos(s,e) = ini%rho_u_ed_pos_0(f) + noise(1)
stt%rho_sgl_mob_edg_neg(s,e) = ini%rho_u_ed_neg_0(f) + noise(1)
stt%rho_sgl_mob_scr_pos(s,e) = ini%rho_u_sc_pos_0(f) + noise(2)
stt%rho_sgl_mob_scr_neg(s,e) = ini%rho_u_sc_neg_0(f) + noise(2)
enddo
stt%rho_dip_edg(from:upto,e) = ini%rhoDipEdge0(f)
stt%rho_dip_scr(from:upto,e) = ini%rhoDipScrew0(f)
stt%rho_dip_edg(from:upto,e) = ini%rho_d_ed_0(f)
stt%rho_dip_scr(from:upto,e) = ini%rho_d_sc_0(f)
enddo
enddo
endif
@ -1732,14 +1732,14 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Tem
!* Effective stress includes non Schmid constributions
!* The derivative only gives absolute values; the correct sign is taken care of in the formula for the derivative of the velocity
tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s)) ! ensure that the effective stress is positive
meanfreepath_P = prm%burgers(s)
jumpWidth_P = prm%burgers(s)
activationLength_P = prm%doublekinkwidth *prm%burgers(s)
activationVolume_P = activationLength_P * jumpWidth_P * prm%burgers(s)
meanfreepath_P = prm%b_sl(s)
jumpWidth_P = prm%b_sl(s)
activationLength_P = prm%w *prm%b_sl(s)
activationVolume_P = activationLength_P * jumpWidth_P * prm%b_sl(s)
criticalStress_P = prm%peierlsStress(s,c)
activationEnergy_P = criticalStress_P * activationVolume_P
tauRel_P = min(1.0_pReal, tauEff / criticalStress_P) ! ensure that the activation probability cannot become greater than one
tPeierls = 1.0_pReal / prm%fattack &
tPeierls = 1.0_pReal / prm%nu_a &
* exp(activationEnergy_P / (kB * Temperature) &
* (1.0_pReal - tauRel_P**prm%p)**prm%q)
if (tauEff < criticalStress_P) then
@ -1752,14 +1752,14 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Tem
!* Contribution from solid solution strengthening
!* The derivative only gives absolute values; the correct sign is taken care of in the formula for the derivative of the velocity
tauEff = abs(tau(s)) - tauThreshold(s)
meanfreepath_S = prm%burgers(s) / sqrt(prm%solidSolutionConcentration)
jumpWidth_S = prm%solidSolutionSize * prm%burgers(s)
activationLength_S = prm%burgers(s) / sqrt(prm%solidSolutionConcentration)
activationVolume_S = activationLength_S * jumpWidth_S * prm%burgers(s)
activationEnergy_S = prm%solidSolutionEnergy
meanfreepath_S = prm%b_sl(s) / sqrt(prm%c_sol)
jumpWidth_S = prm%f_sol * prm%b_sl(s)
activationLength_S = prm%b_sl(s) / sqrt(prm%c_sol)
activationVolume_S = activationLength_S * jumpWidth_S * prm%b_sl(s)
activationEnergy_S = prm%Q_sol
criticalStress_S = activationEnergy_S / activationVolume_S
tauRel_S = min(1.0_pReal, tauEff / criticalStress_S) ! ensure that the activation probability cannot become greater than one
tSolidSolution = 1.0_pReal / prm%fattack &
tSolidSolution = 1.0_pReal / prm%nu_a &
* exp(activationEnergy_S / (kB * Temperature)* (1.0_pReal - tauRel_S**prm%p)**prm%q)
if (tauEff < criticalStress_S) then
dtSolidSolution_dtau = tSolidSolution * prm%p * prm%q * activationVolume_S / (kB * Temperature) &
@ -1770,7 +1770,7 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Tem
!* viscous glide velocity
tauEff = abs(tau(s)) - tauThreshold(s)
mobility = prm%burgers(s) / prm%viscosity
mobility = prm%b_sl(s) / prm%eta
vViscous = mobility * tauEff
!* Mean velocity results from waiting time at peierls barriers and solid solution obstacles with respective meanfreepath of
@ -1805,7 +1805,7 @@ pure function getRho(instance,of,ip,el)
getRho(:,mob) = max(getRho(:,mob),0.0_pReal)
getRho(:,dip) = max(getRho(:,dip),0.0_pReal)
where(abs(getRho) < max(prm%significantN/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%significantRho)) &
where(abs(getRho) < max(prm%rho_min/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%rho_significant)) &
getRho = 0.0_pReal
end associate
@ -1830,7 +1830,7 @@ pure function getRho0(instance,of,ip,el)
getRho0(:,mob) = max(getRho0(:,mob),0.0_pReal)
getRho0(:,dip) = max(getRho0(:,dip),0.0_pReal)
where(abs(getRho0) < max(prm%significantN/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%significantRho)) &
where(abs(getRho0) < max(prm%rho_min/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%rho_significant)) &
getRho0 = 0.0_pReal
end associate

View File

@ -8,28 +8,28 @@ submodule(constitutive:constitutive_plastic) plastic_phenopowerlaw
type :: tParameters
real(pReal) :: &
gdot0_slip = 1.0_pReal, & !< reference shear strain rate for slip
gdot0_twin = 1.0_pReal, & !< reference shear strain rate for twin
n_slip = 1.0_pReal, & !< stress exponent for slip
n_twin = 1.0_pReal, & !< stress exponent for twin
spr = 1.0_pReal, & !< push-up factor for slip saturation due to twinning
c_1 = 1.0_pReal, &
c_2 = 1.0_pReal, &
c_3 = 1.0_pReal, &
c_4 = 1.0_pReal, &
h0_SlipSlip = 1.0_pReal, & !< reference hardening slip - slip
h0_TwinSlip = 1.0_pReal, & !< reference hardening twin - slip
h0_TwinTwin = 1.0_pReal, & !< reference hardening twin - twin
a_slip = 1.0_pReal
dot_gamma_0_sl = 1.0_pReal, & !< reference shear strain rate for slip
dot_gamma_0_tw = 1.0_pReal, & !< reference shear strain rate for twin
n_sl = 1.0_pReal, & !< stress exponent for slip
n_tw = 1.0_pReal, & !< stress exponent for twin
f_sl_sat_tw = 1.0_pReal, & !< push-up factor for slip saturation due to twinning
c_1 = 1.0_pReal, &
c_2 = 1.0_pReal, &
c_3 = 1.0_pReal, &
c_4 = 1.0_pReal, &
h_0_sl_sl = 1.0_pReal, & !< reference hardening slip - slip
h_0_tw_sl = 1.0_pReal, & !< reference hardening twin - slip
h_0_tw_tw = 1.0_pReal, & !< reference hardening twin - twin
a_sl = 1.0_pReal
real(pReal), allocatable, dimension(:) :: &
xi_slip_sat, & !< maximum critical shear stress for slip
H_int, & !< per family hardening activity (optional)
gamma_twin_char !< characteristic shear for twins
xi_inf_sl, & !< maximum critical shear stress for slip
h_int, & !< per family hardening activity (optional)
gamma_tw_char !< characteristic shear for twins
real(pReal), allocatable, dimension(:,:) :: &
interaction_SlipSlip, & !< slip resistance from slip activity
interaction_SlipTwin, & !< slip resistance from twin activity
interaction_TwinSlip, & !< twin resistance from slip activity
interaction_TwinTwin !< twin resistance from twin activity
h_sl_sl, & !< slip resistance from slip activity
h_sl_tw, & !< slip resistance from twin activity
h_tw_sl, & !< twin resistance from slip activity
h_tw_tw !< twin resistance from twin activity
real(pReal), allocatable, dimension(:,:,:) :: &
P_sl, &
P_tw, &
@ -78,8 +78,8 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
integer, dimension(:), allocatable :: &
N_sl, N_tw
real(pReal), dimension(:), allocatable :: &
xi_slip_0, & !< initial critical shear stress for slip
xi_twin_0, & !< initial critical shear stress for twin
xi_0_sl, & !< initial critical shear stress for slip
xi_0_tw, & !< initial critical shear stress for twin
a !< non-Schmid coefficients
character(len=pStringLen) :: &
extmsg = ''
@ -120,7 +120,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
if(phase%get_asString('lattice') == 'bcc') then
a = pl%get_asFloats('nonSchmid_coefficients',defaultVal=emptyRealArray)
a = pl%get_asFloats('a_nonSchmid',defaultVal=emptyRealArray)
if(size(a) > 0) prm%nonSchmidActive = .true.
prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,a,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,a,-1)
@ -128,36 +128,36 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
prm%nonSchmid_pos = prm%P_sl
prm%nonSchmid_neg = prm%P_sl
endif
prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(N_sl, &
pl%get_asFloats('h_sl_sl'), &
phase%get_asString('lattice'))
prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl, &
pl%get_asFloats('h_sl_sl'), &
phase%get_asString('lattice'))
xi_slip_0 = pl%get_asFloats('xi_0_sl', requiredSize=size(N_sl))
prm%xi_slip_sat = pl%get_asFloats('xi_inf_sl', requiredSize=size(N_sl))
prm%H_int = pl%get_asFloats('h_int', requiredSize=size(N_sl), &
defaultVal=[(0.0_pReal,i=1,size(N_sl))])
xi_0_sl = pl%get_asFloats('xi_0_sl', requiredSize=size(N_sl))
prm%xi_inf_sl = pl%get_asFloats('xi_inf_sl', requiredSize=size(N_sl))
prm%h_int = pl%get_asFloats('h_int', requiredSize=size(N_sl), &
defaultVal=[(0.0_pReal,i=1,size(N_sl))])
prm%gdot0_slip = pl%get_asFloat('dot_gamma_0_sl')
prm%n_slip = pl%get_asFloat('n_sl')
prm%a_slip = pl%get_asFloat('a_sl')
prm%h0_SlipSlip = pl%get_asFloat('h_0_sl_sl')
prm%dot_gamma_0_sl = pl%get_asFloat('dot_gamma_0_sl')
prm%n_sl = pl%get_asFloat('n_sl')
prm%a_sl = pl%get_asFloat('a_sl')
prm%h_0_sl_sl = pl%get_asFloat('h_0_sl_sl')
! expand: family => system
xi_slip_0 = math_expand(xi_slip_0, N_sl)
prm%xi_slip_sat = math_expand(prm%xi_slip_sat,N_sl)
prm%H_int = math_expand(prm%H_int, N_sl)
xi_0_sl = math_expand(xi_0_sl, N_sl)
prm%xi_inf_sl = math_expand(prm%xi_inf_sl,N_sl)
prm%h_int = math_expand(prm%h_int, N_sl)
! sanity checks
if ( prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//' dot_gamma_0_sl'
if ( prm%a_slip <= 0.0_pReal) extmsg = trim(extmsg)//' a_sl'
if ( prm%n_slip <= 0.0_pReal) extmsg = trim(extmsg)//' n_sl'
if (any(xi_slip_0 <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_0_sl'
if (any(prm%xi_slip_sat <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_inf_sl'
if ( prm%dot_gamma_0_sl <= 0.0_pReal) extmsg = trim(extmsg)//' dot_gamma_0_sl'
if ( prm%a_sl <= 0.0_pReal) extmsg = trim(extmsg)//' a_sl'
if ( prm%n_sl <= 0.0_pReal) extmsg = trim(extmsg)//' n_sl'
if (any(xi_0_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_0_sl'
if (any(prm%xi_inf_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_inf_sl'
else slipActive
xi_slip_0 = emptyRealArray
allocate(prm%xi_slip_sat,prm%H_int,source=emptyRealArray)
allocate(prm%interaction_SlipSlip(0,0))
xi_0_sl = emptyRealArray
allocate(prm%xi_inf_sl,prm%h_int,source=emptyRealArray)
allocate(prm%h_sl_sl(0,0))
endif slipActive
!--------------------------------------------------------------------------------------------------
@ -165,52 +165,52 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
N_tw = pl%get_asInts('N_tw', defaultVal=emptyIntArray)
prm%sum_N_tw = sum(abs(N_tw))
twinActive: if (prm%sum_N_tw > 0) then
prm%P_tw = lattice_SchmidMatrix_twin(N_tw,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
prm%interaction_TwinTwin = lattice_interaction_TwinByTwin(N_tw,&
pl%get_asFloats('h_tw_tw'), &
phase%get_asString('lattice'))
prm%gamma_twin_char = lattice_characteristicShear_twin(N_tw,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
prm%P_tw = lattice_SchmidMatrix_twin(N_tw,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
prm%h_tw_tw = lattice_interaction_TwinByTwin(N_tw,&
pl%get_asFloats('h_tw_tw'), &
phase%get_asString('lattice'))
prm%gamma_tw_char = lattice_characteristicShear_twin(N_tw,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
xi_twin_0 = pl%get_asFloats('xi_0_tw',requiredSize=size(N_tw))
xi_0_tw = pl%get_asFloats('xi_0_tw',requiredSize=size(N_tw))
prm%c_1 = pl%get_asFloat('c_1',defaultVal=0.0_pReal)
prm%c_2 = pl%get_asFloat('c_2',defaultVal=1.0_pReal)
prm%c_3 = pl%get_asFloat('c_3',defaultVal=0.0_pReal)
prm%c_4 = pl%get_asFloat('c_4',defaultVal=0.0_pReal)
prm%gdot0_twin = pl%get_asFloat('dot_gamma_0_tw')
prm%n_twin = pl%get_asFloat('n_tw')
prm%spr = pl%get_asFloat('f_sl_sat_tw')
prm%h0_TwinTwin = pl%get_asFloat('h_0_tw_tw')
prm%c_1 = pl%get_asFloat('c_1',defaultVal=0.0_pReal)
prm%c_2 = pl%get_asFloat('c_2',defaultVal=1.0_pReal)
prm%c_3 = pl%get_asFloat('c_3',defaultVal=0.0_pReal)
prm%c_4 = pl%get_asFloat('c_4',defaultVal=0.0_pReal)
prm%dot_gamma_0_tw = pl%get_asFloat('dot_gamma_0_tw')
prm%n_tw = pl%get_asFloat('n_tw')
prm%f_sl_sat_tw = pl%get_asFloat('f_sl_sat_tw')
prm%h_0_tw_tw = pl%get_asFloat('h_0_tw_tw')
! expand: family => system
xi_twin_0 = math_expand(xi_twin_0,N_tw)
xi_0_tw = math_expand(xi_0_tw,N_tw)
! sanity checks
if (prm%gdot0_twin <= 0.0_pReal) extmsg = trim(extmsg)//' dot_gamma_0_tw'
if (prm%n_twin <= 0.0_pReal) extmsg = trim(extmsg)//' n_tw'
if (prm%dot_gamma_0_tw <= 0.0_pReal) extmsg = trim(extmsg)//' dot_gamma_0_tw'
if (prm%n_tw <= 0.0_pReal) extmsg = trim(extmsg)//' n_tw'
else twinActive
xi_twin_0 = emptyRealArray
allocate(prm%gamma_twin_char,source=emptyRealArray)
allocate(prm%interaction_TwinTwin(0,0))
xi_0_tw = emptyRealArray
allocate(prm%gamma_tw_char,source=emptyRealArray)
allocate(prm%h_tw_tw(0,0))
endif twinActive
!--------------------------------------------------------------------------------------------------
! slip-twin related parameters
slipAndTwinActive: if (prm%sum_N_sl > 0 .and. prm%sum_N_tw > 0) then
prm%h0_TwinSlip = pl%get_asFloat('h_0_tw_sl')
prm%interaction_SlipTwin = lattice_interaction_SlipByTwin(N_sl,N_tw,&
pl%get_asFloats('h_sl_tw'), &
phase%get_asString('lattice'))
prm%interaction_TwinSlip = lattice_interaction_TwinBySlip(N_tw,N_sl,&
pl%get_asFloats('h_tw_sl'), &
phase%get_asString('lattice'))
prm%h_0_tw_sl = pl%get_asFloat('h_0_tw_sl')
prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,N_tw,&
pl%get_asFloats('h_sl_tw'), &
phase%get_asString('lattice'))
prm%h_tw_sl = lattice_interaction_TwinBySlip(N_tw,N_sl,&
pl%get_asFloats('h_tw_sl'), &
phase%get_asString('lattice'))
else slipAndTwinActive
allocate(prm%interaction_SlipTwin(prm%sum_N_sl,prm%sum_N_tw)) ! at least one dimension is 0
allocate(prm%interaction_TwinSlip(prm%sum_N_tw,prm%sum_N_sl)) ! at least one dimension is 0
prm%h0_TwinSlip = 0.0_pReal
allocate(prm%h_sl_tw(prm%sum_N_sl,prm%sum_N_tw)) ! at least one dimension is 0
allocate(prm%h_tw_sl(prm%sum_N_tw,prm%sum_N_sl)) ! at least one dimension is 0
prm%h_0_tw_sl = 0.0_pReal
endif slipAndTwinActive
!--------------------------------------------------------------------------------------------------
@ -237,7 +237,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
startIndex = 1
endIndex = prm%sum_N_sl
stt%xi_slip => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_slip = spread(xi_slip_0, 2, NipcMyPhase)
stt%xi_slip = spread(xi_0_sl, 2, NipcMyPhase)
dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
@ -245,7 +245,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw
stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_twin = spread(xi_twin_0, 2, NipcMyPhase)
stt%xi_twin = spread(xi_0_tw, 2, NipcMyPhase)
dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
@ -354,20 +354,20 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
associate(prm => param(instance), stt => state(instance), dot => dotState(instance))
sumGamma = sum(stt%gamma_slip(:,of))
sumF = sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)
sumF = sum(stt%gamma_twin(:,of)/prm%gamma_tw_char)
!--------------------------------------------------------------------------------------------------
! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices
c_SlipSlip = prm%h0_slipslip * (1.0_pReal + prm%c_1*sumF** prm%c_2)
c_TwinSlip = prm%h0_TwinSlip * sumGamma**prm%c_3
c_TwinTwin = prm%h0_TwinTwin * sumF**prm%c_4
c_SlipSlip = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2)
c_TwinSlip = prm%h_0_tw_sl * sumGamma**prm%c_3
c_TwinTwin = prm%h_0_tw_tw * sumF**prm%c_4
!--------------------------------------------------------------------------------------------------
! calculate left and right vectors
left_SlipSlip = 1.0_pReal + prm%H_int
xi_slip_sat_offset = prm%spr*sqrt(sumF)
right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) **prm%a_slip &
* sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset))
left_SlipSlip = 1.0_pReal + prm%h_int
xi_slip_sat_offset = prm%f_sl_sat_tw*sqrt(sumF)
right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,of) / (prm%xi_inf_sl+xi_slip_sat_offset)) **prm%a_sl &
* sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,of) / (prm%xi_inf_sl+xi_slip_sat_offset))
!--------------------------------------------------------------------------------------------------
! shear rates
@ -378,11 +378,11 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
!--------------------------------------------------------------------------------------------------
! hardening
dot%xi_slip(:,of) = c_SlipSlip * left_SlipSlip * &
matmul(prm%interaction_SlipSlip,dot%gamma_slip(:,of)*right_SlipSlip) &
+ matmul(prm%interaction_SlipTwin,dot%gamma_twin(:,of))
matmul(prm%h_sl_sl,dot%gamma_slip(:,of)*right_SlipSlip) &
+ matmul(prm%h_sl_tw,dot%gamma_twin(:,of))
dot%xi_twin(:,of) = c_TwinSlip * matmul(prm%interaction_TwinSlip,dot%gamma_slip(:,of)) &
+ c_TwinTwin * matmul(prm%interaction_TwinTwin,dot%gamma_twin(:,of))
dot%xi_twin(:,of) = c_TwinSlip * matmul(prm%h_tw_sl,dot%gamma_slip(:,of)) &
+ c_TwinTwin * matmul(prm%h_tw_tw,dot%gamma_twin(:,of))
end associate
end subroutine plastic_phenopowerlaw_dotState
@ -460,29 +460,29 @@ pure subroutine kinetics_slip(Mp,instance,of, &
enddo
where(dNeq0(tau_slip_pos))
gdot_slip_pos = prm%gdot0_slip * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active
* sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_slip, tau_slip_pos)
gdot_slip_pos = prm%dot_gamma_0_sl * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active
* sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_sl, tau_slip_pos)
else where
gdot_slip_pos = 0.0_pReal
end where
where(dNeq0(tau_slip_neg))
gdot_slip_neg = prm%gdot0_slip * 0.5_pReal & ! only used if non-Schmid active, always 1/2
* sign(abs(tau_slip_neg/stt%xi_slip(:,of))**prm%n_slip, tau_slip_neg)
gdot_slip_neg = prm%dot_gamma_0_sl * 0.5_pReal & ! only used if non-Schmid active, always 1/2
* sign(abs(tau_slip_neg/stt%xi_slip(:,of))**prm%n_sl, tau_slip_neg)
else where
gdot_slip_neg = 0.0_pReal
end where
if (present(dgdot_dtau_slip_pos)) then
where(dNeq0(gdot_slip_pos))
dgdot_dtau_slip_pos = gdot_slip_pos*prm%n_slip/tau_slip_pos
dgdot_dtau_slip_pos = gdot_slip_pos*prm%n_sl/tau_slip_pos
else where
dgdot_dtau_slip_pos = 0.0_pReal
end where
endif
if (present(dgdot_dtau_slip_neg)) then
where(dNeq0(gdot_slip_neg))
dgdot_dtau_slip_neg = gdot_slip_neg*prm%n_slip/tau_slip_neg
dgdot_dtau_slip_neg = gdot_slip_neg*prm%n_sl/tau_slip_neg
else where
dgdot_dtau_slip_neg = 0.0_pReal
end where
@ -524,15 +524,15 @@ pure subroutine kinetics_twin(Mp,instance,of,&
enddo
where(tau_twin > 0.0_pReal)
gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)) & ! only twin in untwinned volume fraction
* prm%gdot0_twin*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_twin
gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,of)/prm%gamma_tw_char)) & ! only twin in untwinned volume fraction
* prm%dot_gamma_0_tw*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_tw
else where
gdot_twin = 0.0_pReal
end where
if (present(dgdot_dtau_twin)) then
where(dNeq0(gdot_twin))
dgdot_dtau_twin = gdot_twin*prm%n_twin/tau_twin
dgdot_dtau_twin = gdot_twin*prm%n_tw/tau_twin
else where
dgdot_dtau_twin = 0.0_pReal
end where

View File

@ -134,7 +134,7 @@ function damage_nonlocal_getDiffusion(ip,el)
damage_nonlocal_getDiffusion = 0.0_pReal
do grain = 1, homogenization_Ngrains(homog)
damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + &
crystallite_push33ToRef(grain,ip,el,lattice_DamageDiffusion(1:3,1:3,material_phaseAt(grain,el)))
crystallite_push33ToRef(grain,ip,el,lattice_D(1:3,1:3,material_phaseAt(grain,el)))
enddo
damage_nonlocal_getDiffusion = &
@ -157,7 +157,7 @@ real(pReal) function damage_nonlocal_getMobility(ip,el)
damage_nonlocal_getMobility = 0.0_pReal
do ipc = 1, homogenization_Ngrains(material_homogenizationAt(el))
damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_DamageMobility(material_phaseAt(ipc,el))
damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_M(material_phaseAt(ipc,el))
enddo
damage_nonlocal_getMobility = damage_nonlocal_getMobility/&

View File

@ -56,7 +56,7 @@ subroutine discretization_grid_init(restart)
myGrid !< domain grid of this process
integer, dimension(:), allocatable :: &
microstructureAt
materialAt
integer :: &
j, &
@ -68,9 +68,9 @@ subroutine discretization_grid_init(restart)
print'(/,a)', ' <<<+- discretization_grid init -+>>>'; flush(IO_STDOUT)
if(index(interface_geomFile,'.vtr') /= 0) then
call readVTR(grid,geomSize,origin,microstructureAt)
call readVTR(grid,geomSize,origin,materialAt)
else
call readGeom(grid,geomSize,origin,microstructureAt)
call readGeom(grid,geomSize,origin,materialAt)
endif
print'(/,a,3(i12 ))', ' grid a b c: ', grid
@ -102,15 +102,14 @@ subroutine discretization_grid_init(restart)
!--------------------------------------------------------------------------------------------------
! general discretization
microstructureAt = microstructureAt(product(grid(1:2))*grid3Offset+1: &
product(grid(1:2))*(grid3Offset+grid3)) ! reallocate/shrink in case of MPI
materialAt = materialAt(product(grid(1:2))*grid3Offset+1:product(grid(1:2))*(grid3Offset+grid3)) ! reallocate/shrink in case of MPI
call discretization_init(microstructureAt, &
call discretization_init(materialAt, &
IPcoordinates0(myGrid,mySize,grid3Offset), &
Nodes0(myGrid,mySize,grid3Offset),&
merge((grid(1)+1) * (grid(2)+1) * (grid3+1),& ! write bottom layer
(grid(1)+1) * (grid(2)+1) * grid3,& ! do not write bottom layer (is top of rank-1)
worldrank<1))
merge((grid(1)+1) * (grid(2)+1) * (grid3+1),& ! write top layer...
(grid(1)+1) * (grid(2)+1) * grid3,& ! ...unless not last process
worldrank+1==worldsize))
FEsolving_execElem = [1,product(myGrid)] ! parallel loop bounds set to comprise all elements
FEsolving_execIP = [1,1] ! parallel loop bounds set to comprise the only IP
@ -147,7 +146,7 @@ end subroutine discretization_grid_init
!> @details important variables have an implicit "save" attribute. Therefore, this function is
! supposed to be called only once!
!--------------------------------------------------------------------------------------------------
subroutine readGeom(grid,geomSize,origin,microstructure)
subroutine readGeom(grid,geomSize,origin,material)
integer, dimension(3), intent(out) :: &
grid ! grid (across all processes!)
@ -155,7 +154,7 @@ subroutine readGeom(grid,geomSize,origin,microstructure)
geomSize, & ! size (across all processes!)
origin ! origin (across all processes!)
integer, dimension(:), intent(out), allocatable :: &
microstructure
material
character(len=:), allocatable :: rawData
character(len=65536) :: line
@ -167,7 +166,7 @@ subroutine readGeom(grid,geomSize,origin,microstructure)
startPos, endPos, &
myStat, &
l, & !< line counter
c, & !< counter for # microstructures in line
c, & !< counter for # materials in line
o, & !< order of "to" packing
e, & !< "element", i.e. spectral collocation point
i, j
@ -266,7 +265,7 @@ subroutine readGeom(grid,geomSize,origin,microstructure)
if(any(geomSize < 0.0_pReal)) &
call IO_error(error_ID = 842, ext_msg='size (readGeom)')
allocate(microstructure(product(grid)), source = -1) ! too large in case of MPI (shrink later, not very elegant)
allocate(material(product(grid)), source = -1) ! too large in case of MPI (shrink later, not very elegant)
!--------------------------------------------------------------------------------------------------
! read and interpret content
@ -281,18 +280,18 @@ subroutine readGeom(grid,geomSize,origin,microstructure)
noCompression: if (chunkPos(1) /= 3) then
c = chunkPos(1)
microstructure(e:e+c-1) = [(IO_intValue(line,chunkPos,i+1), i=0, c-1)]
material(e:e+c-1) = [(IO_intValue(line,chunkPos,i+1), i=0, c-1)]
else noCompression
compression: if (IO_lc(IO_stringValue(line,chunkPos,2)) == 'of') then
c = IO_intValue(line,chunkPos,1)
microstructure(e:e+c-1) = [(IO_intValue(line,chunkPos,3),i = 1,IO_intValue(line,chunkPos,1))]
material(e:e+c-1) = [(IO_intValue(line,chunkPos,3),i = 1,IO_intValue(line,chunkPos,1))]
else if (IO_lc(IO_stringValue(line,chunkPos,2)) == 'to') then compression
c = abs(IO_intValue(line,chunkPos,3) - IO_intValue(line,chunkPos,1)) + 1
o = merge(+1, -1, IO_intValue(line,chunkPos,3) > IO_intValue(line,chunkPos,1))
microstructure(e:e+c-1) = [(i, i = IO_intValue(line,chunkPos,1),IO_intValue(line,chunkPos,3),o)]
material(e:e+c-1) = [(i, i = IO_intValue(line,chunkPos,1),IO_intValue(line,chunkPos,3),o)]
else compression
c = chunkPos(1)
microstructure(e:e+c-1) = [(IO_intValue(line,chunkPos,i+1), i=0, c-1)]
material(e:e+c-1) = [(IO_intValue(line,chunkPos,i+1), i=0, c-1)]
endif compression
endif noCompression
@ -308,7 +307,7 @@ end subroutine readGeom
!> @brief Parse vtk rectilinear grid (.vtr)
!> @details https://vtk.org/Wiki/VTK_XML_Formats
!--------------------------------------------------------------------------------------------------
subroutine readVTR(grid,geomSize,origin,microstructure)
subroutine readVTR(grid,geomSize,origin,material)
integer, dimension(3), intent(out) :: &
grid ! grid (across all processes!)
@ -316,7 +315,7 @@ subroutine readVTR(grid,geomSize,origin,microstructure)
geomSize, & ! size (across all processes!)
origin ! origin (across all processes!)
integer, dimension(:), intent(out), allocatable :: &
microstructure
material
character(len=:), allocatable :: fileContent, dataType, headerType
logical :: inFile,inGrid,gotCoordinates,gotCellData,compressed
@ -364,11 +363,9 @@ subroutine readVTR(grid,geomSize,origin,microstructure)
else
if(index(fileContent(startPos:endPos),'<CellData>',kind=pI64) /= 0_pI64) then
gotCellData = .true.
startPos = endPos + 2_pI64
do while (index(fileContent(startPos:endPos),'</CellData>',kind=pI64) == 0_pI64)
endPos = startPos + index(fileContent(startPos:),IO_EOL,kind=pI64) - 2_pI64
if(index(fileContent(startPos:endPos),'<DataArray',kind=pI64) /= 0_pI64 .and. &
getXMLValue(fileContent(startPos:endPos),'Name') == 'materialpoint' ) then
getXMLValue(fileContent(startPos:endPos),'Name') == 'material' ) then
if(getXMLValue(fileContent(startPos:endPos),'format') /= 'binary') &
call IO_error(error_ID = 844, ext_msg='format (materialpoint)')
@ -377,10 +374,11 @@ subroutine readVTR(grid,geomSize,origin,microstructure)
startPos = endPos + 2_pI64
endPos = startPos + index(fileContent(startPos:),IO_EOL,kind=pI64) - 2_pI64
s = startPos + verify(fileContent(startPos:endPos),IO_WHITESPACE,kind=pI64) -1_pI64 ! start (no leading whitespace)
microstructure = as_Int(fileContent(s:endPos),headerType,compressed,dataType)
material = as_Int(fileContent(s:endPos),headerType,compressed,dataType)
exit
endif
startPos = endPos + 2_pI64
endPos = startPos + index(fileContent(startPos:),IO_EOL,kind=pI64) - 2_pI64
enddo
elseif(index(fileContent(startPos:endPos),'<Coordinates>',kind=pI64) /= 0_pI64) then
gotCoordinates = .true.
@ -415,10 +413,10 @@ subroutine readVTR(grid,geomSize,origin,microstructure)
end do
if(.not. allocated(microstructure)) call IO_error(error_ID = 844, ext_msg='materialpoint not found')
if(size(microstructure) /= product(grid)) call IO_error(error_ID = 844, ext_msg='size(materialpoint)')
if(any(geomSize<=0)) call IO_error(error_ID = 844, ext_msg='size')
if(any(grid<1)) call IO_error(error_ID = 844, ext_msg='grid')
if(.not. allocated(material)) call IO_error(error_ID = 844, ext_msg='material data not found')
if(size(material) /= product(grid)) call IO_error(error_ID = 844, ext_msg='size(material)')
if(any(geomSize<=0)) call IO_error(error_ID = 844, ext_msg='size')
if(any(grid<1)) call IO_error(error_ID = 844, ext_msg='grid')
contains

View File

@ -4,20 +4,20 @@
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Relaxed grain cluster (RGC) homogenization scheme
!> Nconstituents is defined as p x q x r (cluster)
!> N_constituents is defined as p x q x r (cluster)
!--------------------------------------------------------------------------------------------------
submodule(homogenization) homogenization_mech_RGC
use rotations
type :: tParameters
integer, dimension(:), allocatable :: &
Nconstituents
N_constituents
real(pReal) :: &
xiAlpha, &
ciAlpha
xi_alpha, &
c_Alpha
real(pReal), dimension(:), allocatable :: &
dAlpha, &
angles
D_alpha, &
a_g
integer :: &
of_debug = 0
character(len=pStringLen), allocatable, dimension(:) :: &
@ -163,20 +163,20 @@ module subroutine mech_RGC_init(num_homogMech)
prm%output = homogMech%get_asStrings('output',defaultVal=emptyStringArray)
#endif
prm%Nconstituents = homogMech%get_asInts('cluster_size',requiredSize=3)
if (homogenization_Ngrains(h) /= product(prm%Nconstituents)) &
prm%N_constituents = homogMech%get_asInts('cluster_size',requiredSize=3)
if (homogenization_Ngrains(h) /= product(prm%N_constituents)) &
call IO_error(211,ext_msg='clustersize (mech_rgc)')
prm%xiAlpha = homogMech%get_asFloat('xi_alpha')
prm%ciAlpha = homogMech%get_asFloat('c_alpha')
prm%xi_alpha = homogMech%get_asFloat('xi_alpha')
prm%c_alpha = homogMech%get_asFloat('c_alpha')
prm%dAlpha = homogMech%get_asFloats('D_alpha', requiredSize=3)
prm%angles = homogMech%get_asFloats('a_g', requiredSize=3)
prm%D_alpha = homogMech%get_asFloats('D_alpha', requiredSize=3)
prm%a_g = homogMech%get_asFloats('a_g', requiredSize=3)
NofMyHomog = count(material_homogenizationAt == h)
nIntFaceTot = 3*( (prm%Nconstituents(1)-1)*prm%Nconstituents(2)*prm%Nconstituents(3) &
+ prm%Nconstituents(1)*(prm%Nconstituents(2)-1)*prm%Nconstituents(3) &
+ prm%Nconstituents(1)*prm%Nconstituents(2)*(prm%Nconstituents(3)-1))
nIntFaceTot = 3*( (prm%N_constituents(1)-1)*prm%N_constituents(2)*prm%N_constituents(3) &
+ prm%N_constituents(1)*(prm%N_constituents(2)-1)*prm%N_constituents(3) &
+ prm%N_constituents(1)*prm%N_constituents(2)*(prm%N_constituents(3)-1))
sizeState = nIntFaceTot &
+ size(['avg constitutive work ','average penalty energy'])
@ -197,8 +197,8 @@ module subroutine mech_RGC_init(num_homogMech)
!--------------------------------------------------------------------------------------------------
! assigning cluster orientations
dependentState(homogenization_typeInstance(h))%orientation = spread(eu2om(prm%angles*inRad),3,NofMyHomog)
!dst%orientation = spread(eu2om(prm%angles*inRad),3,NofMyHomog) ifort version 18.0.1 crashes (for whatever reason)
dependentState(homogenization_typeInstance(h))%orientation = spread(eu2om(prm%a_g*inRad),3,NofMyHomog)
!dst%orientation = spread(eu2om(prm%a_g*inRad),3,NofMyHomog) ifort version 18.0.1 crashes (for whatever reason)
end associate
@ -229,8 +229,8 @@ module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of)
!--------------------------------------------------------------------------------------------------
! compute the deformation gradient of individual grains due to relaxations
F = 0.0_pReal
do iGrain = 1,product(prm%Nconstituents)
iGrain3 = grain1to3(iGrain,prm%Nconstituents)
do iGrain = 1,product(prm%N_constituents)
iGrain3 = grain1to3(iGrain,prm%N_constituents)
do iFace = 1,6
intFace = getInterface(iFace,iGrain3) ! identifying 6 interfaces of each grain
aVect = relaxationVector(intFace,instance,of) ! get the relaxation vectors for each interface from global relaxation vector array
@ -290,7 +290,7 @@ module procedure mech_RGC_updateState
!--------------------------------------------------------------------------------------------------
! get the dimension of the cluster (grains and interfaces)
nGDim = prm%Nconstituents
nGDim = prm%N_constituents
nGrain = product(nGDim)
nIntFaceTot = (nGDim(1)-1)*nGDim(2)*nGDim(3) &
+ nGDim(1)*(nGDim(2)-1)*nGDim(3) &
@ -324,12 +324,12 @@ module procedure mech_RGC_updateState
!------------------------------------------------------------------------------------------------
! computing the residual stress from the balance of traction at all (interior) interfaces
do iNum = 1,nIntFaceTot
faceID = interface1to4(iNum,param(instance)%Nconstituents) ! identifying the interface ID in local coordinate system (4-dimensional index)
faceID = interface1to4(iNum,param(instance)%N_constituents) ! identifying the interface ID in local coordinate system (4-dimensional index)
!--------------------------------------------------------------------------------------------------
! identify the left/bottom/back grain (-|N)
iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate system (3-dimensional index)
iGrN = grain3to1(iGr3N,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
iGrN = grain3to1(iGr3N,param(instance)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceN = getInterface(2*faceID(1),iGr3N)
normN = interfaceNormal(intFaceN,instance,of)
@ -337,7 +337,7 @@ module procedure mech_RGC_updateState
! identify the right/up/front grain (+|P)
iGr3P = iGr3N
iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate system (3-dimensional index)
iGrP = grain3to1(iGr3P,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
iGrP = grain3to1(iGr3P,param(instance)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceP = getInterface(2*faceID(1)-1,iGr3P)
normP = interfaceNormal(intFaceP,instance,of)
@ -393,7 +393,7 @@ module procedure mech_RGC_updateState
!--------------------------------------------------------------------------------------------------
! compute/update the state for postResult, i.e., all energy densities computed by time-integration
do iGrain = 1,product(prm%Nconstituents)
do iGrain = 1,product(prm%N_constituents)
do i = 1,3;do j = 1,3
stt%work(of) = stt%work(of) &
+ P(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal)
@ -450,18 +450,18 @@ module procedure mech_RGC_updateState
! ... of the constitutive stress tangent, assembled from dPdF or material constitutive model "smatrix"
allocate(smatrix(3*nIntFaceTot,3*nIntFaceTot), source=0.0_pReal)
do iNum = 1,nIntFaceTot
faceID = interface1to4(iNum,param(instance)%Nconstituents) ! assembling of local dPdF into global Jacobian matrix
faceID = interface1to4(iNum,param(instance)%N_constituents) ! assembling of local dPdF into global Jacobian matrix
!--------------------------------------------------------------------------------------------------
! identify the left/bottom/back grain (-|N)
iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate sytem
iGrN = grain3to1(iGr3N,param(instance)%Nconstituents) ! translate into global grain ID
iGrN = grain3to1(iGr3N,param(instance)%N_constituents) ! translate into global grain ID
intFaceN = getInterface(2*faceID(1),iGr3N) ! identifying the connecting interface in local coordinate system
normN = interfaceNormal(intFaceN,instance,of)
do iFace = 1,6
intFaceN = getInterface(iFace,iGr3N) ! identifying all interfaces that influence relaxation of the above interface
mornN = interfaceNormal(intFaceN,instance,of)
iMun = interface4to1(intFaceN,param(instance)%Nconstituents) ! translate the interfaces ID into local 4-dimensional index
iMun = interface4to1(intFaceN,param(instance)%N_constituents) ! translate the interfaces ID into local 4-dimensional index
if (iMun > 0) then ! get the corresponding tangent
do i=1,3; do j=1,3; do k=1,3; do l=1,3
smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) &
@ -476,13 +476,13 @@ module procedure mech_RGC_updateState
! identify the right/up/front grain (+|P)
iGr3P = iGr3N
iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate sytem
iGrP = grain3to1(iGr3P,param(instance)%Nconstituents) ! translate into global grain ID
iGrP = grain3to1(iGr3P,param(instance)%N_constituents) ! translate into global grain ID
intFaceP = getInterface(2*faceID(1)-1,iGr3P) ! identifying the connecting interface in local coordinate system
normP = interfaceNormal(intFaceP,instance,of)
do iFace = 1,6
intFaceP = getInterface(iFace,iGr3P) ! identifying all interfaces that influence relaxation of the above interface
mornP = interfaceNormal(intFaceP,instance,of)
iMun = interface4to1(intFaceP,param(instance)%Nconstituents) ! translate the interfaces ID into local 4-dimensional index
iMun = interface4to1(intFaceP,param(instance)%N_constituents) ! translate the interfaces ID into local 4-dimensional index
if (iMun > 0) then ! get the corresponding tangent
do i=1,3; do j=1,3; do k=1,3; do l=1,3
smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) &
@ -522,12 +522,12 @@ module procedure mech_RGC_updateState
! computing the global stress residual array from the perturbed state
p_resid = 0.0_pReal
do iNum = 1,nIntFaceTot
faceID = interface1to4(iNum,param(instance)%Nconstituents) ! identifying the interface ID in local coordinate system (4-dimensional index)
faceID = interface1to4(iNum,param(instance)%N_constituents) ! identifying the interface ID in local coordinate system (4-dimensional index)
!--------------------------------------------------------------------------------------------------
! identify the left/bottom/back grain (-|N)
iGr3N = faceID(2:4) ! identify the grain ID in local coordinate system (3-dimensional index)
iGrN = grain3to1(iGr3N,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
iGrN = grain3to1(iGr3N,param(instance)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceN = getInterface(2*faceID(1),iGr3N) ! identify the interface ID of the grain
normN = interfaceNormal(intFaceN,instance,of)
@ -535,7 +535,7 @@ module procedure mech_RGC_updateState
! identify the right/up/front grain (+|P)
iGr3P = iGr3N
iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identify the grain ID in local coordinate system (3-dimensional index)
iGrP = grain3to1(iGr3P,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
iGrP = grain3to1(iGr3P,param(instance)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceP = getInterface(2*faceID(1)-1,iGr3P) ! identify the interface ID of the grain
normP = interfaceNormal(intFaceP,instance,of)
@ -664,7 +664,7 @@ module procedure mech_RGC_updateState
real(pReal) :: muGrain,muGNghb,nDefNorm,bgGrain,bgGNghb
real(pReal), parameter :: nDefToler = 1.0e-10_pReal
nGDim = param(instance)%Nconstituents
nGDim = param(instance)%N_constituents
rPen = 0.0_pReal
nMis = 0.0_pReal
@ -685,11 +685,11 @@ module procedure mech_RGC_updateState
!-----------------------------------------------------------------------------------------------
! computing the mismatch and penalty stress tensor of all grains
grainLoop: do iGrain = 1,product(prm%Nconstituents)
grainLoop: do iGrain = 1,product(prm%N_constituents)
Gmoduli = equivalentModuli(iGrain,ip,el)
muGrain = Gmoduli(1) ! collecting the equivalent shear modulus of grain
bgGrain = Gmoduli(2) ! and the lengthh of Burgers vector
iGrain3 = grain1to3(iGrain,prm%Nconstituents) ! get the grain ID in local 3-dimensional index (x,y,z)-position
iGrain3 = grain1to3(iGrain,prm%N_constituents) ! get the grain ID in local 3-dimensional index (x,y,z)-position
interfaceLoop: do iFace = 1,6
intFace = getInterface(iFace,iGrain3) ! get the 4-dimensional index of the interface in local numbering system of the grain
@ -699,7 +699,7 @@ module procedure mech_RGC_updateState
+ int(real(intFace(1),pReal)/real(abs(intFace(1)),pReal))
where(iGNghb3 < 1) iGNghb3 = nGDim
where(iGNghb3 >nGDim) iGNghb3 = 1
iGNghb = grain3to1(iGNghb3,prm%Nconstituents) ! get the ID of the neighboring grain
iGNghb = grain3to1(iGNghb3,prm%N_constituents) ! get the ID of the neighboring grain
Gmoduli = equivalentModuli(iGNghb,ip,el) ! collect the shear modulus and Burgers vector of the neighbor
muGNghb = Gmoduli(1)
bgGNghb = Gmoduli(2)
@ -728,9 +728,9 @@ module procedure mech_RGC_updateState
!-------------------------------------------------------------------------------------------
! compute the stress penalty of all interfaces
do i = 1,3; do j = 1,3; do k = 1,3; do l = 1,3
rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*bgGrain + muGNghb*bgGNghb)*prm%xiAlpha &
*surfCorr(abs(intFace(1)))/prm%dAlpha(abs(intFace(1))) &
*cosh(prm%ciAlpha*nDefNorm) &
rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*bgGrain + muGNghb*bgGNghb)*prm%xi_alpha &
*surfCorr(abs(intFace(1)))/prm%D_alpha(abs(intFace(1))) &
*cosh(prm%c_alpha*nDefNorm) &
*0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_LeviCivita(k,l,j) &
*tanh(nDefNorm/num%xSmoo)
enddo; enddo;enddo; enddo
@ -885,8 +885,8 @@ module procedure mech_RGC_updateState
associate(prm => param(instance))
F = 0.0_pReal
do iGrain = 1,product(prm%Nconstituents)
iGrain3 = grain1to3(iGrain,prm%Nconstituents)
do iGrain = 1,product(prm%N_constituents)
iGrain3 = grain1to3(iGrain,prm%N_constituents)
do iFace = 1,6
intFace = getInterface(iFace,iGrain3)
aVect = relaxationVector(intFace,instance,of)
@ -916,8 +916,8 @@ module subroutine mech_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ins
real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
integer, intent(in) :: instance
avgP = sum(P,3) /real(product(param(instance)%Nconstituents),pReal)
dAvgPdAvgF = sum(dPdF,5)/real(product(param(instance)%Nconstituents),pReal)
avgP = sum(P,3) /real(product(param(instance)%N_constituents),pReal)
dAvgPdAvgF = sum(dPdF,5)/real(product(param(instance)%N_constituents),pReal)
end subroutine mech_RGC_averageStressAndItsTangent
@ -975,7 +975,7 @@ pure function relaxationVector(intFace,instance,of)
!--------------------------------------------------------------------------------------------------
! collect the interface relaxation vector from the global state array
iNum = interface4to1(intFace,param(instance)%Nconstituents) ! identify the position of the interface in global state array
iNum = interface4to1(intFace,param(instance)%N_constituents) ! identify the position of the interface in global state array
if (iNum > 0) then
relaxationVector = state(instance)%relaxationVector((3*iNum-2):(3*iNum),of)
else

View File

@ -13,7 +13,7 @@ submodule(homogenization) homogenization_mech_isostrain
type :: tParameters !< container type for internal constitutive parameters
integer :: &
Nconstituents
N_constituents
integer(kind(average_ID)) :: &
mapping
end type
@ -51,7 +51,7 @@ module subroutine mech_isostrain_init
homogMech => homog%get('mech')
associate(prm => param(homogenization_typeInstance(h)))
prm%Nconstituents = homogMech%get_asInt('N_constituents')
prm%N_constituents = homogMech%get_asInt('N_constituents')
select case(homogMech%get_asString('mapping',defaultVal = 'sum'))
case ('sum')
prm%mapping = parallel_ID
@ -107,8 +107,8 @@ module subroutine mech_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dP
avgP = sum(P,3)
dAvgPdAvgF = sum(dPdF,5)
case (average_ID)
avgP = sum(P,3) /real(prm%Nconstituents,pReal)
dAvgPdAvgF = sum(dPdF,5)/real(prm%Nconstituents,pReal)
avgP = sum(P,3) /real(prm%N_constituents,pReal)
dAvgPdAvgF = sum(dPdF,5)/real(prm%N_constituents,pReal)
end select
end associate

View File

@ -12,10 +12,10 @@ submodule(constitutive:constitutive_damage) kinematics_cleavage_opening
integer :: &
sum_N_cl !< total number of cleavage planes
real(pReal) :: &
sdot0, & !< opening rate of cleavage planes
n !< damage rate sensitivity
dot_o, & !< opening rate of cleavage planes
q !< damage rate sensitivity
real(pReal), dimension(:), allocatable :: &
critLoad
g_crit
real(pReal), dimension(:,:,:,:), allocatable :: &
cleavage_systems
end type tParameters
@ -70,21 +70,21 @@ module function kinematics_cleavage_opening_init(kinematics_length) result(myKin
N_cl = kinematic_type%get_asInts('N_cl')
prm%sum_N_cl = sum(abs(N_cl))
prm%n = kinematic_type%get_asFloat('q')
prm%sdot0 = kinematic_type%get_asFloat('dot_o')
prm%q = kinematic_type%get_asFloat('q')
prm%dot_o = kinematic_type%get_asFloat('dot_o')
prm%critLoad = kinematic_type%get_asFloats('g_crit',requiredSize=size(N_cl))
prm%g_crit = kinematic_type%get_asFloats('g_crit',requiredSize=size(N_cl))
prm%cleavage_systems = lattice_SchmidMatrix_cleavage(N_cl,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
! expand: family => system
prm%critLoad = math_expand(prm%critLoad,N_cl)
prm%g_crit = math_expand(prm%g_crit,N_cl)
! sanity checks
if (prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (prm%sdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' dot_o'
if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' g_crit'
if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (prm%dot_o <= 0.0_pReal) extmsg = trim(extmsg)//' dot_o'
if (any(prm%g_crit < 0.0_pReal)) extmsg = trim(extmsg)//' g_crit'
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
@ -128,13 +128,13 @@ module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S,
dLd_dTstar = 0.0_pReal
associate(prm => param(kinematics_cleavage_opening_instance(material_phaseAt(ipc,el))))
do i = 1,prm%sum_N_cl
traction_crit = prm%critLoad(i)* damage(homog)%p(damageOffset)**2.0_pReal
traction_crit = prm%g_crit(i)* damage(homog)%p(damageOffset)**2.0_pReal
traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i))
if (abs(traction_d) > traction_crit + tol_math_check) then
udotd = sign(1.0_pReal,traction_d)* prm%sdot0 * ((abs(traction_d) - traction_crit)/traction_crit)**prm%n
udotd = sign(1.0_pReal,traction_d)* prm%dot_o * ((abs(traction_d) - traction_crit)/traction_crit)**prm%q
Ld = Ld + udotd*prm%cleavage_systems(1:3,1:3,1,i)
dudotd_dt = sign(1.0_pReal,traction_d)*udotd*prm%n / (abs(traction_d) - traction_crit)
dudotd_dt = sign(1.0_pReal,traction_d)*udotd*prm%q / (abs(traction_d) - traction_crit)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) &
+ dudotd_dt*prm%cleavage_systems(k,l,1,i) * prm%cleavage_systems(m,n,1,i)
@ -142,9 +142,9 @@ module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S,
traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i))
if (abs(traction_t) > traction_crit + tol_math_check) then
udott = sign(1.0_pReal,traction_t)* prm%sdot0 * ((abs(traction_t) - traction_crit)/traction_crit)**prm%n
udott = sign(1.0_pReal,traction_t)* prm%dot_o * ((abs(traction_t) - traction_crit)/traction_crit)**prm%q
Ld = Ld + udott*prm%cleavage_systems(1:3,1:3,2,i)
dudott_dt = sign(1.0_pReal,traction_t)*udott*prm%n / (abs(traction_t) - traction_crit)
dudott_dt = sign(1.0_pReal,traction_t)*udott*prm%q / (abs(traction_t) - traction_crit)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) &
+ dudott_dt*prm%cleavage_systems(k,l,2,i) * prm%cleavage_systems(m,n,2,i)
@ -152,9 +152,9 @@ module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S,
traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i))
if (abs(traction_n) > traction_crit + tol_math_check) then
udotn = sign(1.0_pReal,traction_n)* prm%sdot0 * ((abs(traction_n) - traction_crit)/traction_crit)**prm%n
udotn = sign(1.0_pReal,traction_n)* prm%dot_o * ((abs(traction_n) - traction_crit)/traction_crit)**prm%q
Ld = Ld + udotn*prm%cleavage_systems(1:3,1:3,3,i)
dudotn_dt = sign(1.0_pReal,traction_n)*udotn*prm%n / (abs(traction_n) - traction_crit)
dudotn_dt = sign(1.0_pReal,traction_n)*udotn*prm%q / (abs(traction_n) - traction_crit)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) &
+ dudotn_dt*prm%cleavage_systems(k,l,3,i) * prm%cleavage_systems(m,n,3,i)

View File

@ -12,10 +12,10 @@ submodule(constitutive:constitutive_damage) kinematics_slipplane_opening
integer :: &
sum_N_sl !< total number of cleavage planes
real(pReal) :: &
sdot0, & !< opening rate of cleavage planes
n !< damage rate sensitivity
dot_o, & !< opening rate of cleavage planes
q !< damage rate sensitivity
real(pReal), dimension(:), allocatable :: &
critLoad
g_crit
real(pReal), dimension(:,:,:), allocatable :: &
P_d, &
P_t, &
@ -70,8 +70,8 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi
associate(prm => param(kinematics_slipplane_opening_instance(p)))
kinematic_type => kinematics%get(k)
prm%sdot0 = kinematic_type%get_asFloat('dot_o')
prm%n = kinematic_type%get_asFloat('q')
prm%dot_o = kinematic_type%get_asFloat('dot_o')
prm%q = kinematic_type%get_asFloat('q')
N_sl = pl%get_asInts('N_sl')
prm%sum_N_sl = sum(abs(N_sl))
@ -89,15 +89,15 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi
prm%P_n(1:3,1:3,i) = math_outer(n(1:3,i), n(1:3,i))
enddo
prm%critLoad = kinematic_type%get_asFloats('g_crit',requiredSize=size(N_sl))
prm%g_crit = kinematic_type%get_asFloats('g_crit',requiredSize=size(N_sl))
! expand: family => system
prm%critLoad = math_expand(prm%critLoad,N_sl)
prm%g_crit = math_expand(prm%g_crit,N_sl)
! sanity checks
if (prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' anisoDuctile_n'
if (prm%sdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' anisoDuctile_sdot0'
if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' anisoDuctile_critLoad'
if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' anisoDuctile_n'
if (prm%dot_o <= 0.0_pReal) extmsg = trim(extmsg)//' anisoDuctile_sdot0'
if (any(prm%g_crit < 0.0_pReal)) extmsg = trim(extmsg)//' anisoDuctile_critLoad'
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
@ -150,27 +150,27 @@ module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S
traction_t = math_tensordot(S,prm%P_t(1:3,1:3,i))
traction_n = math_tensordot(S,prm%P_n(1:3,1:3,i))
traction_crit = prm%critLoad(i)* damage(homog)%p(damageOffset) ! degrading critical load carrying capacity by damage
traction_crit = prm%g_crit(i)* damage(homog)%p(damageOffset) ! degrading critical load carrying capacity by damage
udotd = sign(1.0_pReal,traction_d)* prm%sdot0* ( abs(traction_d)/traction_crit &
- abs(traction_d)/prm%critLoad(i))**prm%n
udott = sign(1.0_pReal,traction_t)* prm%sdot0* ( abs(traction_t)/traction_crit &
- abs(traction_t)/prm%critLoad(i))**prm%n
udotn = prm%sdot0* ( max(0.0_pReal,traction_n)/traction_crit &
- max(0.0_pReal,traction_n)/prm%critLoad(i))**prm%n
udotd = sign(1.0_pReal,traction_d)* prm%dot_o* ( abs(traction_d)/traction_crit &
- abs(traction_d)/prm%g_crit(i))**prm%q
udott = sign(1.0_pReal,traction_t)* prm%dot_o* ( abs(traction_t)/traction_crit &
- abs(traction_t)/prm%g_crit(i))**prm%q
udotn = prm%dot_o* ( max(0.0_pReal,traction_n)/traction_crit &
- max(0.0_pReal,traction_n)/prm%g_crit(i))**prm%q
if (dNeq0(traction_d)) then
dudotd_dt = udotd*prm%n/traction_d
dudotd_dt = udotd*prm%q/traction_d
else
dudotd_dt = 0.0_pReal
endif
if (dNeq0(traction_t)) then
dudott_dt = udott*prm%n/traction_t
dudott_dt = udott*prm%q/traction_t
else
dudott_dt = 0.0_pReal
endif
if (dNeq0(traction_n)) then
dudotn_dt = udotn*prm%n/traction_n
dudotn_dt = udotn*prm%q/traction_n
else
dudotn_dt = 0.0_pReal
endif

View File

@ -11,7 +11,7 @@ submodule(constitutive:constitutive_thermal) kinematics_thermal_expansion
real(pReal) :: &
T_ref
real(pReal), dimension(3,3,3) :: &
expansion = 0.0_pReal
A = 0.0_pReal
end type tParameters
type(tParameters), dimension(:), allocatable :: param
@ -64,13 +64,13 @@ module function kinematics_thermal_expansion_init(kinematics_length) result(myKi
! read up to three parameters (constant, linear, quadratic with T)
temp = kinematic_type%get_asFloats('A_11')
prm%expansion(1,1,1:size(temp)) = temp
prm%A(1,1,1:size(temp)) = temp
temp = kinematic_type%get_asFloats('A_22',defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp))
prm%expansion(2,2,1:size(temp)) = temp
prm%A(2,2,1:size(temp)) = temp
temp = kinematic_type%get_asFloats('A_33',defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp))
prm%expansion(3,3,1:size(temp)) = temp
do i=1, size(prm%expansion,3)
prm%expansion(1:3,1:3,i) = lattice_applyLatticeSymmetry33(prm%expansion(1:3,1:3,i),&
prm%A(3,3,1:size(temp)) = temp
do i=1, size(prm%A,3)
prm%A(1:3,1:3,i) = lattice_applyLatticeSymmetry33(prm%A(1:3,1:3,i),&
phase%get_asString('lattice'))
enddo
@ -94,13 +94,13 @@ pure module function kinematics_thermal_expansion_initialStrain(homog,phase,offs
offset
real(pReal), dimension(3,3) :: &
initialStrain !< initial thermal strain (should be small strain, though)
initialStrain !< initial thermal strain (should be small strain, though)
associate(prm => param(kinematics_thermal_expansion_instance(phase)))
initialStrain = &
(temperature(homog)%p(offset) - prm%T_ref)**1 / 1. * prm%expansion(1:3,1:3,1) + & ! constant coefficient
(temperature(homog)%p(offset) - prm%T_ref)**2 / 2. * prm%expansion(1:3,1:3,2) + & ! linear coefficient
(temperature(homog)%p(offset) - prm%T_ref)**3 / 3. * prm%expansion(1:3,1:3,3) ! quadratic coefficient
(temperature(homog)%p(offset) - prm%T_ref)**1 / 1. * prm%A(1:3,1:3,1) + & ! constant coefficient
(temperature(homog)%p(offset) - prm%T_ref)**2 / 2. * prm%A(1:3,1:3,2) + & ! linear coefficient
(temperature(homog)%p(offset) - prm%T_ref)**3 / 3. * prm%A(1:3,1:3,3) ! quadratic coefficient
end associate
end function kinematics_thermal_expansion_initialStrain
@ -133,14 +133,14 @@ module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, i
associate(prm => param(kinematics_thermal_expansion_instance(phase)))
Li = TDot * ( &
prm%expansion(1:3,1:3,1)*(T - prm%T_ref)**0 & ! constant coefficient
+ prm%expansion(1:3,1:3,2)*(T - prm%T_ref)**1 & ! linear coefficient
+ prm%expansion(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient
prm%A(1:3,1:3,1)*(T - prm%T_ref)**0 & ! constant coefficient
+ prm%A(1:3,1:3,2)*(T - prm%T_ref)**1 & ! linear coefficient
+ prm%A(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient
) / &
(1.0_pReal &
+ prm%expansion(1:3,1:3,1)*(T - prm%T_ref)**1 / 1. &
+ prm%expansion(1:3,1:3,2)*(T - prm%T_ref)**2 / 2. &
+ prm%expansion(1:3,1:3,3)*(T - prm%T_ref)**3 / 3. &
+ prm%A(1:3,1:3,1)*(T - prm%T_ref)**1 / 1. &
+ prm%A(1:3,1:3,2)*(T - prm%T_ref)**2 / 2. &
+ prm%A(1:3,1:3,3)*(T - prm%T_ref)**3 / 3. &
)
end associate
dLi_dTstar = 0.0_pReal

View File

@ -394,13 +394,13 @@ module lattice
! SHOULD NOT BE PART OF LATTICE BEGIN
real(pReal), dimension(:), allocatable, public, protected :: &
lattice_mu, lattice_nu, &
lattice_damageMobility, &
lattice_massDensity, &
lattice_specificHeat
lattice_M, &
lattice_rho, &
lattice_c_p
real(pReal), dimension(:,:,:), allocatable, public, protected :: &
lattice_C66, &
lattice_thermalConductivity, &
lattice_damageDiffusion
lattice_K, &
lattice_D
integer(kind(lattice_UNDEFINED_ID)), dimension(:), allocatable, public, protected :: &
lattice_structure
! SHOULD NOT BE PART OF LATTICE END
@ -465,11 +465,11 @@ subroutine lattice_init
allocate(lattice_structure(Nphases),source = lattice_UNDEFINED_ID)
allocate(lattice_C66(6,6,Nphases), source=0.0_pReal)
allocate(lattice_thermalConductivity (3,3,Nphases), source=0.0_pReal)
allocate(lattice_damageDiffusion (3,3,Nphases), source=0.0_pReal)
allocate(lattice_K (3,3,Nphases), source=0.0_pReal)
allocate(lattice_D (3,3,Nphases), source=0.0_pReal)
allocate(lattice_damageMobility,&
lattice_massDensity,lattice_specificHeat, &
allocate(lattice_M,&
lattice_rho,lattice_c_p, &
lattice_mu, lattice_nu,&
source=[(0.0_pReal,i=1,Nphases)])
@ -517,22 +517,22 @@ subroutine lattice_init
! SHOULD NOT BE PART OF LATTICE BEGIN
lattice_thermalConductivity(1,1,p) = phase%get_asFloat('K_11',defaultVal=0.0_pReal)
lattice_thermalConductivity(2,2,p) = phase%get_asFloat('K_22',defaultVal=0.0_pReal)
lattice_thermalConductivity(3,3,p) = phase%get_asFloat('K_33',defaultVal=0.0_pReal)
lattice_thermalConductivity(1:3,1:3,p) = lattice_applyLatticeSymmetry33(lattice_thermalConductivity(1:3,1:3,p), &
phase%get_asString('lattice'))
lattice_K(1,1,p) = phase%get_asFloat('K_11',defaultVal=0.0_pReal)
lattice_K(2,2,p) = phase%get_asFloat('K_22',defaultVal=0.0_pReal)
lattice_K(3,3,p) = phase%get_asFloat('K_33',defaultVal=0.0_pReal)
lattice_K(1:3,1:3,p) = lattice_applyLatticeSymmetry33(lattice_K(1:3,1:3,p), &
phase%get_asString('lattice'))
lattice_specificHeat(p) = phase%get_asFloat('c_p',defaultVal=0.0_pReal)
lattice_massDensity(p) = phase%get_asFloat('rho', defaultVal=0.0_pReal)
lattice_c_p(p) = phase%get_asFloat('c_p', defaultVal=0.0_pReal)
lattice_rho(p) = phase%get_asFloat('rho', defaultVal=0.0_pReal)
lattice_DamageDiffusion(1,1,p) = phase%get_asFloat('D_11',defaultVal=0.0_pReal)
lattice_DamageDiffusion(2,2,p) = phase%get_asFloat('D_22',defaultVal=0.0_pReal)
lattice_DamageDiffusion(3,3,p) = phase%get_asFloat('D_33',defaultVal=0.0_pReal)
lattice_DamageDiffusion(1:3,1:3,p) = lattice_applyLatticeSymmetry33(lattice_DamageDiffusion(1:3,1:3,p), &
phase%get_asString('lattice'))
lattice_D(1,1,p) = phase%get_asFloat('D_11',defaultVal=0.0_pReal)
lattice_D(2,2,p) = phase%get_asFloat('D_22',defaultVal=0.0_pReal)
lattice_D(3,3,p) = phase%get_asFloat('D_33',defaultVal=0.0_pReal)
lattice_D(1:3,1:3,p) = lattice_applyLatticeSymmetry33(lattice_D(1:3,1:3,p), &
phase%get_asString('lattice'))
lattice_DamageMobility(p) = phase%get_asFloat('M',defaultVal=0.0_pReal)
lattice_M(p) = phase%get_asFloat('M',defaultVal=0.0_pReal)
! SHOULD NOT BE PART OF LATTICE END
call selfTest

View File

@ -52,7 +52,7 @@ subroutine discretization_marc_init
type(tElement) :: elem
integer, dimension(:), allocatable :: &
microstructureAt
materialAt
integer:: &
Nnodes, & !< total number of nodes in the mesh
Nelems, & !< total number of elements in the mesh
@ -83,7 +83,7 @@ subroutine discretization_marc_init
mesh_unitlength = num_commercialFEM%get_asFloat('unitlength',defaultVal=1.0_pReal) ! set physical extent of a length unit in mesh
if (mesh_unitlength <= 0.0_pReal) call IO_error(301,ext_msg='unitlength')
call inputRead(elem,node0_elem,connectivity_elem,microstructureAt)
call inputRead(elem,node0_elem,connectivity_elem,materialAt)
nElems = size(connectivity_elem,2)
if (debug_e < 1 .or. debug_e > nElems) call IO_error(602,ext_msg='element')
@ -103,7 +103,7 @@ subroutine discretization_marc_init
call buildIPcoordinates(IP_reshaped,reshape(connectivity_cell,[elem%NcellNodesPerCell,&
elem%nIPs*nElems]),node0_cell)
call discretization_init(microstructureAt,&
call discretization_init(materialAt,&
IP_reshaped,&
node0_cell)
@ -172,7 +172,7 @@ end subroutine writeGeometry
!--------------------------------------------------------------------------------------------------
!> @brief Read mesh from marc input file
!--------------------------------------------------------------------------------------------------
subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt)
subroutine inputRead(elem,node0_elem,connectivity_elem,materialAt)
type(tElement), intent(out) :: elem
real(pReal), dimension(:,:), allocatable, intent(out) :: &
@ -180,7 +180,7 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt)
integer, dimension(:,:), allocatable, intent(out) :: &
connectivity_elem
integer, dimension(:), allocatable, intent(out) :: &
microstructureAt
materialAt
integer :: &
fileFormatVersion, &
@ -226,9 +226,9 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt)
connectivity_elem = inputRead_connectivityElem(nElems,elem%nNodes,inputFile)
call inputRead_microstructure(microstructureAt, &
nElems,elem%nNodes,nameElemSet,mapElemSet,&
initialcondTableStyle,inputFile)
call inputRead_material(materialAt, &
nElems,elem%nNodes,nameElemSet,mapElemSet,&
initialcondTableStyle,inputFile)
end subroutine inputRead
@ -675,13 +675,13 @@ end function inputRead_connectivityElem
!--------------------------------------------------------------------------------------------------
!> @brief Store microstructure ID
!> @brief Store material ID
!--------------------------------------------------------------------------------------------------
subroutine inputRead_microstructure(microstructureAt,&
nElem,nNodes,nameElemSet,mapElemSet,initialcondTableStyle,fileContent)
subroutine inputRead_material(materialAt,&
nElem,nNodes,nameElemSet,mapElemSet,initialcondTableStyle,fileContent)
integer, dimension(:), allocatable, intent(out) :: &
microstructureAt
materialAt
integer, intent(in) :: &
nElem, &
nNodes, & !< number of nodes per element
@ -696,7 +696,7 @@ subroutine inputRead_microstructure(microstructureAt,&
integer :: i,j,t,sv,myVal,e,nNodesAlreadyRead,l,k,m
allocate(microstructureAt(nElem),source=0)
allocate(materialAt(nElem),source=0)
do l = 1, size(fileContent)
chunkPos = IO_stringPos(fileContent(l))
@ -715,7 +715,7 @@ subroutine inputRead_microstructure(microstructureAt,&
contInts = continuousIntValues(fileContent(l+k+m+1:),nElem,nameElemSet,mapElemSet,size(nameElemSet)) ! get affected elements
do i = 1,contInts(1)
e = mesh_FEM2DAMASK_elem(contInts(1+i))
microstructureAt(e) = myVal
materialAt(e) = myVal
enddo
if (initialcondTableStyle == 0) m = m + 1
enddo
@ -723,9 +723,9 @@ subroutine inputRead_microstructure(microstructureAt,&
endif
enddo
if(any(microstructureAt < 1)) call IO_error(180)
if(any(materialAt < 1)) call IO_error(180)
end subroutine inputRead_microstructure
end subroutine inputRead_material
!--------------------------------------------------------------------------------------------------

View File

@ -12,11 +12,11 @@ submodule (constitutive:constitutive_damage) source_damage_anisoBrittle
type :: tParameters !< container type for internal constitutive parameters
real(pReal) :: &
sdot_0, & !< opening rate of cleavage planes
n !< damage rate sensitivity
dot_o, & !< opening rate of cleavage planes
q !< damage rate sensitivity
real(pReal), dimension(:), allocatable :: &
critDisp, & !< critical displacement
critLoad !< critical load
s_crit, & !< critical displacement
g_crit !< critical load
real(pReal), dimension(:,:,:,:), allocatable :: &
cleavage_systems
integer :: &
@ -75,18 +75,18 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources)
N_cl = src%get_asInts('N_cl',defaultVal=emptyIntArray)
prm%sum_N_cl = sum(abs(N_cl))
prm%n = src%get_asFloat('q')
prm%sdot_0 = src%get_asFloat('dot_o')
prm%q = src%get_asFloat('q')
prm%dot_o = src%get_asFloat('dot_o')
prm%critDisp = src%get_asFloats('s_crit', requiredSize=size(N_cl))
prm%critLoad = src%get_asFloats('g_crit', requiredSize=size(N_cl))
prm%s_crit = src%get_asFloats('s_crit', requiredSize=size(N_cl))
prm%g_crit = src%get_asFloats('g_crit', requiredSize=size(N_cl))
prm%cleavage_systems = lattice_SchmidMatrix_cleavage(N_cl,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
! expand: family => system
prm%critDisp = math_expand(prm%critDisp,N_cl)
prm%critLoad = math_expand(prm%critLoad,N_cl)
prm%s_crit = math_expand(prm%s_crit,N_cl)
prm%g_crit = math_expand(prm%g_crit,N_cl)
#if defined (__GFORTRAN__)
prm%output = output_asStrings(src)
@ -95,10 +95,10 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources)
#endif
! sanity checks
if (prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (prm%sdot_0 <= 0.0_pReal) extmsg = trim(extmsg)//' dot_o'
if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' g_crit'
if (any(prm%critDisp < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit'
if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (prm%dot_o <= 0.0_pReal) extmsg = trim(extmsg)//' dot_o'
if (any(prm%g_crit < 0.0_pReal)) extmsg = trim(extmsg)//' g_crit'
if (any(prm%s_crit < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit'
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP
call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0)
@ -152,14 +152,14 @@ module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el)
traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i))
traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i))
traction_crit = prm%critLoad(i)*damage(homog)%p(damageOffset)**2.0_pReal
traction_crit = prm%g_crit(i)*damage(homog)%p(damageOffset)**2.0_pReal
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &
= sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &
+ prm%sdot_0 / prm%critDisp(i) &
* ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%n + &
(max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%n + &
(max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**prm%n)
+ prm%dot_o / prm%s_crit(i) &
* ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%q + &
(max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%q + &
(max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**prm%q)
enddo
end associate

View File

@ -12,9 +12,9 @@ submodule(constitutive:constitutive_damage) source_damage_anisoDuctile
type :: tParameters !< container type for internal constitutive parameters
real(pReal) :: &
n !< damage rate sensitivity
q !< damage rate sensitivity
real(pReal), dimension(:), allocatable :: &
critPlasticStrain !< critical plastic strain per slip system
gamma_crit !< critical plastic strain per slip system
character(len=pStringLen), allocatable, dimension(:) :: &
output
end type tParameters
@ -67,12 +67,12 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources)
associate(prm => param(source_damage_anisoDuctile_instance(p)))
src => sources%get(sourceOffset)
N_sl = pl%get_asInts('N_sl',defaultVal=emptyIntArray)
prm%n = src%get_asFloat('q')
prm%critPlasticStrain = src%get_asFloats('gamma_crit',requiredSize=size(N_sl))
N_sl = pl%get_asInts('N_sl',defaultVal=emptyIntArray)
prm%q = src%get_asFloat('q')
prm%gamma_crit = src%get_asFloats('gamma_crit',requiredSize=size(N_sl))
! expand: family => system
prm%critPlasticStrain = math_expand(prm%critPlasticStrain,N_sl)
prm%gamma_crit = math_expand(prm%gamma_crit,N_sl)
#if defined (__GFORTRAN__)
prm%output = output_asStrings(src)
@ -81,8 +81,8 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources)
#endif
! sanity checks
if (prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (any(prm%critPlasticStrain < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit'
if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit'
NipcMyPhase=count(material_phaseAt==p) * discretization_nIP
call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0)
@ -127,7 +127,7 @@ module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el)
associate(prm => param(source_damage_anisoDuctile_instance(phase)))
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &
= sum(plasticState(phase)%slipRate(:,constituent)/(damage(homog)%p(damageOffset)**prm%n)/prm%critPlasticStrain)
= sum(plasticState(phase)%slipRate(:,constituent)/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit)
end associate
end subroutine source_damage_anisoDuctile_dotState

View File

@ -12,8 +12,8 @@ submodule (constitutive:constitutive_damage) source_damage_isoDuctile
type:: tParameters !< container type for internal constitutive parameters
real(pReal) :: &
critPlasticStrain, & !< critical plastic strain
N
gamma_crit, & !< critical plastic strain
q
character(len=pStringLen), allocatable, dimension(:) :: &
output
end type tParameters
@ -64,8 +64,8 @@ module function source_damage_isoDuctile_init(source_length) result(mySources)
associate(prm => param(source_damage_isoDuctile_instance(p)))
src => sources%get(sourceOffset)
prm%N = src%get_asFloat('q')
prm%critPlasticStrain = src%get_asFloat('gamma_crit')
prm%q = src%get_asFloat('q')
prm%gamma_crit = src%get_asFloat('gamma_crit')
#if defined (__GFORTRAN__)
prm%output = output_asStrings(src)
@ -74,8 +74,8 @@ module function source_damage_isoDuctile_init(source_length) result(mySources)
#endif
! sanity checks
if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (prm%critPlasticStrain <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit'
if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit'
NipcMyPhase=count(material_phaseAt==p) * discretization_nIP
call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0)
@ -120,7 +120,7 @@ module subroutine source_damage_isoDuctile_dotState(ipc, ip, el)
associate(prm => param(source_damage_isoDuctile_instance(phase)))
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = &
sum(plasticState(phase)%slipRate(:,constituent))/(damage(homog)%p(damageOffset)**prm%N)/prm%critPlasticStrain
sum(plasticState(phase)%slipRate(:,constituent))/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit
end associate
end subroutine source_damage_isoDuctile_dotState

View File

@ -13,8 +13,8 @@ submodule(constitutive:constitutive_thermal) source_thermal_externalheat
type :: tParameters !< container type for internal constitutive parameters
real(pReal), dimension(:), allocatable :: &
time, &
heat_rate
t_n, &
f_T
integer :: &
nIntervals
end type tParameters
@ -64,10 +64,10 @@ module function source_thermal_externalheat_init(source_length) result(mySources
associate(prm => param(source_thermal_externalheat_instance(p)))
src => sources%get(sourceOffset)
prm%time = src%get_asFloats('t_n')
prm%nIntervals = size(prm%time) - 1
prm%t_n = src%get_asFloats('t_n')
prm%nIntervals = size(prm%t_n) - 1
prm%heat_rate = src%get_asFloats('f_T',requiredSize = size(prm%time))
prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n))
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP
call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0)
@ -121,13 +121,13 @@ module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_d
associate(prm => param(source_thermal_externalheat_instance(phase)))
do interval = 1, prm%nIntervals ! scan through all rate segments
frac_time = (sourceState(phase)%p(sourceOffset)%state(1,of) - prm%time(interval)) &
/ (prm%time(interval+1) - prm%time(interval)) ! fractional time within segment
frac_time = (sourceState(phase)%p(sourceOffset)%state(1,of) - prm%t_n(interval)) &
/ (prm%t_n(interval+1) - prm%t_n(interval)) ! fractional time within segment
if ( (frac_time < 0.0_pReal .and. interval == 1) &
.or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) &
.or. (frac_time >= 0.0_pReal .and. frac_time < 1.0_pReal) ) &
TDot = prm%heat_rate(interval ) * (1.0_pReal - frac_time) + &
prm%heat_rate(interval+1) * frac_time ! interpolate heat rate between segment boundaries...
TDot = prm%f_T(interval ) * (1.0_pReal - frac_time) + &
prm%f_T(interval+1) * frac_time ! interpolate heat rate between segment boundaries...
! ...or extrapolate if outside of bounds
enddo
dTDot_dT = 0.0

View File

@ -169,7 +169,7 @@ function thermal_adiabatic_getSpecificHeat(ip,el)
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el))
thermal_adiabatic_getSpecificHeat = thermal_adiabatic_getSpecificHeat &
+ lattice_specificHeat(material_phaseAt(grain,el))
+ lattice_c_p(material_phaseAt(grain,el))
enddo
thermal_adiabatic_getSpecificHeat = thermal_adiabatic_getSpecificHeat &
@ -195,7 +195,7 @@ function thermal_adiabatic_getMassDensity(ip,el)
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el))
thermal_adiabatic_getMassDensity = thermal_adiabatic_getMassDensity &
+ lattice_massDensity(material_phaseAt(grain,el))
+ lattice_rho(material_phaseAt(grain,el))
enddo
thermal_adiabatic_getMassDensity = thermal_adiabatic_getMassDensity &

View File

@ -127,7 +127,7 @@ function thermal_conduction_getConductivity(ip,el)
thermal_conduction_getConductivity = 0.0_pReal
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el))
thermal_conduction_getConductivity = thermal_conduction_getConductivity + &
crystallite_push33ToRef(grain,ip,el,lattice_thermalConductivity(:,:,material_phaseAt(grain,el)))
crystallite_push33ToRef(grain,ip,el,lattice_K(:,:,material_phaseAt(grain,el)))
enddo
thermal_conduction_getConductivity = thermal_conduction_getConductivity &
@ -153,7 +153,7 @@ function thermal_conduction_getSpecificHeat(ip,el)
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el))
thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat &
+ lattice_specificHeat(material_phaseAt(grain,el))
+ lattice_c_p(material_phaseAt(grain,el))
enddo
thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat &
@ -180,7 +180,7 @@ function thermal_conduction_getMassDensity(ip,el)
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el))
thermal_conduction_getMassDensity = thermal_conduction_getMassDensity &
+ lattice_massDensity(material_phaseAt(grain,el))
+ lattice_rho(material_phaseAt(grain,el))
enddo
thermal_conduction_getMassDensity = thermal_conduction_getMassDensity &