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

Merge branch 'development' into vectorize_rotation

This commit is contained in:
Martin Diehl 2020-04-24 06:31:14 +02:00
parent 7efe14be35 05fd943cf4
commit 042f64200c
36 changed files with 529 additions and 586 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitlab-ci.yml
View File

@ -203,7 +203,6 @@ Post_OrientationConversion:
stage: postprocessing stage: postprocessing
script: script:
- OrientationConversion/test.py - OrientationConversion/test.py
- OrientationConversion/test2.py
except: except:
- master - master
- release - release

2
VERSION
View File

@ -1 +1 @@
v2.0.3-2303-g2a6132b7 v2.0.3-2364-g62f7363a

2
processing/post/DADF5_postResults.py
View File

@ -33,7 +33,7 @@ for filename in options.filenames:
results = damask.Result(filename) results = damask.Result(filename)
if not results.structured: continue if not results.structured: continue
coords = damask.grid_filters.cell_coord0(results.grid,results.size,results.origin) coords = damask.grid_filters.cell_coord0(results.grid,results.size,results.origin).reshape(-1,3,order='F')
N_digits = int(np.floor(np.log10(int(results.increments[-1][3:]))))+1 N_digits = int(np.floor(np.log10(int(results.increments[-1][3:]))))+1
N_digits = 5 # hack to keep test intact N_digits = 5 # hack to keep test intact

75
processing/post/addCompatibilityMismatch.py
View File

@ -16,8 +16,8 @@ scriptID = ' '.join([scriptName,damask.version])
def volTetrahedron(coords): def volTetrahedron(coords):
""" """
Return the volume of the tetrahedron with given vertices or sides. Return the volume of the tetrahedron with given vertices or sides.
Ifvertices are given they must be in a NumPy array with shape (4,3): the If vertices are given they must be in a NumPy array with shape (4,3): the
position vectors of the 4 vertices in 3 dimensions; if the six sides are position vectors of the 4 vertices in 3 dimensions; if the six sides are
given, they must be an array of length 6. If both are given, the sides given, they must be an array of length 6. If both are given, the sides
will be used in the calculation. will be used in the calculation.
@ -62,19 +62,18 @@ def volTetrahedron(coords):
def volumeMismatch(size,F,nodes): def volumeMismatch(size,F,nodes):
""" """
Calculates the volume mismatch. Calculates the volume mismatch.
volume mismatch is defined as the difference between volume of reconstructed volume mismatch is defined as the difference between volume of reconstructed
(compatible) cube and determinant of deformation gradient at Fourier point. (compatible) cube and determinant of deformation gradient at Fourier point.
""" """
coords = np.empty([8,3]) coords = np.empty([8,3])
vMismatch = np.empty(grid[::-1]) vMismatch = np.empty(F.shape[:3])
volInitial = size.prod()/grid.prod()
#-------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------
# calculate actual volume and volume resulting from deformation gradient # calculate actual volume and volume resulting from deformation gradient
for k in range(grid[2]): for k in range(grid[0]):
for j in range(grid[1]): for j in range(grid[1]):
for i in range(grid[0]): for i in range(grid[2]):
coords[0,0:3] = nodes[k, j, i ,0:3] coords[0,0:3] = nodes[k, j, i ,0:3]
coords[1,0:3] = nodes[k ,j, i+1,0:3] coords[1,0:3] = nodes[k ,j, i+1,0:3]
coords[2,0:3] = nodes[k ,j+1,i+1,0:3] coords[2,0:3] = nodes[k ,j+1,i+1,0:3]
@ -91,47 +90,45 @@ def volumeMismatch(size,F,nodes):
+ abs(volTetrahedron([coords[6,0:3],coords[4,0:3],coords[1,0:3],coords[5,0:3]])) \ + abs(volTetrahedron([coords[6,0:3],coords[4,0:3],coords[1,0:3],coords[5,0:3]])) \
+ abs(volTetrahedron([coords[6,0:3],coords[4,0:3],coords[1,0:3],coords[0,0:3]]))) \ + abs(volTetrahedron([coords[6,0:3],coords[4,0:3],coords[1,0:3],coords[0,0:3]]))) \
/np.linalg.det(F[k,j,i,0:3,0:3]) /np.linalg.det(F[k,j,i,0:3,0:3])
return vMismatch/volInitial return vMismatch/(size.prod()/grid.prod())
def shapeMismatch(size,F,nodes,centres): def shapeMismatch(size,F,nodes,centres):
""" """
Routine to calculate the shape mismatch. Routine to calculate the shape mismatch.
shape mismatch is defined as difference between the vectors from the central point to shape mismatch is defined as difference between the vectors from the central point to
the corners of reconstructed (combatible) volume element and the vectors calculated by deforming the corners of reconstructed (combatible) volume element and the vectors calculated by deforming
the initial volume element with the current deformation gradient. the initial volume element with the current deformation gradient.
""" """
coordsInitial = np.empty([8,3]) sMismatch = np.empty(F.shape[:3])
sMismatch = np.empty(grid[::-1])
#-------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------
# initial positions # initial positions
coordsInitial[0,0:3] = [-size[0]/grid[0],-size[1]/grid[1],-size[2]/grid[2]] delta = size/grid*.5
coordsInitial[1,0:3] = [+size[0]/grid[0],-size[1]/grid[1],-size[2]/grid[2]] coordsInitial = np.vstack((delta * np.array((-1,-1,-1)),
coordsInitial[2,0:3] = [+size[0]/grid[0],+size[1]/grid[1],-size[2]/grid[2]] delta * np.array((+1,-1,-1)),
coordsInitial[3,0:3] = [-size[0]/grid[0],+size[1]/grid[1],-size[2]/grid[2]] delta * np.array((+1,+1,-1)),
coordsInitial[4,0:3] = [-size[0]/grid[0],-size[1]/grid[1],+size[2]/grid[2]] delta * np.array((-1,+1,-1)),
coordsInitial[5,0:3] = [+size[0]/grid[0],-size[1]/grid[1],+size[2]/grid[2]] delta * np.array((-1,-1,+1)),
coordsInitial[6,0:3] = [+size[0]/grid[0],+size[1]/grid[1],+size[2]/grid[2]] delta * np.array((+1,-1,+1)),
coordsInitial[7,0:3] = [-size[0]/grid[0],+size[1]/grid[1],+size[2]/grid[2]] delta * np.array((+1,+1,+1)),
coordsInitial = coordsInitial/2.0 delta * np.array((-1,+1,+1))))
#-------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------
# compare deformed original and deformed positions to actual positions # compare deformed original and deformed positions to actual positions
for k in range(grid[2]): for k in range(grid[0]):
for j in range(grid[1]): for j in range(grid[1]):
for i in range(grid[0]): for i in range(grid[2]):
sMismatch[k,j,i] = \ sMismatch[k,j,i] = \
+ np.linalg.norm(nodes[k, j, i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[0,0:3]))\ + np.linalg.norm(nodes[k, j, i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[0,0:3]))\
+ np.linalg.norm(nodes[k, j, i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[1,0:3]))\ + np.linalg.norm(nodes[k+1,j, i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[1,0:3]))\
+ np.linalg.norm(nodes[k, j+1,i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[2,0:3]))\ + np.linalg.norm(nodes[k+1,j+1,i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[2,0:3]))\
+ np.linalg.norm(nodes[k, j+1,i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[3,0:3]))\ + np.linalg.norm(nodes[k, j+1,i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[3,0:3]))\
+ np.linalg.norm(nodes[k+1,j, i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[4,0:3]))\ + np.linalg.norm(nodes[k, j, i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[4,0:3]))\
+ np.linalg.norm(nodes[k+1,j, i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[5,0:3]))\ + np.linalg.norm(nodes[k+1,j, i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[5,0:3]))\
+ np.linalg.norm(nodes[k+1,j+1,i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[6,0:3]))\ + np.linalg.norm(nodes[k+1,j+1,i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[6,0:3]))\
+ np.linalg.norm(nodes[k+1,j+1,i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[7,0:3])) + np.linalg.norm(nodes[k ,j+1,i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[7,0:3]))
return sMismatch return sMismatch
@ -174,24 +171,24 @@ if filenames == []: filenames = [None]
for name in filenames: for name in filenames:
damask.util.report(scriptName,name) damask.util.report(scriptName,name)
table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
grid,size,origin = damask.grid_filters.cell_coord0_gridSizeOrigin(table.get(options.pos)) grid,size,origin = damask.grid_filters.cell_coord0_gridSizeOrigin(table.get(options.pos))
F = table.get(options.defgrad).reshape(grid[2],grid[1],grid[0],3,3) F = table.get(options.defgrad).reshape(tuple(grid)+(-1,),order='F').reshape(tuple(grid)+(3,3))
nodes = damask.grid_filters.node_coord(size,F) nodes = damask.grid_filters.node_coord(size,F)
if options.shape: if options.shape:
centers = damask.grid_filters.cell_coord(size,F) centers = damask.grid_filters.cell_coord(size,F)
shapeMismatch = shapeMismatch( size,table.get(options.defgrad).reshape(grid[2],grid[1],grid[0],3,3),nodes,centers) shapeMismatch = shapeMismatch(size,F,nodes,centers)
table.add('shapeMismatch(({}))'.format(options.defgrad), table.add('shapeMismatch(({}))'.format(options.defgrad),
shapeMismatch.reshape(-1,1), shapeMismatch.reshape(-1,1,order='F'),
scriptID+' '+' '.join(sys.argv[1:])) scriptID+' '+' '.join(sys.argv[1:]))
if options.volume: if options.volume:
volumeMismatch = volumeMismatch(size,table.get(options.defgrad).reshape(grid[2],grid[1],grid[0],3,3),nodes) volumeMismatch = volumeMismatch(size,F,nodes)
table.add('volMismatch(({}))'.format(options.defgrad), table.add('volMismatch(({}))'.format(options.defgrad),
volumeMismatch.reshape(-1,1), volumeMismatch.reshape(-1,1,order='F'),
scriptID+' '+' '.join(sys.argv[1:])) scriptID+' '+' '.join(sys.argv[1:]))
table.to_ASCII(sys.stdout if name is None else name) table.to_ASCII(sys.stdout if name is None else name)

7
processing/post/addCurl.py
View File

@ -49,9 +49,10 @@ for name in filenames:
for label in options.labels: for label in options.labels:
field = table.get(label) field = table.get(label)
shape = (3,) if np.prod(field.shape)//np.prod(grid) == 3 else (3,3) # vector or tensor shape = (3,) if np.prod(field.shape)//np.prod(grid) == 3 else (3,3) # vector or tensor
field = field.reshape(np.append(grid[::-1],shape)) field = field.reshape(tuple(grid)+(-1,),order='F').reshape(tuple(grid)+shape)
curl = damask.grid_filters.curl(size,field)
table.add('curlFFT({})'.format(label), table.add('curlFFT({})'.format(label),
damask.grid_filters.curl(size[::-1],field).reshape(-1,np.prod(shape)), curl.reshape(tuple(grid)+(-1,)).reshape(-1,np.prod(shape),order='F'),
scriptID+' '+' '.join(sys.argv[1:])) scriptID+' '+' '.join(sys.argv[1:]))
table.to_ASCII(sys.stdout if name is None else name) table.to_ASCII(sys.stdout if name is None else name)

16
processing/post/addDisplacement.py
View File

@ -5,8 +5,6 @@ import sys
from io import StringIO from io import StringIO
from optparse import OptionParser from optparse import OptionParser
import numpy as np
import damask import damask
@ -51,23 +49,23 @@ for name in filenames:
table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
grid,size,origin = damask.grid_filters.cell_coord0_gridSizeOrigin(table.get(options.pos)) grid,size,origin = damask.grid_filters.cell_coord0_gridSizeOrigin(table.get(options.pos))
F = table.get(options.f).reshape(np.append(grid[::-1],(3,3))) F = table.get(options.f).reshape(tuple(grid)+(-1,),order='F').reshape(tuple(grid)+(3,3))
if options.nodal: if options.nodal:
table = damask.Table(damask.grid_filters.node_coord0(grid[::-1],size[::-1]).reshape(-1,3), table = damask.Table(damask.grid_filters.node_coord0(grid,size).reshape(-1,3,order='F'),
{'pos':(3,)}) {'pos':(3,)})
table.add('avg({}).{}'.format(options.f,options.pos), table.add('avg({}).{}'.format(options.f,options.pos),
damask.grid_filters.node_displacement_avg(size[::-1],F).reshape(-1,3), damask.grid_filters.node_displacement_avg(size,F).reshape(-1,3,order='F'),
scriptID+' '+' '.join(sys.argv[1:])) scriptID+' '+' '.join(sys.argv[1:]))
table.add('fluct({}).{}'.format(options.f,options.pos), table.add('fluct({}).{}'.format(options.f,options.pos),
damask.grid_filters.node_displacement_fluct(size[::-1],F).reshape(-1,3), damask.grid_filters.node_displacement_fluct(size,F).reshape(-1,3,order='F'),
scriptID+' '+' '.join(sys.argv[1:])) scriptID+' '+' '.join(sys.argv[1:]))
table.to_ASCII(sys.stdout if name is None else os.path.splitext(name)[0]+'_nodal.txt') table.to_ASCII(sys.stdout if name is None else os.path.splitext(name)[0]+'_nodal.txt')
else: else:
table.add('avg({}).{}'.format(options.f,options.pos), table.add('avg({}).{}'.format(options.f,options.pos),
damask.grid_filters.cell_displacement_avg(size[::-1],F).reshape(-1,3), damask.grid_filters.cell_displacement_avg(size,F).reshape(-1,3,order='F'),
scriptID+' '+' '.join(sys.argv[1:])) scriptID+' '+' '.join(sys.argv[1:]))
table.add('fluct({}).{}'.format(options.f,options.pos), table.add('fluct({}).{}'.format(options.f,options.pos),
damask.grid_filters.cell_displacement_fluct(size[::-1],F).reshape(-1,3), damask.grid_filters.cell_displacement_fluct(size,F).reshape(-1,3,order='F'),
scriptID+' '+' '.join(sys.argv[1:])) scriptID+' '+' '.join(sys.argv[1:]))
table.to_ASCII(sys.stdout if name is None else name) table.to_ASCII(sys.stdout if name is None else name)

7
processing/post/addDivergence.py
View File

