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

merge development into kinematic hardening branch

This commit is contained in:
Zhuowen Zhao 2018-02-07 11:30:57 -05:00
parent 227720dc7d 208c4affa4
commit 7f487bb77b
14 changed files with 1991 additions and 154 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
VERSION
View File

@ -1 +1 @@
v2.0.1-1015-g532d669
v2.0.1-1035-gd80a255

55
lib/damask/test/test.py
View File

@ -109,7 +109,7 @@ class Test():
except Exception as e:
logging.critical('exception during variant execution: "{}"'.format(str(e)))
return variant+1 # return culprit
return variant+1 # return culprit
return 0
def feasible(self):
@ -322,8 +322,10 @@ class Test():
cur1Name = self.fileInCurrent(cur1)
return self.compare_Array(cur0Name,cur1Name)
def compare_Table(self,headings0,file0,headings1,file1,normHeadings='',normType=None,
absoluteTolerance=False,perLine=False,skipLines=[]):
def compare_Table(self,headings0,file0,
headings1,file1,
normHeadings='',normType=None,
absoluteTolerance=False,perLine=False,skipLines=[]):
import numpy as np
logging.info('\n '.join(['comparing ASCII Tables',file0,file1]))
@ -337,7 +339,7 @@ class Test():
data = [[] for i in range(dataLength)]
maxError = [0.0 for i in range(dataLength)]
absTol = [absoluteTolerance for i in range(dataLength)]
column = [[1 for i in range(dataLength)] for j in range(2)]
column = [[1 for i in range(dataLength)] for j in range(2)]
norm = [[] for i in range(dataLength)]
normLength = [1 for i in range(dataLength)]
@ -368,11 +370,11 @@ class Test():
key1 = ('1_' if length[i]>1 else '') + headings1[i]['label']
normKey = ('1_' if normLength[i]>1 else '') + normHeadings[i]['label']
if key0 not in table0.labels(raw = True):
raise Exception('column {} not found in 1. table...\n'.format(key0))
raise Exception('column "{}" not found in first table...\n'.format(key0))
elif key1 not in table1.labels(raw = True):
raise Exception('column {} not found in 2. table...\n'.format(key1))
raise Exception('column "{}" not found in second table...\n'.format(key1))
elif normKey not in table0.labels(raw = True):
raise Exception('column {} not found in 1. table...\n'.format(normKey))
raise Exception('column "{}" not found in first table...\n'.format(normKey))
else:
column[0][i] = table0.label_index(key0)
column[1][i] = table1.label_index(key1)
@ -400,9 +402,9 @@ class Test():
norm[i] = [1.0 for j in range(line0-len(skipLines))]
absTol[i] = True
if perLine:
logging.warning('At least one norm of {} in 1. table is 0.0, using absolute tolerance'.format(headings0[i]['label']))
logging.warning('At least one norm of "{}" in first table is 0.0, using absolute tolerance'.format(headings0[i]['label']))
else:
logging.warning('Maximum norm of {} in 1. table is 0.0, using absolute tolerance'.format(headings0[i]['label']))
logging.warning('Maximum norm of "{}" in first table is 0.0, using absolute tolerance'.format(headings0[i]['label']))
line1 = 0
while table1.data_read(): # read next data line of ASCII table
@ -414,7 +416,7 @@ class Test():
norm[i][line1-len(skipLines)])
line1 +=1
if (line0 != line1): raise Exception('found {} lines in 1. table but {} in 2. table'.format(line0,line1))
if (line0 != line1): raise Exception('found {} lines in first table but {} in second table'.format(line0,line1))
logging.info(' ********')
for i in range(dataLength):
@ -561,25 +563,28 @@ class Test():
return allclose
def compare_TableRefCur(self,headingsRef,ref,headingsCur='',cur='',normHeadings='',normType=None,\
absoluteTolerance=False,perLine=False,skipLines=[]):
def compare_TableRefCur(self,headingsRef,ref,headingsCur='',cur='',
normHeadings='',normType=None,
absoluteTolerance=False,perLine=False,skipLines=[]):
if cur == '': cur = ref
if headingsCur == '': headingsCur = headingsRef
refName = self.fileInReference(ref)
curName = self.fileInCurrent(cur)
return self.compare_Table(headingsRef,refName,headingsCur,curName,normHeadings,normType,
absoluteTolerance,perLine,skipLines)
return self.compare_Table(headingsRef,
self.fileInReference(ref),
headingsRef if headingsCur == '' else headingsCur,
self.fileInCurrent(ref if cur == '' else cur),
normHeadings,normType,
absoluteTolerance,perLine,skipLines)
def compare_TableCurCur(self,headingsCur0,Cur0,Cur1,headingsCur1='',normHeadings='',normType=None,\
absoluteTolerance=False,perLine=False,skipLines=[]):
def compare_TableCurCur(self,headingsCur0,Cur0,Cur1,
headingsCur1='',
normHeadings='',normType=None,
absoluteTolerance=False,perLine=False,skipLines=[]):
if headingsCur1 == '': headingsCur1 = headingsCur0
cur0Name = self.fileInCurrent(Cur0)
cur1Name = self.fileInCurrent(Cur1)
return self.compare_Table(headingsCur0,cur0Name,headingsCur1,cur1Name,normHeadings,normType,
absoluteTolerance,perLine,skipLines)
return self.compare_Table(headingsCur0,
self.fileInCurrent(Cur0),
headingsCur0 if headingsCur1 == '' else headingsCur1,
self.fileInCurrent(Cur1),
normHeadings,normType,absoluteTolerance,perLine,skipLines)
def report_Success(self,culprit):

