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 32_NewStyleNonlocal-3

This commit is contained in:
Martin Diehl 2019-01-31 12:42:35 +01:00
parent 854afb7107 6b66563be7
commit a3e61c82dc
26 changed files with 4929 additions and 7323 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 5ed6a1f60b412eb46ff6820cf03b684095ff1f75 Subproject commit beb9682fff7d4d6c65aba12ffd04c7441dc6ba6b

2
VERSION
View File

@ -1 +1 @@
v2.0.2-1404-g3f40eeac v2.0.2-1634-g370b23d5

2
examples/ConfigFiles/Homogenization_None_Dummy.config
View File

@ -1,3 +1,3 @@
[directSX] [directSX]
type none mech none

4
examples/ConfigFiles/Phase_Isotropic_AluminumIsotropic.config
View File

@ -11,11 +11,11 @@ lattice_structure isotropic
c11 110.9e9 c11 110.9e9
c12 58.34e9 c12 58.34e9
taylorfactor 3 m 3
tau0 31e6 tau0 31e6
gdot0 0.001 gdot0 0.001
n 20 n 20
h0 75e6 h0 75e6
tausat 63e6 tausat 63e6
w0 2.25 a 2.25
atol_resistance 1 atol_resistance 1

2
processing/post/viewTable.py
View File

@ -68,7 +68,7 @@ for name in filenames:
(['data'] if options.data else []) + (['data'] if options.data else []) +
[] []
) )
damask.util.report(scriptName,name + ('' if details == '' else ' -- '+details)) damask.util.report(scriptName,(name if name is not None else '') + ('' if details == '' else ' -- '+details))
# ------------------------------------------ output head --------------------------------------- # ------------------------------------------ output head ---------------------------------------

185
processing/pre/geom_fromVPSC.py
View File