@ -49,9 +49,10 @@ for name in filenames:
for label in options.labels: for label in options.labels:
field = table.get(label) field = table.get(label)
shape = (3,) if np.prod(field.shape)//np.prod(grid) == 3 else (3,3) # vector or tensor shape = (3,) if np.prod(field.shape)//np.prod(grid) == 3 else (3,3) # vector or tensor
field = field.reshape(np.append(grid[::-1],shape)) field = field.reshape(tuple(grid)+(-1,),order='F').reshape(tuple(grid)+shape)
div = damask.grid_filters.divergence(size,field)
table.add('divFFT({})'.format(label), table.add('divFFT({})'.format(label),
damask.grid_filters.divergence(size[::-1],field).reshape(-1,np.prod(shape)//3), div.reshape(tuple(grid)+(-1,)).reshape(-1,np.prod(shape)//3,order='F'),
scriptID+' '+' '.join(sys.argv[1:])) scriptID+' '+' '.join(sys.argv[1:]))
table.to_ASCII(sys.stdout if name is None else name) table.to_ASCII(sys.stdout if name is None else name)

7
processing/post/addGradient.py
View File

@ -49,9 +49,10 @@ for name in filenames:
for label in options.labels: for label in options.labels:
field = table.get(label) field = table.get(label)
shape = (1,) if np.prod(field.shape)//np.prod(grid) == 1 else (3,) # scalar or vector shape = (1,) if np.prod(field.shape)//np.prod(grid) == 1 else (3,) # scalar or vector
field = field.reshape(np.append(grid[::-1],shape)) field = field.reshape(tuple(grid)+(-1,),order='F')
grad = damask.grid_filters.gradient(size,field)
table.add('gradFFT({})'.format(label), table.add('gradFFT({})'.format(label),
damask.grid_filters.gradient(size[::-1],field).reshape(-1,np.prod(shape)*3), grad.reshape(tuple(grid)+(-1,)).reshape(-1,np.prod(shape)*3,order='F'),
scriptID+' '+' '.join(sys.argv[1:])) scriptID+' '+' '.join(sys.argv[1:]))
table.to_ASCII(sys.stdout if name is None else name) table.to_ASCII(sys.stdout if name is None else name)

6
processing/post/averageDown.py
View File

@ -61,7 +61,7 @@ if any(shift != 0): prefix += 'shift{:+}{:+}{:+}_'.format(*shift)
for name in filenames: for name in filenames:
damask.util.report(scriptName,name) damask.util.report(scriptName,name)
table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
if (options.grid is None or options.size is None): if (options.grid is None or options.size is None):
@ -87,11 +87,11 @@ for name in filenames:
origin = list(-(packing//2)) + [0])\ origin = list(-(packing//2)) + [0])\
[::packing[0],::packing[1],::packing[2],:].reshape((packedGrid.prod(),-1),order = 'F') [::packing[0],::packing[1],::packing[2],:].reshape((packedGrid.prod(),-1),order = 'F')
table = damask.Table(averagedDown,table.shapes,table.comments) table = damask.Table(averagedDown,table.shapes,table.comments)
coords = damask.grid_filters.cell_coord0(packedGrid,size,shift/packedGrid*size+origin) coords = damask.grid_filters.cell_coord0(packedGrid,size,shift/packedGrid*size+origin)
table.set(options.pos, coords.reshape(-1,3)) table.set(options.pos, coords.reshape(-1,3,order='F'))
outname = os.path.join(os.path.dirname(name),prefix+os.path.basename(name)) outname = os.path.join(os.path.dirname(name),prefix+os.path.basename(name))

6
processing/post/blowUp.py
View File

@ -59,13 +59,13 @@ for name in filenames:
packing = np.array(options.packing,'i') packing = np.array(options.packing,'i')
outSize = grid*packing outSize = grid*packing
data = table.data.values.reshape(tuple(grid)+(-1,)) data = table.data.values.reshape(tuple(grid)+(-1,),order='F')
blownUp = ndimage.interpolation.zoom(data,tuple(packing)+(1,),order=0,mode='nearest').reshape(outSize.prod(),-1) blownUp = ndimage.interpolation.zoom(data,tuple(packing)+(1,),order=0,mode='nearest').reshape(outSize.prod(),-1,order='F')
table = damask.Table(blownUp,table.shapes,table.comments) table = damask.Table(blownUp,table.shapes,table.comments)
coords = damask.grid_filters.cell_coord0(outSize,size,origin) coords = damask.grid_filters.cell_coord0(outSize,size,origin)
table.set(options.pos,coords.reshape(-1,3)) table.set(options.pos,coords.reshape(-1,3,order='F'))
table.set('elem',np.arange(1,outSize.prod()+1)) table.set('elem',np.arange(1,outSize.prod()+1))
outname = os.path.join(os.path.dirname(name),prefix+os.path.basename(name)) outname = os.path.join(os.path.dirname(name),prefix+os.path.basename(name))

14
processing/post/permuteData.py
View File

@ -60,7 +60,7 @@ for name in filenames:
table.head_read() table.head_read()
# ------------------------------------------ process labels --------------------------------------- # ------------------------------------------ process labels ---------------------------------------
errors = [] errors = []
remarks = [] remarks = []
@ -80,7 +80,7 @@ for name in filenames:
damask.util.croak(errors) damask.util.croak(errors)
table.close(dismiss = True) table.close(dismiss = True)
continue continue
# ------------------------------------------ assemble header --------------------------------------- # ------------------------------------------ assemble header ---------------------------------------
randomSeed = int(os.urandom(4).hex(), 16) if options.randomSeed is None else options.randomSeed # random seed per file randomSeed = int(os.urandom(4).hex(), 16) if options.randomSeed is None else options.randomSeed # random seed per file
@ -97,17 +97,17 @@ for name in filenames:
for col,dim in zip(columns,dims): for col,dim in zip(columns,dims):
if options.unique: if options.unique:
s = set(map(tuple,table.data[:,col:col+dim])) # generate set of (unique) values s = set(map(tuple,table.data[:,col:col+dim])) # generate set of (unique) values
uniques = np.array(map(np.array,s)) # translate set to np.array uniques = np.array(list(map(np.array,s))) # translate set to np.array
shuffler = dict(zip(s,np.random.permutation(len(s)))) # random permutation shuffler = dict(zip(s,np.random.permutation(len(s)))) # random permutation
table.data[:,col:col+dim] = uniques[np.array(map(lambda x: shuffler[tuple(x)], table.data[:,col:col+dim] = uniques[np.array(list(map(lambda x: shuffler[tuple(x)],
table.data[:,col:col+dim]))] # fill table with mapped uniques table.data[:,col:col+dim])))] # fill table with mapped uniques
else: else:
np.random.shuffle(table.data[:,col:col+dim]) # independently shuffle every row np.random.shuffle(table.data[:,col:col+dim]) # independently shuffle every row
# ------------------------------------------ output result ----------------------------------------- # ------------------------------------------ output result -----------------------------------------
table.data_writeArray() table.data_writeArray()
# ------------------------------------------ output finalization ----------------------------------- # ------------------------------------------ output finalization -----------------------------------
table.close() # close ASCII tables table.close() # close ASCII tables

12
processing/pre/geom_fromVoronoiTessellation.py
View File

@ -24,22 +24,22 @@ def findClosestSeed(seeds, weights, point):
def Laguerre_tessellation(grid, size, seeds, weights, origin = np.zeros(3), periodic = True, cpus = 2): def Laguerre_tessellation(grid, size, seeds, weights, origin = np.zeros(3), periodic = True, cpus = 2):
if periodic: if periodic:
weights_p = np.tile(weights,27).flatten(order='F') # Laguerre weights (1,2,3,1,2,3,...,1,2,3) weights_p = np.tile(weights.squeeze(),27) # Laguerre weights (1,2,3,1,2,3,...,1,2,3)
seeds_p = np.vstack((seeds -np.array([size[0],0.,0.]),seeds, seeds +np.array([size[0],0.,0.]))) seeds_p = np.vstack((seeds -np.array([size[0],0.,0.]),seeds, seeds +np.array([size[0],0.,0.])))
seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.]))) seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.])))
seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]]))) seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]])))
coords = damask.grid_filters.cell_coord0(grid*3,size*3,-origin-size).reshape(-1,3,order='F') coords = damask.grid_filters.cell_coord0(grid*3,size*3,-origin-size).reshape(-1,3)
else: else:
weights_p = weights.flatten() weights_p = weights.squeeze()
seeds_p = seeds seeds_p = seeds
coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3,order='F') coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3)
if cpus > 1: if cpus > 1:
pool = multiprocessing.Pool(processes = cpus) pool = multiprocessing.Pool(processes = cpus)
result = pool.map_async(partial(findClosestSeed,seeds_p,weights_p), [coord for coord in coords]) result = pool.map_async(partial(findClosestSeed,seeds_p,weights_p), [coord for coord in coords])
pool.close() pool.close()
pool.join() pool.join()
closest_seed = np.array(result.get()) closest_seed = np.array(result.get()).reshape(-1,3)
else: else:
closest_seed= np.array([findClosestSeed(seeds_p,weights_p,coord) for coord in coords]) closest_seed= np.array([findClosestSeed(seeds_p,weights_p,coord) for coord in coords])
@ -52,7 +52,7 @@ def Laguerre_tessellation(grid, size, seeds, weights, origin = np.zeros(3), peri
def Voronoi_tessellation(grid, size, seeds, origin = np.zeros(3), periodic = True): def Voronoi_tessellation(grid, size, seeds, origin = np.zeros(3), periodic = True):
coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3,order='F') coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3)
KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds) KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds)
devNull,closest_seed = KDTree.query(coords) devNull,closest_seed = KDTree.query(coords)

8
processing/pre/seeds_fromGeom.py
View File

@ -45,7 +45,7 @@ options.blacklist = [int(i) for i in options.blacklist]
for name in filenames: for name in filenames:
damask.util.report(scriptName,name) damask.util.report(scriptName,name)
geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name) geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name)
microstructure = geom.get_microstructure().reshape((-1,1),order='F') microstructure = geom.get_microstructure().reshape((-1,1),order='F')
@ -53,9 +53,9 @@ for name in filenames:
np.full(geom.grid.prod(),True,dtype=bool), np.full(geom.grid.prod(),True,dtype=bool),
np.in1d(microstructure,options.blacklist,invert=True) if options.blacklist else \ np.in1d(microstructure,options.blacklist,invert=True) if options.blacklist else \
np.full(geom.grid.prod(),True,dtype=bool)) np.full(geom.grid.prod(),True,dtype=bool))
seeds = damask.grid_filters.cell_coord0(geom.grid,geom.size).reshape(-1,3) seeds = damask.grid_filters.cell_coord0(geom.grid,geom.size).reshape(-1,3,order='F')
comments = geom.comments \ comments = geom.comments \
+ [scriptID + ' ' + ' '.join(sys.argv[1:]), + [scriptID + ' ' + ' '.join(sys.argv[1:]),
'grid\ta {}\tb {}\tc {}'.format(*geom.grid), 'grid\ta {}\tb {}\tc {}'.format(*geom.grid),

2
processing/pre/seeds_fromRandom.py
View File

@ -128,7 +128,7 @@ for name in filenames:
if not options.selective: if not options.selective:
coords = damask.grid_filters.cell_coord0(grid,size).reshape(-1,3) coords = damask.grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F')
seeds = coords[np.random.choice(np.prod(grid), options.N, replace=False)] \ seeds = coords[np.random.choice(np.prod(grid), options.N, replace=False)] \
+ np.broadcast_to(size/grid,(options.N,3))*(np.random.rand(options.N,3)*.5-.25) # wobble without leaving grid + np.broadcast_to(size/grid,(options.N,3))*(np.random.rand(options.N,3)*.5-.25) # wobble without leaving grid
else: else:

16
python/damask/_geom.py
View File

@ -322,11 +322,10 @@ class Geom:
if i != grid.prod(): if i != grid.prod():
raise TypeError('Invalid file: expected {} entries, found {}'.format(grid.prod(),i)) raise TypeError('Invalid file: expected {} entries, found {}'.format(grid.prod(),i))
microstructure = microstructure.reshape(grid,order='F') if not np.any(np.mod(microstructure,1) != 0.0): # no float present
if not np.any(np.mod(microstructure.flatten(),1) != 0.0): # no float present
microstructure = microstructure.astype('int') microstructure = microstructure.astype('int')
return Geom(microstructure.reshape(grid),size,origin,homogenization,comments) return Geom(microstructure.reshape(grid,order='F'),size,origin,homogenization,comments)
@staticmethod @staticmethod
@ -352,16 +351,15 @@ class Geom:
""" """
if periodic: if periodic:
weights_p = np.tile(weights,27).flatten(order='F') # Laguerre weights (1,2,3,1,2,3,...,1,2,3) weights_p = np.tile(weights,27) # Laguerre weights (1,2,3,1,2,3,...,1,2,3)
seeds_p = np.vstack((seeds -np.array([size[0],0.,0.]),seeds, seeds +np.array([size[0],0.,0.]))) seeds_p = np.vstack((seeds -np.array([size[0],0.,0.]),seeds, seeds +np.array([size[0],0.,0.])))
seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.]))) seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.])))
seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]]))) seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]])))
coords = grid_filters.cell_coord0(grid*3,size*3,-size).reshape(-1,3,order='F') coords = grid_filters.cell_coord0(grid*3,size*3,-size).reshape(-1,3)
else: else:
weights_p = weights.flatten() weights_p = weights
seeds_p = seeds seeds_p = seeds
coords = grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F') coords = grid_filters.cell_coord0(grid,size).reshape(-1,3)
pool = multiprocessing.Pool(processes = int(Environment().options['DAMASK_NUM_THREADS'])) pool = multiprocessing.Pool(processes = int(Environment().options['DAMASK_NUM_THREADS']))
result = pool.map_async(partial(Geom._find_closest_seed,seeds_p,weights_p), [coord for coord in coords]) result = pool.map_async(partial(Geom._find_closest_seed,seeds_p,weights_p), [coord for coord in coords])
@ -396,7 +394,7 @@ class Geom:
perform a periodic tessellation. Defaults to True. perform a periodic tessellation. Defaults to True.
""" """
coords = grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F') coords = grid_filters.cell_coord0(grid,size).reshape(-1,3)
KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds) KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds)
devNull,microstructure = KDTree.query(coords) devNull,microstructure = KDTree.query(coords)

16
python/damask/_result.py
View File