17
lib/damask/util.py Normal file → Executable file
View File

@ -39,7 +39,7 @@ def srepr(arg,glue = '\n'):
if (not hasattr(arg, "strip") and
hasattr(arg, "__getitem__") or
hasattr(arg, "__iter__")):
return glue.join(srepr(x) for x in arg)
return glue.join(str(x) for x in arg)
return arg if isinstance(arg,str) else repr(arg)
# -----------------------------
@ -100,6 +100,18 @@ def execute(cmd,
if process.returncode != 0: raise RuntimeError('{} failed with returncode {}'.format(cmd,process.returncode))
return out,error
def coordGridAndSize(coordinates):
"""Determines grid count and overall physical size along each dimension of an ordered array of coordinates"""
dim = coordinates.shape[1]
coords = [np.unique(coordinates[:,i]) for i in range(dim)]
mincorner = np.array(map(min,coords))
maxcorner = np.array(map(max,coords))
grid = np.array(map(len,coords),'i')
size = grid/np.maximum(np.ones(dim,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1)
size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 equal to smallest among other ones
return grid,size
# -----------------------------
class extendableOption(Option):
"""
@ -130,7 +142,7 @@ class backgroundMessage(threading.Thread):
'hexagon': ['⬢', '⬣'],
'square': ['▖', '▘', '▝', '▗'],
'triangle': ['ᐊ', 'ᐊ', 'ᐃ', 'ᐅ', 'ᐅ', 'ᐃ'],
'amoeba': ['▖', '▏', '▘', '▔', '▝', '▕', '▗', '▂'],
'amoeba': ['▖', '▏', '▘', '▔', '▝', '▕', '▗', '▁'],
'beat': ['▁', '▂', '▃', '▅', '▆', '▇', '▇', '▆', '▅', '▃', '▂'],
'prison': ['ᚋ', 'ᚌ', 'ᚍ', 'ᚏ', 'ᚎ', 'ᚍ', 'ᚌ', 'ᚋ'],
'breath': ['ᚐ', 'ᚑ', 'ᚒ', 'ᚓ', 'ᚔ', 'ᚓ', 'ᚒ', 'ᚑ', 'ᚐ'],
@ -221,6 +233,7 @@ def leastsqBound(func, x0, args=(), bounds=None, Dfun=None, full_output=0,
def _check_func(checker, argname, thefunc, x0, args, numinputs,
output_shape=None):
from numpy import shape
"""The same as that of minpack.py"""
res = np.atleast_1d(thefunc(*((x0[:numinputs],) + args)))
if (output_shape is not None) and (shape(res) != output_shape):

1513
processing/misc/yieldSurfaceFast.py Executable file
View File

File diff suppressed because it is too large Load Diff

9
processing/post/addCurl.py
View File

@ -66,7 +66,7 @@ parser.add_option('-p','--pos','--periodiccellcenter',
dest = 'pos',
type = 'string', metavar = 'string',
help = 'label of coordinates [%default]')
parser.add_option('-d','--data',
parser.add_option('-l','--label',
dest = 'data',
action = 'extend', metavar = '<string LIST>',
help = 'label(s) of field values')
@ -139,12 +139,7 @@ for name in filenames:
table.data_readArray()
coords = [np.unique(table.data[:,coordCol+i]) for i in range(3)]
mincorner = np.array(map(min,coords))
maxcorner = np.array(map(max,coords))
grid = np.array(map(len,coords),'i')
size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1)
size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 equal to smallest among other ones
grid,size = damask.util.coordGridAndSize(table.data[:,table.label_indexrange(options.pos)])
# ------------------------------------------ process value field -----------------------------------

150
processing/post/addDivergence.py
View File

@ -9,36 +9,43 @@ import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
def merge_dicts(*dict_args):
"""Given any number of dicts, shallow copy and merge into a new dict, with precedence going to key value pairs in latter dicts."""
result = {}
for dictionary in dict_args:
result.update(dictionary)
return result
def divFFT(geomdim,field):
shapeFFT = np.array(np.shape(field))[0:3]
grid = np.array(np.shape(field)[2::-1])
N = grid.prod() # field size
n = np.array(np.shape(field)[3:]).prod() # data size
"""Calculate divergence of a vector or tensor field by transforming into Fourier space."""
shapeFFT = np.array(np.shape(field))[0:3]
grid = np.array(np.shape(field)[2::-1])
N = grid.prod() # field size
n = np.array(np.shape(field)[3:]).prod() # data size
if n == 3: dataType = 'vector'
elif n == 9: dataType = 'tensor'
field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT)
div_fourier = np.empty(field_fourier.shape[0:len(np.shape(field))-1],'c16')
field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT)
div_fourier = np.empty(field_fourier.shape[0:len(np.shape(field))-1],'c16')
# differentiation in Fourier space
TWOPIIMG = 2.0j*math.pi
einsums = {
3:'ijkl,ijkl->ijk', # vector, 3 -> 1
9:'ijkm,ijklm->ijkl', # tensor, 3x3 -> 3
}
k_sk = np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2]))/geomdim[0]
if grid[2]%2 == 0: k_sk[grid[2]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011)
# differentiation in Fourier space
TWOPIIMG = 2.0j*math.pi
k_sk = np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2]))/geomdim[0]
if grid[2]%2 == 0: k_sk[grid[2]//2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/geomdim[1]
if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/geomdim[1]
if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011)
k_si = np.arange(grid[0]//2+1)/geomdim[2]
kk, kj, ki = np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij')
k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3).astype('c16')
if dataType == 'tensor': # tensor, 3x3 -> 3
div_fourier = np.einsum('ijklm,ijkm->ijkl',field_fourier,k_s)*TWOPIIMG
elif dataType == 'vector': # vector, 3 -> 1
div_fourier = np.einsum('ijkl,ijkl->ijk',field_fourier,k_s)*TWOPIIMG
return np.fft.irfftn(div_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n/3])
k_si = np.arange(grid[0]//2+1)/geomdim[2]
kk, kj, ki = np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij')
k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3).astype('c16')
div_fourier = np.einsum(einsums[n],k_s,field_fourier)*TWOPIIMG
return np.fft.irfftn(div_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n/3])
# --------------------------------------------------------------------
@ -46,32 +53,38 @@ def divFFT(geomdim,field):
# --------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [ASCIItable(s)]', description = """
Add column(s) containing divergence of requested column(s).
Operates on periodic ordered three-dimensional data sets.
Deals with both vector- and tensor-valued fields.
Add column(s) containing curl of requested column(s).
Operates on periodic ordered three-dimensional data sets
of vector and tensor fields.
""", version = scriptID)
parser.add_option('-p','--pos','--periodiccellcenter',
dest = 'pos',
type = 'string', metavar = 'string',
help = 'label of coordinates [%default]')
parser.add_option('-v','--vector',
dest = 'vector',
parser.add_option('-l','--label',
dest = 'data',
action = 'extend', metavar = '<string LIST>',
help = 'label(s) of vector field values')
parser.add_option('-t','--tensor',
dest = 'tensor',
action = 'extend', metavar = '<string LIST>',
help = 'label(s) of tensor field values')
help = 'label(s) of field values')
parser.set_defaults(pos = 'pos',
)
(options,filenames) = parser.parse_args()
if options.vector is None and options.tensor is None:
parser.error('no data column specified.')
if options.data is None: parser.error('no data column specified.')
# --- define possible data types -------------------------------------------------------------------
datatypes = {
3: {'name': 'vector',
'shape': [3],
},
9: {'name': 'tensor',
'shape': [3,3],
},
}
# --- loop over input files ------------------------------------------------------------------------
@ -82,30 +95,27 @@ for name in filenames:
except: continue
damask.util.report(scriptName,name)
# ------------------------------------------ read header ------------------------------------------
# --- interpret header ----------------------------------------------------------------------------
table.head_read()
# ------------------------------------------ sanity checks ----------------------------------------
items = {
'tensor': {'dim': 9, 'shape': [3,3], 'labels':options.tensor, 'active':[], 'column': []},
'vector': {'dim': 3, 'shape': [3], 'labels':options.vector, 'active':[], 'column': []},
}
errors = []
remarks = []
column = {}
if table.label_dimension(options.pos) != 3: errors.append('coordinates {} are not a vector.'.format(options.pos))
else: colCoord = table.label_index(options.pos)
errors = []
active = []
for type, data in items.iteritems():
for what in (data['labels'] if data['labels'] is not None else []):
dim = table.label_dimension(what)
if dim != data['dim']: remarks.append('column {} is not a {}.'.format(what,type))
else:
items[type]['active'].append(what)
items[type]['column'].append(table.label_index(what))
coordDim = table.label_dimension(options.pos)
if coordDim != 3:
errors.append('coordinates "{}" must be three-dimensional.'.format(options.pos))
else: coordCol = table.label_index(options.pos)
for me in options.data:
dim = table.label_dimension(me)
if dim in datatypes:
active.append(merge_dicts({'label':me},datatypes[dim]))
remarks.append('differentiating {} "{}"...'.format(datatypes[dim]['name'],me))
else:
remarks.append('skipping "{}" of dimension {}...'.format(me,dim) if dim != -1 else \
'"{}" not found...'.format(me) )
if remarks != []: damask.util.croak(remarks)
if errors != []:
@ -116,32 +126,26 @@ for name in filenames:
# ------------------------------------------ assemble header --------------------------------------
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
for type, data in items.iteritems():
for label in data['active']:
table.labels_append(['divFFT({})'.format(label) if type == 'vector' else
'{}_divFFT({})'.format(i+1,label) for i in range(data['dim']//3)]) # extend ASCII header with new labels
for data in active:
table.labels_append(['divFFT({})'.format(data['label']) if data['shape'] == [3] \
else '{}_divFFT({})'.format(i+1,data['label'])
for i in range(np.prod(np.array(data['shape']))//3)]) # extend ASCII header with new labels
table.head_write()
# --------------- figure out size and grid ---------------------------------------------------------
table.data_readArray()
coords = [np.unique(table.data[:,colCoord+i]) for i in range(3)]
mincorner = np.array(map(min,coords))
maxcorner = np.array(map(max,coords))
grid = np.array(map(len,coords),'i')
size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1)
size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 equal to smallest among other ones
grid,size = damask.util.coordGridAndSize(table.data[:,table.label_indexrange(options.pos)])
# ------------------------------------------ process value field -----------------------------------
stack = [table.data]
for type, data in items.iteritems():
for i,label in enumerate(data['active']):
# we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation
stack.append(divFFT(size[::-1],
table.data[:,data['column'][i]:data['column'][i]+data['dim']].
reshape(grid[::-1].tolist()+data['shape'])))
for data in active:
# we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation
stack.append(divFFT(size[::-1],
table.data[:,table.label_indexrange(data['label'])].
reshape(grid[::-1].tolist()+data['shape'])))
# ------------------------------------------ output result -----------------------------------------

14
processing/post/addGaussian.py
View File

@ -18,7 +18,7 @@ scriptID = ' '.join([scriptName,damask.version])
parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [ASCIItable(s)]', description = """
Add column(s) containing Gaussian filtered values of requested column(s).
Operates on periodic and non-periodic ordered three-dimensional data sets.
For Details see scipy.ndimage documentation.
For details see scipy.ndimage documentation.
""", version = scriptID)
@ -43,15 +43,14 @@ parser.add_option('--sigma',
parser.add_option('--periodic',
dest = 'periodic',
action = 'store_true',
help = 'assume periodic grain structure'
)
help = 'assume periodic grain structure')
parser.set_defaults(pos = 'pos',
order = 0,
sigma = 1,
periodic = False
periodic = False,
)
(options,filenames) = parser.parse_args()
@ -110,12 +109,7 @@ for name in filenames:
table.data_readArray()
coords = [np.unique(table.data[:,colCoord+i]) for i in range(3)]
mincorner = np.array(map(min,coords))
maxcorner = np.array(map(max,coords))
grid = np.array(map(len,coords),'i')
size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1)
size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 equal to smallest among other ones
grid,size = damask.util.coordGridAndSize(table.data[:,table.label_indexrange(options.pos)])
# ------------------------------------------ process value field -----------------------------------

9
processing/post/addGradient.py
View File

@ -63,7 +63,7 @@ parser.add_option('-p','--pos','--periodiccellcenter',
dest = 'pos',
type = 'string', metavar = 'string',
help = 'label of coordinates [%default]')
parser.add_option('-d','--data',
parser.add_option('-l','--label',
dest = 'data',
action = 'extend', metavar = '<string LIST>',
help = 'label(s) of field values')
@ -135,12 +135,7 @@ for name in filenames:
table.data_readArray()
coords = [np.unique(table.data[:,coordCol+i]) for i in range(3)]
mincorner = np.array(map(min,coords))
maxcorner = np.array(map(max,coords))
grid = np.array(map(len,coords),'i')
size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1)
size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 equal to smallest among other ones
grid,size = damask.util.coordGridAndSize(table.data[:,table.label_indexrange(options.pos)])
# ------------------------------------------ process value field -----------------------------------

115
src/DAMASK_spectral.f90 Normal file → Executable file
View File

@ -78,7 +78,8 @@ program DAMASK_spectral
FIELD_UNDEFINED_ID, &
FIELD_MECH_ID, &
FIELD_THERMAL_ID, &
FIELD_DAMAGE_ID
FIELD_DAMAGE_ID, &
utilities_calcPlasticity
use spectral_mech_Basic
use spectral_mech_AL
use spectral_mech_Polarisation
@ -156,6 +157,19 @@ program DAMASK_spectral
MPI_finalize, &
MPI_allreduce, &
PETScFinalize
!--------------------------------------------------------------------------------------------------
! variables related to stop criterion for yielding
real(pReal) :: plasticWorkOld, plasticWorkNew, & ! plastic work
eqTotalStrainOld, eqTotalStrainNew, & ! total equivalent strain
eqPlasticStrainOld, eqPlasticStrainNew, & ! total equivalent plastic strain
eqStressOld, eqStressNew , & ! equivalent stress
yieldStopValue
real(pReal), dimension(3,3) :: yieldStress,yieldStressOld,yieldStressNew, &
plasticStrainOld, plasticStrainNew, plasticStrainRate
integer(pInt) :: yieldResUnit = 0_pInt
integer(pInt) :: stressstrainUnit = 0_pInt
character(len=13) :: stopFlag
logical :: yieldStop, yieldStopSatisfied
!--------------------------------------------------------------------------------------------------
! init DAMASK (all modules)
@ -213,6 +227,8 @@ program DAMASK_spectral
!--------------------------------------------------------------------------------------------------
! reading the load case and assign values to the allocated data structure
yieldStop = .False.
yieldStopSatisfied = .False.
rewind(FILEUNIT)
do
line = IO_read(FILEUNIT)
@ -287,10 +303,30 @@ program DAMASK_spectral
temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j)
enddo
loadCases(currentLoadCase)%rotation = math_plain9to33(temp_valueVector)
case('totalstrain')
yieldStop = .True.
stopFlag = 'totalStrain'
yieldStopValue = IO_floatValue(line,chunkPos,i+1_pInt)
case('plasticstrain')
yieldStop = .True.
stopFlag = 'plasticStrain'
yieldStopValue = IO_floatValue(line,chunkPos,i+1_pInt)
case('plasticwork')
yieldStop = .True.
stopFlag = 'plasticWork'
yieldStopValue = IO_floatValue(line,chunkPos,i+1_pInt)
end select
enddo; enddo
close(FILEUNIT)
if(yieldStop) then ! initialize variables related to yield stop
yieldStressNew = 0.0_pReal
plasticStrainNew = 0.0_pReal
eqStressNew = 0.0_pReal
eqTotalStrainNew = 0.0_pReal
eqPlasticStrainNew = 0.0_pReal
plasticWorkNew = 0.0_pReal
endif
!--------------------------------------------------------------------------------------------------
! consistency checks and output of load case
loadCases(1)%followFormerTrajectory = .false. ! cannot guess along trajectory for first inc of first currentLoadCase
@ -704,11 +740,84 @@ program DAMASK_spectral
restartWrite = .true. ! set restart parameter for FEsolving
lastRestartWritten = inc ! QUESTION: first call to CPFEM_general will write?
endif
<<<<<<< HEAD
endif skipping
=======
else forwarding
time = time + timeinc
guess = .true.
endif forwarding
yieldCheck: if(yieldStop) then ! check if it yields or satisfies the certain stop condition
yieldStressOld = yieldStressNew
plasticStrainOld = plasticStrainNew
eqStressOld = eqStressNew
eqTotalStrainOld = eqTotalStrainNew
eqPlasticStrainOld = eqPlasticStrainNew
plasticWorkOld = plasticWorkNew
call utilities_calcPlasticity(yieldStressNew, plasticStrainNew, eqStressNew, eqTotalStrainNew, &
eqPlasticStrainNew, plasticWorkNew, loadCases(currentLoadCase)%rotation)
if (worldrank == 0) then ! output the stress-strain curve to file if yield stop criterion is used
if ((currentLoadCase == 1_pInt) .and. (inc == 1_pInt)) then
open(newunit=stressstrainUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//&
'.stressstrain',form='FORMATTED',status='REPLACE')
write(stressstrainUnit,*) 0.0_pReal, 0.0_pReal
write(stressstrainUnit,*) eqTotalStrainNew, eqStressNew
close(stressstrainUnit)
else
open(newunit=stressstrainUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//&
'.stressstrain',form='FORMATTED', position='APPEND', status='OLD')
write(stressstrainUnit,*) eqTotalStrainNew, eqStressNew
close(stressstrainUnit)
endif
endif
if(stopFlag == 'totalStrain') then
if(eqTotalStrainNew > yieldStopValue) then
yieldStress = yieldStressOld * (eqTotalStrainNew - yieldStopValue)/(eqTotalStrainNew - eqTotalStrainOld) & ! linear interpolation of stress values
+ yieldStressNew * (yieldStopValue - eqTotalStrainOld)/(eqTotalStrainNew - eqTotalStrainOld)
plasticStrainRate = (plasticStrainNew - plasticStrainOld)/(time - time0) ! calculate plastic strain rate
yieldStopSatisfied = .True.
endif
elseif(stopFlag == 'plasticStrain') then
if(eqPlasticStrainNew > yieldStopValue) then
yieldStress = yieldStressOld * (eqPlasticStrainNew - yieldStopValue)/(eqPlasticStrainNew - eqPlasticStrainOld) &
+ yieldStressNew * (yieldStopValue - eqPlasticStrainOld)/(eqPlasticStrainNew - eqPlasticStrainOld)
plasticStrainRate = (plasticStrainNew - plasticStrainOld)/(time - time0)
yieldStopSatisfied = .True.
endif
elseif(stopFlag == 'plasticWork') then
if(plasticWorkNew > yieldStopValue) then
yieldStress = yieldStressOld * (plasticWorkNew - yieldStopValue)/(plasticWorkNew - plasticWorkOld) &
+ yieldStressNew * (yieldStopValue - plasticWorkOld)/(plasticWorkNew - plasticWorkOld)
plasticStrainRate = (plasticStrainNew - plasticStrainOld)/(time - time0)
yieldStopSatisfied = .True.
endif
endif
endif yieldCheck
if (yieldStopSatisfied) then ! when yield, write the yield stress and strain rate to file and quit the job
if (worldrank == 0) then
open(newunit=yieldResUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//&
'.yield',form='FORMATTED',status='REPLACE')
do i = 1_pInt,3_pInt
write(yieldResUnit,*) (yieldStress(i,j), j=1,3)
enddo
do i = 1_pInt,3_pInt
write(yieldResUnit,*) (plasticStrainRate(i,j), j=1,3)
enddo
close(yieldResUnit)
call quit(0_pInt)
endif
endif
>>>>>>> development
enddo incLooping
enddo loadCaseLooping
!--------------------------------------------------------------------------------------------------
! report summary of whole calculation
write(6,'(/,a)') ' ###########################################################################'

6
src/IO.f90
View File

@ -1268,7 +1268,7 @@ integer(pInt) function IO_countNumericalDataLines(fileUnit)
line = IO_read(fileUnit)
chunkPos = IO_stringPos(line)
tmp = IO_lc(IO_stringValue(line,chunkPos,1_pInt))
if (scan(tmp,"abcdefghijklmnopqrstuvwxyz")/=0) then ! found keyword
if (verify(trim(tmp) ,"0123456789")/=0) then ! found keyword
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
else
@ -1839,12 +1839,12 @@ integer(pInt) function IO_verifyIntValue (string,validChars,myName)
invalidWhere = verify(string,validChars)
if (invalidWhere == 0_pInt) then
read(UNIT=string,iostat=readStatus,FMT=*) IO_verifyIntValue ! no offending chars found
if (readStatus /= 0_pInt) & ! error during string to float conversion
if (readStatus /= 0_pInt) & ! error during string to integer conversion
call IO_warning(203_pInt,ext_msg=myName//'"'//string//'"')
else
call IO_warning(202_pInt,ext_msg=myName//'"'//string//'"') ! complain about offending characters
read(UNIT=string(1_pInt:invalidWhere-1_pInt),iostat=readStatus,FMT=*) IO_verifyIntValue ! interpret remaining string
if (readStatus /= 0_pInt) & ! error during string to float conversion
if (readStatus /= 0_pInt) & ! error during string to integer conversion
call IO_warning(203_pInt,ext_msg=myName//'"'//string(1_pInt:invalidWhere-1_pInt)//'"')
endif

4
src/crystallite.f90 Normal file → Executable file
View File

@ -55,14 +55,14 @@ module crystallite
crystallite_Li0, & !< intermediate velocitiy grad at start of FE inc
crystallite_partionedLi0,& !< intermediate velocity grad at start of homog inc
crystallite_Fe, & !< current "elastic" def grad (end of converged time step)
crystallite_P !< 1st Piola-Kirchhoff stress per grain
crystallite_P, & !< 1st Piola-Kirchhoff stress per grain
crystallite_subF !< def grad to be reached at end of crystallite inc
real(pReal), dimension(:,:,:,:,:), allocatable, private :: &
crystallite_subFe0,& !< "elastic" def grad at start of crystallite inc
crystallite_invFp, & !< inverse of current plastic def grad (end of converged time step)
crystallite_subFp0,& !< plastic def grad at start of crystallite inc
crystallite_invFi, & !< inverse of current intermediate def grad (end of converged time step)
crystallite_subFi0,& !< intermediate def grad at start of crystallite inc
crystallite_subF, & !< def grad to be reached at end of crystallite inc
crystallite_subF0, & !< def grad at start of crystallite inc
crystallite_subLp0,& !< plastic velocity grad at start of crystallite inc
crystallite_subLi0,& !< intermediate velocity grad at start of crystallite inc

65
src/math.f90 Normal file → Executable file
View File

@ -145,6 +145,7 @@ module math
math_sampleGaussVar, &
math_symmetricEulers, &
math_eigenvectorBasisSym33, &
math_eigenvectorBasisSym33_log, &
math_eigenvectorBasisSym, &
math_eigenValuesVectorsSym33, &
math_eigenValuesVectorsSym, &
@ -2120,6 +2121,70 @@ function math_eigenvectorBasisSym33(m)
end function math_eigenvectorBasisSym33
!--------------------------------------------------------------------------------------------------
!> @brief logarithm eigenvector basis of symmetric 33 matrix m
!--------------------------------------------------------------------------------------------------
function math_eigenvectorBasisSym33_log(m)
implicit none
real(pReal), dimension(3,3) :: math_eigenvectorBasisSym33_log
real(pReal), dimension(3) :: invariants, values
real(pReal), dimension(3,3), intent(in) :: m
real(pReal) :: P, Q, rho, phi
real(pReal), parameter :: TOL=1.e-14_pReal
real(pReal), dimension(3,3,3) :: N, EB
invariants = math_invariantsSym33(m)
EB = 0.0_pReal
P = invariants(2)-invariants(1)**2.0_pReal/3.0_pReal
Q = -2.0_pReal/27.0_pReal*invariants(1)**3.0_pReal+product(invariants(1:2))/3.0_pReal-invariants(3)
threeSimilarEigenvalues: if(all(abs([P,Q]) < TOL)) then
values = invariants(1)/3.0_pReal
! this is not really correct, but at least the basis is correct
EB(1,1,1)=1.0_pReal
EB(2,2,2)=1.0_pReal
EB(3,3,3)=1.0_pReal
else threeSimilarEigenvalues
rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal
phi=acos(math_limit(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal))
values = 2.0_pReal*rho**(1.0_pReal/3.0_pReal)* &
[cos(phi/3.0_pReal), &
cos((phi+2.0_pReal*PI)/3.0_pReal), &
cos((phi+4.0_pReal*PI)/3.0_pReal) &
] + invariants(1)/3.0_pReal
N(1:3,1:3,1) = m-values(1)*math_I3
N(1:3,1:3,2) = m-values(2)*math_I3
N(1:3,1:3,3) = m-values(3)*math_I3
twoSimilarEigenvalues: if(abs(values(1)-values(2)) < TOL) then
EB(1:3,1:3,3)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,2))/ &
((values(3)-values(1))*(values(3)-values(2)))
EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,3)
elseif(abs(values(2)-values(3)) < TOL) then twoSimilarEigenvalues
EB(1:3,1:3,1)=math_mul33x33(N(1:3,1:3,2),N(1:3,1:3,3))/ &
((values(1)-values(2))*(values(1)-values(3)))
EB(1:3,1:3,2)=math_I3-EB(1:3,1:3,1)
elseif(abs(values(3)-values(1)) < TOL) then twoSimilarEigenvalues
EB(1:3,1:3,2)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,3))/ &
((values(2)-values(1))*(values(2)-values(3)))
EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,2)
else twoSimilarEigenvalues
EB(1:3,1:3,1)=math_mul33x33(N(1:3,1:3,2),N(1:3,1:3,3))/ &
((values(1)-values(2))*(values(1)-values(3)))
EB(1:3,1:3,2)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,3))/ &
((values(2)-values(1))*(values(2)-values(3)))
EB(1:3,1:3,3)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,2))/ &
((values(3)-values(1))*(values(3)-values(2)))
endif twoSimilarEigenvalues
endif threeSimilarEigenvalues
math_eigenvectorBasisSym33_log = log(sqrt(values(1))) * EB(1:3,1:3,1) &
+ log(sqrt(values(2))) * EB(1:3,1:3,2) &
+ log(sqrt(values(3))) * EB(1:3,1:3,3)
end function math_eigenvectorBasisSym33_log
!--------------------------------------------------------------------------------------------------
!> @brief rotational part from polar decomposition of 33 tensor m
!--------------------------------------------------------------------------------------------------

64
src/mesh.f90
View File

@ -1642,11 +1642,12 @@ subroutine mesh_marc_get_matNumber(fileUnit)
if (len(trim(line))/=0_pInt) then
chunkPos = IO_stringPos(line)
data_blocks = IO_intValue(line,chunkPos,1_pInt)
endif
endif
allocate(Marc_matNumber(data_blocks))
do i=1_pInt,data_blocks ! read all data blocks
read (fileUnit,610,END=620) line
chunkPos = IO_stringPos(line)
Marc_matNumber = (/Marc_matNumber, IO_intValue(line,chunkPos,1_pInt)/)
Marc_matNumber(i) = IO_intValue(line,chunkPos,1_pInt)
do j=1_pint,2_pInt + hypoelasticTableStyle ! read 2 or 3 remaining lines of data block
read (fileUnit,610,END=620) line
enddo
@ -1809,12 +1810,11 @@ subroutine mesh_marc_count_cpElements(fileUnit)
do i=1_pInt,3_pInt+hypoelasticTableStyle ! Skip 3 or 4 lines
read (fileUnit,610,END=620) line
enddo
mesh_NcpElems = mesh_NcpElems + IO_countContinuousIntValues(fileUnit) ! why not simply mesh_NcpElems = IO_countContinuousIntValues(fileUnit)? keyword hypoelastic might appear several times so the next line probably should not be there
mesh_NcpElems = mesh_NcpElems + IO_countContinuousIntValues(fileUnit) ! why not simply mesh_NcpElems = IO_countContinuousIntValues(fileUnit)? not fully correct as hypoelastic can have multiple data fields, needs update
exit
endif
enddo
else ! Marc2017 and later
call IO_error(error_ID=701_pInt)
else ! Marc2017 and later
do
read (fileUnit,610,END=620) line
chunkPos = IO_stringPos(line)
@ -1839,6 +1839,7 @@ subroutine mesh_marc_map_elements(fileUnit)
use math, only: math_qsort
use IO, only: IO_lc, &
IO_intValue, &
IO_stringValue, &
IO_stringPos, &
IO_continuousIntValues
@ -1847,7 +1848,8 @@ subroutine mesh_marc_map_elements(fileUnit)
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos
character(len=300) :: line
character(len=300) :: line, &
tmp
integer(pInt), dimension (1_pInt+mesh_NcpElems) :: contInts
integer(pInt) :: i,cpElem = 0_pInt
@ -1856,25 +1858,47 @@ subroutine mesh_marc_map_elements(fileUnit)
610 FORMAT(A300)
contInts = 0_pInt
rewind(fileUnit)
do
read (fileUnit,610,END=660) line
chunkPos = IO_stringPos(line)
if( IO_lc(IO_stringValue(line,chunkPos,1_pInt)) == 'hypoelastic' ) then
do i=1_pInt,3_pInt+hypoelasticTableStyle ! skip three (or four if new table style!) lines
read (fileUnit,610,END=660) line
enddo
contInts = IO_continuousIntValues(fileUnit,mesh_NcpElems,mesh_nameElemSet,&
if (MarcVersion < 13) then ! Marc 2016 or earlier
if( IO_lc(IO_stringValue(line,chunkPos,1_pInt)) == 'hypoelastic' ) then
do i=1_pInt,3_pInt+hypoelasticTableStyle ! skip three (or four if new table style!) lines
read (fileUnit,610,END=660) line
enddo
contInts = IO_continuousIntValues(fileUnit,mesh_NcpElems,mesh_nameElemSet,&
mesh_mapElemSet,mesh_NelemSets)
do i = 1_pInt,contInts(1)
cpElem = cpElem+1_pInt
mesh_mapFEtoCPelem(1,cpElem) = contInts(1_pInt+i)
mesh_mapFEtoCPelem(2,cpElem) = cpElem
enddo
endif
enddo
660 call math_qsort(mesh_mapFEtoCPelem,1_pInt,int(size(mesh_mapFEtoCPelem,2_pInt),pInt)) ! should be mesh_NcpElems
exit
endif
else ! Marc2017 and later
if ( IO_lc(IO_stringValue(line,chunkPos,1_pInt)) == 'connectivity') then
read (fileUnit,610,END=660) line
chunkPos = IO_stringPos(line)
if(any(Marc_matNumber==IO_intValue(line,chunkPos,6_pInt))) then
do
read (fileUnit,610,END=660) line
chunkPos = IO_stringPos(line)
tmp = IO_lc(IO_stringValue(line,chunkPos,1_pInt))
if (verify(trim(tmp),"0123456789")/=0) then ! found keyword
exit
else
contInts(1) = contInts(1) + 1_pInt
read (tmp,*) contInts(contInts(1)+1)
endif
enddo
endif
endif
endif
enddo
660 do i = 1_pInt,contInts(1)
cpElem = cpElem+1_pInt
mesh_mapFEtoCPelem(1,cpElem) = contInts(1_pInt+i)
mesh_mapFEtoCPelem(2,cpElem) = cpElem
enddo
call math_qsort(mesh_mapFEtoCPelem,1_pInt,int(size(mesh_mapFEtoCPelem,2_pInt),pInt)) ! should be mesh_NcpElems
end subroutine mesh_marc_map_elements

122
src/spectral_utilities.f90 Normal file → Executable file
View File

@ -145,7 +145,8 @@ module spectral_utilities
FIELD_UNDEFINED_ID, &
FIELD_MECH_ID, &
FIELD_THERMAL_ID, &
FIELD_DAMAGE_ID
FIELD_DAMAGE_ID, &
utilities_calcPlasticity
private :: &
utilities_getFreqDerivative
@ -1041,6 +1042,125 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
end subroutine utilities_constitutiveResponse
!--------------------------------------------------------------------------------------------------
!> @brief calculates yield stress, plastic strain, total strain and their equivalent values
!--------------------------------------------------------------------------------------------------
subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTotalStrain, &
eqPlasticStrain, plasticWork, rotation_BC)
use crystallite, only: &
crystallite_Fe, &
crystallite_P, &
crystallite_subF
use material, only: &
homogenization_maxNgrains
use mesh, only: &
mesh_maxNips,&
mesh_NcpElems
use math, only: &
math_det33, &
math_inv33, &
math_mul33x33, &
math_trace33, &
math_transpose33, &
math_equivStrain33, &
math_equivStress33, &
math_rotate_forward33, &
math_identity2nd, &
math_crossproduct, &
math_eigenvectorBasisSym, &
math_eigenvectorBasisSym33, &
math_eigenvectorBasisSym33_log, &
math_eigenValuesVectorsSym33
implicit none
real(pReal), intent(inout) :: eqStress, eqPlasticStrain, plasticWork
real(pReal), intent(out) :: eqTotalStrain
real(pReal), dimension(3,3),intent(out) :: yieldStress, plasticStrain
real(pReal), intent(in), dimension(3,3) :: rotation_BC !< rotation of load frame
real(pReal), dimension(3,3) :: cauchy, P_av, F_av, Ve_av !< average
real(pReal), dimension(3) :: Values, S
real(pReal), dimension(3,3) :: Vectors, diag
real(pReal), dimension(3,3) :: &
Vp, F_temp, U, VT, R, V, V_total
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
Be, Ve, Fe
real(pReal), dimension(15) :: WORK !< previous deformation gradient
integer(pInt) :: INFO, i, j, k, l, ierr
real(pReal) :: wgtm
real(pReal) :: eqStressOld, eqPlasticStrainOld, plasticWorkOld
external :: dgesvd
eqStressOld = eqStress
eqPlasticStrainOld = eqPlasticStrain
plasticWorkOld = plasticWork
wgtm = 1.0_pReal/real(mesh_NcpElems*mesh_maxNips*homogenization_maxNgrains,pReal)
diag = 0.0_pReal
P_av = sum(sum(sum(crystallite_P,dim=5),dim=4),dim=3) * wgtm
call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
P_av = math_rotate_forward33(P_av,rotation_BC)
F_av = sum(sum(sum(crystallite_subF,dim=5),dim=4),dim=3) * wgtm
call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
F_av = math_rotate_forward33(F_av,rotation_BC)
cauchy = 1.0_pReal/math_det33(F_av)*math_mul33x33(P_av,transpose(F_av))
yieldStress = cauchy
eqStress = math_equivStress33(cauchy)
F_temp = F_av
call dgesvd ('A', 'A', 3, 3, F_temp, 3, S, U, 3, VT, 3, WORK, 15, INFO) ! singular value decomposition
R = math_mul33x33(U, VT) ! rotation of polar decomposition
V = math_mul33x33(F_av,math_inv33(R))
call math_eigenValuesVectorsSym33(V,Values,Vectors)
do l = 1_pInt, 3_pInt
if (Values(l) < 0.0_pReal) then
Values(l) = -Values(l)
Vectors(1:3, l) = -Vectors(1:3, l)
endif
Values(l) = log(Values(l))
diag(l,l) = Values(l)
enddo
if (dot_product(Vectors(1:3,1),Vectors(1:3,2)) /= 0) then
Vectors(1:3,2) = math_crossproduct(Vectors(1:3,3), Vectors(1:3,1))
Vectors(1:3,2) = Vectors(1:3,2)/sqrt(dot_product(Vectors(1:3,2),Vectors(1:3,2)))
endif
if (dot_product(Vectors(1:3,2),Vectors(1:3,3)) /= 0) then
Vectors(1:3,3) = math_crossproduct(Vectors(1:3,1), Vectors(1:3,2))
Vectors(1:3,3) = Vectors(1:3,3)/sqrt(dot_product(Vectors(1:3,3),Vectors(1:3,3)))
endif
if (dot_product(Vectors(1:3,3),Vectors(1:3,1)) /= 0) then
Vectors(1:3,1) = math_crossproduct(Vectors(1:3,2), Vectors(1:3,3))
Vectors(1:3,1) = Vectors(1:3,1)/sqrt(dot_product(Vectors(1:3,1),Vectors(1:3,1)))
endif
V_total = REAL(math_mul33x33(Vectors, math_mul33x33(diag, transpose(Vectors))))
eqTotalStrain = math_equivStrain33(V_total)
do k = 1_pInt, mesh_NcpElems; do j = 1_pInt, mesh_maxNips; do i = 1_pInt,homogenization_maxNgrains
Fe(1:3,1:3,i,j,k) = crystallite_Fe(1:3,1:3,i,j,k)
Fe(1:3,1:3,i,j,k) = math_rotate_forward33(Fe(1:3,1:3,i,j,k),rotation_BC)
Be(1:3,1:3,i,j,k) = math_mul33x33(Fe(1:3,1:3,i,j,k),math_transpose33(Fe(1:3,1:3,i,j,k))) ! elastic part of left Cauchy–Green deformation tensor
Ve(1:3,1:3,i,j,k) = math_eigenvectorBasisSym33_log(Be(1:3,1:3,i,j,k))
enddo; enddo; enddo
Ve_av = sum(sum(sum(Ve,dim=5),dim=4),dim=3) * wgtm
call MPI_Allreduce(MPI_IN_PLACE,Ve_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
Vp = V_total - Ve_av
eqPlasticStrain = math_equivStrain33(Vp)
plasticStrain = Vp
plasticWork = plasticWorkOld + 0.5*(eqStressOld + eqStress) * (eqPlasticStrain - eqPlasticStrainOld)
end subroutine utilities_calcPlasticity
!--------------------------------------------------------------------------------------------------
!> @brief calculates forward rate, either guessing or just add delta/timeinc