@ -1,185 +0,0 @@
#!/usr/bin/env python2.7
# -*- coding: UTF-8 no BOM -*-
import os,sys,math
import numpy as np
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
#--------------------------------------------------------------------------------------------------
# MAIN
#--------------------------------------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
Generate geometry description and material configuration from input files used by R.A. Lebensohn.
""", version = scriptID)
parser.add_option('--column', dest='column', type='int', metavar = 'int',
help='data column to discriminate between both phases [%default]')
parser.add_option('-t','--threshold', dest='threshold', type='float', metavar = 'float',
help='threshold value for phase discrimination [%default]')
parser.add_option('--homogenization', dest='homogenization', type='int', metavar = 'int',
help='homogenization index for <microstructure> configuration [%default]')
parser.add_option('--phase', dest='phase', type='int', nargs = 2, metavar = 'int int',
help='phase indices for <microstructure> configuration %default')
parser.add_option('--crystallite', dest='crystallite', type='int', metavar = 'int',
help='crystallite index for <microstructure> configuration [%default]')
parser.add_option('--compress', dest='compress', action='store_true',
help='lump identical microstructure and texture information [%default]')
parser.add_option('-p', '--precision', dest='precision', choices=['0','1','2','3'], metavar = 'int',
help = 'euler angles decimal places for output format and compressing {0,1,2,3} [2]')
parser.set_defaults(column = 7)
parser.set_defaults(threshold = 1.0)
parser.set_defaults(homogenization = 1)
parser.set_defaults(phase = [1,2])
parser.set_defaults(crystallite = 1)
parser.set_defaults(config = False)
parser.set_defaults(compress = False)
parser.set_defaults(precision = '2')
(options,filenames) = parser.parse_args()
if filenames == []: filenames = [None]
for name in filenames:
try:
table = damask.ASCIItable(name = name,
outname = os.path.splitext(name)[-2]+'.geom' if name else name,
buffered = False,
labeled = False)
except: continue
damask.util.report(scriptName,name)
info = {
'grid': np.zeros(3,'i'),
'size': np.zeros(3,'d'),
'origin': np.zeros(3,'d'),
'microstructures': 0,
'homogenization': options.homogenization
}
coords = [{},{},{}]
pos = {'min':[ float("inf"), float("inf"), float("inf")],
'max':[-float("inf"),-float("inf"),-float("inf")]}
phase = []
eulerangles = []
outputAlive = True
# ------------------------------------------ process data ------------------------------------------
while outputAlive and table.data_read():
if table.data != []:
currPos = table.data[3:6]
for i in range(3):
coords[i][currPos[i]] = True
currPos = map(float,currPos)
for i in range(3):
pos['min'][i] = min(pos['min'][i],currPos[i])
pos['max'][i] = max(pos['max'][i],currPos[i])
eulerangles.append(map(math.degrees,map(float,table.data[:3])))
phase.append(options.phase[int(float(table.data[options.column-1]) > options.threshold)])
# --------------- determine size and grid ---------------------------------------------------------
info['grid'] = np.array(map(len,coords),'i')
info['size'] = info['grid']/np.maximum(np.ones(3,'d'),info['grid']-1.0)* \
np.array([pos['max'][0]-pos['min'][0],
pos['max'][1]-pos['min'][1],
pos['max'][2]-pos['min'][2]],'d')
eulerangles = np.array(eulerangles,dtype='f').reshape(info['grid'].prod(),3)
phase = np.array(phase,dtype='i').reshape(info['grid'].prod())
limits = [360,180,360]
if any([np.any(eulerangles[:,i]>=limits[i]) for i in [0,1,2]]):
damask.util.croak.write('Error: euler angles out of bound. Ang file might contain unidexed poins.\n')
for i,angle in enumerate(['phi1','PHI','phi2']):
for n in np.nditer(np.where(eulerangles[:,i]>=limits[i]),['zerosize_ok']):
damask.util.croak.write('%s in line %i (%4.2f %4.2f %4.2f)\n'
%(angle,n,eulerangles[n,0],eulerangles[n,1],eulerangles[n,2]))
continue
eulerangles=np.around(eulerangles,int(options.precision)) # round to desired precision
# ensure, that rounded euler angles are not out of bounds (modulo by limits)
for i,angle in enumerate(['phi1','PHI','phi2']):
eulerangles[:,i]%=limits[i]
# scale angles by desired precision and convert to int. create unique integer key from three euler angles by
# concatenating the string representation with leading zeros and store as integer and search unique euler angle keys.
# Texture IDs are the indices of the first occurrence, the inverse is used to construct the microstructure
# create a microstructure (texture/phase pair) for each point using unique texture IDs.
# Use longInt (64bit, i8) because the keys might be long
if options.compress:
formatString='{0:0>'+str(int(options.precision)+3)+'}'
euleranglesRadInt = (eulerangles*10**int(options.precision)).astype('int')
eulerKeys = np.array([int(''.join(map(formatString.format,euleranglesRadInt[i,:]))) \
for i in range(info['grid'].prod())])
devNull, texture, eulerKeys_idx = np.unique(eulerKeys, return_index = True, return_inverse=True)
msFull = np.array([[eulerKeys_idx[i],phase[i]] for i in range(info['grid'].prod())],'i8')
devNull,msUnique,matPoints = np.unique(msFull.view('c16'),True,True)
matPoints+=1
microstructure = np.array([msFull[i] for i in msUnique]) # pick only unique microstructures
else:
texture = np.arange(info['grid'].prod())
microstructure = np.hstack( zip(texture,phase) ).reshape(info['grid'].prod(),2) # create texture/phase pairs
formatOut = 1+int(math.log10(len(texture)))
config_header = []
formatwidth = 1+int(math.log10(len(microstructure)))
config_header += ['<microstructure>']
for i in range(len(microstructure)):
config_header += ['[Grain%s]'%str(i+1).zfill(formatwidth),
'crystallite\t%i'%options.crystallite,
'(constituent)\tphase %i\ttexture %i\tfraction 1.0'%(microstructure[i,1],microstructure[i,0]+1)
]
config_header += ['<texture>']
eulerFormatOut='%%%i.%if'%(int(options.precision)+4,int(options.precision))
outStringAngles='(gauss) phi1 '+eulerFormatOut+' Phi '+eulerFormatOut+' phi2 '+eulerFormatOut+' scatter 0.0 fraction 1.0'
for i in range(len(texture)):
config_header += ['[Texture%s]'%str(i+1).zfill(formatOut),
outStringAngles%tuple(eulerangles[texture[i],...])
]
table.labels_clear()
table.info_clear()
info['microstructures'] = len(microstructure)
#--- report ---------------------------------------------------------------------------------------
damask.util.croak('grid a b c: %s\n'%(' x '.join(map(str,info['grid']))) +
'size x y z: %s\n'%(' x '.join(map(str,info['size']))) +
'origin x y z: %s\n'%(' : '.join(map(str,info['origin']))) +
'homogenization: %i\n'%info['homogenization'] +
'microstructures: %i\n\n'%info['microstructures'])
if np.any(info['grid'] < 1):
damask.util.croak('invalid grid a b c.\n')
continue
if np.any(info['size'] <= 0.0):
damask.util.croak('invalid size x y z.\n')
continue
#--- write data -----------------------------------------------------------------------------------
table.info_append([' '.join([scriptID] + sys.argv[1:]),
"grid\ta %i\tb %i\tc %i"%(info['grid'][0],info['grid'][1],info['grid'][2],),
"size\tx %f\ty %f\tz %f"%(info['size'][0],info['size'][1],info['size'][2],),
"origin\tx %f\ty %f\tz %f"%(info['origin'][0],info['origin'][1],info['origin'][2],),
"microstructures\t%i"%info['microstructures'],
"homogenization\t%i"%info['homogenization'],
config_header
])
table.head_write()
if options.compress:
table.data = matPoints.reshape(info['grid'][1]*info['grid'][2],info['grid'][0])
table.data_writeArray('%%%ii'%(formatwidth),delimiter=' ')
else:
table.data = ["1 to %i\n"%(info['microstructures'])]
# ------------------------------------------ output finalization -----------------------------------
table.close()

52
src/CPFEM.f90
View File

@ -155,7 +155,6 @@ subroutine CPFEM_init
crystallite_Lp0, & crystallite_Lp0, &
crystallite_Fi0, & crystallite_Fi0, &
crystallite_Li0, & crystallite_Li0, &
crystallite_dPdF0, &
crystallite_Tstar0_v crystallite_Tstar0_v
implicit none implicit none
@ -207,9 +206,6 @@ subroutine CPFEM_init
read (777,rec=1) crystallite_Li0 read (777,rec=1) crystallite_Li0
close (777) close (777)
call IO_read_realFile(777,'convergeddPdF'//trim(rankStr),modelName,size(crystallite_dPdF0))
read (777,rec=1) crystallite_dPdF0
close (777)
call IO_read_realFile(777,'convergedTstar'//trim(rankStr),modelName,size(crystallite_Tstar0_v)) call IO_read_realFile(777,'convergedTstar'//trim(rankStr),modelName,size(crystallite_Tstar0_v))
read (777,rec=1) crystallite_Tstar0_v read (777,rec=1) crystallite_Tstar0_v
@ -286,12 +282,11 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
math_identity2nd, & math_identity2nd, &
math_mul33x33, & math_mul33x33, &
math_det33, & math_det33, &
math_transpose33, & math_delta, &
math_I3, & math_sym3333to66, &
math_Mandel3333to66, & math_66toSym3333, &
math_Mandel66to3333, & math_sym33to6, &
math_Mandel33to6, & math_6toSym33
math_Mandel6to33
use mesh, only: & use mesh, only: &
mesh_FEasCP, & mesh_FEasCP, &
mesh_NcpElems, & mesh_NcpElems, &
@ -326,7 +321,6 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
crystallite_Lp, & crystallite_Lp, &
crystallite_Li0, & crystallite_Li0, &
crystallite_Li, & crystallite_Li, &
crystallite_dPdF0, &
crystallite_dPdF, & crystallite_dPdF, &
crystallite_Tstar0_v, & crystallite_Tstar0_v, &
crystallite_Tstar_v crystallite_Tstar_v
@ -353,8 +347,8 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
integer(pInt), intent(in) :: mode !< computation mode 1: regular computation plus aging of results integer(pInt), intent(in) :: mode !< computation mode 1: regular computation plus aging of results
real(pReal), intent(in) :: temperature_inp !< temperature real(pReal), intent(in) :: temperature_inp !< temperature
logical, intent(in) :: parallelExecution !< flag indicating parallel computation of requested IPs logical, intent(in) :: parallelExecution !< flag indicating parallel computation of requested IPs
real(pReal), dimension(6), intent(out) :: cauchyStress !< stress vector in Mandel notation real(pReal), dimension(6), intent(out) :: cauchyStress !< stress as 6 vector
real(pReal), dimension(6,6), intent(out) :: jacobian !< jacobian in Mandel notation (Consistent tangent dcs/dE) real(pReal), dimension(6,6), intent(out) :: jacobian !< jacobian as 66 tensor (Consistent tangent dcs/dE)
real(pReal) J_inverse, & ! inverse of Jacobian real(pReal) J_inverse, & ! inverse of Jacobian
rnd rnd
@ -398,7 +392,6 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity
crystallite_Fi0 = crystallite_Fi ! crystallite intermediate deformation crystallite_Fi0 = crystallite_Fi ! crystallite intermediate deformation
crystallite_Li0 = crystallite_Li ! crystallite intermediate velocity crystallite_Li0 = crystallite_Li ! crystallite intermediate velocity
crystallite_dPdF0 = crystallite_dPdF ! crystallite stiffness
crystallite_Tstar0_v = crystallite_Tstar_v ! crystallite 2nd Piola Kirchhoff stress crystallite_Tstar0_v = crystallite_Tstar_v ! crystallite 2nd Piola Kirchhoff stress
forall ( i = 1:size(plasticState )) plasticState(i)%state0 = plasticState(i)%state ! copy state in this lenghty way because: A component cannot be an array if the encompassing structure is an array forall ( i = 1:size(plasticState )) plasticState(i)%state0 = plasticState(i)%state ! copy state in this lenghty way because: A component cannot be an array if the encompassing structure is an array
@ -454,10 +447,6 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
write (777,rec=1) crystallite_Li0 write (777,rec=1) crystallite_Li0
close (777) close (777)
call IO_write_jobRealFile(777,'convergeddPdF'//trim(rankStr),size(crystallite_dPdF0))
write (777,rec=1) crystallite_dPdF0
close (777)
call IO_write_jobRealFile(777,'convergedTstar'//trim(rankStr),size(crystallite_Tstar0_v)) call IO_write_jobRealFile(777,'convergedTstar'//trim(rankStr),size(crystallite_Tstar0_v))
write (777,rec=1) crystallite_Tstar0_v write (777,rec=1) crystallite_Tstar0_v
close (777) close (777)
@ -534,8 +523,8 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) then
write(6,'(a,1x,i8,1x,i2)') '<< CPFEM >> OUTDATED at elFE ip',elFE,ip write(6,'(a,1x,i8,1x,i2)') '<< CPFEM >> OUTDATED at elFE ip',elFE,ip
write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 old:',& write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 old:',&
math_transpose33(materialpoint_F(1:3,1:3,ip,elCP)) transpose(materialpoint_F(1:3,1:3,ip,elCP))
write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 now:',math_transpose33(ffn1) write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 now:',transpose(ffn1)
endif endif
outdatedFFN1 = .true. outdatedFFN1 = .true.
endif endif
@ -593,26 +582,25 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
endif endif
! translate from P to CS ! translate from P to CS
Kirchhoff = math_mul33x33(materialpoint_P(1:3,1:3,ip,elCP), math_transpose33(materialpoint_F(1:3,1:3,ip,elCP))) Kirchhoff = math_mul33x33(materialpoint_P(1:3,1:3,ip,elCP), transpose(materialpoint_F(1:3,1:3,ip,elCP)))
J_inverse = 1.0_pReal / math_det33(materialpoint_F(1:3,1:3,ip,elCP)) J_inverse = 1.0_pReal / math_det33(materialpoint_F(1:3,1:3,ip,elCP))
CPFEM_cs(1:6,ip,elCP) = math_Mandel33to6(J_inverse * Kirchhoff) CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.)
! translate from dP/dF to dCS/dE ! translate from dP/dF to dCS/dE
H = 0.0_pReal H = 0.0_pReal
do i=1,3; do j=1,3; do k=1,3; do l=1,3; do m=1,3; do n=1,3 do i=1,3; do j=1,3; do k=1,3; do l=1,3; do m=1,3; do n=1,3
H(i,j,k,l) = H(i,j,k,l) + & H(i,j,k,l) = H(i,j,k,l) &
materialpoint_F(j,m,ip,elCP) * & + materialpoint_F(j,m,ip,elCP) * materialpoint_F(l,n,ip,elCP) &
materialpoint_F(l,n,ip,elCP) * & * materialpoint_dPdF(i,m,k,n,ip,elCP) &
materialpoint_dPdF(i,m,k,n,ip,elCP) - & - math_delta(j,l) * materialpoint_F(i,m,ip,elCP) * materialpoint_P(k,m,ip,elCP) &
math_I3(j,l) * materialpoint_F(i,m,ip,elCP) * materialpoint_P(k,m,ip,elCP) + & + 0.5_pReal * ( Kirchhoff(j,l)*math_delta(i,k) + Kirchhoff(i,k)*math_delta(j,l) &
0.5_pReal * (math_I3(i,k) * Kirchhoff(j,l) + math_I3(j,l) * Kirchhoff(i,k) + & + Kirchhoff(j,k)*math_delta(i,l) + Kirchhoff(i,l)*math_delta(j,k))
math_I3(i,l) * Kirchhoff(j,k) + math_I3(j,k) * Kirchhoff(i,l))
enddo; enddo; enddo; enddo; enddo; enddo enddo; enddo; enddo; enddo; enddo; enddo
forall(i=1:3, j=1:3,k=1:3,l=1:3) & forall(i=1:3, j=1:3,k=1:3,l=1:3) &
H_sym(i,j,k,l) = 0.25_pReal * (H(i,j,k,l) + H(j,i,k,l) + H(i,j,l,k) + H(j,i,l,k)) H_sym(i,j,k,l) = 0.25_pReal * (H(i,j,k,l) + H(j,i,k,l) + H(i,j,l,k) + H(j,i,l,k))
CPFEM_dcsde(1:6,1:6,ip,elCP) = math_Mandel3333to66(J_inverse * H_sym) CPFEM_dcsde(1:6,1:6,ip,elCP) = math_sym3333to66(J_inverse * H_sym,weighted=.false.)
endif terminalIllness endif terminalIllness
endif validCalculation endif validCalculation
@ -639,7 +627,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
!*** remember extreme values of stress ... !*** remember extreme values of stress ...
cauchyStress33 = math_Mandel6to33(CPFEM_cs(1:6,ip,elCP)) cauchyStress33 = math_6toSym33(CPFEM_cs(1:6,ip,elCP),weighted=.false.)
if (maxval(cauchyStress33) > debug_stressMax) then if (maxval(cauchyStress33) > debug_stressMax) then
debug_stressMaxLocation = [elCP, ip] debug_stressMaxLocation = [elCP, ip]
debug_stressMax = maxval(cauchyStress33) debug_stressMax = maxval(cauchyStress33)
@ -649,7 +637,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
debug_stressMin = minval(cauchyStress33) debug_stressMin = minval(cauchyStress33)
endif endif
!*** ... and Jacobian !*** ... and Jacobian
jacobian3333 = math_Mandel66to3333(CPFEM_dcsdE(1:6,1:6,ip,elCP)) jacobian3333 = math_66toSym3333(CPFEM_dcsdE(1:6,1:6,ip,elCP),weighted=.false.)
if (maxval(jacobian3333) > debug_jacobianMax) then if (maxval(jacobian3333) > debug_jacobianMax) then
debug_jacobianMaxLocation = [elCP, ip] debug_jacobianMaxLocation = [elCP, ip]
debug_jacobianMax = maxval(jacobian3333) debug_jacobianMax = maxval(jacobian3333)

5
src/CPFEM2.f90
View File

@ -121,7 +121,6 @@ subroutine CPFEM_init
crystallite_Lp0, & crystallite_Lp0, &
crystallite_Fi0, & crystallite_Fi0, &
crystallite_Li0, & crystallite_Li0, &
crystallite_dPdF0, &
crystallite_Tstar0_v crystallite_Tstar0_v
use hdf5 use hdf5
use HDF5_utilities, only: & use HDF5_utilities, only: &
@ -160,7 +159,6 @@ subroutine CPFEM_init
call HDF5_read(fileHandle,crystallite_Fi0, 'convergedFi') call HDF5_read(fileHandle,crystallite_Fi0, 'convergedFi')
call HDF5_read(fileHandle,crystallite_Lp0, 'convergedLp') call HDF5_read(fileHandle,crystallite_Lp0, 'convergedLp')
call HDF5_read(fileHandle,crystallite_Li0, 'convergedLi') call HDF5_read(fileHandle,crystallite_Li0, 'convergedLi')
call HDF5_read(fileHandle,crystallite_dPdF0, 'convergeddPdF')
call HDF5_read(fileHandle,crystallite_Tstar0_v,'convergedTstar') call HDF5_read(fileHandle,crystallite_Tstar0_v,'convergedTstar')
groupPlasticID = HDF5_openGroup(fileHandle,'PlasticPhases') groupPlasticID = HDF5_openGroup(fileHandle,'PlasticPhases')
@ -224,7 +222,6 @@ subroutine CPFEM_age()
crystallite_Lp, & crystallite_Lp, &
crystallite_Li0, & crystallite_Li0, &
crystallite_Li, & crystallite_Li, &
crystallite_dPdF0, &
crystallite_dPdF, & crystallite_dPdF, &
crystallite_Tstar0_v, & crystallite_Tstar0_v, &
crystallite_Tstar_v crystallite_Tstar_v
@ -254,7 +251,6 @@ subroutine CPFEM_age()
crystallite_Lp0 = crystallite_Lp crystallite_Lp0 = crystallite_Lp
crystallite_Fi0 = crystallite_Fi crystallite_Fi0 = crystallite_Fi
crystallite_Li0 = crystallite_Li crystallite_Li0 = crystallite_Li
crystallite_dPdF0 = crystallite_dPdF
crystallite_Tstar0_v = crystallite_Tstar_v crystallite_Tstar0_v = crystallite_Tstar_v
forall (i = 1:size(plasticState)) plasticState(i)%state0 = plasticState(i)%state ! copy state in this lengthy way because: A component cannot be an array if the encompassing structure is an array forall (i = 1:size(plasticState)) plasticState(i)%state0 = plasticState(i)%state ! copy state in this lengthy way because: A component cannot be an array if the encompassing structure is an array
@ -283,7 +279,6 @@ subroutine CPFEM_age()
call HDF5_write(fileHandle,crystallite_Fi0, 'convergedFi') call HDF5_write(fileHandle,crystallite_Fi0, 'convergedFi')
call HDF5_write(fileHandle,crystallite_Lp0, 'convergedLp') call HDF5_write(fileHandle,crystallite_Lp0, 'convergedLp')
call HDF5_write(fileHandle,crystallite_Li0, 'convergedLi') call HDF5_write(fileHandle,crystallite_Li0, 'convergedLi')
call HDF5_write(fileHandle,crystallite_dPdF0, 'convergeddPdF')
call HDF5_write(fileHandle,crystallite_Tstar0_v,'convergedTstar') call HDF5_write(fileHandle,crystallite_Tstar0_v,'convergedTstar')
groupPlastic = HDF5_addGroup(fileHandle,'PlasticPhases') groupPlastic = HDF5_addGroup(fileHandle,'PlasticPhases')

8
src/DAMASK_abaqus.f
View File

@ -102,8 +102,6 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
calcMode, & calcMode, &
terminallyIll, & terminallyIll, &
symmetricSolver symmetricSolver
use math, only: &
invnrmMandel
use debug, only: & use debug, only: &
debug_info, & debug_info, &
debug_reset, & debug_reset, &
@ -305,9 +303,9 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
! ABAQUS implicit: 11, 22, 33, 12, 13, 23 ! ABAQUS implicit: 11, 22, 33, 12, 13, 23
! ABAQUS implicit: 11, 22, 33, 12 ! ABAQUS implicit: 11, 22, 33, 12
forall(i=1:ntens) ddsdde(1:ntens,i) = invnrmMandel(i)*ddsdde_h(1:ntens,i)*invnrmMandel(1:ntens) ddsdde = ddsdde_h(1:ntens,1:ntens)
stress(1:ntens) = stress_h(1:ntens)*invnrmMandel(1:ntens) stress = stress_h(1:ntens)
if(symmetricSolver) ddsdde(1:ntens,1:ntens) = 0.5_pReal*(ddsdde(1:ntens,1:ntens) + transpose(ddsdde(1:ntens,1:ntens))) if(symmetricSolver) ddsdde = 0.5_pReal*(ddsdde + transpose(ddsdde))
if(ntens == 6) then if(ntens == 6) then
stress_h = stress stress_h = stress
stress(5) = stress_h(6) stress(5) = stress_h(6)

11
src/DAMASK_marc.f90
View File

@ -127,9 +127,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
calcMode, & calcMode, &
terminallyIll, & terminallyIll, &
symmetricSolver symmetricSolver
use math, only: &
math_transpose33,&
invnrmMandel
use debug, only: & use debug, only: &
debug_level, & debug_level, &
debug_LEVELBASIC, & debug_LEVELBASIC, &
@ -235,9 +232,9 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
write(6,'(a,i12)') ' Nodes: ', nnode write(6,'(a,i12)') ' Nodes: ', nnode
write(6,'(a,i1)') ' Deformation gradient: ', itel write(6,'(a,i1)') ' Deformation gradient: ', itel
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n:', & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n:', &
math_transpose33(ffn) transpose(ffn)
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n+1:', & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n+1:', &
math_transpose33(ffn1) transpose(ffn1)
endif endif
!$ defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc !$ defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc
@ -357,8 +354,8 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
! Marc: 11, 22, 33, 12, 23, 13 ! Marc: 11, 22, 33, 12, 23, 13
! Marc: 11, 22, 33, 12 ! Marc: 11, 22, 33, 12
forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*ddsdde(1:ngens,i)*invnrmMandel(1:ngens) d = ddsdde(1:ngens,1:ngens)
s(1:ndi+nshear) = stress(1:ndi+nshear)*invnrmMandel(1:ndi+nshear) s = stress(1:ndi+nshear)
g = 0.0_pReal g = 0.0_pReal
if(symmetricSolver) d = 0.5_pReal*(d+transpose(d)) if(symmetricSolver) d = 0.5_pReal*(d+transpose(d))

53
src/IO.f90
View File

@ -186,11 +186,10 @@ recursive function IO_recursiveRead(fileName,cnt) result(fileContent)
fileUnit, & fileUnit, &
startPos, endPos, & startPos, endPos, &
myTotalLines, & !< # lines read from file without include statements myTotalLines, & !< # lines read from file without include statements
includedLines, & !< # lines included from other file(s)
missingLines, & !< # lines missing from current file
l,i, & l,i, &
myStat myStat
logical :: warned
if (present(cnt)) then if (present(cnt)) then
if (cnt>10_pInt) call IO_error(106_pInt,ext_msg=trim(fileName)) if (cnt>10_pInt) call IO_error(106_pInt,ext_msg=trim(fileName))
endif endif
@ -207,37 +206,39 @@ recursive function IO_recursiveRead(fileName,cnt) result(fileContent)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! count lines to allocate string array ! count lines to allocate string array
myTotalLines = 0_pInt myTotalLines = 1_pInt
do l=1_pInt, len(rawData) do l=1_pInt, len(rawData)
if (rawData(l:l) == new_line('') .or. l==len(rawData)) myTotalLines = myTotalLines+1 ! end of line or end of file without new line if (rawData(l:l) == new_line('')) myTotalLines = myTotalLines+1
enddo enddo
allocate(fileContent(myTotalLines)) allocate(fileContent(myTotalLines))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! split raw data at end of line and handle includes ! split raw data at end of line and handle includes
warned = .false.
startPos = 1_pInt startPos = 1_pInt
endPos = 0_pInt l = 1_pInt
do while (l <= myTotalLines)
endPos = merge(startPos + scan(rawData(startPos:),new_line('')) - 2_pInt,len(rawData),l /= myTotalLines)
if (endPos - startPos > 255_pInt) then
line = rawData(startPos:startPos+255_pInt)
if (.not. warned) then
call IO_warning(207_pInt,ext_msg=trim(fileName),el=l)
warned = .true.
endif
else
line = rawData(startPos:endpos)
endif
startPos = endPos + 2_pInt ! jump to next line start
includedLines=0_pInt recursion: if (scan(trim(adjustl(line)),'{') == 1 .and. scan(trim(line),'}') > 2) then
l=0_pInt
do while (startPos <= len(rawData))
l = l + 1_pInt
endPos = endPos + scan(rawData(startPos:),new_line(''))
if(endPos < startPos) endPos = len(rawData) ! end of file without end of line
if(endPos - startPos >256) call IO_error(107_pInt,ext_msg=trim(fileName))
line = rawData(startPos:endPos-1_pInt)
startPos = endPos + 1_pInt
recursion: if(scan(trim(line),'{') < scan(trim(line),'}')) then
myTotalLines = myTotalLines - 1_pInt
includedContent = IO_recursiveRead(trim(line(scan(line,'{')+1_pInt:scan(line,'}')-1_pInt)), & includedContent = IO_recursiveRead(trim(line(scan(line,'{')+1_pInt:scan(line,'}')-1_pInt)), &
merge(cnt,1_pInt,present(cnt))) ! to track recursion depth merge(cnt,1_pInt,present(cnt))) ! to track recursion depth
includedLines = includedLines + size(includedContent) fileContent = [ fileContent(1:l-1_pInt), includedContent, [(dummy,i=1,myTotalLines-l)] ] ! add content and grow array
missingLines = myTotalLines + includedLines - size(fileContent(1:l-1)) -size(includedContent) myTotalLines = myTotalLines - 1_pInt + size(includedContent)
fileContent = [ fileContent(1:l-1_pInt), includedContent, [(dummy,i=1,missingLines)] ] ! add content and grow array l = l - 1_pInt + size(includedContent)
l = l - 1_pInt + size(includedContent)
else recursion else recursion
fileContent(l) = line fileContent(l) = line
l = l + 1_pInt
endif recursion endif recursion
enddo enddo
@ -1236,6 +1237,10 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
msg = 'zero entry on stiffness diagonal' msg = 'zero entry on stiffness diagonal'
case (136_pInt) case (136_pInt)
msg = 'zero entry on stiffness diagonal for transformed phase' msg = 'zero entry on stiffness diagonal for transformed phase'
case (137_pInt)
msg = 'not defined for lattice structure'
case (138_pInt)
msg = 'not enough interaction parameters given'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! errors related to the parsing of material.config ! errors related to the parsing of material.config
@ -1494,6 +1499,8 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg)
msg = 'invalid character in string chunk' msg = 'invalid character in string chunk'
case (203_pInt) case (203_pInt)
msg = 'interpretation of string chunk failed' msg = 'interpretation of string chunk failed'
case (207_pInt)
msg = 'line truncated'
case (600_pInt) case (600_pInt)
msg = 'crystallite responds elastically' msg = 'crystallite responds elastically'
case (601_pInt) case (601_pInt)

11
src/config.f90
View File

@ -550,7 +550,7 @@ end function getString
!> @details for cumulative keys, "()", values from all occurrences are return. Otherwise only all !> @details for cumulative keys, "()", values from all occurrences are return. Otherwise only all
!! values from the last occurrence. If key is not found exits with error unless default is given. !! values from the last occurrence. If key is not found exits with error unless default is given.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function getFloats(this,key,defaultVal,requiredShape,requiredSize) function getFloats(this,key,defaultVal,requiredSize)
use IO, only: & use IO, only: &
IO_error, & IO_error, &
IO_stringValue, & IO_stringValue, &
@ -561,7 +561,6 @@ function getFloats(this,key,defaultVal,requiredShape,requiredSize)
class(tPartitionedStringList), target, intent(in) :: this class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key character(len=*), intent(in) :: key
real(pReal), dimension(:), intent(in), optional :: defaultVal real(pReal), dimension(:), intent(in), optional :: defaultVal
integer(pInt), dimension(:), intent(in), optional :: requiredShape ! not useful (is always 1D array)
integer(pInt), intent(in), optional :: requiredSize integer(pInt), intent(in), optional :: requiredSize
type(tPartitionedStringList), pointer :: item type(tPartitionedStringList), pointer :: item
integer(pInt) :: i integer(pInt) :: i
@ -601,7 +600,7 @@ end function getFloats
!> @details for cumulative keys, "()", values from all occurrences are return. Otherwise only all !> @details for cumulative keys, "()", values from all occurrences are return. Otherwise only all
!! values from the last occurrence. If key is not found exits with error unless default is given. !! values from the last occurrence. If key is not found exits with error unless default is given.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function getInts(this,key,defaultVal,requiredShape,requiredSize) function getInts(this,key,defaultVal,requiredSize)
use IO, only: & use IO, only: &
IO_error, & IO_error, &
IO_stringValue, & IO_stringValue, &
@ -611,8 +610,7 @@ function getInts(this,key,defaultVal,requiredShape,requiredSize)
integer(pInt), dimension(:), allocatable :: getInts integer(pInt), dimension(:), allocatable :: getInts
class(tPartitionedStringList), target, intent(in) :: this class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key character(len=*), intent(in) :: key
integer(pInt), dimension(:), intent(in), optional :: defaultVal, & integer(pInt), dimension(:), intent(in), optional :: defaultVal
requiredShape ! not useful (is always 1D array)
integer(pInt), intent(in), optional :: requiredSize integer(pInt), intent(in), optional :: requiredSize
type(tPartitionedStringList), pointer :: item type(tPartitionedStringList), pointer :: item
integer(pInt) :: i integer(pInt) :: i
@ -653,7 +651,7 @@ end function getInts
!! values from the last occurrence. If key is not found exits with error unless default is given. !! values from the last occurrence. If key is not found exits with error unless default is given.
!! If raw is true, the the complete string is returned, otherwise the individual chunks are returned !! If raw is true, the the complete string is returned, otherwise the individual chunks are returned
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function getStrings(this,key,defaultVal,requiredShape,raw) function getStrings(this,key,defaultVal,raw)
use IO, only: & use IO, only: &
IO_error, & IO_error, &
IO_StringValue IO_StringValue
@ -663,7 +661,6 @@ function getStrings(this,key,defaultVal,requiredShape,raw)
class(tPartitionedStringList), target, intent(in) :: this class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key character(len=*), intent(in) :: key
character(len=65536),dimension(:), intent(in), optional :: defaultVal character(len=65536),dimension(:), intent(in), optional :: defaultVal
integer(pInt), dimension(:), intent(in), optional :: requiredShape
logical, intent(in), optional :: raw logical, intent(in), optional :: raw
type(tPartitionedStringList), pointer :: item type(tPartitionedStringList), pointer :: item
character(len=65536) :: str character(len=65536) :: str

60
src/constitutive.f90
View File

@ -151,7 +151,7 @@ subroutine constitutive_init()
if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID)) call plastic_isotropic_init if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID)) call plastic_isotropic_init
if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init
if (any(phase_plasticity == PLASTICITY_KINEHARDENING_ID)) call plastic_kinehardening_init if (any(phase_plasticity == PLASTICITY_KINEHARDENING_ID)) call plastic_kinehardening_init
if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID)) call plastic_dislotwin_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID)) call plastic_dislotwin_init
if (any(phase_plasticity == PLASTICITY_DISLOUCLA_ID)) call plastic_disloucla_init if (any(phase_plasticity == PLASTICITY_DISLOUCLA_ID)) call plastic_disloucla_init
if (any(phase_plasticity == PLASTICITY_NONLOCAL_ID)) then if (any(phase_plasticity == PLASTICITY_NONLOCAL_ID)) then
call plastic_nonlocal_init(FILEUNIT) call plastic_nonlocal_init(FILEUNIT)
@ -365,7 +365,7 @@ subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el)
use plastic_nonlocal, only: & use plastic_nonlocal, only: &
plastic_nonlocal_microstructure plastic_nonlocal_microstructure
use plastic_dislotwin, only: & use plastic_dislotwin, only: &
plastic_dislotwin_microstructure plastic_dislotwin_dependentState
use plastic_disloUCLA, only: & use plastic_disloUCLA, only: &
plastic_disloUCLA_dependentState plastic_disloUCLA_dependentState
@ -389,7 +389,9 @@ subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_DISLOTWIN_ID) plasticityType case (PLASTICITY_DISLOTWIN_ID) plasticityType
call plastic_dislotwin_microstructure(temperature(ho)%p(tme),ipc,ip,el) of = phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
call plastic_dislotwin_dependentState(temperature(ho)%p(tme),instance,of)
case (PLASTICITY_DISLOUCLA_ID) plasticityType case (PLASTICITY_DISLOUCLA_ID) plasticityType
of = phasememberAt(ipc,ip,el) of = phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el)) instance = phase_plasticityInstance(material_phase(ipc,ip,el))
@ -409,9 +411,9 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
pReal pReal
use math, only: & use math, only: &
math_mul33x33, & math_mul33x33, &
math_Mandel6to33, & math_6toSym33, &
math_Mandel33to6, & math_sym33to6, &
math_Plain99to3333 math_99to3333
use material, only: & use material, only: &
phasememberAt, & phasememberAt, &
phase_plasticity, & phase_plasticity, &
@ -470,7 +472,7 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
ho = material_homogenizationAt(el) ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el) tme = thermalMapping(ho)%p(ip,el)
S = math_Mandel6to33(S6) S = math_6toSym33(S6)
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S) Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
@ -495,9 +497,9 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMp, Mp,instance,of) call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMp, Mp,instance,of)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp), & call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMp99, math_sym33to6(Mp), &
temperature(ho)%p(tme),ip,el) temperature(ho)%p(tme),ip,el)
dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget dLp_dMp = math_99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
case (PLASTICITY_DISLOTWIN_ID) plasticityType case (PLASTICITY_DISLOTWIN_ID) plasticityType
of = phasememberAt(ipc,ip,el) of = phasememberAt(ipc,ip,el)
@ -540,7 +542,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, e
math_inv33, & math_inv33, &
math_det33, & math_det33, &
math_mul33x33, & math_mul33x33, &
math_Mandel6to33 math_6toSym33
use material, only: & use material, only: &
phasememberAt, & phasememberAt, &
phase_plasticity, & phase_plasticity, &
@ -597,7 +599,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, e
case (PLASTICITY_isotropic_ID) plasticityType case (PLASTICITY_isotropic_ID) plasticityType
of = phasememberAt(ipc,ip,el) of = phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el)) instance = phase_plasticityInstance(material_phase(ipc,ip,el))
call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, math_Mandel6to33(S6),instance,of) call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, math_6toSym33(S6),instance,of)
case default plasticityType case default plasticityType
my_Li = 0.0_pReal my_Li = 0.0_pReal
my_dLi_dS = 0.0_pReal my_dLi_dS = 0.0_pReal
@ -716,7 +718,7 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip
use math, only : & use math, only : &
math_mul33x33, & math_mul33x33, &
math_mul3333xx33, & math_mul3333xx33, &
math_Mandel66to3333, & math_66toSym3333, &
math_I3 math_I3
use material, only: & use material, only: &
material_phase, & material_phase, &
@ -749,7 +751,7 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip
i, j i, j
ho = material_homogenizationAt(el) ho = material_homogenizationAt(el)
C = math_Mandel66to3333(constitutive_homogenizedC(ipc,ip,el)) C = math_66toSym3333(constitutive_homogenizedC(ipc,ip,el))
DegradationLoop: do d = 1_pInt, phase_NstiffnessDegradations(material_phase(ipc,ip,el)) DegradationLoop: do d = 1_pInt, phase_NstiffnessDegradations(material_phase(ipc,ip,el))
degradationType: select case(phase_stiffnessDegradation(d,material_phase(ipc,ip,el))) degradationType: select case(phase_stiffnessDegradation(d,material_phase(ipc,ip,el)))
@ -784,8 +786,8 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac
debug_levelBasic debug_levelBasic
use math, only: & use math, only: &
math_mul33x33, & math_mul33x33, &
math_Mandel6to33, & math_6toSym33, &
math_Mandel33to6, & math_sym33to6, &
math_mul33x33 math_mul33x33
use mesh, only: & use mesh, only: &
mesh_NcpElems, & mesh_NcpElems, &
@ -854,13 +856,13 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac
integer(pInt) :: & integer(pInt) :: &
ho, & !< homogenization ho, & !< homogenization
tme, & !< thermal member position tme, & !< thermal member position
s, & !< counter in source loop s, & !< counter in source loop
instance, of instance, of
ho = material_homogenizationAt(el) ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el) tme = thermalMapping(ho)%p(ip,el)
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6)) Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_6toSym33(S6))
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
@ -890,7 +892,7 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac
call plastic_disloucla_dotState (Mp,temperature(ho)%p(tme),instance,of) call plastic_disloucla_dotState (Mp,temperature(ho)%p(tme),instance,of)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_dotState (math_Mandel33to6(Mp),FeArray,FpArray,temperature(ho)%p(tme), & call plastic_nonlocal_dotState (math_sym33to6(Mp),FeArray,FpArray,temperature(ho)%p(tme), &
subdt,subfracArray,ip,el) subdt,subfracArray,ip,el)
end select plasticityType end select plasticityType
@ -920,7 +922,7 @@ end subroutine 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(S6, Fe, Fi, ipc, ip, el) subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
use prec, only: & use prec, only: &
pReal, & pReal, &
pLongInt pLongInt
@ -929,8 +931,7 @@ subroutine constitutive_collectDeltaState(S6, Fe, Fi, ipc, ip, el)
debug_constitutive, & debug_constitutive, &
debug_levelBasic debug_levelBasic
use math, only: & use math, only: &
math_Mandel6to33, & math_sym33to6, &
math_Mandel33to6, &
math_mul33x33 math_mul33x33
use material, only: & use material, only: &
phasememberAt, & phasememberAt, &
@ -954,18 +955,17 @@ subroutine constitutive_collectDeltaState(S6, Fe, Fi, ipc, ip, el)
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
real(pReal), intent(in), dimension(6) :: &
S6 !< 2nd Piola Kirchhoff stress (vector notation)
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
S, & !< 2nd Piola Kirchhoff stress
Fe, & !< elastic deformation gradient Fe, & !< elastic deformation gradient
Fi !< intermediate deformation gradient Fi !< intermediate deformation gradient
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
Mp Mp
integer(pInt) :: & integer(pInt) :: &
s, & !< counter in source loop i, &
instance, of instance, of
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6)) Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
@ -975,13 +975,13 @@ subroutine constitutive_collectDeltaState(S6, Fe, Fi, ipc, ip, el)
call plastic_kinehardening_deltaState(Mp,instance,of) call plastic_kinehardening_deltaState(Mp,instance,of)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_deltaState(math_Mandel33to6(Mp),ip,el) call plastic_nonlocal_deltaState(math_sym33to6(Mp),ip,el)
end select plasticityType end select plasticityType
sourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el)) sourceLoop: do i = 1_pInt, phase_Nsources(material_phase(ipc,ip,el))
sourceType: select case (phase_source(s,material_phase(ipc,ip,el))) sourceType: select case (phase_source(i,material_phase(ipc,ip,el)))
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, &
@ -1001,7 +1001,7 @@ function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el)
use prec, only: & use prec, only: &
pReal pReal
use math, only: & use math, only: &
math_Mandel6to33, & math_6toSym33, &
math_mul33x33 math_mul33x33
use mesh, only: & use mesh, only: &
mesh_NcpElems, & mesh_NcpElems, &
@ -1076,7 +1076,7 @@ function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el)
constitutive_postResults = 0.0_pReal constitutive_postResults = 0.0_pReal
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6)) Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_6toSym33(S6))
ho = material_homogenizationAt(el) ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el) tme = thermalMapping(ho)%p(ip,el)

4878
src/crystallite.f90
View File

File diff suppressed because it is too large Load Diff

180
src/homogenization.f90
View File

@ -45,10 +45,10 @@ module homogenization
materialpoint_stressAndItsTangent, & materialpoint_stressAndItsTangent, &
materialpoint_postResults materialpoint_postResults
private :: & private :: &
homogenization_partitionDeformation, & partitionDeformation, &
homogenization_updateState, & updateState, &
homogenization_averageStressAndItsTangent, & averageStressAndItsTangent, &
homogenization_postResults postResults
contains contains
@ -118,12 +118,9 @@ subroutine homogenization_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! parse homogenization from config file ! parse homogenization from config file
if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) & if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call homogenization_none_init
call homogenization_none_init() if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call homogenization_isostrain_init
if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) & if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call homogenization_RGC_init
call homogenization_isostrain_init(FILEUNIT)
if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) &
call homogenization_RGC_init(FILEUNIT)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! parse thermal from config file ! parse thermal from config file
@ -156,17 +153,14 @@ subroutine homogenization_init
select case(homogenization_type(p)) ! split per homogenization type select case(homogenization_type(p)) ! split per homogenization type
case (HOMOGENIZATION_NONE_ID) case (HOMOGENIZATION_NONE_ID)
outputName = HOMOGENIZATION_NONE_label outputName = HOMOGENIZATION_NONE_label
thisNoutput => null()
thisOutput => null() thisOutput => null()
thisSize => null() thisSize => null()
case (HOMOGENIZATION_ISOSTRAIN_ID) case (HOMOGENIZATION_ISOSTRAIN_ID)
outputName = HOMOGENIZATION_ISOSTRAIN_label outputName = HOMOGENIZATION_ISOSTRAIN_label
thisNoutput => homogenization_isostrain_Noutput thisOutput => null()
thisOutput => homogenization_isostrain_output thisSize => null()
thisSize => homogenization_isostrain_sizePostResult
case (HOMOGENIZATION_RGC_ID) case (HOMOGENIZATION_RGC_ID)
outputName = HOMOGENIZATION_RGC_label outputName = HOMOGENIZATION_RGC_label
thisNoutput => homogenization_RGC_Noutput
thisOutput => homogenization_RGC_output thisOutput => homogenization_RGC_output
thisSize => homogenization_RGC_sizePostResult thisSize => homogenization_RGC_sizePostResult
case default case default
@ -176,8 +170,9 @@ subroutine homogenization_init
if (valid) then if (valid) then
write(FILEUNIT,'(a)') '(type)'//char(9)//trim(outputName) write(FILEUNIT,'(a)') '(type)'//char(9)//trim(outputName)
write(FILEUNIT,'(a,i4)') '(ngrains)'//char(9),homogenization_Ngrains(p) write(FILEUNIT,'(a,i4)') '(ngrains)'//char(9),homogenization_Ngrains(p)
if (homogenization_type(p) /= HOMOGENIZATION_NONE_ID) then if (homogenization_type(p) /= HOMOGENIZATION_NONE_ID .and. &
do e = 1,thisNoutput(i) homogenization_type(p) /= HOMOGENIZATION_ISOSTRAIN_ID) then
do e = 1,size(thisOutput(:,i))
write(FILEUNIT,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i) write(FILEUNIT,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i)
enddo enddo
endif endif
@ -352,7 +347,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
crystallite_Li0, & crystallite_Li0, &
crystallite_Li, & crystallite_Li, &
crystallite_dPdF, & crystallite_dPdF, &
crystallite_dPdF0, &
crystallite_Tstar0_v, & crystallite_Tstar0_v, &
crystallite_Tstar_v, & crystallite_Tstar_v, &
crystallite_partionedF0, & crystallite_partionedF0, &
@ -361,12 +355,11 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
crystallite_partionedLp0, & crystallite_partionedLp0, &
crystallite_partionedFi0, & crystallite_partionedFi0, &
crystallite_partionedLi0, & crystallite_partionedLi0, &
crystallite_partioneddPdF0, &
crystallite_partionedTstar0_v, & crystallite_partionedTstar0_v, &
crystallite_dt, & crystallite_dt, &
crystallite_requested, & crystallite_requested, &
crystallite_converged, & crystallite_stress, &
crystallite_stressAndItsTangent, & crystallite_stressTangent, &
crystallite_orientations crystallite_orientations
#ifdef DEBUG #ifdef DEBUG
use debug, only: & use debug, only: &
@ -419,7 +412,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
crystallite_partionedLp0(1:3,1:3,g,i,e) = crystallite_Lp0(1:3,1:3,g,i,e) ! ...plastic velocity grads crystallite_partionedLp0(1:3,1:3,g,i,e) = crystallite_Lp0(1:3,1:3,g,i,e) ! ...plastic velocity grads
crystallite_partionedFi0(1:3,1:3,g,i,e) = crystallite_Fi0(1:3,1:3,g,i,e) ! ...intermediate def grads crystallite_partionedFi0(1:3,1:3,g,i,e) = crystallite_Fi0(1:3,1:3,g,i,e) ! ...intermediate def grads
crystallite_partionedLi0(1:3,1:3,g,i,e) = crystallite_Li0(1:3,1:3,g,i,e) ! ...intermediate velocity grads crystallite_partionedLi0(1:3,1:3,g,i,e) = crystallite_Li0(1:3,1:3,g,i,e) ! ...intermediate velocity grads
crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,g,i,e) = crystallite_dPdF0(1:3,1:3,1:3,1:3,g,i,e) ! ...stiffness
crystallite_partionedF0(1:3,1:3,g,i,e) = crystallite_F0(1:3,1:3,g,i,e) ! ...def grads crystallite_partionedF0(1:3,1:3,g,i,e) = crystallite_F0(1:3,1:3,g,i,e) ! ...def grads
crystallite_partionedTstar0_v(1:6,g,i,e) = crystallite_Tstar0_v(1:6,g,i,e) ! ...2nd PK stress crystallite_partionedTstar0_v(1:6,g,i,e) = crystallite_Tstar0_v(1:6,g,i,e) ! ...2nd PK stress
@ -489,9 +481,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
crystallite_partionedLi0(1:3,1:3,1:myNgrains,i,e) = & crystallite_partionedLi0(1:3,1:3,1:myNgrains,i,e) = &
crystallite_Li(1:3,1:3,1:myNgrains,i,e) ! ...intermediate velocity grads crystallite_Li(1:3,1:3,1:myNgrains,i,e) ! ...intermediate velocity grads
crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,1:myNgrains,i,e) = &
crystallite_dPdF(1:3,1:3,1:3,1:3,1:myNgrains,i,e) ! ...stiffness
crystallite_partionedTstar0_v(1:6,1:myNgrains,i,e) = & crystallite_partionedTstar0_v(1:6,1:myNgrains,i,e) = &
crystallite_Tstar_v(1:6,1:myNgrains,i,e) ! ...2nd PK stress crystallite_Tstar_v(1:6,1:myNgrains,i,e) ! ...2nd PK stress
@ -555,8 +544,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
crystallite_partionedFi0(1:3,1:3,1:myNgrains,i,e) ! ...intermediate def grads crystallite_partionedFi0(1:3,1:3,1:myNgrains,i,e) ! ...intermediate def grads
crystallite_Li(1:3,1:3,1:myNgrains,i,e) = & crystallite_Li(1:3,1:3,1:myNgrains,i,e) = &
crystallite_partionedLi0(1:3,1:3,1:myNgrains,i,e) ! ...intermediate velocity grads crystallite_partionedLi0(1:3,1:3,1:myNgrains,i,e) ! ...intermediate velocity grads
crystallite_dPdF(1:3,1:3,1:3,1:3,1:myNgrains,i,e) = &
crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,1:myNgrains,i,e) ! ...stiffness
crystallite_Tstar_v(1:6,1:myNgrains,i,e) = & crystallite_Tstar_v(1:6,1:myNgrains,i,e) = &
crystallite_partionedTstar0_v(1:6,1:myNgrains,i,e) ! ...2nd PK stress crystallite_partionedTstar0_v(1:6,1:myNgrains,i,e) ! ...2nd PK stress
do g = 1, myNgrains do g = 1, myNgrains
@ -613,7 +600,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
IpLooping2: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) IpLooping2: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
if ( materialpoint_requested(i,e) .and. & ! process requested but... if ( materialpoint_requested(i,e) .and. & ! process requested but...
.not. materialpoint_doneAndHappy(1,i,e)) then ! ...not yet done material points .not. materialpoint_doneAndHappy(1,i,e)) then ! ...not yet done material points
call homogenization_partitionDeformation(i,e) ! partition deformation onto constituents call partitionDeformation(i,e) ! partition deformation onto constituents
crystallite_dt(1:myNgrains,i,e) = materialpoint_subdt(i,e) ! propagate materialpoint dt to grains crystallite_dt(1:myNgrains,i,e) = materialpoint_subdt(i,e) ! propagate materialpoint dt to grains
crystallite_requested(1:myNgrains,i,e) = .true. ! request calculation for constituents crystallite_requested(1:myNgrains,i,e) = .true. ! request calculation for constituents
else else
@ -627,7 +614,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
! crystallite integration ! crystallite integration
! based on crystallite_partionedF0,.._partionedF ! based on crystallite_partionedF0,.._partionedF
! incrementing by crystallite_dt ! incrementing by crystallite_dt
call crystallite_stressAndItsTangent(updateJaco) ! request stress and tangent calculation for constituent grains materialpoint_converged = crystallite_stress() !ToDo: MD not sure if that is the best logic
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state update ! state update
@ -636,11 +623,10 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
IpLooping3: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) IpLooping3: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
if ( materialpoint_requested(i,e) .and. & if ( materialpoint_requested(i,e) .and. &
.not. materialpoint_doneAndHappy(1,i,e)) then .not. materialpoint_doneAndHappy(1,i,e)) then
if (.not. all(crystallite_converged(:,i,e))) then if (.not. materialpoint_converged(i,e)) then
materialpoint_doneAndHappy(1:2,i,e) = [.true.,.false.] materialpoint_doneAndHappy(1:2,i,e) = [.true.,.false.]
materialpoint_converged(i,e) = .false.
else else
materialpoint_doneAndHappy(1:2,i,e) = homogenization_updateState(i,e) materialpoint_doneAndHappy(1:2,i,e) = updateState(i,e)
materialpoint_converged(i,e) = all(materialpoint_doneAndHappy(1:2,i,e)) ! converged if done and happy materialpoint_converged(i,e) = all(materialpoint_doneAndHappy(1:2,i,e)) ! converged if done and happy
endif endif
endif endif
@ -653,13 +639,15 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
NiterationHomog = NiterationHomog + 1_pInt NiterationHomog = NiterationHomog + 1_pInt
enddo cutBackLooping enddo cutBackLooping
if(updateJaco) call crystallite_stressTangent
if (.not. terminallyIll ) then if (.not. terminallyIll ) then
call crystallite_orientations() ! calculate crystal orientations call crystallite_orientations() ! calculate crystal orientations
!$OMP PARALLEL DO !$OMP PARALLEL DO
elementLooping4: do e = FEsolving_execElem(1),FEsolving_execElem(2) elementLooping4: do e = FEsolving_execElem(1),FEsolving_execElem(2)
IpLooping4: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) IpLooping4: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
call homogenization_averageStressAndItsTangent(i,e) call averageStressAndItsTangent(i,e)
enddo IpLooping4 enddo IpLooping4
enddo elementLooping4 enddo elementLooping4
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
@ -717,7 +705,7 @@ subroutine materialpoint_postResults
thePos = thePos + 1_pInt thePos = thePos + 1_pInt
if (theSize > 0_pInt) then ! any homogenization results to mention? if (theSize > 0_pInt) then ! any homogenization results to mention?
materialpoint_results(thePos+1:thePos+theSize,i,e) = homogenization_postResults(i,e) ! tell homogenization results materialpoint_results(thePos+1:thePos+theSize,i,e) = postResults(i,e) ! tell homogenization results
thePos = thePos + theSize thePos = thePos + theSize
endif endif
@ -741,12 +729,12 @@ end subroutine materialpoint_postResults
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief partition material point def grad onto constituents !> @brief partition material point def grad onto constituents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine homogenization_partitionDeformation(ip,el) subroutine partitionDeformation(ip,el)
use mesh, only: & use mesh, only: &
mesh_element mesh_element
use material, only: & use material, only: &
homogenization_type, & homogenization_type, &
homogenization_maxNgrains, & homogenization_Ngrains, &
HOMOGENIZATION_NONE_ID, & HOMOGENIZATION_NONE_ID, &
HOMOGENIZATION_ISOSTRAIN_ID, & HOMOGENIZATION_ISOSTRAIN_ID, &
HOMOGENIZATION_RGC_ID HOMOGENIZATION_RGC_ID
@ -765,38 +753,36 @@ subroutine homogenization_partitionDeformation(ip,el)
chosenHomogenization: select case(homogenization_type(mesh_element(3,el))) chosenHomogenization: select case(homogenization_type(mesh_element(3,el)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization case (HOMOGENIZATION_NONE_ID) chosenHomogenization
crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el) = 0.0_pReal crystallite_partionedF(1:3,1:3,1,ip,el) = materialpoint_subF(1:3,1:3,ip,el)
crystallite_partionedF(1:3,1:3,1:1,ip,el) = &
spread(materialpoint_subF(1:3,1:3,ip,el),3,1)
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
call homogenization_isostrain_partitionDeformation(& call homogenization_isostrain_partitionDeformation(&
crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), & crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), &
materialpoint_subF(1:3,1:3,ip,el),& materialpoint_subF(1:3,1:3,ip,el))
el)
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
call homogenization_RGC_partitionDeformation(& call homogenization_RGC_partitionDeformation(&
crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), & crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), &
materialpoint_subF(1:3,1:3,ip,el),& materialpoint_subF(1:3,1:3,ip,el),&
ip, & ip, &
el) el)
end select chosenHomogenization end select chosenHomogenization
end subroutine homogenization_partitionDeformation end subroutine partitionDeformation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief update the internal state of the homogenization scheme and tell whether "done" and !> @brief update the internal state of the homogenization scheme and tell whether "done" and
!> "happy" with result !> "happy" with result
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function homogenization_updateState(ip,el) function updateState(ip,el)
use mesh, only: & use mesh, only: &
mesh_element mesh_element
use material, only: & use material, only: &
homogenization_type, & homogenization_type, &
thermal_type, & thermal_type, &
damage_type, & damage_type, &
homogenization_maxNgrains, & homogenization_Ngrains, &
HOMOGENIZATION_RGC_ID, & HOMOGENIZATION_RGC_ID, &
THERMAL_adiabatic_ID, & THERMAL_adiabatic_ID, &
DAMAGE_local_ID DAMAGE_local_ID
@ -816,27 +802,27 @@ function homogenization_updateState(ip,el)
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ip, & !< integration point ip, & !< integration point
el !< element number el !< element number
logical, dimension(2) :: homogenization_updateState logical, dimension(2) :: updateState
homogenization_updateState = .true. updateState = .true.
chosenHomogenization: select case(homogenization_type(mesh_element(3,el))) chosenHomogenization: select case(homogenization_type(mesh_element(3,el)))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
homogenization_updateState = & updateState = &
homogenization_updateState .and. & updateState .and. &
homogenization_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_maxNgrains,ip,el), & homogenization_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), &
crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), & crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), &
crystallite_partionedF0(1:3,1:3,1:homogenization_maxNgrains,ip,el),& crystallite_partionedF0(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el),&
materialpoint_subF(1:3,1:3,ip,el),& materialpoint_subF(1:3,1:3,ip,el),&
materialpoint_subdt(ip,el), & materialpoint_subdt(ip,el), &
crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,ip,el), & crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), &
ip, & ip, &
el) el)
end select chosenHomogenization end select chosenHomogenization
chosenThermal: select case (thermal_type(mesh_element(3,el))) chosenThermal: select case (thermal_type(mesh_element(3,el)))
case (THERMAL_adiabatic_ID) chosenThermal case (THERMAL_adiabatic_ID) chosenThermal
homogenization_updateState = & updateState = &
homogenization_updateState .and. & updateState .and. &
thermal_adiabatic_updateState(materialpoint_subdt(ip,el), & thermal_adiabatic_updateState(materialpoint_subdt(ip,el), &
ip, & ip, &
el) el)
@ -844,25 +830,26 @@ function homogenization_updateState(ip,el)
chosenDamage: select case (damage_type(mesh_element(3,el))) chosenDamage: select case (damage_type(mesh_element(3,el)))
case (DAMAGE_local_ID) chosenDamage case (DAMAGE_local_ID) chosenDamage
homogenization_updateState = & updateState = &
homogenization_updateState .and. & updateState .and. &
damage_local_updateState(materialpoint_subdt(ip,el), & damage_local_updateState(materialpoint_subdt(ip,el), &
ip, & ip, &
el) el)
end select chosenDamage end select chosenDamage
end function homogenization_updateState end function updateState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief derive average stress and stiffness from constituent quantities !> @brief derive average stress and stiffness from constituent quantities
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine homogenization_averageStressAndItsTangent(ip,el) subroutine averageStressAndItsTangent(ip,el)
use mesh, only: & use mesh, only: &
mesh_element mesh_element
use material, only: & use material, only: &
homogenization_type, & homogenization_type, &
homogenization_maxNgrains, & homogenization_typeInstance, &
homogenization_Ngrains, &
HOMOGENIZATION_NONE_ID, & HOMOGENIZATION_NONE_ID, &
HOMOGENIZATION_ISOSTRAIN_ID, & HOMOGENIZATION_ISOSTRAIN_ID, &
HOMOGENIZATION_RGC_ID HOMOGENIZATION_RGC_ID
@ -880,36 +867,39 @@ subroutine homogenization_averageStressAndItsTangent(ip,el)
chosenHomogenization: select case(homogenization_type(mesh_element(3,el))) chosenHomogenization: select case(homogenization_type(mesh_element(3,el)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization case (HOMOGENIZATION_NONE_ID) chosenHomogenization
materialpoint_P(1:3,1:3,ip,el) = sum(crystallite_P(1:3,1:3,1:1,ip,el),3) materialpoint_P(1:3,1:3,ip,el) = crystallite_P(1:3,1:3,1,ip,el)
materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el) & materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el) = crystallite_dPdF(1:3,1:3,1:3,1:3,1,ip,el)
= sum(crystallite_dPdF(1:3,1:3,1:3,1:3,1:1,ip,el),5)
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
call homogenization_isostrain_averageStressAndItsTangent(& call homogenization_isostrain_averageStressAndItsTangent(&
materialpoint_P(1:3,1:3,ip,el), & materialpoint_P(1:3,1:3,ip,el), &
materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),& materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),&
crystallite_P(1:3,1:3,1:homogenization_maxNgrains,ip,el), & crystallite_P(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), &
crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,ip,el), & crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), &
el) homogenization_typeInstance(mesh_element(3,el)))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
call homogenization_RGC_averageStressAndItsTangent(& call homogenization_RGC_averageStressAndItsTangent(&
materialpoint_P(1:3,1:3,ip,el), & materialpoint_P(1:3,1:3,ip,el), &
materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),& materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),&
crystallite_P(1:3,1:3,1:homogenization_maxNgrains,ip,el), & crystallite_P(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), &
crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,ip,el), & crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), &
el) homogenization_typeInstance(mesh_element(3,el)))
end select chosenHomogenization end select chosenHomogenization
end subroutine homogenization_averageStressAndItsTangent end subroutine averageStressAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return array of homogenization results for post file inclusion. call only, !> @brief return array of homogenization results for post file inclusion. call only,
!> if homogenization_sizePostResults(i,e) > 0 !! !> if homogenization_sizePostResults(i,e) > 0 !!
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function homogenization_postResults(ip,el) function postResults(ip,el)
use mesh, only: & use mesh, only: &
mesh_element mesh_element
use material, only: & use material, only: &
material_homogenizationAt, &
homogenization_typeInstance,&
mappingHomogenization, & mappingHomogenization, &
homogState, & homogState, &
thermalState, & thermalState, &
@ -926,8 +916,6 @@ function homogenization_postResults(ip,el)
DAMAGE_none_ID, & DAMAGE_none_ID, &
DAMAGE_local_ID, & DAMAGE_local_ID, &
DAMAGE_nonlocal_ID DAMAGE_nonlocal_ID
use homogenization_isostrain, only: &
homogenization_isostrain_postResults
use homogenization_RGC, only: & use homogenization_RGC, only: &
homogenization_RGC_postResults homogenization_RGC_postResults
use thermal_adiabatic, only: & use thermal_adiabatic, only: &
@ -946,60 +934,46 @@ function homogenization_postResults(ip,el)
real(pReal), dimension( homogState (mappingHomogenization(2,ip,el))%sizePostResults & real(pReal), dimension( homogState (mappingHomogenization(2,ip,el))%sizePostResults &
+ thermalState (mappingHomogenization(2,ip,el))%sizePostResults & + thermalState (mappingHomogenization(2,ip,el))%sizePostResults &
+ damageState (mappingHomogenization(2,ip,el))%sizePostResults) :: & + damageState (mappingHomogenization(2,ip,el))%sizePostResults) :: &
homogenization_postResults postResults
integer(pInt) :: & integer(pInt) :: &
startPos, endPos startPos, endPos ,&
of, instance
homogenization_postResults = 0.0_pReal
postResults = 0.0_pReal
startPos = 1_pInt startPos = 1_pInt
endPos = homogState(mappingHomogenization(2,ip,el))%sizePostResults endPos = homogState(mappingHomogenization(2,ip,el))%sizePostResults
chosenHomogenization: select case (homogenization_type(mesh_element(3,el))) chosenHomogenization: select case (homogenization_type(mesh_element(3,el)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
homogenization_postResults(startPos:endPos) = &
homogenization_isostrain_postResults(&
ip, &
el, &
materialpoint_P(1:3,1:3,ip,el), &
materialpoint_F(1:3,1:3,ip,el))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
homogenization_postResults(startPos:endPos) = & instance = homogenization_typeInstance(material_homogenizationAt(el))
homogenization_RGC_postResults(& of = mappingHomogenization(1,ip,el)
ip, & postResults(startPos:endPos) = homogenization_RGC_postResults(instance,of)
el, &
materialpoint_P(1:3,1:3,ip,el), &
materialpoint_F(1:3,1:3,ip,el))
end select chosenHomogenization end select chosenHomogenization
startPos = endPos + 1_pInt startPos = endPos + 1_pInt
endPos = endPos + thermalState(mappingHomogenization(2,ip,el))%sizePostResults endPos = endPos + thermalState(mappingHomogenization(2,ip,el))%sizePostResults
chosenThermal: select case (thermal_type(mesh_element(3,el))) chosenThermal: select case (thermal_type(mesh_element(3,el)))
case (THERMAL_isothermal_ID) chosenThermal
case (THERMAL_adiabatic_ID) chosenThermal case (THERMAL_adiabatic_ID) chosenThermal
homogenization_postResults(startPos:endPos) = & postResults(startPos:endPos) = thermal_adiabatic_postResults(ip, el)
thermal_adiabatic_postResults(ip, el)
case (THERMAL_conduction_ID) chosenThermal case (THERMAL_conduction_ID) chosenThermal
homogenization_postResults(startPos:endPos) = & postResults(startPos:endPos) = thermal_conduction_postResults(ip, el)
thermal_conduction_postResults(ip, el)
end select chosenThermal end select chosenThermal
startPos = endPos + 1_pInt startPos = endPos + 1_pInt
endPos = endPos + damageState(mappingHomogenization(2,ip,el))%sizePostResults endPos = endPos + damageState(mappingHomogenization(2,ip,el))%sizePostResults
chosenDamage: select case (damage_type(mesh_element(3,el))) chosenDamage: select case (damage_type(mesh_element(3,el)))
case (DAMAGE_none_ID) chosenDamage
case (DAMAGE_local_ID) chosenDamage case (DAMAGE_local_ID) chosenDamage
homogenization_postResults(startPos:endPos) = & postResults(startPos:endPos) = damage_local_postResults(ip, el)
damage_local_postResults(ip, el)
case (DAMAGE_nonlocal_ID) chosenDamage case (DAMAGE_nonlocal_ID) chosenDamage
homogenization_postResults(startPos:endPos) = & postResults(startPos:endPos) = damage_nonlocal_postResults(ip, el)
damage_nonlocal_postResults(ip, el)
end select chosenDamage end select chosenDamage
end function homogenization_postResults end function postResults
end module homogenization end module homogenization

1789
src/homogenization_RGC.f90
View File

File diff suppressed because it is too large Load Diff

321
src/homogenization_isostrain.f90
View File

@ -1,4 +1,5 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Isostrain (full constraint Taylor assuption) homogenization scheme !> @brief Isostrain (full constraint Taylor assuption) homogenization scheme
@ -6,220 +7,119 @@
module homogenization_isostrain module homogenization_isostrain
use prec, only: & use prec, only: &
pInt pInt
implicit none implicit none
private private
integer(pInt), dimension(:), allocatable, public, protected :: &
homogenization_isostrain_sizePostResults
integer(pInt), dimension(:,:), allocatable, target, public :: &
homogenization_isostrain_sizePostResult
character(len=64), dimension(:,:), allocatable, target, public :: &
homogenization_isostrain_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
homogenization_isostrain_Noutput !< number of outputs per homog instance
integer(pInt), dimension(:), allocatable, private :: &
homogenization_isostrain_Ngrains
enum, bind(c) enum, bind(c)
enumerator :: undefined_ID, & enumerator :: &
nconstituents_ID, & parallel_ID, &
ipcoords_ID, & average_ID
avgdefgrad_ID, &
avgfirstpiola_ID
end enum end enum
enum, bind(c)
enumerator :: parallel_ID, &
average_ID
end enum
integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: &
homogenization_isostrain_outputID !< ID of each post result output
integer(kind(average_ID)), dimension(:), allocatable, private :: &
homogenization_isostrain_mapping !< mapping type
type, private :: tParameters !< container type for internal constitutive parameters
integer(pInt) :: &
Nconstituents
integer(kind(average_ID)) :: &
mapping
end type
type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance)
public :: & public :: &
homogenization_isostrain_init, & homogenization_isostrain_init, &
homogenization_isostrain_partitionDeformation, & homogenization_isostrain_partitionDeformation, &
homogenization_isostrain_averageStressAndItsTangent, & homogenization_isostrain_averageStressAndItsTangent
homogenization_isostrain_postResults
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields, reads information from material configuration file !> @brief allocates all neccessary fields, reads information from material configuration file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine homogenization_isostrain_init(fileUnit) subroutine homogenization_isostrain_init()
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 #if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: & use, intrinsic :: iso_fortran_env, only: &
compiler_version, & compiler_version, &
compiler_options compiler_options
#endif #endif
use prec, only: &
pReal
use debug, only: & use debug, only: &
debug_HOMOGENIZATION, & debug_HOMOGENIZATION, &
debug_level, & debug_level, &
debug_levelBasic debug_levelBasic
use IO use IO, only: &
use material IO_timeStamp, &
use config IO_error
use material, only: &
homogenization_type, &
material_homog, &
homogState, &
HOMOGENIZATION_ISOSTRAIN_ID, &
HOMOGENIZATION_ISOSTRAIN_LABEL, &
homogenization_typeInstance
use config, only: &
config_homogenization
implicit none implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: & integer(pInt) :: &
section = 0_pInt, i, mySize, o Ninstance, &
integer :: & h, &
maxNinstance, & NofMyHomog
homog, &
instance
integer :: &
NofMyHomog ! no pInt (stores a system dependen value from 'count'
character(len=65536) :: & character(len=65536) :: &
tag = '', & tag = ''
line = ''
write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_ISOSTRAIN_label//' init -+>>>' write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_ISOSTRAIN_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
maxNinstance = count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID) Ninstance = int(count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID),pInt)
if (maxNinstance == 0) return
if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0_pInt) & if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(homogenization_isostrain_sizePostResults(maxNinstance), source=0_pInt)
allocate(homogenization_isostrain_sizePostResult(maxval(homogenization_Noutput),maxNinstance), & allocate(param(Ninstance)) ! one container of parameters per instance
source=0_pInt)
allocate(homogenization_isostrain_Noutput(maxNinstance), source=0_pInt) do h = 1_pInt, size(homogenization_type)
allocate(homogenization_isostrain_Ngrains(maxNinstance), source=0_pInt) if (homogenization_type(h) /= HOMOGENIZATION_ISOSTRAIN_ID) cycle
allocate(homogenization_isostrain_mapping(maxNinstance), source=average_ID)
allocate(homogenization_isostrain_output(maxval(homogenization_Noutput),maxNinstance)) associate(prm => param(homogenization_typeInstance(h)),&
homogenization_isostrain_output = '' config => config_homogenization(h))
allocate(homogenization_isostrain_outputID(maxval(homogenization_Noutput),maxNinstance), &
source=undefined_ID) prm%Nconstituents = config_homogenization(h)%getInt('nconstituents')
tag = 'sum'
select case(trim(config%getString('mapping',defaultVal = tag)))
case ('sum')
prm%mapping = parallel_ID
case ('avg')
prm%mapping = average_ID
case default
call IO_error(211_pInt,ext_msg=trim(tag)//' ('//HOMOGENIZATION_isostrain_label//')')
end select
NofMyHomog = count(material_homog == h)
homogState(h)%sizeState = 0_pInt
homogState(h)%sizePostResults = 0_pInt
allocate(homogState(h)%state0 (0_pInt,NofMyHomog))
allocate(homogState(h)%subState0(0_pInt,NofMyHomog))
allocate(homogState(h)%state (0_pInt,NofMyHomog))
end associate
rewind(fileUnit)
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization)! wind forward to <homogenization>
line = IO_read(fileUnit)
enddo enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of homogenization part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1_pInt
cycle
endif
if (section > 0_pInt ) then ! do not short-circuit here (.and. with next if-statement). It's not safe in Fortran
if (homogenization_type(section) == HOMOGENIZATION_ISOSTRAIN_ID) then ! one of my sections
i = homogenization_typeInstance(section) ! which instance of my type is present homogenization
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('(output)')
select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
case('nconstituents','ngrains')
homogenization_isostrain_Noutput(i) = homogenization_isostrain_Noutput(i) + 1_pInt
homogenization_isostrain_outputID(homogenization_isostrain_Noutput(i),i) = nconstituents_ID
homogenization_isostrain_output(homogenization_isostrain_Noutput(i),i) = &
IO_lc(IO_stringValue(line,chunkPos,2_pInt))
case('ipcoords')
homogenization_isostrain_Noutput(i) = homogenization_isostrain_Noutput(i) + 1_pInt
homogenization_isostrain_outputID(homogenization_isostrain_Noutput(i),i) = ipcoords_ID
homogenization_isostrain_output(homogenization_isostrain_Noutput(i),i) = &
IO_lc(IO_stringValue(line,chunkPos,2_pInt))
case('avgdefgrad','avgf')
homogenization_isostrain_Noutput(i) = homogenization_isostrain_Noutput(i) + 1_pInt
homogenization_isostrain_outputID(homogenization_isostrain_Noutput(i),i) = avgdefgrad_ID
homogenization_isostrain_output(homogenization_isostrain_Noutput(i),i) = &
IO_lc(IO_stringValue(line,chunkPos,2_pInt))
case('avgp','avgfirstpiola','avg1stpiola')
homogenization_isostrain_Noutput(i) = homogenization_isostrain_Noutput(i) + 1_pInt
homogenization_isostrain_outputID(homogenization_isostrain_Noutput(i),i) = avgfirstpiola_ID
homogenization_isostrain_output(homogenization_isostrain_Noutput(i),i) = &
IO_lc(IO_stringValue(line,chunkPos,2_pInt))
end select
case ('nconstituents','ngrains')
homogenization_isostrain_Ngrains(i) = IO_intValue(line,chunkPos,2_pInt)
case ('mapping')
select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
case ('parallel','sum')
homogenization_isostrain_mapping(i) = parallel_ID
case ('average','mean','avg')
homogenization_isostrain_mapping(i) = average_ID
case default
call IO_error(211_pInt,ext_msg=trim(tag)//' ('//HOMOGENIZATION_isostrain_label//')')
end select
end select
endif
endif
enddo parsingFile
initializeInstances: do homog = 1_pInt, material_Nhomogenization
myHomog: if (homogenization_type(homog) == HOMOGENIZATION_ISOSTRAIN_ID) then
NofMyHomog = count(material_homog == homog)
instance = homogenization_typeInstance(homog)
! * Determine size of postResults array
outputsLoop: do o = 1_pInt, homogenization_isostrain_Noutput(instance)
select case(homogenization_isostrain_outputID(o,instance))
case(nconstituents_ID)
mySize = 1_pInt
case(ipcoords_ID)
mySize = 3_pInt
case(avgdefgrad_ID, avgfirstpiola_ID)
mySize = 9_pInt
case default
mySize = 0_pInt
end select
outputFound: if (mySize > 0_pInt) then
homogenization_isostrain_sizePostResult(o,instance) = mySize
homogenization_isostrain_sizePostResults(instance) = &
homogenization_isostrain_sizePostResults(instance) + mySize
endif outputFound
enddo outputsLoop
! allocate state arrays
homogState(homog)%sizeState = 0_pInt
homogState(homog)%sizePostResults = homogenization_isostrain_sizePostResults(instance)
allocate(homogState(homog)%state0 (0_pInt,NofMyHomog), source=0.0_pReal)
allocate(homogState(homog)%subState0(0_pInt,NofMyHomog), source=0.0_pReal)
allocate(homogState(homog)%state (0_pInt,NofMyHomog), source=0.0_pReal)
endif myHomog
enddo initializeInstances
end subroutine homogenization_isostrain_init end subroutine homogenization_isostrain_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief partitions the deformation gradient onto the constituents !> @brief partitions the deformation gradient onto the constituents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine homogenization_isostrain_partitionDeformation(F,avgF,el) subroutine homogenization_isostrain_partitionDeformation(F,avgF)
use prec, only: & use prec, only: &
pReal pReal
use mesh, only: &
mesh_element
use material, only: &
homogenization_maxNgrains, &
homogenization_Ngrains
implicit none implicit none
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F !< partioned def grad per grain real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient
real(pReal), dimension (3,3), intent(in) :: avgF !< my average def grad
integer(pInt), intent(in) :: & real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point
el !< element number
F=0.0_pReal F = spread(avgF,3,size(F,3))
F(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)))= &
spread(avgF,3,homogenization_Ngrains(mesh_element(3,el)))
end subroutine homogenization_isostrain_partitionDeformation end subroutine homogenization_isostrain_partitionDeformation
@ -227,90 +127,31 @@ end subroutine homogenization_isostrain_partitionDeformation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief derive average stress and stiffness from constituent quantities !> @brief derive average stress and stiffness from constituent quantities
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine homogenization_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,el) subroutine homogenization_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance)
use prec, only: & use prec, only: &
pReal pReal
use mesh, only: &
mesh_element
use material, only: &
homogenization_maxNgrains, &
homogenization_Ngrains, &
homogenization_typeInstance
implicit none implicit none
real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: P !< array of current grain stresses
real(pReal), dimension (3,3,3,3,homogenization_maxNgrains), intent(in) :: dPdF !< array of current grain stiffnesses
integer(pInt), intent(in) :: el !< element number
integer(pInt) :: &
homID, &
Ngrains
homID = homogenization_typeInstance(mesh_element(3,el)) real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses
Ngrains = homogenization_Ngrains(mesh_element(3,el)) real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
integer(pInt), intent(in) :: instance
select case (homogenization_isostrain_mapping(homID)) associate(prm => param(instance))
select case (prm%mapping)
case (parallel_ID) case (parallel_ID)
avgP = sum(P,3) avgP = sum(P,3)
dAvgPdAvgF = sum(dPdF,5) dAvgPdAvgF = sum(dPdF,5)
case (average_ID) case (average_ID)
avgP = sum(P,3) /real(Ngrains,pReal) avgP = sum(P,3) /real(prm%Nconstituents,pReal)
dAvgPdAvgF = sum(dPdF,5)/real(Ngrains,pReal) dAvgPdAvgF = sum(dPdF,5)/real(prm%Nconstituents,pReal)
end select end select
end associate
end subroutine homogenization_isostrain_averageStressAndItsTangent end subroutine homogenization_isostrain_averageStressAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief return array of homogenization results for post file inclusion
!--------------------------------------------------------------------------------------------------
pure function homogenization_isostrain_postResults(ip,el,avgP,avgF)
use prec, only: &
pReal
use mesh, only: &
mesh_element, &
mesh_ipCoordinates
use material, only: &
homogenization_typeInstance, &
homogenization_Noutput
implicit none
integer(pInt), intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), dimension(3,3), intent(in) :: &
avgP, & !< average stress at material point
avgF !< average deformation gradient at material point
real(pReal), dimension(homogenization_isostrain_sizePostResults &
(homogenization_typeInstance(mesh_element(3,el)))) :: &
homogenization_isostrain_postResults
integer(pInt) :: &
homID, &
o, c
c = 0_pInt
homID = homogenization_typeInstance(mesh_element(3,el))
homogenization_isostrain_postResults = 0.0_pReal
do o = 1_pInt,homogenization_Noutput(mesh_element(3,el))
select case(homogenization_isostrain_outputID(o,homID))
case (nconstituents_ID)
homogenization_isostrain_postResults(c+1_pInt) = real(homogenization_isostrain_Ngrains(homID),pReal)
c = c + 1_pInt
case (avgdefgrad_ID)
homogenization_isostrain_postResults(c+1_pInt:c+9_pInt) = reshape(transpose(avgF),[9])
c = c + 9_pInt
case (avgfirstpiola_ID)
homogenization_isostrain_postResults(c+1_pInt:c+9_pInt) = reshape(transpose(avgP),[9])
c = c + 9_pInt
case (ipcoords_ID)
homogenization_isostrain_postResults(c+1_pInt:c+3_pInt) = mesh_ipCoordinates(1:3,ip,el) ! current ip coordinates
c = c + 3_pInt
end select
enddo
end function homogenization_isostrain_postResults
end module homogenization_isostrain end module homogenization_isostrain

47
src/homogenization_none.f90
View File

@ -2,7 +2,7 @@
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief dummy homogenization homogenization scheme !> @brief dummy homogenization homogenization scheme for 1 constituent per material point
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module homogenization_none module homogenization_none
@ -24,35 +24,46 @@ subroutine homogenization_none_init()
compiler_options compiler_options
#endif #endif
use prec, only: & use prec, only: &
pReal, & pInt
pInt use debug, only: &
debug_HOMOGENIZATION, &
debug_level, &
debug_levelBasic
use IO, only: & use IO, only: &
IO_timeStamp IO_timeStamp
use material
use config use material, only: &
homogenization_type, &
material_homog, &
homogState, &
HOMOGENIZATION_NONE_LABEL, &
HOMOGENIZATION_NONE_ID
implicit none implicit none
integer(pInt) :: & integer(pInt) :: &
homog, & Ninstance, &
h, &
NofMyHomog NofMyHomog
write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_NONE_label//' init -+>>>' write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_NONE_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
initializeInstances: do homog = 1_pInt, material_Nhomogenization Ninstance = int(count(homogenization_type == HOMOGENIZATION_NONE_ID),pInt)
if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
do h = 1_pInt, size(homogenization_type)
if (homogenization_type(h) /= HOMOGENIZATION_NONE_ID) cycle
myhomog: if (homogenization_type(homog) == HOMOGENIZATION_none_ID) then NofMyHomog = count(material_homog == h)
NofMyHomog = count(material_homog == homog) homogState(h)%sizeState = 0_pInt
homogState(homog)%sizeState = 0_pInt homogState(h)%sizePostResults = 0_pInt
homogState(homog)%sizePostResults = 0_pInt allocate(homogState(h)%state0 (0_pInt,NofMyHomog))
allocate(homogState(homog)%state0 (0_pInt,NofMyHomog), source=0.0_pReal) allocate(homogState(h)%subState0(0_pInt,NofMyHomog))
allocate(homogState(homog)%subState0(0_pInt,NofMyHomog), source=0.0_pReal) allocate(homogState(h)%state (0_pInt,NofMyHomog))
allocate(homogState(homog)%state (0_pInt,NofMyHomog), source=0.0_pReal)
endif myhomog
enddo initializeInstances
enddo
end subroutine homogenization_none_init end subroutine homogenization_none_init

2371
src/lattice.f90
View File

File diff suppressed because it is too large Load Diff

16
src/material.f90
View File

@ -327,19 +327,19 @@ subroutine material_init()
#include "compilation_info.f90" #include "compilation_info.f90"
call material_parsePhase() call material_parsePhase()
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Phase parsed'; flush(6) if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Phase parsed'; flush(6)
call material_parseMicrostructure() call material_parseMicrostructure()
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Microstructure parsed'; flush(6) if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Microstructure parsed'; flush(6)
call material_parseCrystallite() call material_parseCrystallite()
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Crystallite parsed'; flush(6) if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Crystallite parsed'; flush(6)
call material_parseHomogenization() call material_parseHomogenization()
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Homogenization parsed'; flush(6) if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Homogenization parsed'; flush(6)
call material_parseTexture() call material_parseTexture()
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Texture parsed'; flush(6) if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Texture parsed'; flush(6)
allocate(plasticState (size(config_phase))) allocate(plasticState (size(config_phase)))
allocate(sourceState (size(config_phase))) allocate(sourceState (size(config_phase)))
@ -918,7 +918,8 @@ end subroutine material_parseTexture
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates the plastic state of a phase !> @brief allocates the plastic state of a phase
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine material_allocatePlasticState(phase,NofMyPhase,sizeState,sizeDotState,sizeDeltaState,& subroutine material_allocatePlasticState(phase,NofMyPhase,&
sizeState,sizeDotState,sizeDeltaState,&
Nslip,Ntwin,Ntrans) Nslip,Ntwin,Ntrans)
use numerics, only: & use numerics, only: &
numerics_integrator2 => numerics_integrator ! compatibility hack numerics_integrator2 => numerics_integrator ! compatibility hack
@ -936,9 +937,10 @@ subroutine material_allocatePlasticState(phase,NofMyPhase,sizeState,sizeDotState
integer(pInt) :: numerics_integrator ! compatibility hack integer(pInt) :: numerics_integrator ! compatibility hack
numerics_integrator = numerics_integrator2(1) ! compatibility hack numerics_integrator = numerics_integrator2(1) ! compatibility hack
plasticState(phase)%sizeState = sizeState plasticState(phase)%sizeState = sizeState
plasticState(phase)%sizeDotState = sizeDotState plasticState(phase)%sizeDotState = sizeDotState
plasticState(phase)%sizeDeltaState = sizeDeltaState plasticState(phase)%sizeDeltaState = sizeDeltaState
plasticState(phase)%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition
plasticState(phase)%Nslip = Nslip plasticState(phase)%Nslip = Nslip
plasticState(phase)%Ntwin = Ntwin plasticState(phase)%Ntwin = Ntwin
plasticState(phase)%Ntrans= Ntrans plasticState(phase)%Ntrans= Ntrans

461
src/math.f90
View File

@ -24,33 +24,25 @@ module math
0.0_pReal,0.0_pReal,1.0_pReal & 0.0_pReal,0.0_pReal,1.0_pReal &
],[3,3]) !< 3x3 Identity ],[3,3]) !< 3x3 Identity
! ToDo MD: Our naming scheme is a little bit odd: We use essentially the re-ordering according to Nye real(pReal), dimension(6), parameter, private :: &
! (convenient because Abaqus and Marc want to have 12 on position 4) nrmMandel = [&
! but weight the shear components according to Mandel (convenient for matrix multiplications) 1.0_pReal, 1.0_pReal, 1.0_pReal, &
! I suggest to keep Voigt3333to66 (required for reading in elasticity matrices) but rename sqrt(2.0_pReal), sqrt(2.0_pReal), sqrt(2.0_pReal) ] !< weighting for Mandel notation (forward)
! mapMandel to mapNye, math_MandelXtoY to math_XtoY and math_PlainXtoY to math_XtoY.
! It is then clear that math_33to9 just reorders and math_33to6 does the "DAMASK conversion" real(pReal), dimension(6), parameter , private :: &
! without leaving the impression that it follows any established convention invnrmMandel = [&
1.0_pReal, 1.0_pReal, 1.0_pReal, &
1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal) ] !< weighting for Mandel notation (backward)
integer(pInt), dimension (2,6), parameter, private :: & integer(pInt), dimension (2,6), parameter, private :: &
mapMandel = reshape([& mapNye = reshape([&
1_pInt,1_pInt, & 1_pInt,1_pInt, &
2_pInt,2_pInt, & 2_pInt,2_pInt, &
3_pInt,3_pInt, & 3_pInt,3_pInt, &
1_pInt,2_pInt, & 1_pInt,2_pInt, &
2_pInt,3_pInt, & 2_pInt,3_pInt, &
1_pInt,3_pInt & 1_pInt,3_pInt &
],[2,6]) !< arrangement in Mandel notation. Differs from https://en.wikipedia.org/wiki/Voigt_notation#Mandel_notation ],[2,6]) !< arrangement in Nye notation.
real(pReal), dimension(6), parameter, private :: &
nrmMandel = [&
1.0_pReal, 1.0_pReal, 1.0_pReal, &
sqrt(2.0_pReal), sqrt(2.0_pReal), sqrt(2.0_pReal) ] !< weighting for Mandel notation (forward)
real(pReal), dimension(6), parameter , public :: &
invnrmMandel = [&
1.0_pReal, 1.0_pReal, 1.0_pReal, &
1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal) ] !< weighting for Mandel notation (backward)
integer(pInt), dimension (2,6), parameter, private :: & integer(pInt), dimension (2,6), parameter, private :: &
mapVoigt = reshape([& mapVoigt = reshape([&
@ -62,10 +54,6 @@ module math
1_pInt,2_pInt & 1_pInt,2_pInt &
],[2,6]) !< arrangement in Voigt notation ],[2,6]) !< arrangement in Voigt notation
real(pReal), dimension(6), parameter, private :: &
nrmVoigt = 1.0_pReal, & !< weighting for Voigt notation (forward)
invnrmVoigt = 1.0_pReal !< weighting for Voigt notation (backward)
integer(pInt), dimension (2,9), parameter, private :: & integer(pInt), dimension (2,9), parameter, private :: &
mapPlain = reshape([& mapPlain = reshape([&
1_pInt,1_pInt, & 1_pInt,1_pInt, &
@ -78,6 +66,56 @@ module math
3_pInt,2_pInt, & 3_pInt,2_pInt, &
3_pInt,3_pInt & 3_pInt,3_pInt &
],[2,9]) !< arrangement in Plain notation ],[2,9]) !< arrangement in Plain notation
!--------------------------------------------------------------------------------------------------
! Provide deprecated names for compatibility
! ToDo MD: Our naming scheme was a little bit odd: We use essentially the re-ordering according to Nye
! (convenient because Abaqus and Marc want to have 12 on position 4)
! but weight the shear components according to Mandel (convenient for matrix multiplications)
interface math_Plain33to9
module procedure math_33to9
end interface math_Plain33to9
interface math_Plain9to33
module procedure math_9to33
end interface math_Plain9to33
interface math_Mandel33to6
module procedure math_sym33to6
end interface math_Mandel33to6
interface math_Mandel6to33
module procedure math_6toSym33
end interface math_Mandel6to33
interface math_Plain3333to99
module procedure math_3333to99
end interface math_Plain3333to99
interface math_Plain99to3333
module procedure math_99to3333
end interface math_Plain99to3333
interface math_Mandel3333to66
module procedure math_sym3333to66
end interface math_Mandel3333to66
interface math_Mandel66to3333
module procedure math_66toSym3333
end interface math_Mandel66to3333
public :: &
math_Plain33to9, &
math_Plain9to33, &
math_Mandel33to6, &
math_Mandel6to33, &
math_Plain3333to99, &
math_Plain99to3333, &
math_Mandel3333to66, &
math_Mandel66to3333
!---------------------------------------------------------------------------------------------------
public :: & public :: &
math_init, & math_init, &
@ -107,6 +145,7 @@ module math
math_invert33, & math_invert33, &
math_invSym3333, & math_invSym3333, &
math_invert, & math_invert, &
math_invert2, &
math_symmetric33, & math_symmetric33, &
math_symmetric66, & math_symmetric66, &
math_skew33, & math_skew33, &
@ -116,16 +155,14 @@ module math
math_equivStress33, & math_equivStress33, &
math_trace33, & math_trace33, &
math_det33, & math_det33, &
math_Plain33to9, & math_33to9, &
math_Plain9to33, & math_9to33, &
math_Mandel33to6, & math_sym33to6, &
math_Mandel6to33, & math_6toSym33, &
math_Plain3333to99, & math_3333to99, &
math_Plain99to3333, & math_99to3333, &
math_Mandel66toPlain66, & math_sym3333to66, &
math_Plain66toMandel66, & math_66toSym3333, &
math_Mandel3333to66, &
math_Mandel66to3333, &
math_Voigt66to3333, & math_Voigt66to3333, &
math_qRand, & math_qRand, &
math_qMul, & math_qMul, &
@ -423,7 +460,7 @@ pure function math_identity2nd(dimen)
real(pReal), dimension(dimen,dimen) :: math_identity2nd real(pReal), dimension(dimen,dimen) :: math_identity2nd
math_identity2nd = 0.0_pReal math_identity2nd = 0.0_pReal
forall (i=1_pInt:dimen) math_identity2nd(i,i) = 1.0_pReal forall(i=1_pInt:dimen) math_identity2nd(i,i) = 1.0_pReal
end function math_identity2nd end function math_identity2nd
@ -437,9 +474,11 @@ pure function math_identity4th(dimen)
integer(pInt), intent(in) :: dimen !< tensor dimension integer(pInt), intent(in) :: dimen !< tensor dimension
integer(pInt) :: i,j,k,l integer(pInt) :: i,j,k,l
real(pReal), dimension(dimen,dimen,dimen,dimen) :: math_identity4th real(pReal), dimension(dimen,dimen,dimen,dimen) :: math_identity4th
real(pReal), dimension(dimen,dimen) :: identity2nd
forall (i=1_pInt:dimen,j=1_pInt:dimen,k=1_pInt:dimen,l=1_pInt:dimen) math_identity4th(i,j,k,l) = & identity2nd = math_identity2nd(dimen)
0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k)) forall(i=1_pInt:dimen,j=1_pInt:dimen,k=1_pInt:dimen,l=1_pInt:dimen) &
math_identity4th(i,j,k,l) = 0.5_pReal*(identity2nd(i,k)*identity2nd(j,l)+identity2nd(i,l)*identity2nd(j,k))
end function math_identity4th end function math_identity4th
@ -508,7 +547,7 @@ pure function math_tensorproduct(A,B)
real(pReal), dimension(size(A,1),size(B,1)) :: math_tensorproduct real(pReal), dimension(size(A,1),size(B,1)) :: math_tensorproduct
integer(pInt) :: i,j integer(pInt) :: i,j
forall (i=1_pInt:size(A,1),j=1_pInt:size(B,1)) math_tensorproduct(i,j) = A(i)*B(j) forall(i=1_pInt:size(A,1),j=1_pInt:size(B,1)) math_tensorproduct(i,j) = A(i)*B(j)
end function math_tensorproduct end function math_tensorproduct
@ -523,7 +562,7 @@ pure function math_tensorproduct33(A,B)
real(pReal), dimension(3), intent(in) :: A,B real(pReal), dimension(3), intent(in) :: A,B
integer(pInt) :: i,j integer(pInt) :: i,j
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) math_tensorproduct33(i,j) = A(i)*B(j) forall(i=1_pInt:3_pInt,j=1_pInt:3_pInt) math_tensorproduct33(i,j) = A(i)*B(j)
end function math_tensorproduct33 end function math_tensorproduct33
@ -564,7 +603,7 @@ real(pReal) pure function math_mul33xx33(A,B)
integer(pInt) :: i,j integer(pInt) :: i,j
real(pReal), dimension(3,3) :: C real(pReal), dimension(3,3) :: C
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) C(i,j) = A(i,j) * B(i,j) forall(i=1_pInt:3_pInt,j=1_pInt:3_pInt) C(i,j) = A(i,j) * B(i,j)
math_mul33xx33 = sum(C) math_mul33xx33 = sum(C)
end function math_mul33xx33 end function math_mul33xx33
@ -581,9 +620,8 @@ pure function math_mul3333xx33(A,B)
real(pReal), dimension(3,3), intent(in) :: B real(pReal), dimension(3,3), intent(in) :: B
integer(pInt) :: i,j integer(pInt) :: i,j
forall(i = 1_pInt:3_pInt,j = 1_pInt:3_pInt) & forall(i = 1_pInt:3_pInt,j = 1_pInt:3_pInt) math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
end function math_mul3333xx33 end function math_mul3333xx33
@ -614,8 +652,7 @@ pure function math_mul33x33(A,B)
real(pReal), dimension(3,3), intent(in) :: A,B real(pReal), dimension(3,3), intent(in) :: A,B
integer(pInt) :: i,j integer(pInt) :: i,j
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) & forall(i=1_pInt:3_pInt,j=1_pInt:3_pInt) math_mul33x33(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j)
math_mul33x33(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j)
end function math_mul33x33 end function math_mul33x33
@ -630,9 +667,9 @@ pure function math_mul66x66(A,B)
real(pReal), dimension(6,6), intent(in) :: A,B real(pReal), dimension(6,6), intent(in) :: A,B
integer(pInt) :: i,j integer(pInt) :: i,j
forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt) math_mul66x66(i,j) = & forall(i=1_pInt:6_pInt,j=1_pInt:6_pInt) &
A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) + & math_mul66x66(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) &
A(i,4)*B(4,j) + A(i,5)*B(5,j) + A(i,6)*B(6,j) + A(i,4)*B(4,j) + A(i,5)*B(5,j) + A(i,6)*B(6,j)
end function math_mul66x66 end function math_mul66x66
@ -647,10 +684,10 @@ pure function math_mul99x99(A,B)
real(pReal), dimension(9,9), intent(in) :: A,B real(pReal), dimension(9,9), intent(in) :: A,B
integer(pInt) i,j integer(pInt) i,j
forall (i=1_pInt:9_pInt,j=1_pInt:9_pInt) math_mul99x99(i,j) = & forall(i=1_pInt:9_pInt,j=1_pInt:9_pInt) &
A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) + & math_mul99x99(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) &
A(i,4)*B(4,j) + A(i,5)*B(5,j) + A(i,6)*B(6,j) + & + A(i,4)*B(4,j) + A(i,5)*B(5,j) + A(i,6)*B(6,j) &
A(i,7)*B(7,j) + A(i,8)*B(8,j) + A(i,9)*B(9,j) + A(i,7)*B(7,j) + A(i,8)*B(8,j) + A(i,9)*B(9,j)
end function math_mul99x99 end function math_mul99x99
@ -698,9 +735,8 @@ pure function math_mul66x6(A,B)
real(pReal), dimension(6), intent(in) :: B real(pReal), dimension(6), intent(in) :: B
integer(pInt) :: i integer(pInt) :: i
forall (i=1_pInt:6_pInt) math_mul66x6(i) = & forall (i=1_pInt:6_pInt) math_mul66x6(i) = A(i,1)*B(1) + A(i,2)*B(2) + A(i,3)*B(3) &
A(i,1)*B(1) + A(i,2)*B(2) + A(i,3)*B(3) + & + A(i,4)*B(4) + A(i,5)*B(5) + A(i,6)*B(6)
A(i,4)*B(4) + A(i,5)*B(5) + A(i,6)*B(6)
end function math_mul66x6 end function math_mul66x6
@ -747,8 +783,8 @@ end function math_transpose33
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Cramer inversion of 33 matrix (function) !> @brief Cramer inversion of 33 matrix (function)
! direct Cramer inversion of matrix A. !> @details Direct Cramer inversion of matrix A. Returns all zeroes if not possible, i.e.
! returns all zeroes if not possible, i.e. if det close to zero ! if determinant is close to zero
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function math_inv33(A) pure function math_inv33(A)
use prec, only: & use prec, only: &
@ -784,9 +820,9 @@ end function math_inv33
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Cramer inversion of 33 matrix (subroutine) !> @brief Cramer inversion of 33 matrix (subroutine)
! direct Cramer inversion of matrix A. !> @details Direct Cramer inversion of matrix A. Also returns determinant
! also returns determinant ! Returns an error if not possible, i.e. if determinant is close to zero
! returns error if not possible, i.e. if det close to zero ! ToDo: Output arguments should be first
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine math_invert33(A, InvA, DetA, error) pure subroutine math_invert33(A, InvA, DetA, error)
use prec, only: & use prec, only: &
@ -843,11 +879,11 @@ function math_invSym3333(A)
dgetrf, & dgetrf, &
dgetri dgetri
temp66_real = math_Mandel3333to66(A) temp66_real = math_sym3333to66(A)
call dgetrf(6,6,temp66_real,6,ipiv6,ierr) call dgetrf(6,6,temp66_real,6,ipiv6,ierr)
call dgetri(6,temp66_real,6,ipiv6,work6,6,ierr) call dgetri(6,temp66_real,6,ipiv6,work6,6,ierr)
if (ierr == 0_pInt) then if (ierr == 0_pInt) then
math_invSym3333 = math_Mandel66to3333(temp66_real) math_invSym3333 = math_66toSym3333(temp66_real)
else else
call IO_error(400_pInt, ext_msg = 'math_invSym3333') call IO_error(400_pInt, ext_msg = 'math_invSym3333')
endif endif
@ -855,8 +891,27 @@ function math_invSym3333(A)
end function math_invSym3333 end function math_invSym3333
!--------------------------------------------------------------------------------------------------
!> @brief invert quadratic matrix of arbitrary dimension
! ToDo: replaces math_invert
!--------------------------------------------------------------------------------------------------
subroutine math_invert2(InvA, error, A)
implicit none
real(pReal), dimension(:,:), intent(in) :: A
real(pReal), dimension(size(A,1),size(A,1)), intent(out) :: invA
logical, intent(out) :: error
call math_invert(size(A,1), A, InvA, error)
end subroutine math_invert2
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief invert matrix of arbitrary dimension !> @brief invert matrix of arbitrary dimension
! ToDo: Wrong order of arguments and superfluous myDim argument.
! Use math_invert2 instead
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine math_invert(myDim,A, InvA, error) subroutine math_invert(myDim,A, InvA, error)
@ -961,15 +1016,14 @@ pure function math_equivStrain33(m)
real(pReal), dimension(3,3), intent(in) :: m real(pReal), dimension(3,3), intent(in) :: m
real(pReal), dimension(3) :: e,s real(pReal), dimension(3) :: e,s
real(pReal) :: math_equivStrain33 real(pReal) :: math_equivStrain33
real(pReal), parameter :: TWOTHIRD = 2.0_pReal/3.0_pReal
e = [2.0_pReal*m(1,1)-m(2,2)-m(3,3), & e = [2.0_pReal*m(1,1)-m(2,2)-m(3,3), &
2.0_pReal*m(2,2)-m(3,3)-m(1,1), & 2.0_pReal*m(2,2)-m(3,3)-m(1,1), &
2.0_pReal*m(3,3)-m(1,1)-m(2,2)]/3.0_pReal 2.0_pReal*m(3,3)-m(1,1)-m(2,2)]/3.0_pReal
s = [m(1,2),m(2,3),m(1,3)]*2.0_pReal s = [m(1,2),m(2,3),m(1,3)]*2.0_pReal
math_equivStrain33 = TWOTHIRD*(1.50_pReal*(sum(e**2.0_pReal)) + & math_equivStrain33 = 2.0_pReal/3.0_pReal &
0.75_pReal*(sum(s**2.0_pReal)))**(0.5_pReal) * (1.50_pReal*(sum(e**2.0_pReal))+ 0.75_pReal*(sum(s**2.0_pReal)))**(0.5_pReal)
end function math_equivStrain33 end function math_equivStrain33
@ -1041,172 +1095,188 @@ end function math_detSym33
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief convert 33 matrix into vector 9 !> @brief convert 33 matrix into vector 9
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function math_Plain33to9(m33) pure function math_33to9(m33)
implicit none implicit none
real(pReal), dimension(9) :: math_Plain33to9 real(pReal), dimension(9) :: math_33to9
real(pReal), dimension(3,3), intent(in) :: m33
integer(pInt) :: i
forall (i=1_pInt:9_pInt) math_Plain33to9(i) = m33(mapPlain(1,i),mapPlain(2,i))
end function math_Plain33to9
!--------------------------------------------------------------------------------------------------
!> @brief convert Plain 9 back to 33 matrix
!--------------------------------------------------------------------------------------------------
pure function math_Plain9to33(v9)
implicit none
real(pReal), dimension(3,3) :: math_Plain9to33
real(pReal), dimension(9), intent(in) :: v9
integer(pInt) :: i
forall (i=1_pInt:9_pInt) math_Plain9to33(mapPlain(1,i),mapPlain(2,i)) = v9(i)
end function math_Plain9to33
!--------------------------------------------------------------------------------------------------
!> @brief convert symmetric 33 matrix into Mandel vector 6
!--------------------------------------------------------------------------------------------------
pure function math_Mandel33to6(m33)
implicit none
real(pReal), dimension(6) :: math_Mandel33to6
real(pReal), dimension(3,3), intent(in) :: m33 real(pReal), dimension(3,3), intent(in) :: m33
integer(pInt) :: i integer(pInt) :: i
forall (i=1_pInt:6_pInt) math_Mandel33to6(i) = nrmMandel(i)*m33(mapMandel(1,i),mapMandel(2,i)) forall(i=1_pInt:9_pInt) math_33to9(i) = m33(mapPlain(1,i),mapPlain(2,i))
end function math_Mandel33to6 end function math_33to9
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief convert Mandel 6 back to symmetric 33 matrix !> @brief convert 9 vector into 33 matrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function math_Mandel6to33(v6) pure function math_9to33(v9)
implicit none implicit none
real(pReal), dimension(6), intent(in) :: v6 real(pReal), dimension(3,3) :: math_9to33
real(pReal), dimension(3,3) :: math_Mandel6to33 real(pReal), dimension(9), intent(in) :: v9
integer(pInt) :: i integer(pInt) :: i
forall (i=1_pInt:6_pInt) forall(i=1_pInt:9_pInt) math_9to33(mapPlain(1,i),mapPlain(2,i)) = v9(i)
math_Mandel6to33(mapMandel(1,i),mapMandel(2,i)) = invnrmMandel(i)*v6(i)
math_Mandel6to33(mapMandel(2,i),mapMandel(1,i)) = invnrmMandel(i)*v6(i)
end forall
end function math_Mandel6to33 end function math_9to33
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief convert 3333 tensor into plain matrix 99 !> @brief convert symmetric 33 matrix into 6 vector
!> @details Weighted conversion (default) rearranges according to Nye and weights shear
! components according to Mandel. Advisable for matrix operations.
! Unweighted conversion only changes order according to Nye
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function math_Plain3333to99(m3333) pure function math_sym33to6(m33,weighted)
implicit none implicit none
real(pReal), dimension(6) :: math_sym33to6
real(pReal), dimension(3,3), intent(in) :: m33
logical, optional, intent(in) :: weighted
real(pReal), dimension(6) :: w
integer(pInt) :: i
if(present(weighted)) then
w = merge(nrmMandel,1.0_pReal,weighted)
else
w = nrmMandel
endif
forall(i=1_pInt:6_pInt) math_sym33to6(i) = w(i)*m33(mapNye(1,i),mapNye(2,i))
end function math_sym33to6
!--------------------------------------------------------------------------------------------------
!> @brief convert 6 vector into symmetric 33 matrix
!> @details Weighted conversion (default) rearranges according to Nye and weights shear
! components according to Mandel. Advisable for matrix operations.
! Unweighted conversion only changes order according to Nye
!--------------------------------------------------------------------------------------------------
pure function math_6toSym33(v6,weighted)
implicit none
real(pReal), dimension(3,3) :: math_6toSym33
real(pReal), dimension(6), intent(in) :: v6
logical, optional, intent(in) :: weighted
real(pReal), dimension(6) :: w
integer(pInt) :: i
if(present(weighted)) then
w = merge(invnrmMandel,1.0_pReal,weighted)
else
w = invnrmMandel
endif
do i=1_pInt,6_pInt
math_6toSym33(mapNye(1,i),mapNye(2,i)) = w(i)*v6(i)
math_6toSym33(mapNye(2,i),mapNye(1,i)) = w(i)*v6(i)
enddo
end function math_6toSym33
!--------------------------------------------------------------------------------------------------
!> @brief convert 3333 matrix into 99 matrix
!--------------------------------------------------------------------------------------------------
pure function math_3333to99(m3333)
implicit none
real(pReal), dimension(9,9) :: math_3333to99
real(pReal), dimension(3,3,3,3), intent(in) :: m3333 real(pReal), dimension(3,3,3,3), intent(in) :: m3333
real(pReal), dimension(9,9) :: math_Plain3333to99
integer(pInt) :: i,j integer(pInt) :: i,j
forall (i=1_pInt:9_pInt,j=1_pInt:9_pInt) math_Plain3333to99(i,j) = & forall(i=1_pInt:9_pInt,j=1_pInt:9_pInt) &
m3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j)) math_3333to99(i,j) = m3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j))
end function math_3333to99
end function math_Plain3333to99
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief plain matrix 99 into 3333 tensor !> @brief convert 99 matrix into 3333 matrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function math_Plain99to3333(m99) pure function math_99to3333(m99)
implicit none implicit none
real(pReal), dimension(3,3,3,3) :: math_99to3333
real(pReal), dimension(9,9), intent(in) :: m99 real(pReal), dimension(9,9), intent(in) :: m99
real(pReal), dimension(3,3,3,3) :: math_Plain99to3333
integer(pInt) :: i,j integer(pInt) :: i,j
forall (i=1_pInt:9_pInt,j=1_pInt:9_pInt) math_Plain99to3333(mapPlain(1,i),mapPlain(2,i),& forall(i=1_pInt:9_pInt,j=1_pInt:9_pInt) &
mapPlain(1,j),mapPlain(2,j)) = m99(i,j) math_99to3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j)) = m99(i,j)
end function math_Plain99to3333 end function math_99to3333
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief convert Mandel matrix 66 into Plain matrix 66 !> @brief convert symmetric 3333 matrix into 66 matrix
!> @details Weighted conversion (default) rearranges according to Nye and weights shear
! components according to Mandel. Advisable for matrix operations.
! Unweighted conversion only changes order according to Nye
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function math_Mandel66toPlain66(m66) pure function math_sym3333to66(m3333,weighted)
implicit none
real(pReal), dimension(6,6), intent(in) :: m66
real(pReal), dimension(6,6) :: math_Mandel66toPlain66
integer(pInt) :: i,j
forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt) &
math_Mandel66toPlain66(i,j) = invnrmMandel(i) * invnrmMandel(j) * m66(i,j)
end function math_Mandel66toPlain66
!--------------------------------------------------------------------------------------------------
!> @brief convert Plain matrix 66 into Mandel matrix 66
!--------------------------------------------------------------------------------------------------
pure function math_Plain66toMandel66(m66)
implicit none
real(pReal), dimension(6,6), intent(in) :: m66
real(pReal), dimension(6,6) :: math_Plain66toMandel66
integer(pInt) :: i,j
forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt) &
math_Plain66toMandel66(i,j) = nrmMandel(i) * nrmMandel(j) * m66(i,j)
end function math_Plain66toMandel66
!--------------------------------------------------------------------------------------------------
!> @brief convert symmetric 3333 tensor into Mandel matrix 66
!--------------------------------------------------------------------------------------------------
pure function math_Mandel3333to66(m3333)
implicit none implicit none
real(pReal), dimension(6,6) :: math_sym3333to66
real(pReal), dimension(3,3,3,3), intent(in) :: m3333 real(pReal), dimension(3,3,3,3), intent(in) :: m3333
real(pReal), dimension(6,6) :: math_Mandel3333to66 logical, optional, intent(in) :: weighted
real(pReal), dimension(6) :: w
integer(pInt) :: i,j integer(pInt) :: i,j
if(present(weighted)) then
w = merge(nrmMandel,1.0_pReal,weighted)
else
w = nrmMandel
endif
forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt) math_Mandel3333to66(i,j) = & forall(i=1_pInt:6_pInt,j=1_pInt:6_pInt) &
nrmMandel(i)*nrmMandel(j)*m3333(mapMandel(1,i),mapMandel(2,i),mapMandel(1,j),mapMandel(2,j)) math_sym3333to66(i,j) = w(i)*w(j)*m3333(mapNye(1,i),mapNye(2,i),mapNye(1,j),mapNye(2,j))
end function math_Mandel3333to66 end function math_sym3333to66
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief convert Mandel matrix 66 back to symmetric 3333 tensor !> @brief convert 66 matrix into symmetric 3333 matrix
!> @details Weighted conversion (default) rearranges according to Nye and weights shear
! components according to Mandel. Advisable for matrix operations.
! Unweighted conversion only changes order according to Nye
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function math_Mandel66to3333(m66) pure function math_66toSym3333(m66,weighted)
implicit none implicit none
real(pReal), dimension(3,3,3,3) :: math_Mandel66to3333 real(pReal), dimension(3,3,3,3) :: math_66toSym3333
real(pReal), dimension(6,6), intent(in) :: m66 real(pReal), dimension(6,6), intent(in) :: m66
logical, optional, intent(in) :: weighted
real(pReal), dimension(6) :: w
integer(pInt) :: i,j integer(pInt) :: i,j
if(present(weighted)) then
w = merge(invnrmMandel,1.0_pReal,weighted)
else
w = invnrmMandel
endif
forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt) do i=1_pInt,6_pInt; do j=1_pInt, 6_pInt
math_Mandel66to3333(mapMandel(1,i),mapMandel(2,i),mapMandel(1,j),mapMandel(2,j)) = & math_66toSym3333(mapNye(1,i),mapNye(2,i),mapNye(1,j),mapNye(2,j)) = w(i)*w(j)*m66(i,j)
invnrmMandel(i)*invnrmMandel(j)*m66(i,j) math_66toSym3333(mapNye(2,i),mapNye(1,i),mapNye(1,j),mapNye(2,j)) = w(i)*w(j)*m66(i,j)
math_Mandel66to3333(mapMandel(2,i),mapMandel(1,i),mapMandel(1,j),mapMandel(2,j)) = & math_66toSym3333(mapNye(1,i),mapNye(2,i),mapNye(2,j),mapNye(1,j)) = w(i)*w(j)*m66(i,j)
invnrmMandel(i)*invnrmMandel(j)*m66(i,j) math_66toSym3333(mapNye(2,i),mapNye(1,i),mapNye(2,j),mapNye(1,j)) = w(i)*w(j)*m66(i,j)
math_Mandel66to3333(mapMandel(1,i),mapMandel(2,i),mapMandel(2,j),mapMandel(1,j)) = & enddo; enddo
invnrmMandel(i)*invnrmMandel(j)*m66(i,j)
math_Mandel66to3333(mapMandel(2,i),mapMandel(1,i),mapMandel(2,j),mapMandel(1,j)) = &
invnrmMandel(i)*invnrmMandel(j)*m66(i,j)
end forall
end function math_Mandel66to3333 end function math_66toSym3333
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief convert Voigt matrix 66 back to symmetric 3333 tensor !> @brief convert 66 Voigt matrix into symmetric 3333 matrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function math_Voigt66to3333(m66) pure function math_Voigt66to3333(m66)
@ -1215,16 +1285,12 @@ pure function math_Voigt66to3333(m66)
real(pReal), dimension(6,6), intent(in) :: m66 real(pReal), dimension(6,6), intent(in) :: m66
integer(pInt) :: i,j integer(pInt) :: i,j
forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt) do i=1_pInt,6_pInt; do j=1_pInt, 6_pInt
math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(1,j),mapVoigt(2,j)) = & math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(1,j),mapVoigt(2,j)) = m66(i,j)
invnrmVoigt(i)*invnrmVoigt(j)*m66(i,j) math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(1,j),mapVoigt(2,j)) = m66(i,j)
math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(1,j),mapVoigt(2,j)) = & math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(2,j),mapVoigt(1,j)) = m66(i,j)
invnrmVoigt(i)*invnrmVoigt(j)*m66(i,j) math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(2,j),mapVoigt(1,j)) = m66(i,j)
math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(2,j),mapVoigt(1,j)) = & enddo; enddo
invnrmVoigt(i)*invnrmVoigt(j)*m66(i,j)
math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(2,j),mapVoigt(1,j)) = &
invnrmVoigt(i)*invnrmVoigt(j)*m66(i,j)
end forall
end function math_Voigt66to3333 end function math_Voigt66to3333
@ -1632,8 +1698,7 @@ pure function math_qToR(q)
real(pReal), dimension(3,3) :: math_qToR, T,S real(pReal), dimension(3,3) :: math_qToR, T,S
integer(pInt) :: i, j integer(pInt) :: i, j
forall (i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) & forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) T(i,j) = q(i+1_pInt) * q(j+1_pInt)
T(i,j) = q(i+1_pInt) * q(j+1_pInt)
S = reshape( [0.0_pReal, -q(4), q(3), & S = reshape( [0.0_pReal, -q(4), q(3), &
q(4), 0.0_pReal, -q(2), & q(4), 0.0_pReal, -q(2), &
@ -1933,6 +1998,7 @@ end function math_symmetricEulers
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief eigenvalues and eigenvectors of symmetric matrix m !> @brief eigenvalues and eigenvectors of symmetric matrix m
! ToDo: has wrong oder of arguments
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine math_eigenValuesVectorsSym(m,values,vectors,error) subroutine math_eigenValuesVectorsSym(m,values,vectors,error)
@ -1956,9 +2022,10 @@ end subroutine math_eigenValuesVectorsSym
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief eigenvalues and eigenvectors of symmetric 33 matrix m using an analytical expression !> @brief eigenvalues and eigenvectors of symmetric 33 matrix m using an analytical expression
!> and the general LAPACK powered version for arbritrary sized matrices as fallback !> and the general LAPACK powered version for arbritrary sized matrices as fallback
!> @author Joachim Kopp, Max–Planck–Institut für Kernphysik, Heidelberg (Copyright (C) 2006) !> @author Joachim Kopp, Max-Planck-Institut für Kernphysik, Heidelberg (Copyright (C) 2006)
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @details See http://arxiv.org/abs/physics/0610206 (DSYEVH3) !> @details See http://arxiv.org/abs/physics/0610206 (DSYEVH3)
! ToDo: has wrong oder of arguments
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine math_eigenValuesVectorsSym33(m,values,vectors) subroutine math_eigenValuesVectorsSym33(m,values,vectors)
@ -2038,7 +2105,7 @@ end function math_eigenvectorBasisSym
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief eigenvector basis of symmetric 33 matrix m !> @brief eigenvector basis of symmetric 33 matrix m
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function math_eigenvectorBasisSym33(m) pure function math_eigenvectorBasisSym33(m)
implicit none implicit none
real(pReal), dimension(3,3) :: math_eigenvectorBasisSym33 real(pReal), dimension(3,3) :: math_eigenvectorBasisSym33
@ -2103,7 +2170,7 @@ end function math_eigenvectorBasisSym33
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief logarithm eigenvector basis of symmetric 33 matrix m !> @brief logarithm eigenvector basis of symmetric 33 matrix m
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function math_eigenvectorBasisSym33_log(m) pure function math_eigenvectorBasisSym33_log(m)
implicit none implicit none
real(pReal), dimension(3,3) :: math_eigenvectorBasisSym33_log real(pReal), dimension(3,3) :: math_eigenvectorBasisSym33_log
@ -2159,11 +2226,12 @@ function math_eigenvectorBasisSym33_log(m)
endif threeSimilarEigenvalues endif threeSimilarEigenvalues
math_eigenvectorBasisSym33_log = log(sqrt(values(1))) * EB(1:3,1:3,1) & 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(2))) * EB(1:3,1:3,2) &
+ log(sqrt(values(3))) * EB(1:3,1:3,3) + log(sqrt(values(3))) * EB(1:3,1:3,3)
end function math_eigenvectorBasisSym33_log end function math_eigenvectorBasisSym33_log
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief rotational part from polar decomposition of 33 tensor m !> @brief rotational part from polar decomposition of 33 tensor m
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -2608,13 +2676,12 @@ pure function math_rotate_forward3333(tensor,rot_tensor)
real(pReal), dimension(3,3,3,3), intent(in) :: tensor real(pReal), dimension(3,3,3,3), intent(in) :: tensor
integer(pInt) :: i,j,k,l,m,n,o,p integer(pInt) :: i,j,k,l,m,n,o,p
math_rotate_forward3333= 0.0_pReal math_rotate_forward3333 = 0.0_pReal
do i = 1_pInt,3_pInt;do j = 1_pInt,3_pInt;do k = 1_pInt,3_pInt;do l = 1_pInt,3_pInt
do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt; do k = 1_pInt,3_pInt; do l = 1_pInt,3_pInt do m = 1_pInt,3_pInt;do n = 1_pInt,3_pInt;do o = 1_pInt,3_pInt;do p = 1_pInt,3_pInt
do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt; do o = 1_pInt,3_pInt; do p = 1_pInt,3_pInt math_rotate_forward3333(i,j,k,l) &
math_rotate_forward3333(i,j,k,l) = math_rotate_forward3333(i,j,k,l) & = math_rotate_forward3333(i,j,k,l) &
+ rot_tensor(i,m) * rot_tensor(j,n) & + rot_tensor(i,m) * rot_tensor(j,n) * rot_tensor(k,o) * rot_tensor(l,p) * tensor(m,n,o,p)
* rot_tensor(k,o) * rot_tensor(l,p) * tensor(m,n,o,p)
enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo
end function math_rotate_forward3333 end function math_rotate_forward3333

5
src/numerics.f90
View File

@ -276,8 +276,6 @@ subroutine numerics_init
numerics_integrator = IO_intValue(line,chunkPos,2_pInt) numerics_integrator = IO_intValue(line,chunkPos,2_pInt)
case ('usepingpong') case ('usepingpong')
usepingpong = IO_intValue(line,chunkPos,2_pInt) > 0_pInt usepingpong = IO_intValue(line,chunkPos,2_pInt) > 0_pInt
case ('timesyncing')
numerics_timeSyncing = IO_intValue(line,chunkPos,2_pInt) > 0_pInt
case ('unitlength') case ('unitlength')
numerics_unitlength = IO_floatValue(line,chunkPos,2_pInt) numerics_unitlength = IO_floatValue(line,chunkPos,2_pInt)
@ -454,8 +452,6 @@ subroutine numerics_init
end select end select
#endif #endif
numerics_timeSyncing = numerics_timeSyncing .and. all(numerics_integrator==2_pInt) ! timeSyncing only allowed for explicit Euler integrator
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! writing parameters to output ! writing parameters to output
write(6,'(a24,1x,es8.1)') ' relevantStrain: ',relevantStrain write(6,'(a24,1x,es8.1)') ' relevantStrain: ',relevantStrain
@ -476,7 +472,6 @@ subroutine numerics_init
write(6,'(a24,1x,es8.1)') ' rTol_crystalliteStress: ',rTol_crystalliteStress write(6,'(a24,1x,es8.1)') ' rTol_crystalliteStress: ',rTol_crystalliteStress
write(6,'(a24,1x,es8.1)') ' aTol_crystalliteStress: ',aTol_crystalliteStress write(6,'(a24,1x,es8.1)') ' aTol_crystalliteStress: ',aTol_crystalliteStress
write(6,'(a24,2(1x,i8))') ' integrator: ',numerics_integrator write(6,'(a24,2(1x,i8))') ' integrator: ',numerics_integrator
write(6,'(a24,1x,L8)') ' timeSyncing: ',numerics_timeSyncing
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

96
src/plastic_disloUCLA.f90
View File

@ -28,8 +28,7 @@ module plastic_disloUCLA
shearrate_ID, & shearrate_ID, &
accumulatedshear_ID, & accumulatedshear_ID, &
mfp_ID, & mfp_ID, &
thresholdstress_ID, & thresholdstress_ID
dipoledistance_ID
end enum end enum
type, private :: tParameters type, private :: tParameters
@ -73,7 +72,7 @@ module plastic_disloUCLA
integer(kind(undefined_ID)), allocatable, dimension(:) :: & integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID !< ID of each post result output outputID !< ID of each post result output
logical :: & logical :: &
dipoleformation dipoleFormation !< flag indicating consideration of dipole formation
end type !< container type for internal constitutive parameters end type !< container type for internal constitutive parameters
type, private :: tDisloUCLAState type, private :: tDisloUCLAState
@ -93,7 +92,7 @@ module plastic_disloUCLA
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! containers for parameters and state ! containers for parameters and state
type(tParameters), allocatable, dimension(:), private :: param type(tParameters), allocatable, dimension(:), private :: param
type(tDisloUCLAState ), allocatable, dimension(:), private :: & type(tDisloUCLAState), allocatable, dimension(:), private :: &
dotState, & dotState, &
state state
type(tDisloUCLAdependentState), allocatable, dimension(:), private :: dependentState type(tDisloUCLAdependentState), allocatable, dimension(:), private :: dependentState
@ -127,7 +126,6 @@ subroutine plastic_disloUCLA_init()
debug_constitutive,& debug_constitutive,&
debug_levelBasic debug_levelBasic
use math, only: & use math, only: &
math_mul3x3, &
math_expand math_expand
use IO, only: & use IO, only: &
IO_error, & IO_error, &
@ -148,8 +146,6 @@ subroutine plastic_disloUCLA_init()
implicit none implicit none
integer(pInt) :: & integer(pInt) :: &
index_myFamily, index_otherFamily, &
f,j,k,o, &
Ninstance, & Ninstance, &
p, i, & p, i, &
NipcMyPhase, & NipcMyPhase, &
@ -164,7 +160,6 @@ subroutine plastic_disloUCLA_init()
outputID outputID
character(len=pStringLen) :: & character(len=pStringLen) :: &
structure = '',&
extmsg = '' extmsg = ''
character(len=65536), dimension(:), allocatable :: & character(len=65536), dimension(:), allocatable :: &
outputs outputs
@ -197,8 +192,6 @@ subroutine plastic_disloUCLA_init()
dst => dependentState(phase_plasticityInstance(p)), & dst => dependentState(phase_plasticityInstance(p)), &
config => config_phase(p)) config => config_phase(p))
structure = config%getString('lattice_structure')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! optional parameters that need to be defined ! optional parameters that need to be defined
prm%mu = lattice_mu(p) prm%mu = lattice_mu(p)
@ -213,36 +206,41 @@ subroutine plastic_disloUCLA_init()
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%totalNslip = sum(prm%Nslip) prm%totalNslip = sum(prm%Nslip)
slipActive: if (prm%totalNslip > 0_pInt) then slipActive: if (prm%totalNslip > 0_pInt) then
prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,structure(1:3),& prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal)) config%getFloat('c/a',defaultVal=0.0_pReal))
if(structure=='bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',& if(trim(config%getString('lattice_structure')) == 'bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
defaultVal = emptyRealArray) defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt) prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt) prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt)
else else
prm%nonSchmid_pos = prm%Schmid prm%nonSchmid_pos = prm%Schmid
prm%nonSchmid_neg = prm%Schmid prm%nonSchmid_neg = prm%Schmid
endif endif
prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, & prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, &
config%getFloats('interaction_slipslip'), & config%getFloats('interaction_slipslip'), &
structure(1:3)) config%getString('lattice_structure'))
prm%rho0 = config%getFloats('rhoedge0', requiredShape=shape(prm%Nslip)) prm%forestProjectionEdge = lattice_forestProjection(prm%Nslip,config%getString('lattice_structure'),&
prm%rhoDip0 = config%getFloats('rhoedgedip0', requiredShape=shape(prm%Nslip)) config%getFloat('c/a',defaultVal=0.0_pReal))
prm%v0 = config%getFloats('v0', requiredShape=shape(prm%Nslip))
prm%burgers = config%getFloats('slipburgers', requiredShape=shape(prm%Nslip))
prm%H0kp = config%getFloats('qedge', requiredShape=shape(prm%Nslip))
prm%clambda = config%getFloats('clambdaslip', requiredShape=shape(prm%Nslip)) prm%rho0 = config%getFloats('rhoedge0', requiredSize=size(prm%Nslip))
prm%tau_Peierls = config%getFloats('tau_peierls', requiredShape=shape(prm%Nslip)) ! ToDo: Deprecated prm%rhoDip0 = config%getFloats('rhoedgedip0', requiredSize=size(prm%Nslip))
prm%p = config%getFloats('p_slip', requiredShape=shape(prm%Nslip), & prm%v0 = config%getFloats('v0', requiredSize=size(prm%Nslip))
prm%burgers = config%getFloats('slipburgers', requiredSize=size(prm%Nslip))
prm%H0kp = config%getFloats('qedge', requiredSize=size(prm%Nslip))
prm%clambda = config%getFloats('clambdaslip', requiredSize=size(prm%Nslip))
prm%tau_Peierls = config%getFloats('tau_peierls', requiredSize=size(prm%Nslip)) ! ToDo: Deprecated
prm%p = config%getFloats('p_slip', requiredSize=size(prm%Nslip), &
defaultVal=[(1.0_pReal,i=1_pInt,size(prm%Nslip))]) defaultVal=[(1.0_pReal,i=1_pInt,size(prm%Nslip))])
prm%q = config%getFloats('q_slip', requiredShape=shape(prm%Nslip), & prm%q = config%getFloats('q_slip', requiredSize=size(prm%Nslip), &
defaultVal=[(1.0_pReal,i=1_pInt,size(prm%Nslip))]) defaultVal=[(1.0_pReal,i=1_pInt,size(prm%Nslip))])
prm%kink_height = config%getFloats('kink_height', requiredShape=shape(prm%Nslip)) prm%kink_height = config%getFloats('kink_height', requiredSize=size(prm%Nslip))
prm%w = config%getFloats('kink_width', requiredShape=shape(prm%Nslip)) prm%w = config%getFloats('kink_width', requiredSize=size(prm%Nslip))
prm%omega = config%getFloats('omega', requiredShape=shape(prm%Nslip)) prm%omega = config%getFloats('omega', requiredSize=size(prm%Nslip))
prm%B = config%getFloats('friction_coeff', requiredShape=shape(prm%Nslip)) prm%B = config%getFloats('friction_coeff', requiredSize=size(prm%Nslip))
prm%SolidSolutionStrength = config%getFloat('solidsolutionstrength') ! ToDo: Deprecated prm%SolidSolutionStrength = config%getFloat('solidsolutionstrength') ! ToDo: Deprecated
prm%grainSize = config%getFloat('grainsize') prm%grainSize = config%getFloat('grainsize')
@ -250,7 +248,7 @@ subroutine plastic_disloUCLA_init()
prm%Qsd = config%getFloat('qsd') prm%Qsd = config%getFloat('qsd')
prm%atomicVolume = config%getFloat('catomicvolume') * prm%burgers**3.0_pReal prm%atomicVolume = config%getFloat('catomicvolume') * prm%burgers**3.0_pReal
prm%minDipDistance = config%getFloat('cedgedipmindistance') * prm%burgers prm%minDipDistance = config%getFloat('cedgedipmindistance') * prm%burgers
prm%dipoleformation = config%getFloat('dipoleformationfactor') > 0.0_pReal !should be on by default, ToDo: change to /key/-key prm%dipoleformation = config%getFloat('dipoleformationfactor') > 0.0_pReal !should be on by default, ToDo: change to /key/-type key
! expand: family => system ! expand: family => system
prm%rho0 = math_expand(prm%rho0, prm%Nslip) prm%rho0 = math_expand(prm%rho0, prm%Nslip)
@ -313,8 +311,6 @@ subroutine plastic_disloUCLA_init()
outputID = merge(mfp_ID,undefined_ID,prm%totalNslip>0_pInt) outputID = merge(mfp_ID,undefined_ID,prm%totalNslip>0_pInt)
case ('threshold_stress','threshold_stress_slip') case ('threshold_stress','threshold_stress_slip')
outputID = merge(thresholdstress_ID,undefined_ID,prm%totalNslip>0_pInt) outputID = merge(thresholdstress_ID,undefined_ID,prm%totalNslip>0_pInt)
case ('edge_dipole_distance')
outputID = merge(dipoleDistance_ID,undefined_ID,prm%totalNslip>0_pInt)
end select end select
@ -336,24 +332,6 @@ subroutine plastic_disloUCLA_init()
prm%totalNslip,0_pInt,0_pInt) prm%totalNslip,0_pInt,0_pInt)
plasticState(p)%sizePostResults = sum(plastic_disloUCLA_sizePostResult(:,phase_plasticityInstance(p))) plasticState(p)%sizePostResults = sum(plastic_disloUCLA_sizePostResult(:,phase_plasticityInstance(p)))
allocate(prm%forestProjectionEdge(prm%totalNslip,prm%totalNslip),source = 0.0_pReal)
i = 0_pInt
mySlipFamilies: do f = 1_pInt,size(prm%Nslip,1)
index_myFamily = sum(prm%Nslip(1:f-1_pInt))
slipSystemsLoop: do j = 1_pInt,prm%Nslip(f)
i = i + 1_pInt
do o = 1_pInt, size(prm%Nslip,1)
index_otherFamily = sum(prm%Nslip(1:o-1_pInt))
do k = 1_pInt,prm%Nslip(o) ! loop over (active) systems in other family (slip)
prm%forestProjectionEdge(index_myFamily+j,index_otherFamily+k) = &
abs(math_mul3x3(lattice_sn(:,sum(lattice_NslipSystem(1:f-1,p))+j,p), &
lattice_st(:,sum(lattice_NslipSystem(1:o-1,p))+k,p)))
enddo; enddo
enddo slipSystemsLoop
enddo mySlipFamilies
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and aTolState ! locally defined state aliases and initialization of state0 and aTolState
startIndex = 1_pInt startIndex = 1_pInt
@ -374,7 +352,7 @@ subroutine plastic_disloUCLA_init()
endIndex = endIndex + prm%totalNslip endIndex = endIndex + prm%totalNslip
stt%accshear=>plasticState(p)%state(startIndex:endIndex,:) stt%accshear=>plasticState(p)%state(startIndex:endIndex,:)
dot%accshear=>plasticState(p)%dotState(startIndex:endIndex,:) dot%accshear=>plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = 1e6_pReal !ToDo: better make optional parameter plasticState(p)%aTolState(startIndex:endIndex) = 1.0e6_pReal !ToDo: better make optional parameter
! global alias ! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:) plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:)
@ -579,16 +557,6 @@ function plastic_disloUCLA_postResults(Mp,Temperature,instance,of) result(postRe
postResults(c+1_pInt:c+prm%totalNslip) = dst%mfp(1_pInt:prm%totalNslip, of) postResults(c+1_pInt:c+prm%totalNslip) = dst%mfp(1_pInt:prm%totalNslip, of)
case (thresholdstress_ID) case (thresholdstress_ID)
postResults(c+1_pInt:c+prm%totalNslip) = dst%threshold_stress(1_pInt:prm%totalNslip,of) postResults(c+1_pInt:c+prm%totalNslip) = dst%threshold_stress(1_pInt:prm%totalNslip,of)
case (dipoleDistance_ID) ! ToDo: Discuss required changes with Franz
do i = 1_pInt, prm%totalNslip
if (dNeq0(abs(math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i))))) then
postResults(c+i) = (3.0_pReal*prm%mu*prm%burgers(i)) &
/ (16.0_pReal*pi*abs(math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i))))
else
postResults(c+i) = huge(1.0_pReal)
endif
postResults(c+i)=min(postResults(c+i),dst%mfp(i,of))
enddo
end select end select

1602
src/plastic_dislotwin.f90
View File

File diff suppressed because it is too large Load Diff

39
src/plastic_kinematichardening.f90
View File

@ -151,7 +151,6 @@ subroutine plastic_kinehardening_init
outputID outputID
character(len=pStringLen) :: & character(len=pStringLen) :: &
structure = '',&
extmsg = '' extmsg = ''
character(len=65536), dimension(:), allocatable :: & character(len=65536), dimension(:), allocatable :: &
outputs outputs
@ -187,8 +186,6 @@ subroutine plastic_kinehardening_init
endif endif
#endif #endif
structure = config%getString('lattice_structure')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! optional parameters that need to be defined ! optional parameters that need to be defined
prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal) prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal)
@ -203,28 +200,29 @@ subroutine plastic_kinehardening_init
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%totalNslip = sum(prm%Nslip) prm%totalNslip = sum(prm%Nslip)
slipActive: if (prm%totalNslip > 0_pInt) then slipActive: if (prm%totalNslip > 0_pInt) then
prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,structure(1:3),& prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal)) config%getFloat('c/a',defaultVal=0.0_pReal))
if(structure=='bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',& if(trim(config%getString('lattice_structure')) == 'bcc') then
defaultVal = emptyRealArray) prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt) defaultVal = emptyRealArray)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt) prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt)
else else
prm%nonSchmid_pos = prm%Schmid prm%nonSchmid_pos = prm%Schmid
prm%nonSchmid_neg = prm%Schmid prm%nonSchmid_neg = prm%Schmid
endif endif
prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, & prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, &
config%getFloats('interaction_slipslip'), & config%getFloats('interaction_slipslip'), &
structure(1:3)) config%getString('lattice_structure'))
prm%crss0 = config%getFloats('crss0', requiredShape=shape(prm%Nslip)) prm%crss0 = config%getFloats('crss0', requiredSize=size(prm%Nslip))
prm%tau1 = config%getFloats('tau1', requiredShape=shape(prm%Nslip)) prm%tau1 = config%getFloats('tau1', requiredSize=size(prm%Nslip))
prm%tau1_b = config%getFloats('tau1_b', requiredShape=shape(prm%Nslip)) prm%tau1_b = config%getFloats('tau1_b', requiredSize=size(prm%Nslip))
prm%theta0 = config%getFloats('theta0', requiredShape=shape(prm%Nslip)) prm%theta0 = config%getFloats('theta0', requiredSize=size(prm%Nslip))
prm%theta1 = config%getFloats('theta1', requiredShape=shape(prm%Nslip)) prm%theta1 = config%getFloats('theta1', requiredSize=size(prm%Nslip))
prm%theta0_b = config%getFloats('theta0_b', requiredShape=shape(prm%Nslip)) prm%theta0_b = config%getFloats('theta0_b', requiredSize=size(prm%Nslip))
prm%theta1_b = config%getFloats('theta1_b', requiredShape=shape(prm%Nslip)) prm%theta1_b = config%getFloats('theta1_b', requiredSize=size(prm%Nslip))
prm%gdot0 = config%getFloat('gdot0') prm%gdot0 = config%getFloat('gdot0')
prm%n = config%getFloat('n_slip') prm%n = config%getFloat('n_slip')
@ -302,7 +300,6 @@ subroutine plastic_kinehardening_init
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState, & call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState, &
prm%totalNslip,0_pInt,0_pInt) prm%totalNslip,0_pInt,0_pInt)
plasticState(p)%sizePostResults = sum(plastic_kinehardening_sizePostResult(:,phase_plasticityInstance(p))) plasticState(p)%sizePostResults = sum(plastic_kinehardening_sizePostResult(:,phase_plasticityInstance(p)))
plasticState(p)%offsetDeltaState = sizeDotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and aTolState ! locally defined state aliases and initialization of state0 and aTolState

50
src/plastic_phenopowerlaw.f90
View File

@ -153,7 +153,6 @@ subroutine plastic_phenopowerlaw_init
outputID outputID
character(len=pStringLen) :: & character(len=pStringLen) :: &
structure = '',&
extmsg = '' extmsg = ''
character(len=65536), dimension(:), allocatable :: & character(len=65536), dimension(:), allocatable :: &
outputs outputs
@ -181,8 +180,6 @@ subroutine plastic_phenopowerlaw_init
stt => state(phase_plasticityInstance(p)), & stt => state(phase_plasticityInstance(p)), &
config => config_phase(p)) config => config_phase(p))
structure = config%getString('lattice_structure')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! optional parameters that need to be defined ! optional parameters that need to be defined
prm%twinB = config%getFloat('twin_b',defaultVal=1.0_pReal) prm%twinB = config%getFloat('twin_b',defaultVal=1.0_pReal)
@ -204,30 +201,31 @@ subroutine plastic_phenopowerlaw_init
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%totalNslip = sum(prm%Nslip) prm%totalNslip = sum(prm%Nslip)
slipActive: if (prm%totalNslip > 0_pInt) then slipActive: if (prm%totalNslip > 0_pInt) then
prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,structure(1:3),& prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal)) config%getFloat('c/a',defaultVal=0.0_pReal))
if(structure=='bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',& if(trim(config%getString('lattice_structure')) == 'bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
defaultVal = emptyRealArray) defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt) prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt) prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt)
else else
prm%nonSchmid_pos = prm%Schmid_slip prm%nonSchmid_pos = prm%Schmid_slip
prm%nonSchmid_neg = prm%Schmid_slip prm%nonSchmid_neg = prm%Schmid_slip
endif endif
prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, & prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, &
config%getFloats('interaction_slipslip'), & config%getFloats('interaction_slipslip'), &
structure(1:3)) config%getString('lattice_structure'))
prm%xi_slip_0 = config%getFloats('tau0_slip', requiredSize=size(prm%Nslip)) prm%xi_slip_0 = config%getFloats('tau0_slip', requiredSize=size(prm%Nslip))
prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(prm%Nslip)) prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(prm%Nslip))
prm%H_int = config%getFloats('h_int', requiredSize=size(prm%Nslip), & prm%H_int = config%getFloats('h_int', requiredSize=size(prm%Nslip), &
defaultVal=[(0.0_pReal,i=1_pInt,size(prm%Nslip))]) defaultVal=[(0.0_pReal,i=1_pInt,size(prm%Nslip))])
prm%gdot0_slip = config%getFloat('gdot0_slip') prm%gdot0_slip = config%getFloat('gdot0_slip')
prm%n_slip = config%getFloat('n_slip') prm%n_slip = config%getFloat('n_slip')
prm%a_slip = config%getFloat('a_slip') prm%a_slip = config%getFloat('a_slip')
prm%h0_SlipSlip = config%getFloat('h0_slipslip') prm%h0_SlipSlip = config%getFloat('h0_slipslip')
! expand: family => system ! expand: family => system
prm%xi_slip_0 = math_expand(prm%xi_slip_0, prm%Nslip) prm%xi_slip_0 = math_expand(prm%xi_slip_0, prm%Nslip)
@ -239,7 +237,7 @@ subroutine plastic_phenopowerlaw_init
if ( prm%a_slip <= 0.0_pReal) extmsg = trim(extmsg)//' a_slip' if ( prm%a_slip <= 0.0_pReal) extmsg = trim(extmsg)//' a_slip'
if ( prm%n_slip <= 0.0_pReal) extmsg = trim(extmsg)//' n_slip' if ( prm%n_slip <= 0.0_pReal) extmsg = trim(extmsg)//' n_slip'
if (any(prm%xi_slip_0 <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_slip_0' if (any(prm%xi_slip_0 <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_slip_0'
if (any(prm%xi_slip_sat < prm%xi_slip_0)) extmsg = trim(extmsg)//' xi_slip_sat' if (any(prm%xi_slip_sat <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_slip_sat'
else slipActive else slipActive
allocate(prm%interaction_SlipSlip(0,0)) allocate(prm%interaction_SlipSlip(0,0))
allocate(prm%xi_slip_0(0)) allocate(prm%xi_slip_0(0))
@ -250,12 +248,12 @@ subroutine plastic_phenopowerlaw_init
prm%Ntwin = config%getInts('ntwin', defaultVal=emptyIntArray) prm%Ntwin = config%getInts('ntwin', defaultVal=emptyIntArray)
prm%totalNtwin = sum(prm%Ntwin) prm%totalNtwin = sum(prm%Ntwin)
twinActive: if (prm%totalNtwin > 0_pInt) then twinActive: if (prm%totalNtwin > 0_pInt) then
prm%Schmid_twin = lattice_SchmidMatrix_twin(prm%Ntwin,structure(1:3),& prm%Schmid_twin = lattice_SchmidMatrix_twin(prm%Ntwin,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal)) config%getFloat('c/a',defaultVal=0.0_pReal))
prm%interaction_TwinTwin = lattice_interaction_TwinTwin(prm%Ntwin,& prm%interaction_TwinTwin = lattice_interaction_TwinTwin(prm%Ntwin,&
config%getFloats('interaction_twintwin'), & config%getFloats('interaction_twintwin'), &
structure(1:3)) config%getString('lattice_structure'))
prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,structure(1:3),& prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,config%getString('lattice_structure'),&
config%getFloat('c/a')) config%getFloat('c/a'))
prm%xi_twin_0 = config%getFloats('tau0_twin',requiredSize=size(prm%Ntwin)) prm%xi_twin_0 = config%getFloats('tau0_twin',requiredSize=size(prm%Ntwin))
@ -282,10 +280,10 @@ subroutine plastic_phenopowerlaw_init
slipAndTwinActive: if (prm%totalNslip > 0_pInt .and. prm%totalNtwin > 0_pInt) then slipAndTwinActive: if (prm%totalNslip > 0_pInt .and. prm%totalNtwin > 0_pInt) then
prm%interaction_SlipTwin = lattice_interaction_SlipTwin(prm%Nslip,prm%Ntwin,& prm%interaction_SlipTwin = lattice_interaction_SlipTwin(prm%Nslip,prm%Ntwin,&
config%getFloats('interaction_sliptwin'), & config%getFloats('interaction_sliptwin'), &
structure(1:3)) config%getString('lattice_structure'))
prm%interaction_TwinSlip = lattice_interaction_TwinSlip(prm%Ntwin,prm%Nslip,& prm%interaction_TwinSlip = lattice_interaction_TwinSlip(prm%Ntwin,prm%Nslip,&
config%getFloats('interaction_twinslip'), & config%getFloats('interaction_twinslip'), &
structure(1:3)) config%getString('lattice_structure'))
else slipAndTwinActive else slipAndTwinActive
allocate(prm%interaction_SlipTwin(prm%totalNslip,prm%TotalNtwin)) ! at least one dimension is 0 allocate(prm%interaction_SlipTwin(prm%totalNslip,prm%TotalNtwin)) ! at least one dimension is 0
allocate(prm%interaction_TwinSlip(prm%totalNtwin,prm%TotalNslip)) ! at least one dimension is 0 allocate(prm%interaction_TwinSlip(prm%totalNtwin,prm%TotalNslip)) ! at least one dimension is 0