@ -111,7 +111,7 @@ class Result:
select from 'set', 'add', and 'del' select from 'set', 'add', and 'del'
what : str what : str
attribute to change (must be from self.selection) attribute to change (must be from self.selection)
datasets : list of str or Boolean datasets : list of str or bool
name of datasets as list, supports ? and * wildcards. name of datasets as list, supports ? and * wildcards.
True is equivalent to [*], False is equivalent to [] True is equivalent to [*], False is equivalent to []
@ -203,7 +203,7 @@ class Result:
---------- ----------
what : str what : str
attribute to change (must be from self.selection) attribute to change (must be from self.selection)
datasets : list of str or Boolean datasets : list of str or bool
name of datasets as list, supports ? and * wildcards. name of datasets as list, supports ? and * wildcards.
True is equivalent to [*], False is equivalent to [] True is equivalent to [*], False is equivalent to []
@ -219,7 +219,7 @@ class Result:
---------- ----------
what : str what : str
attribute to change (must be from self.selection) attribute to change (must be from self.selection)
datasets : list of str or Boolean datasets : list of str or bool
name of datasets as list, supports ? and * wildcards. name of datasets as list, supports ? and * wildcards.
True is equivalent to [*], False is equivalent to [] True is equivalent to [*], False is equivalent to []
@ -235,7 +235,7 @@ class Result:
---------- ----------
what : str what : str
attribute to change (must be from self.selection) attribute to change (must be from self.selection)
datasets : list of str or Boolean datasets : list of str or bool
name of datasets as list, supports ? and * wildcards. name of datasets as list, supports ? and * wildcards.
True is equivalent to [*], False is equivalent to [] True is equivalent to [*], False is equivalent to []
@ -262,10 +262,10 @@ class Result:
datasets : iterable or str datasets : iterable or str
component : int component : int
homogenization component to consider for constituent data homogenization component to consider for constituent data
tagged : Boolean tagged : bool
tag Table.column name with '#component' tag Table.column name with '#component'
defaults to False defaults to False
split : Boolean split : bool
split Table by increment and return dictionary of Tables split Table by increment and return dictionary of Tables
defaults to True defaults to True
@ -326,7 +326,7 @@ class Result:
Parameters Parameters
---------- ----------
datasets : iterable or str or Boolean datasets : iterable or str or bool
Examples Examples
-------- --------
@ -460,7 +460,7 @@ class Result:
def cell_coordinates(self): def cell_coordinates(self):
"""Return initial coordinates of the cell centers.""" """Return initial coordinates of the cell centers."""
if self.structured: if self.structured:
return grid_filters.cell_coord0(self.grid,self.size,self.origin).reshape(-1,3) return grid_filters.cell_coord0(self.grid,self.size,self.origin).reshape(-1,3,order='F')
else: else:
with h5py.File(self.fname,'r') as f: with h5py.File(self.fname,'r') as f:
return f['geometry/x_c'][()] return f['geometry/x_c'][()]

130
python/damask/grid_filters.py
View File

@ -1,3 +1,17 @@
"""
Filters for operations on regular grids.
Notes
-----
The grids are defined as (x,y,z,...) where x is fastest and z is slowest.
This convention is consistent with the geom file format.
When converting to/from a plain list (e.g. storage in ASCII table),
the following operations are required for tensorial data:
D3 = D1.reshape(grid+(-1,),order='F').reshape(grid+(3,3))
D1 = D3.reshape(grid+(-1,)).reshape(-1,9,order='F')
"""
from scipy import spatial as _spatial from scipy import spatial as _spatial
import numpy as _np import numpy as _np
@ -7,8 +21,12 @@ def _ks(size,grid,first_order=False):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray of shape (3)
physical size of the periodic field. physical size of the periodic field.
grid : numpy.ndarray of shape (3)
number of grid points.
first_order : bool, optional
correction for first order derivatives, defaults to False.
""" """
k_sk = _np.where(_np.arange(grid[0])>grid[0]//2,_np.arange(grid[0])-grid[0],_np.arange(grid[0]))/size[0] k_sk = _np.where(_np.arange(grid[0])>grid[0]//2,_np.arange(grid[0])-grid[0],_np.arange(grid[0]))/size[0]
@ -19,8 +37,7 @@ def _ks(size,grid,first_order=False):
k_si = _np.arange(grid[2]//2+1)/size[2] k_si = _np.arange(grid[2]//2+1)/size[2]
kk, kj, ki = _np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij') return _np.stack(_np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij'), axis=-1)
return _np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3)
def curl(size,field): def curl(size,field):
@ -29,8 +46,10 @@ def curl(size,field):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray of shape (3)
physical size of the periodic field. physical size of the periodic field.
field : numpy.ndarray of shape (:,:,:,3) or (:,:,:,3,3)
periodic field of which the curl is calculated.
""" """
n = _np.prod(field.shape[3:]) n = _np.prod(field.shape[3:])
@ -53,8 +72,10 @@ def divergence(size,field):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray of shape (3)
physical size of the periodic field. physical size of the periodic field.
field : numpy.ndarray of shape (:,:,:,3) or (:,:,:,3,3)
periodic field of which the divergence is calculated.
""" """
n = _np.prod(field.shape[3:]) n = _np.prod(field.shape[3:])
@ -69,12 +90,14 @@ def divergence(size,field):
def gradient(size,field): def gradient(size,field):
""" """
Calculate gradient of a vector or scalar field in Fourier space. Calculate gradient of a scalar or vector field in Fourier space.
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray of shape (3)
physical size of the periodic field. physical size of the periodic field.
field : numpy.ndarray of shape (:,:,:,1) or (:,:,:,3)
periodic field of which the gradient is calculated.
""" """
n = _np.prod(field.shape[3:]) n = _np.prod(field.shape[3:])
@ -93,9 +116,9 @@ def cell_coord0(grid,size,origin=_np.zeros(3)):
Parameters Parameters
---------- ----------
grid : numpy.ndarray grid : numpy.ndarray of shape (3)
number of grid points. number of grid points.
size : numpy.ndarray size : numpy.ndarray of shape (3)
physical size of the periodic field. physical size of the periodic field.
origin : numpy.ndarray, optional origin : numpy.ndarray, optional
physical origin of the periodic field. Defaults to [0.0,0.0,0.0]. physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
@ -103,7 +126,11 @@ def cell_coord0(grid,size,origin=_np.zeros(3)):
""" """
start = origin + size/grid*.5 start = origin + size/grid*.5
end = origin + size - size/grid*.5 end = origin + size - size/grid*.5
return _np.mgrid[start[0]:end[0]:grid[0]*1j,start[1]:end[1]:grid[1]*1j,start[2]:end[2]:grid[2]*1j].T
return _np.stack(_np.meshgrid(_np.linspace(start[0],end[0],grid[0]),
_np.linspace(start[1],end[1],grid[1]),
_np.linspace(start[2],end[2],grid[2]),indexing = 'ij'),
axis = -1)
def cell_displacement_fluct(size,F): def cell_displacement_fluct(size,F):
@ -112,7 +139,7 @@ def cell_displacement_fluct(size,F):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray of shape (3)
physical size of the periodic field. physical size of the periodic field.
F : numpy.ndarray F : numpy.ndarray
deformation gradient field. deformation gradient field.
@ -139,14 +166,14 @@ def cell_displacement_avg(size,F):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray of shape (3)
physical size of the periodic field. physical size of the periodic field.
F : numpy.ndarray F : numpy.ndarray
deformation gradient field. deformation gradient field.
""" """
F_avg = _np.average(F,axis=(0,1,2)) F_avg = _np.average(F,axis=(0,1,2))
return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),cell_coord0(F.shape[:3][::-1],size)) return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),cell_coord0(F.shape[:3],size))
def cell_displacement(size,F): def cell_displacement(size,F):
@ -155,7 +182,7 @@ def cell_displacement(size,F):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray of shape (3)
physical size of the periodic field. physical size of the periodic field.
F : numpy.ndarray F : numpy.ndarray
deformation gradient field. deformation gradient field.
@ -170,25 +197,25 @@ def cell_coord(size,F,origin=_np.zeros(3)):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray of shape (3)
physical size of the periodic field. physical size of the periodic field.
F : numpy.ndarray F : numpy.ndarray
deformation gradient field. deformation gradient field.
origin : numpy.ndarray, optional origin : numpy.ndarray of shape (3), optional
physical origin of the periodic field. Defaults to [0.0,0.0,0.0]. physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
""" """
return cell_coord0(F.shape[:3][::-1],size,origin) + cell_displacement(size,F) return cell_coord0(F.shape[:3],size,origin) + cell_displacement(size,F)
def cell_coord0_gridSizeOrigin(coord0,ordered=True): def cell_coord0_gridSizeOrigin(coord0,ordered=True):
""" """
Return grid 'DNA', i.e. grid, size, and origin from array of cell positions. Return grid 'DNA', i.e. grid, size, and origin from 1D array of cell positions.
Parameters Parameters
---------- ----------
coord0 : numpy.ndarray coord0 : numpy.ndarray of shape (:,3)
array of undeformed cell coordinates. undeformed cell coordinates.
ordered : bool, optional ordered : bool, optional
expect coord0 data to be ordered (x fast, z slow). expect coord0 data to be ordered (x fast, z slow).
@ -211,13 +238,13 @@ def cell_coord0_gridSizeOrigin(coord0,ordered=True):
start = origin + delta*.5 start = origin + delta*.5
end = origin - delta*.5 + size end = origin - delta*.5 + size
if not _np.allclose(coords[0],_np.linspace(start[0],end[0],grid[0])) and \ if not (_np.allclose(coords[0],_np.linspace(start[0],end[0],grid[0])) and \
_np.allclose(coords[1],_np.linspace(start[1],end[1],grid[1])) and \ _np.allclose(coords[1],_np.linspace(start[1],end[1],grid[1])) and \
_np.allclose(coords[2],_np.linspace(start[2],end[2],grid[2])): _np.allclose(coords[2],_np.linspace(start[2],end[2],grid[2]))):
raise ValueError('Regular grid spacing violated.') raise ValueError('Regular grid spacing violated.')
if ordered and not _np.allclose(coord0.reshape(tuple(grid[::-1])+(3,)),cell_coord0(grid,size,origin)): if ordered and not _np.allclose(coord0.reshape(tuple(grid)+(3,),order='F'),cell_coord0(grid,size,origin)):
raise ValueError('Input data is not a regular grid.') raise ValueError('Input data is not ordered (x fast, z slow).')
return (grid,size,origin) return (grid,size,origin)
@ -241,17 +268,18 @@ def node_coord0(grid,size,origin=_np.zeros(3)):
Parameters Parameters
---------- ----------
grid : numpy.ndarray grid : numpy.ndarray of shape (3)
number of grid points. number of grid points.
size : numpy.ndarray size : numpy.ndarray of shape (3)
physical size of the periodic field. physical size of the periodic field.
origin : numpy.ndarray, optional origin : numpy.ndarray of shape (3), optional
physical origin of the periodic field. Defaults to [0.0,0.0,0.0]. physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
""" """
return _np.mgrid[origin[0]:size[0]+origin[0]:(grid[0]+1)*1j, return _np.stack(_np.meshgrid(_np.linspace(origin[0],size[0]+origin[0],grid[0]+1),
origin[1]:size[1]+origin[1]:(grid[1]+1)*1j, _np.linspace(origin[1],size[1]+origin[1],grid[1]+1),
origin[2]:size[2]+origin[2]:(grid[2]+1)*1j].T _np.linspace(origin[2],size[2]+origin[2],grid[2]+1),indexing = 'ij'),
axis = -1)
def node_displacement_fluct(size,F): def node_displacement_fluct(size,F):
@ -260,7 +288,7 @@ def node_displacement_fluct(size,F):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray of shape (3)
physical size of the periodic field. physical size of the periodic field.
F : numpy.ndarray F : numpy.ndarray
deformation gradient field. deformation gradient field.
@ -275,14 +303,14 @@ def node_displacement_avg(size,F):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray of shape (3)
physical size of the periodic field. physical size of the periodic field.
F : numpy.ndarray F : numpy.ndarray
deformation gradient field. deformation gradient field.
""" """
F_avg = _np.average(F,axis=(0,1,2)) F_avg = _np.average(F,axis=(0,1,2))
return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),node_coord0(F.shape[:3][::-1],size)) return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),node_coord0(F.shape[:3],size))
def node_displacement(size,F): def node_displacement(size,F):
@ -291,7 +319,7 @@ def node_displacement(size,F):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray of shape (3)
physical size of the periodic field. physical size of the periodic field.
F : numpy.ndarray F : numpy.ndarray
deformation gradient field. deformation gradient field.
@ -306,15 +334,15 @@ def node_coord(size,F,origin=_np.zeros(3)):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray of shape (3)
physical size of the periodic field. physical size of the periodic field.
F : numpy.ndarray F : numpy.ndarray
deformation gradient field. deformation gradient field.
origin : numpy.ndarray, optional origin : numpy.ndarray of shape (3), optional
physical origin of the periodic field. Defaults to [0.0,0.0,0.0]. physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
""" """
return node_coord0(F.shape[:3][::-1],size,origin) + node_displacement(size,F) return node_coord0(F.shape[:3],size,origin) + node_displacement(size,F)
def cell_2_node(cell_data): def cell_2_node(cell_data):
@ -335,14 +363,14 @@ def node_2_cell(node_data):
return c[:-1,:-1,:-1] return c[:-1,:-1,:-1]
def node_coord0_gridSizeOrigin(coord0,ordered=False): def node_coord0_gridSizeOrigin(coord0,ordered=True):
""" """
Return grid 'DNA', i.e. grid, size, and origin from array of nodal positions. Return grid 'DNA', i.e. grid, size, and origin from 1D array of nodal positions.
Parameters Parameters
---------- ----------
coord0 : numpy.ndarray coord0 : numpy.ndarray of shape (:,3)
array of undeformed nodal coordinates. undeformed nodal coordinates.
ordered : bool, optional ordered : bool, optional
expect coord0 data to be ordered (x fast, z slow). expect coord0 data to be ordered (x fast, z slow).
@ -357,13 +385,13 @@ def node_coord0_gridSizeOrigin(coord0,ordered=False):
if (grid+1).prod() != len(coord0): if (grid+1).prod() != len(coord0):
raise ValueError('Data count {} does not match grid {}.'.format(len(coord0),grid)) raise ValueError('Data count {} does not match grid {}.'.format(len(coord0),grid))
if not _np.allclose(coords[0],_np.linspace(mincorner[0],maxcorner[0],grid[0]+1)) and \ if not (_np.allclose(coords[0],_np.linspace(mincorner[0],maxcorner[0],grid[0]+1)) and \
_np.allclose(coords[1],_np.linspace(mincorner[1],maxcorner[1],grid[1]+1)) and \ _np.allclose(coords[1],_np.linspace(mincorner[1],maxcorner[1],grid[1]+1)) and \
_np.allclose(coords[2],_np.linspace(mincorner[2],maxcorner[2],grid[2]+1)): _np.allclose(coords[2],_np.linspace(mincorner[2],maxcorner[2],grid[2]+1))):
raise ValueError('Regular grid spacing violated.') raise ValueError('Regular grid spacing violated.')
if ordered and not _np.allclose(coord0.reshape(tuple((grid+1)[::-1])+(3,)),node_coord0(grid,size,origin)): if ordered and not _np.allclose(coord0.reshape(tuple(grid+1)+(3,),order='F'),node_coord0(grid,size,origin)):
raise ValueError('Input data is not a regular grid.') raise ValueError('Input data is not ordered (x fast, z slow).')
return (grid,size,origin) return (grid,size,origin)
@ -374,15 +402,15 @@ def regrid(size,F,new_grid):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray of shape (3)
physical size physical size
F : numpy.ndarray F : numpy.ndarray of shape (:,:,:,3,3)
deformation gradient field deformation gradient field
new_grid : numpy.ndarray new_grid : numpy.ndarray of shape (3)
new grid for undeformed coordinates new grid for undeformed coordinates
""" """
c = cell_coord0(F.shape[:3][::-1],size) \ c = cell_coord0(F.shape[:3],size) \
+ cell_displacement_avg(size,F) \ + cell_displacement_avg(size,F) \
+ cell_displacement_fluct(size,F) + cell_displacement_fluct(size,F)

49
python/tests/test_grid_filters.py
View File

@ -4,18 +4,18 @@ import numpy as np
from damask import grid_filters from damask import grid_filters
class TestGridFilters: class TestGridFilters:
def test_cell_coord0(self): def test_cell_coord0(self):
size = np.random.random(3) size = np.random.random(3)
grid = np.random.randint(8,32,(3)) grid = np.random.randint(8,32,(3))
coord = grid_filters.cell_coord0(grid,size) coord = grid_filters.cell_coord0(grid,size)
assert np.allclose(coord[0,0,0],size/grid*.5) and coord.shape == tuple(grid[::-1]) + (3,) assert np.allclose(coord[0,0,0],size/grid*.5) and coord.shape == tuple(grid) + (3,)
def test_node_coord0(self): def test_node_coord0(self):
size = np.random.random(3) size = np.random.random(3)
grid = np.random.randint(8,32,(3)) grid = np.random.randint(8,32,(3))
coord = grid_filters.node_coord0(grid,size) coord = grid_filters.node_coord0(grid,size)
assert np.allclose(coord[-1,-1,-1],size) and coord.shape == tuple(grid[::-1]+1) + (3,) assert np.allclose(coord[-1,-1,-1],size) and coord.shape == tuple(grid+1) + (3,)
def test_coord0(self): def test_coord0(self):
size = np.random.random(3) size = np.random.random(3)
@ -31,7 +31,7 @@ class TestGridFilters:
size = np.random.random(3) size = np.random.random(3)
origin = np.random.random(3) origin = np.random.random(3)
coord0 = eval('grid_filters.{}_coord0(grid,size,origin)'.format(mode)) # noqa coord0 = eval('grid_filters.{}_coord0(grid,size,origin)'.format(mode)) # noqa
_grid,_size,_origin = eval('grid_filters.{}_coord0_gridSizeOrigin(coord0.reshape(-1,3))'.format(mode)) _grid,_size,_origin = eval('grid_filters.{}_coord0_gridSizeOrigin(coord0.reshape(-1,3,order="F"))'.format(mode))
assert np.allclose(grid,_grid) and np.allclose(size,_size) and np.allclose(origin,_origin) assert np.allclose(grid,_grid) and np.allclose(size,_size) and np.allclose(origin,_origin)
def test_displacement_fluct_equivalence(self): def test_displacement_fluct_equivalence(self):
@ -57,9 +57,9 @@ class TestGridFilters:
shifted = eval('grid_filters.{}_coord0(grid,size,origin)'.format(mode)) shifted = eval('grid_filters.{}_coord0(grid,size,origin)'.format(mode))
unshifted = eval('grid_filters.{}_coord0(grid,size)'.format(mode)) unshifted = eval('grid_filters.{}_coord0(grid,size)'.format(mode))
if mode == 'cell': if mode == 'cell':
assert np.allclose(shifted,unshifted+np.broadcast_to(origin,tuple(grid[::-1]) +(3,))) assert np.allclose(shifted,unshifted+np.broadcast_to(origin,tuple(grid) +(3,)))
elif mode == 'node': elif mode == 'node':
assert np.allclose(shifted,unshifted+np.broadcast_to(origin,tuple(grid[::-1]+1)+(3,))) assert np.allclose(shifted,unshifted+np.broadcast_to(origin,tuple(grid+1)+(3,)))
@pytest.mark.parametrize('function',[grid_filters.cell_displacement_avg, @pytest.mark.parametrize('function',[grid_filters.cell_displacement_avg,
grid_filters.node_displacement_avg]) grid_filters.node_displacement_avg])
@ -80,8 +80,43 @@ class TestGridFilters:
F = np.broadcast_to(np.random.random((3,3)), tuple(grid)+(3,3)) F = np.broadcast_to(np.random.random((3,3)), tuple(grid)+(3,3))
assert np.allclose(function(size,F),0.0) assert np.allclose(function(size,F),0.0)
@pytest.mark.parametrize('function',[grid_filters.coord0_check,
grid_filters.node_coord0_gridSizeOrigin,
grid_filters.cell_coord0_gridSizeOrigin])
def test_invalid_coordinates(self,function):
invalid_coordinates = np.random.random((np.random.randint(12,52),3))
with pytest.raises(ValueError):
function(invalid_coordinates)
@pytest.mark.parametrize('function',[grid_filters.node_coord0_gridSizeOrigin,
grid_filters.cell_coord0_gridSizeOrigin])
def test_uneven_spaced_coordinates(self,function):
start = np.random.random(3)
end = np.random.random(3)*10. + start
grid = np.random.randint(8,32,(3))
uneven = np.stack(np.meshgrid(np.logspace(start[0],end[0],grid[0]),
np.logspace(start[1],end[1],grid[1]),
np.logspace(start[2],end[2],grid[2]),indexing = 'ij'),
axis = -1).reshape((grid.prod(),3),order='F')
with pytest.raises(ValueError):
function(uneven)
@pytest.mark.parametrize('mode',[True,False])
@pytest.mark.parametrize('function',[grid_filters.node_coord0_gridSizeOrigin,
grid_filters.cell_coord0_gridSizeOrigin])
def test_unordered_coordinates(self,function,mode):
origin = np.random.random(3)
size = np.random.random(3)*10.+origin
grid = np.random.randint(8,32,(3))
unordered = grid_filters.node_coord0(grid,size,origin).reshape(-1,3)
if mode:
with pytest.raises(ValueError):
function(unordered,mode)
else:
function(unordered,mode)
def test_regrid(self): def test_regrid(self):
size = np.random.random(3) size = np.random.random(3)
grid = np.random.randint(8,32,(3)) grid = np.random.randint(8,32,(3))
F = np.broadcast_to(np.eye(3), tuple(grid[::-1])+(3,3)) F = np.broadcast_to(np.eye(3), tuple(grid)+(3,3))
assert all(grid_filters.regrid(size,F,grid) == np.arange(grid.prod())) assert all(grid_filters.regrid(size,F,grid) == np.arange(grid.prod()))

74
src/constitutive.f90
View File

@ -327,7 +327,7 @@ module constitutive
constitutive_initialFi, & constitutive_initialFi, &
constitutive_SandItsTangents, & constitutive_SandItsTangents, &
constitutive_collectDotState, & constitutive_collectDotState, &
constitutive_collectDeltaState, & constitutive_deltaState, &
constitutive_results constitutive_results
contains contains
@ -709,12 +709,14 @@ end subroutine constitutive_hooke_SandItsTangents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure !> @brief contains the constitutive equation for calculating the rate of change of microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el) function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken)
integer, intent(in) :: & integer, intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el, & !< element
phase, &
of
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
subdt !< timestep subdt !< timestep
real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
@ -730,16 +732,16 @@ subroutine constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip,
ho, & !< homogenization ho, & !< homogenization
tme, & !< thermal member position tme, & !< thermal member position
i, & !< counter in source loop i, & !< counter in source loop
instance, of instance
logical :: broken
ho = material_homogenizationAt(el) ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el) tme = thermalMapping(ho)%p(ip,el)
of = material_phasememberAt(ipc,ip,el) instance = phase_plasticityInstance(phase)
instance = phase_plasticityInstance(material_phaseAt(ipc,el))
Mp = matmul(matmul(transpose(Fi),Fi),S) Mp = matmul(matmul(transpose(Fi),Fi),S)
plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) plasticityType: select case (phase_plasticity(phase))
case (PLASTICITY_ISOTROPIC_ID) plasticityType case (PLASTICITY_ISOTROPIC_ID) plasticityType
call plastic_isotropic_dotState (Mp,instance,of) call plastic_isotropic_dotState (Mp,instance,of)
@ -760,10 +762,11 @@ subroutine constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip,
call plastic_nonlocal_dotState (Mp,FArray,FpArray,temperature(ho)%p(tme),subdt, & call plastic_nonlocal_dotState (Mp,FArray,FpArray,temperature(ho)%p(tme),subdt, &
instance,of,ip,el) instance,of,ip,el)
end select plasticityType end select plasticityType
broken = any(IEEE_is_NaN(plasticState(phase)%dotState(:,of)))
SourceLoop: do i = 1, phase_Nsources(material_phaseAt(ipc,el)) SourceLoop: do i = 1, phase_Nsources(phase)
sourceType: select case (phase_source(i,material_phaseAt(ipc,el))) sourceType: select case (phase_source(i,phase))
case (SOURCE_damage_anisoBrittle_ID) sourceType case (SOURCE_damage_anisoBrittle_ID) sourceType
call source_damage_anisoBrittle_dotState (S, ipc, ip, el) !< correct stress? call source_damage_anisoBrittle_dotState (S, ipc, ip, el) !< correct stress?
@ -775,25 +778,29 @@ subroutine constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip,
call source_damage_anisoDuctile_dotState ( ipc, ip, el) call source_damage_anisoDuctile_dotState ( ipc, ip, el)
case (SOURCE_thermal_externalheat_ID) sourceType case (SOURCE_thermal_externalheat_ID) sourceType
call source_thermal_externalheat_dotState(material_phaseAt(ipc,el),of) call source_thermal_externalheat_dotState(phase,of)
end select sourceType end select sourceType
broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%dotState(:,of)))
enddo SourceLoop enddo SourceLoop
end subroutine constitutive_collectDotState end function constitutive_collectDotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief for constitutive models having an instantaneous change of state !> @brief for constitutive models having an instantaneous change of state
!> will return false if delta state is not needed/supported by the constitutive model !> will return false if delta state is not needed/supported by the constitutive model
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el) function constitutive_deltaState(S, Fe, Fi, ipc, ip, el, phase, of) result(broken)
integer, intent(in) :: & integer, intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el, & !< element
phase, &
of
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
S, & !< 2nd Piola Kirchhoff stress S, & !< 2nd Piola Kirchhoff stress
Fe, & !< elastic deformation gradient Fe, & !< elastic deformation gradient
@ -802,35 +809,62 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
Mp Mp
integer :: & integer :: &
i, & i, &
instance, of instance, &
myOffset, &
mySize
logical :: &
broken
Mp = matmul(matmul(transpose(Fi),Fi),S) Mp = matmul(matmul(transpose(Fi),Fi),S)
of = material_phasememberAt(ipc,ip,el) instance = phase_plasticityInstance(phase)
instance = phase_plasticityInstance(material_phaseAt(ipc,el))
plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) plasticityType: select case (phase_plasticity(phase))
case (PLASTICITY_KINEHARDENING_ID) plasticityType case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_deltaState(Mp,instance,of) call plastic_kinehardening_deltaState(Mp,instance,of)
broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of)))
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_deltaState(Mp,instance,of,ip,el) call plastic_nonlocal_deltaState(Mp,instance,of,ip,el)
broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of)))
case default
broken = .false.
end select plasticityType end select plasticityType
sourceLoop: do i = 1, phase_Nsources(material_phaseAt(ipc,el)) if(.not. broken) then
select case(phase_plasticity(phase))
case (PLASTICITY_NONLOCAL_ID,PLASTICITY_KINEHARDENING_ID)
sourceType: select case (phase_source(i,material_phaseAt(ipc,el))) myOffset = plasticState(phase)%offsetDeltaState
mySize = plasticState(phase)%sizeDeltaState
plasticState(phase)%state(myOffset + 1:myOffset + mySize,of) = &
plasticState(phase)%state(myOffset + 1:myOffset + mySize,of) + plasticState(phase)%deltaState(1:mySize,of)
end select
endif
sourceLoop: do i = 1, phase_Nsources(phase)
sourceType: select case (phase_source(i,phase))
case (SOURCE_damage_isoBrittle_ID) sourceType case (SOURCE_damage_isoBrittle_ID) sourceType
call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(ipc,ip,el), Fe, & call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(ipc,ip,el), Fe, &
ipc, ip, el) ipc, ip, el)
broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%deltaState(:,of)))
if(.not. broken) then
myOffset = sourceState(phase)%p(i)%offsetDeltaState
mySize = sourceState(phase)%p(i)%sizeDeltaState
sourceState(phase)%p(i)%state(myOffset + 1: myOffset + mySize,of) = &
sourceState(phase)%p(i)%state(myOffset + 1: myOffset + mySize,of) + sourceState(phase)%p(i)%deltaState(1:mySize,of)
endif
end select sourceType end select sourceType
enddo SourceLoop enddo SourceLoop
end subroutine constitutive_collectDeltaState end function constitutive_deltaState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

2
src/constitutive_plastic_disloUCLA.f90
View File

@ -209,7 +209,7 @@ module subroutine plastic_disloUCLA_init
sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl
sizeState = sizeDotState sizeState = sizeDotState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization

2
src/constitutive_plastic_dislotwin.f90
View File

@ -399,7 +399,7 @@ module subroutine plastic_dislotwin_init
+ size(['f_tr']) * prm%sum_N_tr + size(['f_tr']) * prm%sum_N_tr
sizeState = sizeDotState sizeState = sizeDotState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and atol ! locally defined state aliases and initialization of state0 and atol

2
src/constitutive_plastic_isotropic.f90
View File

@ -117,7 +117,7 @@ module subroutine plastic_isotropic_init
sizeDotState = size(['xi ','accumulated_shear']) sizeDotState = size(['xi ','accumulated_shear'])
sizeState = sizeDotState sizeState = sizeDotState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization

2
src/constitutive_plastic_kinehardening.f90
View File

@ -164,7 +164,7 @@ module subroutine plastic_kinehardening_init
sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl
sizeState = sizeDotState + sizeDeltaState sizeState = sizeDotState + sizeDeltaState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,sizeDeltaState)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization

2
src/constitutive_plastic_none.f90
View File

@ -29,7 +29,7 @@ module subroutine plastic_none_init
if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
call material_allocatePlasticState(p,NipcMyPhase,0,0,0) call material_allocateState(plasticState(p),NipcMyPhase,0,0,0)
enddo enddo

5
src/constitutive_plastic_nonlocal.f90
View File

@ -320,6 +320,7 @@ module subroutine plastic_nonlocal_init
prm%fEdgeMultiplication = config%getFloat('edgemultiplication') prm%fEdgeMultiplication = config%getFloat('edgemultiplication')
prm%shortRangeStressCorrection = config%keyExists('/shortrangestresscorrection/') prm%shortRangeStressCorrection = config%keyExists('/shortrangestresscorrection/')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! sanity checks ! sanity checks
if (any(prm%burgers < 0.0_pReal)) extmsg = trim(extmsg)//' burgers' if (any(prm%burgers < 0.0_pReal)) extmsg = trim(extmsg)//' burgers'
@ -384,9 +385,9 @@ module subroutine plastic_nonlocal_init
'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure 'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure
sizeDeltaState = sizeDotState sizeDeltaState = sizeDotState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,sizeDeltaState)
plasticState(p)%nonlocal = .true. plasticState(p)%nonlocal = config%KeyExists('/nonlocal/')
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,:) st0%rho => plasticState(p)%state0 (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)

2
src/constitutive_plastic_phenopowerlaw.f90
View File

@ -213,7 +213,7 @@ module subroutine plastic_phenopowerlaw_init
+ size(['xi_tw ','gamma_tw']) * prm%sum_N_tw + size(['xi_tw ','gamma_tw']) * prm%sum_N_tw
sizeState = sizeDotState sizeState = sizeDotState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization

544
src/crystallite.f90
View File

@ -15,7 +15,6 @@ module crystallite
use DAMASK_interface use DAMASK_interface
use config use config
use debug use debug
use numerics
use rotations use rotations
use math use math
use FEsolving use FEsolving
@ -70,9 +69,7 @@ module crystallite
logical, dimension(:,:,:), allocatable, public :: & logical, dimension(:,:,:), allocatable, public :: &
crystallite_requested !< used by upper level (homogenization) to request crystallite calculation crystallite_requested !< used by upper level (homogenization) to request crystallite calculation
logical, dimension(:,:,:), allocatable :: & logical, dimension(:,:,:), allocatable :: &
crystallite_converged, & !< convergence flag crystallite_converged !< convergence flag
crystallite_todo, & !< flag to indicate need for further computation
crystallite_localPlasticity !< indicates this grain to have purely local constitutive law
type :: tOutput !< new requested output (per phase) type :: tOutput !< new requested output (per phase)
character(len=pStringLen), allocatable, dimension(:) :: & character(len=pStringLen), allocatable, dimension(:) :: &
@ -84,7 +81,8 @@ module crystallite
integer :: & integer :: &
iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp
nState, & !< state loop limit nState, & !< state loop limit
nStress !< stress loop limit nStress, & !< stress loop limit
integrator !< integration scheme (ToDo: better use a string)
real(pReal) :: & real(pReal) :: &
subStepMinCryst, & !< minimum (relative) size of sub-step allowed during cutback subStepMinCryst, & !< minimum (relative) size of sub-step allowed during cutback
subStepSizeCryst, & !< size of first substep when cutback subStepSizeCryst, & !< size of first substep when cutback
@ -98,7 +96,7 @@ module crystallite
type(tNumerics) :: num ! numerics parameters. Better name? type(tNumerics) :: num ! numerics parameters. Better name?
procedure(), pointer :: integrateState procedure(integrateStateFPI), pointer :: integrateState
public :: & public :: &
crystallite_init, & crystallite_init, &
@ -159,9 +157,7 @@ subroutine crystallite_init
allocate(crystallite_orientation(cMax,iMax,eMax)) allocate(crystallite_orientation(cMax,iMax,eMax))
allocate(crystallite_localPlasticity(cMax,iMax,eMax), source=.true.)
allocate(crystallite_requested(cMax,iMax,eMax), source=.false.) allocate(crystallite_requested(cMax,iMax,eMax), source=.false.)
allocate(crystallite_todo(cMax,iMax,eMax), source=.false.)
allocate(crystallite_converged(cMax,iMax,eMax), source=.true.) allocate(crystallite_converged(cMax,iMax,eMax), source=.true.)
num%subStepMinCryst = config_numerics%getFloat('substepmincryst', defaultVal=1.0e-3_pReal) num%subStepMinCryst = config_numerics%getFloat('substepmincryst', defaultVal=1.0e-3_pReal)
@ -177,6 +173,8 @@ subroutine crystallite_init
num%iJacoLpresiduum = config_numerics%getInt ('ijacolpresiduum', defaultVal=1) num%iJacoLpresiduum = config_numerics%getInt ('ijacolpresiduum', defaultVal=1)
num%integrator = config_numerics%getInt ('integrator', defaultVal=1)
num%nState = config_numerics%getInt ('nstate', defaultVal=20) num%nState = config_numerics%getInt ('nstate', defaultVal=20)
num%nStress = config_numerics%getInt ('nstress', defaultVal=40) num%nStress = config_numerics%getInt ('nstress', defaultVal=40)
@ -193,10 +191,14 @@ subroutine crystallite_init
if(num%iJacoLpresiduum < 1) call IO_error(301,ext_msg='iJacoLpresiduum') if(num%iJacoLpresiduum < 1) call IO_error(301,ext_msg='iJacoLpresiduum')
if(num%integrator < 1 .or. num%integrator > 5) &
call IO_error(301,ext_msg='integrator')
if(num%nState < 1) call IO_error(301,ext_msg='nState') if(num%nState < 1) call IO_error(301,ext_msg='nState')
if(num%nStress< 1) call IO_error(301,ext_msg='nStress') if(num%nStress< 1) call IO_error(301,ext_msg='nStress')
select case(numerics_integrator)
select case(num%integrator)
case(1) case(1)
integrateState => integrateStateFPI integrateState => integrateStateFPI
case(2) case(2)
@ -234,7 +236,6 @@ subroutine crystallite_init
/ math_det33(crystallite_Fp0(1:3,1:3,c,i,e))**(1.0_pReal/3.0_pReal) / math_det33(crystallite_Fp0(1:3,1:3,c,i,e))**(1.0_pReal/3.0_pReal)
crystallite_Fi0(1:3,1:3,c,i,e) = constitutive_initialFi(c,i,e) crystallite_Fi0(1:3,1:3,c,i,e) = constitutive_initialFi(c,i,e)
crystallite_F0(1:3,1:3,c,i,e) = math_I3 crystallite_F0(1:3,1:3,c,i,e) = math_I3
crystallite_localPlasticity(c,i,e) = phase_localPlasticity(material_phaseAt(c,e))
crystallite_Fe(1:3,1:3,c,i,e) = math_inv33(matmul(crystallite_Fi0(1:3,1:3,c,i,e), & crystallite_Fe(1:3,1:3,c,i,e) = math_inv33(matmul(crystallite_Fi0(1:3,1:3,c,i,e), &
crystallite_Fp0(1:3,1:3,c,i,e))) ! assuming that euler angles are given in internal strain free configuration crystallite_Fp0(1:3,1:3,c,i,e))) ! assuming that euler angles are given in internal strain free configuration
crystallite_Fp(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e) crystallite_Fp(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e)
@ -244,7 +245,7 @@ subroutine crystallite_init
enddo enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
if(any(.not. crystallite_localPlasticity) .and. .not. usePingPong) call IO_error(601) ! exit if nonlocal but no ping-pong ToDo: Why not check earlier? or in nonlocal? if(any(plasticState%nonlocal) .and. .not. usePingPong) call IO_error(601) ! exit if nonlocal but no ping-pong ToDo: Why not check earlier? or in nonlocal?
crystallite_partionedFp0 = crystallite_Fp0 crystallite_partionedFp0 = crystallite_Fp0
crystallite_partionedFi0 = crystallite_Fi0 crystallite_partionedFi0 = crystallite_Fi0
@ -271,9 +272,8 @@ subroutine crystallite_init
#ifdef DEBUG #ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) then if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) then
write(6,'(a42,1x,i10)') ' # of elements: ', eMax write(6,'(a42,1x,i10)') ' # of elements: ', eMax
write(6,'(a42,1x,i10)') 'max # of integration points/element: ', iMax write(6,'(a42,1x,i10)') ' # of integration points/element: ', iMax
write(6,'(a42,1x,i10)') 'max # of constituents/integration point: ', cMax write(6,'(a42,1x,i10)') 'max # of constituents/integration point: ', cMax
write(6,'(a42,1x,i10)') ' # of nonlocal constituents: ',count(.not. crystallite_localPlasticity)
flush(6) flush(6)
endif endif
@ -301,6 +301,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
e, & !< counter in element loop e, & !< counter in element loop
startIP, endIP, & startIP, endIP, &
s s
logical, dimension(homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: todo !ToDo: need to set some values to false for different Ngrains
#ifdef DEBUG #ifdef DEBUG
if (iand(debug_level(debug_crystallite),debug_levelSelective) /= 0 & if (iand(debug_level(debug_crystallite),debug_levelSelective) /= 0 &
@ -344,7 +345,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
crystallite_subF0(1:3,1:3,c,i,e) = crystallite_partionedF0(1:3,1:3,c,i,e) crystallite_subF0(1:3,1:3,c,i,e) = crystallite_partionedF0(1:3,1:3,c,i,e)
crystallite_subFrac(c,i,e) = 0.0_pReal crystallite_subFrac(c,i,e) = 0.0_pReal
crystallite_subStep(c,i,e) = 1.0_pReal/num%subStepSizeCryst crystallite_subStep(c,i,e) = 1.0_pReal/num%subStepSizeCryst
crystallite_todo(c,i,e) = .true. todo(c,i,e) = .true.
crystallite_converged(c,i,e) = .false. ! pretend failed step of 1/subStepSizeCryst crystallite_converged(c,i,e) = .false. ! pretend failed step of 1/subStepSizeCryst
endif homogenizationRequestsCalculation endif homogenizationRequestsCalculation
enddo; enddo enddo; enddo
@ -361,7 +362,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
endif singleRun endif singleRun
NiterationCrystallite = 0 NiterationCrystallite = 0
cutbackLooping: do while (any(crystallite_todo(:,startIP:endIP,FEsolving_execELem(1):FEsolving_execElem(2)))) cutbackLooping: do while (any(todo(:,startIP:endIP,FEsolving_execELem(1):FEsolving_execElem(2))))
NiterationCrystallite = NiterationCrystallite + 1 NiterationCrystallite = NiterationCrystallite + 1
#ifdef DEBUG #ifdef DEBUG
@ -380,8 +381,8 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
crystallite_subStep(c,i,e) = min(1.0_pReal - crystallite_subFrac(c,i,e), & crystallite_subStep(c,i,e) = min(1.0_pReal - crystallite_subFrac(c,i,e), &
num%stepIncreaseCryst * crystallite_subStep(c,i,e)) num%stepIncreaseCryst * crystallite_subStep(c,i,e))
crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > 0.0_pReal ! still time left to integrate on? todo(c,i,e) = crystallite_subStep(c,i,e) > 0.0_pReal ! still time left to integrate on?
if (crystallite_todo(c,i,e)) then if (todo(c,i,e)) then
crystallite_subF0 (1:3,1:3,c,i,e) = crystallite_subF(1:3,1:3,c,i,e) crystallite_subF0 (1:3,1:3,c,i,e) = crystallite_subF(1:3,1:3,c,i,e)
crystallite_subLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e) crystallite_subLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e)
crystallite_subLi0(1:3,1:3,c,i,e) = crystallite_Li (1:3,1:3,c,i,e) crystallite_subLi0(1:3,1:3,c,i,e) = crystallite_Li (1:3,1:3,c,i,e)
@ -415,12 +416,12 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
enddo enddo
! cant restore dotState here, since not yet calculated in first cutback after initialization ! cant restore dotState here, since not yet calculated in first cutback after initialization
crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > num%subStepMinCryst ! still on track or already done (beyond repair) todo(c,i,e) = crystallite_subStep(c,i,e) > num%subStepMinCryst ! still on track or already done (beyond repair)
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! prepare for integration ! prepare for integration
if (crystallite_todo(c,i,e)) then if (todo(c,i,e)) then
crystallite_subF(1:3,1:3,c,i,e) = crystallite_subF0(1:3,1:3,c,i,e) & crystallite_subF(1:3,1:3,c,i,e) = crystallite_subF0(1:3,1:3,c,i,e) &
+ crystallite_subStep(c,i,e) *( crystallite_partionedF (1:3,1:3,c,i,e) & + crystallite_subStep(c,i,e) *( crystallite_partionedF (1:3,1:3,c,i,e) &
-crystallite_partionedF0(1:3,1:3,c,i,e)) -crystallite_partionedF0(1:3,1:3,c,i,e))
@ -438,9 +439,9 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! integrate --- requires fully defined state array (basic + dependent state) ! integrate --- requires fully defined state array (basic + dependent state)
if (any(crystallite_todo)) call integrateState ! TODO: unroll into proper elementloop to avoid N^2 for single point evaluation if (any(todo)) call integrateState(todo) ! TODO: unroll into proper elementloop to avoid N^2 for single point evaluation
where(.not. crystallite_converged .and. crystallite_subStep > num%subStepMinCryst) & ! do not try non-converged but fully cutbacked any further where(.not. crystallite_converged .and. crystallite_subStep > num%subStepMinCryst) & ! do not try non-converged but fully cutbacked any further
crystallite_todo = .true. ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation todo = .true. ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation
enddo cutbackLooping enddo cutbackLooping
@ -610,14 +611,16 @@ subroutine crystallite_orientations
enddo; enddo; enddo enddo; enddo; enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
nonlocalPresent: if (any(plasticState%nonLocal)) then nonlocalPresent: if (any(plasticState%nonlocal)) then
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2) if (plasticState(material_phaseAt(1,e))%nonlocal) then
if (plasticState(material_phaseAt(1,e))%nonLocal) & do i = FEsolving_execIP(1),FEsolving_execIP(2)
call plastic_nonlocal_updateCompatibility(crystallite_orientation, & call plastic_nonlocal_updateCompatibility(crystallite_orientation, &
phase_plasticityInstance(material_phaseAt(i,e)),i,e) phase_plasticityInstance(material_phaseAt(i,e)),i,e)
enddo; enddo enddo
endif
enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
endif nonlocalPresent endif nonlocalPresent
@ -777,7 +780,7 @@ end subroutine crystallite_results
!> @brief calculation of stress (P) with time integration based on a residuum in Lp and !> @brief calculation of stress (P) with time integration based on a residuum in Lp and
!> intermediate acceleration of the Newton-Raphson correction !> intermediate acceleration of the Newton-Raphson correction
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
logical function integrateStress(ipc,ip,el,timeFraction) function integrateStress(ipc,ip,el,timeFraction) result(broken)
integer, intent(in):: el, & ! element index integer, intent(in):: el, & ! element index
ip, & ! integration point index ip, & ! integration point index
@ -834,9 +837,9 @@ logical function integrateStress(ipc,ip,el,timeFraction)
p, & p, &
jacoCounterLp, & jacoCounterLp, &
jacoCounterLi ! counters to check for Jacobian update jacoCounterLi ! counters to check for Jacobian update
logical :: error logical :: error,broken
integrateStress = .false. broken = .true.
if (present(timeFraction)) then if (present(timeFraction)) then
dt = crystallite_subdt(ipc,ip,el) * timeFraction dt = crystallite_subdt(ipc,ip,el) * timeFraction
@ -847,6 +850,9 @@ logical function integrateStress(ipc,ip,el,timeFraction)
F = crystallite_subF(1:3,1:3,ipc,ip,el) F = crystallite_subF(1:3,1:3,ipc,ip,el)
endif endif
call constitutive_dependentState(crystallite_partionedF(1:3,1:3,ipc,ip,el), &
crystallite_Fp(1:3,1:3,ipc,ip,el),ipc,ip,el)
Lpguess = crystallite_Lp(1:3,1:3,ipc,ip,el) ! take as first guess Lpguess = crystallite_Lp(1:3,1:3,ipc,ip,el) ! take as first guess
Liguess = crystallite_Li(1:3,1:3,ipc,ip,el) ! take as first guess Liguess = crystallite_Li(1:3,1:3,ipc,ip,el) ! take as first guess
@ -977,7 +983,6 @@ logical function integrateStress(ipc,ip,el,timeFraction)
call math_invert33(Fp_new,devNull,error,invFp_new) call math_invert33(Fp_new,devNull,error,invFp_new)
if (error) return ! error if (error) return ! error
integrateStress = .true.
crystallite_P (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) crystallite_P (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new)))
crystallite_S (1:3,1:3,ipc,ip,el) = S crystallite_S (1:3,1:3,ipc,ip,el) = S
crystallite_Lp (1:3,1:3,ipc,ip,el) = Lpguess crystallite_Lp (1:3,1:3,ipc,ip,el) = Lpguess
@ -985,6 +990,7 @@ logical function integrateStress(ipc,ip,el,timeFraction)
crystallite_Fp (1:3,1:3,ipc,ip,el) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize crystallite_Fp (1:3,1:3,ipc,ip,el) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize
crystallite_Fi (1:3,1:3,ipc,ip,el) = Fi_new crystallite_Fi (1:3,1:3,ipc,ip,el) = Fi_new
crystallite_Fe (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),invFi_new) crystallite_Fe (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),invFi_new)
broken = .false.
end function integrateStress end function integrateStress
@ -993,8 +999,9 @@ end function integrateStress
!> @brief integrate stress, state with adaptive 1st order explicit Euler method !> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize !> using Fixed Point Iteration to adapt the stepsize
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine integrateStateFPI subroutine integrateStateFPI(todo)
logical, dimension(:,:,:), intent(in) :: todo
integer :: & integer :: &
NiterationState, & !< number of iterations in state loop NiterationState, & !< number of iterations in state loop
e, & !< element index in element loop e, & !< element index in element loop
@ -1003,118 +1010,107 @@ subroutine integrateStateFPI
p, & p, &
c, & c, &
s, & s, &
sizeDotState size_pl
integer, dimension(maxval(phase_Nsources)) :: &
size_so
real(pReal) :: & real(pReal) :: &
zeta zeta
real(pReal), dimension(max(constitutive_plasticity_maxSizeDotState,constitutive_source_maxSizeDotState)) :: & real(pReal), dimension(max(constitutive_plasticity_maxSizeDotState,constitutive_source_maxSizeDotState)) :: &
r ! state residuum r ! state residuum
real(pReal), dimension(:), allocatable :: plastic_dotState_p1, plastic_dotState_p2 real(pReal), dimension(constitutive_plasticity_maxSizeDotState,2) :: &
plastic_dotState
real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState
logical :: & logical :: &
nonlocalBroken nonlocalBroken, broken
nonlocalBroken = .false. nonlocalBroken = .false.
!$OMP PARALLEL DO PRIVATE(sizeDotState,r,zeta,p,c,plastic_dotState_p1, plastic_dotState_p2,source_dotState) !$OMP PARALLEL DO PRIVATE(size_pl,size_so,r,zeta,p,c,plastic_dotState,source_dotState,broken)
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2) do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then p = material_phaseAt(g,e)
if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) c = material_phaseMemberAt(g,i,e)
call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, & crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, & crystallite_partionedFp0, &
crystallite_subdt(g,i,e), g,i,e) crystallite_subdt(g,i,e), g,i,e,p,c)
crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true.
do s = 1, phase_Nsources(p) if(broken) cycle
crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c)))
enddo
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) cycle
sizeDotState = plasticState(p)%sizeDotState size_pl = plasticState(p)%sizeDotState
plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & plasticState(p)%state(1:size_pl,c) = plasticState(p)%subState0(1:size_pl,c) &
+ plasticState(p)%dotState (1:sizeDotState,c) & + plasticState(p)%dotState (1:size_pl,c) &
* crystallite_subdt(g,i,e) * crystallite_subdt(g,i,e)
plastic_dotState_p2 = 0.0_pReal * plasticState(p)%dotState (1:sizeDotState,c) ! ToDo can be done smarter/clearer plastic_dotState(1:size_pl,2) = 0.0_pReal
do s = 1, phase_Nsources(p) do s = 1, phase_Nsources(p)
sizeDotState = sourceState(p)%p(s)%sizeDotState size_so(s) = sourceState(p)%p(s)%sizeDotState
sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%subState0(1:size_so(s),c) &
+ sourceState(p)%p(s)%dotState (1:sizeDotState,c) & + sourceState(p)%p(s)%dotState (1:size_so(s),c) &
* crystallite_subdt(g,i,e) * crystallite_subdt(g,i,e)
source_dotState(1:sizeDotState,2,s) = 0.0_pReal source_dotState(1:size_so(s),2,s) = 0.0_pReal
enddo enddo
iteration: do NiterationState = 1, num%nState iteration: do NiterationState = 1, num%nState
if(nIterationState > 1) plastic_dotState_p2 = plastic_dotState_p1 if(nIterationState > 1) plastic_dotState(1:size_pl,2) = plastic_dotState(1:size_pl,1)
plastic_dotState_p1 = plasticState(p)%dotState(:,c) plastic_dotState(1:size_pl,1) = plasticState(p)%dotState(:,c)
do s = 1, phase_Nsources(p) do s = 1, phase_Nsources(p)
sizeDotState = sourceState(p)%p(s)%sizeDotState if(nIterationState > 1) source_dotState(1:size_so(s),2,s) = source_dotState(1:size_so(s),1,s)
if(nIterationState > 1) source_dotState(1:sizeDotState,2,s) = source_dotState(1:sizeDotState,1,s) source_dotState(1:size_so(s),1,s) = sourceState(p)%p(s)%dotState(:,c)
source_dotState(1:sizeDotState,1,s) = sourceState(p)%p(s)%dotState(:,c)
enddo enddo
call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & broken = integrateStress(g,i,e)
crystallite_Fp(1:3,1:3,g,i,e), & if(broken) exit iteration
g, i, e)
crystallite_todo(g,i,e) = integrateStress(g,i,e) broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
if(.not. crystallite_todo(g,i,e)) exit iteration crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, &
crystallite_subdt(g,i,e), g,i,e,p,c)
if(broken) exit iteration
call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & zeta = damper(plasticState(p)%dotState(:,c),plastic_dotState(1:size_pl,1),&
crystallite_partionedF0, & plastic_dotState(1:size_pl,2))
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, &
crystallite_subdt(g,i,e), g,i,e)
crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do s = 1, phase_Nsources(p)
crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c)))
enddo
if(.not. crystallite_todo(g,i,e)) exit iteration
sizeDotState = plasticState(p)%sizeDotState
zeta = damper(plasticState(p)%dotState(:,c),plastic_dotState_p1,plastic_dotState_p2)
plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta & plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta &
+ plastic_dotState_p1 * (1.0_pReal - zeta) + plastic_dotState(1:size_pl,1) * (1.0_pReal - zeta)
r(1:SizeDotState) = plasticState(p)%state (1:sizeDotState,c) & r(1:size_pl) = plasticState(p)%state (1:size_pl,c) &
- plasticState(p)%subState0(1:sizeDotState,c) & - plasticState(p)%subState0(1:size_pl,c) &
- plasticState(p)%dotState (1:sizeDotState,c) * crystallite_subdt(g,i,e) - plasticState(p)%dotState (1:size_pl,c) * crystallite_subdt(g,i,e)
plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%state(1:sizeDotState,c) & plasticState(p)%state(1:size_pl,c) = plasticState(p)%state(1:size_pl,c) &
- r(1:sizeDotState) - r(1:size_pl)
crystallite_converged(g,i,e) = converged(r(1:sizeDotState), & crystallite_converged(g,i,e) = converged(r(1:size_pl), &
plasticState(p)%state(1:sizeDotState,c), & plasticState(p)%state(1:size_pl,c), &
plasticState(p)%atol(1:sizeDotState)) plasticState(p)%atol(1:size_pl))
do s = 1, phase_Nsources(p) do s = 1, phase_Nsources(p)
sizeDotState = sourceState(p)%p(s)%sizeDotState
zeta = damper(sourceState(p)%p(s)%dotState(:,c), & zeta = damper(sourceState(p)%p(s)%dotState(:,c), &
source_dotState(1:sizeDotState,1,s),& source_dotState(1:size_so(s),1,s),&
source_dotState(1:sizeDotState,2,s)) source_dotState(1:size_so(s),2,s))
sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta & sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta &
+ source_dotState(1:sizeDotState,1,s)* (1.0_pReal - zeta) + source_dotState(1:size_so(s),1,s)* (1.0_pReal - zeta)
r(1:sizeDotState) = sourceState(p)%p(s)%state (1:sizeDotState,c) & r(1:size_so(s)) = sourceState(p)%p(s)%state (1:size_so(s),c) &
- sourceState(p)%p(s)%subState0(1:sizeDotState,c) & - sourceState(p)%p(s)%subState0(1:size_so(s),c) &
- sourceState(p)%p(s)%dotState (1:sizeDotState,c) * crystallite_subdt(g,i,e) - sourceState(p)%p(s)%dotState (1:size_so(s),c) * crystallite_subdt(g,i,e)
sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%state(1:sizeDotState,c) & sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%state(1:size_so(s),c) &
- r(1:sizeDotState) - r(1:size_so(s))
crystallite_converged(g,i,e) = & crystallite_converged(g,i,e) = &
crystallite_converged(g,i,e) .and. converged(r(1:sizeDotState), & crystallite_converged(g,i,e) .and. converged(r(1:size_so(s)), &
sourceState(p)%p(s)%state(1:sizeDotState,c), & sourceState(p)%p(s)%state(1:size_so(s),c), &
sourceState(p)%p(s)%atol(1:sizeDotState)) sourceState(p)%p(s)%atol(1:size_so(s)))
enddo enddo
if(crystallite_converged(g,i,e)) then if(crystallite_converged(g,i,e)) then
crystallite_todo(g,i,e) = stateJump(g,i,e) broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c)
exit iteration exit iteration
endif endif
enddo iteration enddo iteration
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true.
nonlocalBroken = .true.
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
@ -1149,7 +1145,9 @@ end subroutine integrateStateFPI
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief integrate state with 1st order explicit Euler method !> @brief integrate state with 1st order explicit Euler method
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine integrateStateEuler subroutine integrateStateEuler(todo)
logical, dimension(:,:,:), intent(in) :: todo
integer :: & integer :: &
e, & !< element index in element loop e, & !< element index in element loop
@ -1160,29 +1158,25 @@ subroutine integrateStateEuler
s, & s, &
sizeDotState sizeDotState
logical :: & logical :: &
nonlocalBroken nonlocalBroken, broken
nonlocalBroken = .false. nonlocalBroken = .false.
!$OMP PARALLEL DO PRIVATE (sizeDotState,p,c) !$OMP PARALLEL DO PRIVATE (sizeDotState,p,c,broken)
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2) do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then p = material_phaseAt(g,e)
if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) c = material_phaseMemberAt(g,i,e)
call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, & crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, & crystallite_partionedFp0, &
crystallite_subdt(g,i,e), g,i,e) crystallite_subdt(g,i,e), g,i,e,p,c)
crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true.
do s = 1, phase_Nsources(p) if(broken) cycle
crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c)))
enddo
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) cycle
sizeDotState = plasticState(p)%sizeDotState sizeDotState = plasticState(p)%sizeDotState
plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) &
@ -1195,21 +1189,15 @@ subroutine integrateStateEuler
* crystallite_subdt(g,i,e) * crystallite_subdt(g,i,e)
enddo enddo
crystallite_todo(g,i,e) = stateJump(g,i,e) broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), &
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & crystallite_Fe(1:3,1:3,g,i,e), &
nonlocalBroken = .true. crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c)
if(.not. crystallite_todo(g,i,e)) cycle if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true.
if(broken) cycle
call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), &
g, i, e)
crystallite_todo(g,i,e) = integrateStress(g,i,e)
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
crystallite_converged(g,i,e) = crystallite_todo(g,i,e)
broken = integrateStress(g,i,e)
if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true.
crystallite_converged(g,i,e) = .not. broken
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
@ -1222,7 +1210,9 @@ end subroutine integrateStateEuler
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with 1st order Euler method with adaptive step size !> @brief integrate stress, state with 1st order Euler method with adaptive step size
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine integrateStateAdaptiveEuler subroutine integrateStateAdaptiveEuler(todo)
logical, dimension(:,:,:), intent(in) :: todo
integer :: & integer :: &
e, & ! element index in element loop e, & ! element index in element loop
@ -1233,32 +1223,28 @@ subroutine integrateStateAdaptiveEuler
s, & s, &
sizeDotState sizeDotState
logical :: & logical :: &
nonlocalBroken nonlocalBroken, broken
real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: residuum_plastic real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: residuum_plastic
real(pReal), dimension(constitutive_source_maxSizeDotState,maxval(phase_Nsources)) :: residuum_source real(pReal), dimension(constitutive_source_maxSizeDotState,maxval(phase_Nsources)) :: residuum_source
nonlocalBroken = .false. nonlocalBroken = .false.
!$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,residuum_plastic,residuum_source) !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,residuum_plastic,residuum_source,broken)
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2) do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then broken = .false.
p = material_phaseAt(g,e)
if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) c = material_phaseMemberAt(g,i,e)
call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, & crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, & crystallite_partionedFp0, &
crystallite_subdt(g,i,e), g,i,e) crystallite_subdt(g,i,e), g,i,e,p,c)
crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) if(broken) cycle
do s = 1, phase_Nsources(p)
crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c)))
enddo
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) cycle
sizeDotState = plasticState(p)%sizeDotState sizeDotState = plasticState(p)%sizeDotState
@ -1274,36 +1260,23 @@ subroutine integrateStateAdaptiveEuler
+ sourceState(p)%p(s)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) + sourceState(p)%p(s)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e)
enddo enddo
crystallite_todo(g,i,e) = stateJump(g,i,e) broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), &
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & crystallite_Fe(1:3,1:3,g,i,e), &
nonlocalBroken = .true. crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c)
if(.not. crystallite_todo(g,i,e)) cycle if(broken) cycle
call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & broken = integrateStress(g,i,e)
crystallite_Fp(1:3,1:3,g,i,e), & if(broken) cycle
g, i, e)
crystallite_todo(g,i,e) = integrateStress(g,i,e) broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) cycle
call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, & crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, & crystallite_partionedFp0, &
crystallite_subdt(g,i,e), g,i,e) crystallite_subdt(g,i,e), g,i,e,p,c)
crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) if(broken) cycle
do s = 1, phase_Nsources(p)
crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c)))
enddo
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) cycle
sizeDotState = plasticState(p)%sizeDotState sizeDotState = plasticState(p)%sizeDotState
crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState) & crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState) &
+ 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e), & + 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e), &
plasticState(p)%state(1:sizeDotState,c), & plasticState(p)%state(1:sizeDotState,c), &
@ -1311,7 +1284,6 @@ subroutine integrateStateAdaptiveEuler
do s = 1, phase_Nsources(p) do s = 1, phase_Nsources(p)
sizeDotState = sourceState(p)%p(s)%sizeDotState sizeDotState = sourceState(p)%p(s)%sizeDotState
crystallite_converged(g,i,e) = & crystallite_converged(g,i,e) = &
crystallite_converged(g,i,e) .and. converged(residuum_source(1:sizeDotState,s) & crystallite_converged(g,i,e) .and. converged(residuum_source(1:sizeDotState,s) &
+ 0.5_pReal*sourceState(p)%p(s)%dotState(:,c)*crystallite_subdt(g,i,e), & + 0.5_pReal*sourceState(p)%p(s)%dotState(:,c)*crystallite_subdt(g,i,e), &
@ -1320,6 +1292,7 @@ subroutine integrateStateAdaptiveEuler
enddo enddo
endif endif
if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true.
enddo; enddo; enddo enddo; enddo; enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
@ -1331,18 +1304,20 @@ end subroutine integrateStateAdaptiveEuler
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with 4th order explicit Runge Kutta method !> @brief integrate stress, state with 4th order explicit Runge Kutta method
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine integrateStateRK4 subroutine integrateStateRK4(todo)
real(pReal), dimension(3,3), parameter :: & logical, dimension(:,:,:), intent(in) :: todo
A = reshape([&
real(pReal), dimension(3,3), parameter :: &
A = reshape([&
0.5_pReal, 0.0_pReal, 0.0_pReal, & 0.5_pReal, 0.0_pReal, 0.0_pReal, &
0.0_pReal, 0.5_pReal, 0.0_pReal, & 0.0_pReal, 0.5_pReal, 0.0_pReal, &
0.0_pReal, 0.0_pReal, 1.0_pReal], & 0.0_pReal, 0.0_pReal, 1.0_pReal], &
[3,3]) [3,3])
real(pReal), dimension(3), parameter :: & real(pReal), dimension(3), parameter :: &
CC = [0.5_pReal, 0.5_pReal, 1.0_pReal] ! factor giving the fraction of the original timestep used for Runge Kutta Integration CC = [0.5_pReal, 0.5_pReal, 1.0_pReal] ! factor giving the fraction of the original timestep used for Runge Kutta Integration
real(pReal), dimension(4), parameter :: & real(pReal), dimension(4), parameter :: &
B = [1.0_pReal/6.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/6.0_pReal] ! weight of slope used for Runge Kutta integration (final weight divided by 6) B = [1.0_pReal/6.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/6.0_pReal] ! weight of slope used for Runge Kutta integration (final weight divided by 6)
integer :: & integer :: &
e, & ! element index in element loop e, & ! element index in element loop
@ -1355,31 +1330,28 @@ subroutine integrateStateRK4
s, & s, &
sizeDotState sizeDotState
logical :: & logical :: &
nonlocalBroken nonlocalBroken, broken
real(pReal), dimension(constitutive_plasticity_maxSizeDotState,4) :: plastic_RK4dotState
real(pReal), dimension(constitutive_source_maxSizeDotState,4,maxval(phase_Nsources)) :: source_RK4dotState real(pReal), dimension(constitutive_source_maxSizeDotState,4,maxval(phase_Nsources)) :: source_RK4dotState
real(pReal), dimension(constitutive_plasticity_maxSizeDotState,4) :: plastic_RK4dotState
nonlocalBroken = .false. nonlocalBroken = .false.
!$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,plastic_RK4dotState,source_RK4dotState) !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,source_RK4dotState,plastic_RK4dotState,broken)
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2) do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then broken = .false.
p = material_phaseAt(g,e)
if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) c = material_phaseMemberAt(g,i,e)
call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, & crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, & crystallite_partionedFp0, &
crystallite_subdt(g,i,e), g,i,e) crystallite_subdt(g,i,e), g,i,e,p,c)
crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) if(broken) cycle
do s = 1, phase_Nsources(p)
crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c)))
enddo
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) cycle
do stage = 1,3 do stage = 1,3
sizeDotState = plasticState(p)%sizeDotState sizeDotState = plasticState(p)%sizeDotState
@ -1413,31 +1385,18 @@ subroutine integrateStateRK4
* crystallite_subdt(g,i,e) * crystallite_subdt(g,i,e)
enddo enddo
call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & broken = integrateStress(g,i,e,CC(stage))
crystallite_Fp(1:3,1:3,g,i,e), & if(broken) exit
g, i, e)
crystallite_todo(g,i,e) = integrateStress(g,i,e,CC(stage)) broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) exit
call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, & crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, & crystallite_partionedFp0, &
crystallite_subdt(g,i,e)*CC(stage), g,i,e) crystallite_subdt(g,i,e)*CC(stage), g,i,e,p,c)
crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) if(broken) exit
do s = 1, phase_Nsources(p)
crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c)))
enddo
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) exit
enddo enddo
if(broken) cycle
if(.not. crystallite_todo(g,i,e)) cycle
sizeDotState = plasticState(p)%sizeDotState sizeDotState = plasticState(p)%sizeDotState
@ -1459,25 +1418,16 @@ subroutine integrateStateRK4
* crystallite_subdt(g,i,e) * crystallite_subdt(g,i,e)
enddo enddo
crystallite_todo(g,i,e) = stateJump(g,i,e) broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), &
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & crystallite_Fe(1:3,1:3,g,i,e), &
nonlocalBroken = .true. crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c)
if(.not. crystallite_todo(g,i,e)) cycle if(broken) cycle
call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & broken = integrateStress(g,i,e)
crystallite_Fp(1:3,1:3,g,i,e), & crystallite_converged(g,i,e) = .not. broken
g, i, e)
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) cycle
crystallite_todo(g,i,e) = integrateStress(g,i,e)
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
crystallite_converged(g,i,e) = crystallite_todo(g,i,e) ! consider converged if not broken
endif endif
if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true.
enddo; enddo; enddo enddo; enddo; enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
@ -1490,7 +1440,9 @@ end subroutine integrateStateRK4
!> @brief integrate stress, state with 5th order Runge-Kutta Cash-Karp method with !> @brief integrate stress, state with 5th order Runge-Kutta Cash-Karp method with
!> adaptive step size (use 5th order solution to advance = "local extrapolation") !> adaptive step size (use 5th order solution to advance = "local extrapolation")
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine integrateStateRKCK45 subroutine integrateStateRKCK45(todo)
logical, dimension(:,:,:), intent(in) :: todo
real(pReal), dimension(5,5), parameter :: & real(pReal), dimension(5,5), parameter :: &
A = reshape([& A = reshape([&
@ -1523,31 +1475,27 @@ subroutine integrateStateRKCK45
s, & s, &
sizeDotState sizeDotState
logical :: & logical :: &
nonlocalBroken nonlocalBroken, broken
real(pReal), dimension(constitutive_plasticity_maxSizeDotState,6) :: plastic_RKdotState
real(pReal), dimension(constitutive_source_maxSizeDotState,6,maxval(phase_Nsources)) :: source_RKdotState real(pReal), dimension(constitutive_source_maxSizeDotState,6,maxval(phase_Nsources)) :: source_RKdotState
real(pReal), dimension(constitutive_plasticity_maxSizeDotState,6) :: plastic_RKdotState
nonlocalBroken = .false. nonlocalBroken = .false.
!$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,plastic_RKdotState,source_RKdotState) !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,plastic_RKdotState,source_RKdotState,broken)
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2) do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then broken = .false.
p = material_phaseAt(g,e)
if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) c = material_phaseMemberAt(g,i,e)
call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, & crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, & crystallite_partionedFp0, &
crystallite_subdt(g,i,e), g,i,e) crystallite_subdt(g,i,e), g,i,e,p,c)
crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) if(broken) cycle
do s = 1, phase_Nsources(p)
crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c)))
enddo
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) cycle
do stage = 1,5 do stage = 1,5
sizeDotState = plasticState(p)%sizeDotState sizeDotState = plasticState(p)%sizeDotState
@ -1581,31 +1529,18 @@ subroutine integrateStateRKCK45
* crystallite_subdt(g,i,e) * crystallite_subdt(g,i,e)
enddo enddo
call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & broken = integrateStress(g,i,e,CC(stage))
crystallite_Fp(1:3,1:3,g,i,e), & if(broken) exit
g, i, e)
crystallite_todo(g,i,e) = integrateStress(g,i,e,CC(stage)) broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) exit
call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, & crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, & crystallite_partionedFp0, &
crystallite_subdt(g,i,e)*CC(stage), g,i,e) crystallite_subdt(g,i,e)*CC(stage), g,i,e,p,c)
crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) if(broken) exit
do s = 1, phase_Nsources(p)
crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c)))
enddo
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) exit
enddo enddo
if(broken) cycle
if(.not. crystallite_todo(g,i,e)) cycle
sizeDotState = plasticState(p)%sizeDotState sizeDotState = plasticState(p)%sizeDotState
@ -1614,7 +1549,7 @@ subroutine integrateStateRKCK45
plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) &
+ plasticState(p)%dotState (1:sizeDotState,c) & + plasticState(p)%dotState (1:sizeDotState,c) &
* crystallite_subdt(g,i,e) * crystallite_subdt(g,i,e)
crystallite_todo(g,i,e) = converged( matmul(plastic_RKdotState(1:sizeDotState,1:6),DB) & broken = .not. converged( matmul(plastic_RKdotState(1:sizeDotState,1:6),DB) &
* crystallite_subdt(g,i,e), & * crystallite_subdt(g,i,e), &
plasticState(p)%state(1:sizeDotState,c), & plasticState(p)%state(1:sizeDotState,c), &
plasticState(p)%atol(1:sizeDotState)) plasticState(p)%atol(1:sizeDotState))
@ -1627,35 +1562,28 @@ subroutine integrateStateRKCK45
sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) &
+ sourceState(p)%p(s)%dotState (1:sizeDotState,c) & + sourceState(p)%p(s)%dotState (1:sizeDotState,c) &
* crystallite_subdt(g,i,e) * crystallite_subdt(g,i,e)
crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. & broken = broken .and. .not. &
converged(matmul(source_RKdotState(1:sizeDotState,1:6,s),DB) & converged(matmul(source_RKdotState(1:sizeDotState,1:6,s),DB) &
* crystallite_subdt(g,i,e), & * crystallite_subdt(g,i,e), &
sourceState(p)%p(s)%state(1:sizeDotState,c), & sourceState(p)%p(s)%state(1:sizeDotState,c), &
sourceState(p)%p(s)%atol(1:sizeDotState)) sourceState(p)%p(s)%atol(1:sizeDotState))
enddo enddo
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & if(broken) cycle
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) cycle
crystallite_todo(g,i,e) = stateJump(g,i,e) broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), &
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & crystallite_Fe(1:3,1:3,g,i,e), &
nonlocalBroken = .true. crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c)
if(.not. crystallite_todo(g,i,e)) cycle if(broken) cycle
call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & broken = integrateStress(g,i,e)
crystallite_Fp(1:3,1:3,g,i,e), & crystallite_converged(g,i,e) = .not. broken
g, i, e)
crystallite_todo(g,i,e) = integrateStress(g,i,e)
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
crystallite_converged(g,i,e) = crystallite_todo(g,i,e) ! consider converged if not broken
endif endif
if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true.
enddo; enddo; enddo enddo; enddo; enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
if (nonlocalBroken) call nonlocalConvergenceCheck if(nonlocalBroken) call nonlocalConvergenceCheck
end subroutine integrateStateRKCK45 end subroutine integrateStateRKCK45
@ -1666,7 +1594,16 @@ end subroutine integrateStateRKCK45
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine nonlocalConvergenceCheck subroutine nonlocalConvergenceCheck
where( .not. crystallite_localPlasticity) crystallite_converged = .false. integer :: e,i,p
!$OMP PARALLEL DO PRIVATE(p)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
p = material_phaseAt(1,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
if(plasticState(p)%nonlocal) crystallite_converged(1,i,e) = .false.
enddo
enddo
!$OMP END PARALLEL DO
end subroutine nonlocalConvergenceCheck end subroutine nonlocalConvergenceCheck
@ -1688,59 +1625,6 @@ logical pure function converged(residuum,state,atol)
end function converged end function converged
!--------------------------------------------------------------------------------------------------
!> @brief calculates a jump in the state according to the current state and the current stress
!> returns true, if state jump was successfull or not needed. false indicates NaN in delta state
!--------------------------------------------------------------------------------------------------
logical function stateJump(ipc,ip,el)
integer, intent(in):: &
el, & ! element index
ip, & ! integration point index
ipc ! grain index
integer :: &
c, &
p, &
mySource, &
myOffset, &
mySize
c = material_phaseMemberAt(ipc,ip,el)
p = material_phaseAt(ipc,el)
call constitutive_collectDeltaState(crystallite_S(1:3,1:3,ipc,ip,el), &
crystallite_Fe(1:3,1:3,ipc,ip,el), &
crystallite_Fi(1:3,1:3,ipc,ip,el), &
ipc,ip,el)
myOffset = plasticState(p)%offsetDeltaState
mySize = plasticState(p)%sizeDeltaState
if( any(IEEE_is_NaN(plasticState(p)%deltaState(1:mySize,c)))) then
stateJump = .false.
return
endif
plasticState(p)%state(myOffset + 1:myOffset + mySize,c) = &
plasticState(p)%state(myOffset + 1:myOffset + mySize,c) + plasticState(p)%deltaState(1:mySize,c)
do mySource = 1, phase_Nsources(p)
myOffset = sourceState(p)%p(mySource)%offsetDeltaState
mySize = sourceState(p)%p(mySource)%sizeDeltaState
if (any(IEEE_is_NaN(sourceState(p)%p(mySource)%deltaState(1:mySize,c)))) then
stateJump = .false.
return
endif
sourceState(p)%p(mySource)%state(myOffset + 1: myOffset + mySize,c) = &
sourceState(p)%p(mySource)%state(myOffset + 1: myOffset + mySize,c) + sourceState(p)%p(mySource)%deltaState(1:mySize,c)
enddo
stateJump = .true.
end function stateJump
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Write current restart information (Field and constitutive data) to file. !> @brief Write current restart information (Field and constitutive data) to file.
! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively ! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively

66
src/material.f90
View File

@ -11,7 +11,6 @@ module material
use results use results
use IO use IO
use debug use debug
use numerics
use rotations use rotations
use discretization use discretization
@ -174,8 +173,7 @@ module material
public :: & public :: &
material_init, & material_init, &
material_allocatePlasticState, & material_allocateState, &
material_allocateSourceState, &
ELASTICITY_HOOKE_ID ,& ELASTICITY_HOOKE_ID ,&
PLASTICITY_NONE_ID, & PLASTICITY_NONE_ID, &
PLASTICITY_ISOTROPIC_ID, & PLASTICITY_ISOTROPIC_ID, &
@ -700,63 +698,35 @@ end subroutine material_parseTexture
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates the plastic state of a phase !> @brief Allocate the components of the state structure for a given phase
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine material_allocatePlasticState(phase,NipcMyPhase,& subroutine material_allocateState(state, &
sizeState,sizeDotState,sizeDeltaState) NipcMyPhase,sizeState,sizeDotState,sizeDeltaState)
class(tState), intent(out) :: &
state
integer, intent(in) :: & integer, intent(in) :: &
phase, &
NipcMyPhase, & NipcMyPhase, &
sizeState, & sizeState, &
sizeDotState, & sizeDotState, &
sizeDeltaState sizeDeltaState
plasticState(phase)%sizeState = sizeState state%sizeState = sizeState
plasticState(phase)%sizeDotState = sizeDotState state%sizeDotState = sizeDotState
plasticState(phase)%sizeDeltaState = sizeDeltaState state%sizeDeltaState = sizeDeltaState
plasticState(phase)%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition state%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition
allocate(plasticState(phase)%atol (sizeState), source=0.0_pReal) allocate(state%atol (sizeState), source=0.0_pReal)
allocate(plasticState(phase)%state0 (sizeState,NipcMyPhase), source=0.0_pReal) allocate(state%state0 (sizeState,NipcMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%partionedState0 (sizeState,NipcMyPhase), source=0.0_pReal) allocate(state%partionedState0(sizeState,NipcMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%subState0 (sizeState,NipcMyPhase), source=0.0_pReal) allocate(state%subState0 (sizeState,NipcMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%state (sizeState,NipcMyPhase), source=0.0_pReal) allocate(state%state (sizeState,NipcMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) allocate(state%dotState (sizeDotState,NipcMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%deltaState (sizeDeltaState,NipcMyPhase),source=0.0_pReal) allocate(state%deltaState(sizeDeltaState,NipcMyPhase), source=0.0_pReal)
end subroutine material_allocatePlasticState end subroutine material_allocateState
!--------------------------------------------------------------------------------------------------
!> @brief allocates the source state of a phase
!--------------------------------------------------------------------------------------------------
subroutine material_allocateSourceState(phase,of,NipcMyPhase,&
sizeState,sizeDotState,sizeDeltaState)
integer, intent(in) :: &
phase, &
of, &
NipcMyPhase, &
sizeState, sizeDotState,sizeDeltaState
sourceState(phase)%p(of)%sizeState = sizeState
sourceState(phase)%p(of)%sizeDotState = sizeDotState
sourceState(phase)%p(of)%sizeDeltaState = sizeDeltaState
sourceState(phase)%p(of)%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition
allocate(sourceState(phase)%p(of)%atol (sizeState), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%state0 (sizeState,NipcMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%partionedState0 (sizeState,NipcMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%subState0 (sizeState,NipcMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%state (sizeState,NipcMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal)
allocate(sourceState(phase)%p(of)%deltaState (sizeDeltaState,NipcMyPhase),source=0.0_pReal)
end subroutine material_allocateSourceState
end module material end module material

23
src/numerics.f90
View File

@ -20,8 +20,7 @@ module numerics
iJacoStiffness = 1, & !< frequency of stiffness update iJacoStiffness = 1, & !< frequency of stiffness update
randomSeed = 0, & !< fixed seeding for pseudo-random number generator, Default 0: use random seed randomSeed = 0, & !< fixed seeding for pseudo-random number generator, Default 0: use random seed
worldrank = 0, & !< MPI worldrank (/=0 for MPI simulations only) worldrank = 0, & !< MPI worldrank (/=0 for MPI simulations only)
worldsize = 1, & !< MPI worldsize (/=1 for MPI simulations only) worldsize = 1 !< MPI worldsize (/=1 for MPI simulations only)
numerics_integrator = 1 !< method used for state integration Default 1: fix-point iteration
integer(4), protected, public :: & integer(4), protected, public :: &
DAMASK_NumThreadsInt = 0 !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive DAMASK_NumThreadsInt = 0 !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive
real(pReal), protected, public :: & real(pReal), protected, public :: &
@ -134,8 +133,6 @@ subroutine numerics_init
defgradTolerance = IO_floatValue(line,chunkPos,2) defgradTolerance = IO_floatValue(line,chunkPos,2)
case ('ijacostiffness') case ('ijacostiffness')
iJacoStiffness = IO_intValue(line,chunkPos,2) iJacoStiffness = IO_intValue(line,chunkPos,2)
case ('integrator')
numerics_integrator = IO_intValue(line,chunkPos,2)
case ('usepingpong') case ('usepingpong')
usepingpong = IO_intValue(line,chunkPos,2) > 0 usepingpong = IO_intValue(line,chunkPos,2) > 0
case ('unitlength') case ('unitlength')
@ -176,6 +173,11 @@ subroutine numerics_init
case ('maxstaggerediter') case ('maxstaggerediter')
stagItMax = IO_intValue(line,chunkPos,2) stagItMax = IO_intValue(line,chunkPos,2)
#ifdef PETSC
case ('petsc_options')
petsc_options = trim(line(chunkPos(4):))
#endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! spectral parameters ! spectral parameters
#ifdef Grid #ifdef Grid
@ -187,8 +189,6 @@ subroutine numerics_init
err_stress_tolrel = IO_floatValue(line,chunkPos,2) err_stress_tolrel = IO_floatValue(line,chunkPos,2)
case ('err_stress_tolabs') case ('err_stress_tolabs')
err_stress_tolabs = IO_floatValue(line,chunkPos,2) err_stress_tolabs = IO_floatValue(line,chunkPos,2)
case ('petsc_options')
petsc_options = trim(line(chunkPos(4):))
case ('err_curl_tolabs') case ('err_curl_tolabs')
err_curl_tolAbs = IO_floatValue(line,chunkPos,2) err_curl_tolAbs = IO_floatValue(line,chunkPos,2)
case ('err_curl_tolrel') case ('err_curl_tolrel')
@ -206,8 +206,6 @@ subroutine numerics_init
integrationorder = IO_intValue(line,chunkPos,2) integrationorder = IO_intValue(line,chunkPos,2)
case ('structorder') case ('structorder')
structorder = IO_intValue(line,chunkPos,2) structorder = IO_intValue(line,chunkPos,2)
case ('petsc_options')
petsc_options = trim(line(chunkPos(4):))
case ('bbarstabilisation') case ('bbarstabilisation')
BBarStabilisation = IO_intValue(line,chunkPos,2) > 0 BBarStabilisation = IO_intValue(line,chunkPos,2) > 0
#endif #endif
@ -223,7 +221,6 @@ subroutine numerics_init
! writing parameters to output ! writing parameters to output
write(6,'(a24,1x,es8.1)') ' defgradTolerance: ',defgradTolerance write(6,'(a24,1x,es8.1)') ' defgradTolerance: ',defgradTolerance
write(6,'(a24,1x,i8)') ' iJacoStiffness: ',iJacoStiffness write(6,'(a24,1x,i8)') ' iJacoStiffness: ',iJacoStiffness
write(6,'(a24,1x,i8)') ' integrator: ',numerics_integrator
write(6,'(a24,1x,L8)') ' use ping pong scheme: ',usepingpong write(6,'(a24,1x,L8)') ' use ping pong scheme: ',usepingpong
write(6,'(a24,1x,es8.1,/)')' unitlength: ',numerics_unitlength write(6,'(a24,1x,es8.1,/)')' unitlength: ',numerics_unitlength
@ -266,7 +263,6 @@ subroutine numerics_init
write(6,'(a24,1x,es8.1)') ' err_curl_tolRel: ',err_curl_tolRel write(6,'(a24,1x,es8.1)') ' err_curl_tolRel: ',err_curl_tolRel
write(6,'(a24,1x,es8.1)') ' polarAlpha: ',polarAlpha write(6,'(a24,1x,es8.1)') ' polarAlpha: ',polarAlpha
write(6,'(a24,1x,es8.1)') ' polarBeta: ',polarBeta write(6,'(a24,1x,es8.1)') ' polarBeta: ',polarBeta
write(6,'(a24,1x,a)') ' PETSc_options: ',trim(petsc_options)
#endif #endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -274,16 +270,17 @@ subroutine numerics_init
#ifdef FEM #ifdef FEM
write(6,'(a24,1x,i8)') ' integrationOrder: ',integrationOrder write(6,'(a24,1x,i8)') ' integrationOrder: ',integrationOrder
write(6,'(a24,1x,i8)') ' structOrder: ',structOrder write(6,'(a24,1x,i8)') ' structOrder: ',structOrder
write(6,'(a24,1x,a)') ' PETSc_options: ',trim(petsc_options)
write(6,'(a24,1x,L8)') ' B-Bar stabilisation: ',BBarStabilisation write(6,'(a24,1x,L8)') ' B-Bar stabilisation: ',BBarStabilisation
#endif #endif
#ifdef PETSC
write(6,'(a24,1x,a)') ' PETSc_options: ',trim(petsc_options)
#endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! sanity checks ! sanity checks
if (defgradTolerance <= 0.0_pReal) call IO_error(301,ext_msg='defgradTolerance') if (defgradTolerance <= 0.0_pReal) call IO_error(301,ext_msg='defgradTolerance')
if (iJacoStiffness < 1) call IO_error(301,ext_msg='iJacoStiffness') if (iJacoStiffness < 1) call IO_error(301,ext_msg='iJacoStiffness')
if (numerics_integrator <= 0 .or. numerics_integrator >= 6) &
call IO_error(301,ext_msg='integrator')
if (numerics_unitlength <= 0.0_pReal) call IO_error(301,ext_msg='unitlength') if (numerics_unitlength <= 0.0_pReal) call IO_error(301,ext_msg='unitlength')
if (residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness') if (residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness')
if (itmax <= 1) call IO_error(301,ext_msg='itmax') if (itmax <= 1) call IO_error(301,ext_msg='itmax')

3
src/prec.f90
View File

@ -53,8 +53,7 @@ module prec
logical :: & logical :: &
nonlocal = .false. nonlocal = .false.
real(pReal), pointer, dimension(:,:) :: & real(pReal), pointer, dimension(:,:) :: &
slipRate, & !< slip rate slipRate !< slip rate
accumulatedSlip !< accumulated plastic slip
end type end type
type :: tSourceState type :: tSourceState

2
src/source_damage_anisoBrittle.f90
View File

@ -107,7 +107,7 @@ subroutine source_damage_anisoBrittle_init
if (any(prm%critDisp < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_critDisp' if (any(prm%critDisp < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_critDisp'
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP NipcMyPhase = count(material_phaseAt==p) * discretization_nIP
call material_allocateSourceState(p,sourceOffset,NipcMyPhase,1,1,0) call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0)
sourceState(p)%p(sourceOffset)%atol = config%getFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal) sourceState(p)%p(sourceOffset)%atol = config%getFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal)
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol' if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol'

2
src/source_damage_anisoDuctile.f90
View File

@ -89,7 +89,7 @@ subroutine source_damage_anisoDuctile_init
if (any(prm%critPlasticStrain < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_criticalplasticstrain' if (any(prm%critPlasticStrain < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_criticalplasticstrain'
NipcMyPhase=count(material_phaseAt==p) * discretization_nIP NipcMyPhase=count(material_phaseAt==p) * discretization_nIP
call material_allocateSourceState(p,sourceOffset,NipcMyPhase,1,1,0) call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0)
sourceState(p)%p(sourceOffset)%atol = config%getFloat('anisoductile_atol',defaultVal=1.0e-3_pReal) sourceState(p)%p(sourceOffset)%atol = config%getFloat('anisoductile_atol',defaultVal=1.0e-3_pReal)
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol'

2
src/source_damage_isoBrittle.f90
View File

@ -83,7 +83,7 @@ subroutine source_damage_isoBrittle_init
if (prm%critStrainEnergy <= 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_criticalstrainenergy' if (prm%critStrainEnergy <= 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_criticalstrainenergy'
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP NipcMyPhase = count(material_phaseAt==p) * discretization_nIP
call material_allocateSourceState(p,sourceOffset,NipcMyPhase,1,1,1) call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,1)
sourceState(p)%p(sourceOffset)%atol = config%getFloat('isobrittle_atol',defaultVal=1.0e-3_pReal) sourceState(p)%p(sourceOffset)%atol = config%getFloat('isobrittle_atol',defaultVal=1.0e-3_pReal)
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol'

2
src/source_damage_isoDuctile.f90
View File

@ -82,7 +82,7 @@ subroutine source_damage_isoDuctile_init
if (prm%critPlasticStrain <= 0.0_pReal) extmsg = trim(extmsg)//' isoductile_criticalplasticstrain' if (prm%critPlasticStrain <= 0.0_pReal) extmsg = trim(extmsg)//' isoductile_criticalplasticstrain'
NipcMyPhase=count(material_phaseAt==p) * discretization_nIP NipcMyPhase=count(material_phaseAt==p) * discretization_nIP
call material_allocateSourceState(p,sourceOffset,NipcMyPhase,1,1,0) call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0)
sourceState(p)%p(sourceOffset)%atol = config%getFloat('isoductile_atol',defaultVal=1.0e-3_pReal) sourceState(p)%p(sourceOffset)%atol = config%getFloat('isoductile_atol',defaultVal=1.0e-3_pReal)
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol'

2
src/source_thermal_dissipation.f90
View File

@ -67,7 +67,7 @@ subroutine source_thermal_dissipation_init
prm%kappa = config%getFloat('dissipation_coldworkcoeff') prm%kappa = config%getFloat('dissipation_coldworkcoeff')
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP NipcMyPhase = count(material_phaseAt==p) * discretization_nIP
call material_allocateSourceState(p,sourceOffset,NipcMyPhase,0,0,0) call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,0,0,0)
end associate end associate
enddo enddo

2
src/source_thermal_externalheat.f90
View File

@ -74,7 +74,7 @@ subroutine source_thermal_externalheat_init
prm%heat_rate = config%getFloats('externalheat_rate',requiredSize = size(prm%time)) prm%heat_rate = config%getFloats('externalheat_rate',requiredSize = size(prm%time))
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP NipcMyPhase = count(material_phaseAt==p) * discretization_nIP
call material_allocateSourceState(p,sourceOffset,NipcMyPhase,1,1,0) call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0)
end associate end associate
enddo enddo