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

Merge remote-tracking branch 'origin/development' into 42-new-coding-style-for-homogenization-NEW

This commit is contained in:
Martin Diehl 2019-01-27 08:39:11 +01:00
parent 9058587a2b 012aa4c697
commit 7484849b22
23 changed files with 3602 additions and 5004 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 5ed6a1f60b412eb46ff6820cf03b684095ff1f75
Subproject commit 683bf0074f3fa079989b51f5a67aa593b7577f0b

2
VERSION
View File

@ -1 +1 @@
v2.0.2-1395-g4848e600
v2.0.2-1516-gffd29bdc

2
examples/ConfigFiles/Homogenization_None_Dummy.config
View File

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

4
examples/ConfigFiles/Phase_Isotropic_AluminumIsotropic.config
View File

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

2
processing/post/viewTable.py
View File

@ -68,7 +68,7 @@ for name in filenames:
(['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 ---------------------------------------

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_Fi0, &
crystallite_Li0, &
crystallite_dPdF0, &
crystallite_Tstar0_v
implicit none
@ -207,9 +206,6 @@ subroutine CPFEM_init
read (777,rec=1) crystallite_Li0
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))
read (777,rec=1) crystallite_Tstar0_v
@ -286,12 +282,11 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
math_identity2nd, &
math_mul33x33, &
math_det33, &
math_transpose33, &
math_I3, &
math_Mandel3333to66, &
math_Mandel66to3333, &
math_Mandel33to6, &
math_Mandel6to33
math_delta, &
math_sym3333to66, &
math_66toSym3333, &
math_sym33to6, &
math_6toSym33
use mesh, only: &
mesh_FEasCP, &
mesh_NcpElems, &
@ -326,7 +321,6 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
crystallite_Lp, &
crystallite_Li0, &
crystallite_Li, &
crystallite_dPdF0, &
crystallite_dPdF, &
crystallite_Tstar0_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
real(pReal), intent(in) :: temperature_inp !< temperature
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,6), intent(out) :: jacobian !< jacobian in Mandel notation (Consistent tangent dcs/dE)
real(pReal), dimension(6), intent(out) :: cauchyStress !< stress as 6 vector
real(pReal), dimension(6,6), intent(out) :: jacobian !< jacobian as 66 tensor (Consistent tangent dcs/dE)
real(pReal) J_inverse, & ! inverse of Jacobian
rnd
@ -398,7 +392,6 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity
crystallite_Fi0 = crystallite_Fi ! crystallite intermediate deformation
crystallite_Li0 = crystallite_Li ! crystallite intermediate velocity
crystallite_dPdF0 = crystallite_dPdF ! crystallite stiffness
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
@ -454,10 +447,6 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
write (777,rec=1) crystallite_Li0
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))
write (777,rec=1) crystallite_Tstar0_v
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
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:',&
math_transpose33(materialpoint_F(1:3,1:3,ip,elCP))
write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 now:',math_transpose33(ffn1)
transpose(materialpoint_F(1:3,1:3,ip,elCP))
write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 now:',transpose(ffn1)
endif
outdatedFFN1 = .true.
endif
@ -593,26 +582,25 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
endif
! 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))
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
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
H(i,j,k,l) = H(i,j,k,l) + &
materialpoint_F(j,m,ip,elCP) * &
materialpoint_F(l,n,ip,elCP) * &
materialpoint_dPdF(i,m,k,n,ip,elCP) - &
math_I3(j,l) * materialpoint_F(i,m,ip,elCP) * materialpoint_P(k,m,ip,elCP) + &
0.5_pReal * (math_I3(i,k) * Kirchhoff(j,l) + math_I3(j,l) * Kirchhoff(i,k) + &
math_I3(i,l) * Kirchhoff(j,k) + math_I3(j,k) * Kirchhoff(i,l))
H(i,j,k,l) = H(i,j,k,l) &
+ materialpoint_F(j,m,ip,elCP) * materialpoint_F(l,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) &
+ 0.5_pReal * ( Kirchhoff(j,l)*math_delta(i,k) + Kirchhoff(i,k)*math_delta(j,l) &
+ Kirchhoff(j,k)*math_delta(i,l) + Kirchhoff(i,l)*math_delta(j,k))
enddo; enddo; enddo; enddo; enddo; enddo
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))
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 validCalculation
@ -639,7 +627,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
!*** 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
debug_stressMaxLocation = [elCP, ip]
debug_stressMax = maxval(cauchyStress33)
@ -649,7 +637,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
debug_stressMin = minval(cauchyStress33)
endif
!*** ... 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
debug_jacobianMaxLocation = [elCP, ip]
debug_jacobianMax = maxval(jacobian3333)

5
src/CPFEM2.f90
View File

@ -121,7 +121,6 @@ subroutine CPFEM_init
crystallite_Lp0, &
crystallite_Fi0, &
crystallite_Li0, &
crystallite_dPdF0, &
crystallite_Tstar0_v
use hdf5
use HDF5_utilities, only: &
@ -160,7 +159,6 @@ subroutine CPFEM_init
call HDF5_read(fileHandle,crystallite_Fi0, 'convergedFi')
call HDF5_read(fileHandle,crystallite_Lp0, 'convergedLp')
call HDF5_read(fileHandle,crystallite_Li0, 'convergedLi')
call HDF5_read(fileHandle,crystallite_dPdF0, 'convergeddPdF')
call HDF5_read(fileHandle,crystallite_Tstar0_v,'convergedTstar')
groupPlasticID = HDF5_openGroup(fileHandle,'PlasticPhases')
@ -224,7 +222,6 @@ subroutine CPFEM_age()
crystallite_Lp, &
crystallite_Li0, &
crystallite_Li, &
crystallite_dPdF0, &
crystallite_dPdF, &
crystallite_Tstar0_v, &
crystallite_Tstar_v
@ -254,7 +251,6 @@ subroutine CPFEM_age()
crystallite_Lp0 = crystallite_Lp
crystallite_Fi0 = crystallite_Fi
crystallite_Li0 = crystallite_Li
crystallite_dPdF0 = crystallite_dPdF
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
@ -283,7 +279,6 @@ subroutine CPFEM_age()
call HDF5_write(fileHandle,crystallite_Fi0, 'convergedFi')
call HDF5_write(fileHandle,crystallite_Lp0, 'convergedLp')
call HDF5_write(fileHandle,crystallite_Li0, 'convergedLi')
call HDF5_write(fileHandle,crystallite_dPdF0, 'convergeddPdF')
call HDF5_write(fileHandle,crystallite_Tstar0_v,'convergedTstar')
groupPlastic = HDF5_addGroup(fileHandle,'PlasticPhases')

14
src/DAMASK_abaqus.f
View File

@ -102,8 +102,6 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
calcMode, &
terminallyIll, &
symmetricSolver
use math, only: &
invnrmMandel
use debug, only: &
debug_info, &
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
forall(i=1:ntens) ddsdde(1:ntens,i) = invnrmMandel(i)*ddsdde_h(1:ntens,i)*invnrmMandel(1:ntens)
stress(1:ntens) = stress_h(1:ntens)*invnrmMandel(1:ntens)
if(symmetricSolver) ddsdde(1:ntens,1:ntens) = 0.5_pReal*(ddsdde(1:ntens,1:ntens) + transpose(ddsdde(1:ntens,1:ntens)))
ddsdde = ddsdde_h(1:ntens,1:ntens)
stress = stress_h(1:ntens)
if(symmetricSolver) ddsdde = 0.5_pReal*(ddsdde + transpose(ddsdde))
if(ntens == 6) then
stress_h = stress
stress(5) = stress_h(6)
@ -322,7 +320,7 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
statev = materialpoint_results(1:min(nstatv,materialpoint_sizeResults),npt,mesh_FEasCP('elem', noel))
if ( terminallyIll ) pnewdt = 0.5_pReal ! force cutback directly ?
if (terminallyIll) pnewdt = 0.5_pReal ! force cutback directly ?
!$ call omp_set_num_threads(defaultNumThreadsInt) ! reset number of threads to stored default value
end subroutine UMAT
@ -331,12 +329,12 @@ end subroutine UMAT
!--------------------------------------------------------------------------------------------------
!> @brief calls the exit function of Abaqus/Standard
!--------------------------------------------------------------------------------------------------
subroutine quit(mpie_error)
subroutine quit(DAMASK_error)
use prec, only: &
pInt
implicit none
integer(pInt) :: mpie_error
integer(pInt) :: DAMASK_error
flush(6)
call xit

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, &
terminallyIll, &
symmetricSolver
use math, only: &
math_transpose33,&
invnrmMandel
use debug, only: &
debug_level, &
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,i1)') ' Deformation gradient: ', itel
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:', &
math_transpose33(ffn1)
transpose(ffn1)
endif
!$ 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
forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*ddsdde(1:ngens,i)*invnrmMandel(1:ngens)
s(1:ndi+nshear) = stress(1:ndi+nshear)*invnrmMandel(1:ndi+nshear)
d = ddsdde(1:ngens,1:ngens)
s = stress(1:ndi+nshear)
g = 0.0_pReal
if(symmetricSolver) d = 0.5_pReal*(d+transpose(d))

4
src/IO.f90
View File

@ -1236,6 +1236,10 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
msg = 'zero entry on stiffness diagonal'
case (136_pInt)
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

20
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_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_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_NONLOCAL_ID)) then
call plastic_nonlocal_init(FILEUNIT)
@ -920,7 +920,7 @@ end subroutine constitutive_collectDotState
!> @brief for constitutive models having an instantaneous change of state
!> 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: &
pReal, &
pLongInt
@ -929,8 +929,7 @@ subroutine constitutive_collectDeltaState(S6, Fe, Fi, ipc, ip, el)
debug_constitutive, &
debug_levelBasic
use math, only: &
math_Mandel6to33, &
math_Mandel33to6, &
math_sym33to6, &
math_mul33x33
use material, only: &
phasememberAt, &
@ -954,18 +953,17 @@ subroutine constitutive_collectDeltaState(S6, Fe, Fi, ipc, ip, el)
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(6) :: &
S6 !< 2nd Piola Kirchhoff stress (vector notation)
real(pReal), intent(in), dimension(3,3) :: &
S, & !< 2nd Piola Kirchhoff stress
Fe, & !< elastic deformation gradient
Fi !< intermediate deformation gradient
real(pReal), dimension(3,3) :: &
Mp
integer(pInt) :: &
s, & !< counter in source loop
i, &
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)))
@ -975,13 +973,13 @@ subroutine constitutive_collectDeltaState(S6, Fe, Fi, ipc, ip, el)
call plastic_kinehardening_deltaState(Mp,instance,of)
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
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
call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(ipc,ip,el), Fe, &

4624
src/crystallite.f90
View File

File diff suppressed because it is too large Load Diff

19
src/homogenization.f90
View File

@ -347,7 +347,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
crystallite_Li0, &
crystallite_Li, &
crystallite_dPdF, &
crystallite_dPdF0, &
crystallite_Tstar0_v, &
crystallite_Tstar_v, &
crystallite_partionedF0, &
@ -356,12 +355,11 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
crystallite_partionedLp0, &
crystallite_partionedFi0, &
crystallite_partionedLi0, &
crystallite_partioneddPdF0, &
crystallite_partionedTstar0_v, &
crystallite_dt, &
crystallite_requested, &
crystallite_converged, &
crystallite_stressAndItsTangent, &
crystallite_stress, &
crystallite_stressTangent, &
crystallite_orientations
#ifdef DEBUG
use debug, only: &
@ -414,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_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_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_partionedTstar0_v(1:6,g,i,e) = crystallite_Tstar0_v(1:6,g,i,e) ! ...2nd PK stress
@ -484,9 +481,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
crystallite_partionedLi0(1:3,1:3,1:myNgrains,i,e) = &
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_Tstar_v(1:6,1:myNgrains,i,e) ! ...2nd PK stress
@ -550,8 +544,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
crystallite_partionedFi0(1:3,1:3,1:myNgrains,i,e) ! ...intermediate def grads
crystallite_Li(1:3,1:3,1:myNgrains,i,e) = &
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_partionedTstar0_v(1:6,1:myNgrains,i,e) ! ...2nd PK stress
do g = 1, myNgrains
@ -622,7 +614,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
! crystallite integration
! based on crystallite_partionedF0,.._partionedF
! 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
@ -631,9 +623,8 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
IpLooping3: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
if ( materialpoint_requested(i,e) .and. &
.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_converged(i,e) = .false.
else
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
@ -649,6 +640,8 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
enddo cutBackLooping
if(updateJaco) call crystallite_stressTangent
if (.not. terminallyIll ) then
call crystallite_orientations() ! calculate crystal orientations
!$OMP PARALLEL DO

2311
src/lattice.f90
View File

File diff suppressed because it is too large Load Diff

4
src/material.f90
View File

@ -918,7 +918,8 @@ end subroutine material_parseTexture
!--------------------------------------------------------------------------------------------------
!> @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)
use numerics, only: &
numerics_integrator2 => numerics_integrator ! compatibility hack
@ -939,6 +940,7 @@ subroutine material_allocatePlasticState(phase,NofMyPhase,sizeState,sizeDotState
plasticState(phase)%sizeState = sizeState
plasticState(phase)%sizeDotState = sizeDotState
plasticState(phase)%sizeDeltaState = sizeDeltaState
plasticState(phase)%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition
plasticState(phase)%Nslip = Nslip
plasticState(phase)%Ntwin = Ntwin
plasticState(phase)%Ntrans= Ntrans

413
src/math.f90
View File

@ -24,25 +24,25 @@ module math
0.0_pReal,0.0_pReal,1.0_pReal &
],[3,3]) !< 3x3 Identity
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 , private :: &
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 :: &
mapMandel = reshape([&
mapNye = reshape([&
1_pInt,1_pInt, &
2_pInt,2_pInt, &
3_pInt,3_pInt, &
1_pInt,2_pInt, &
2_pInt,3_pInt, &
1_pInt,3_pInt &
],[2,6]) !< arrangement in Mandel 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)
],[2,6]) !< arrangement in Nye notation.
integer(pInt), dimension (2,6), parameter, private :: &
mapVoigt = reshape([&
@ -54,10 +54,6 @@ module math
1_pInt,2_pInt &
],[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 :: &
mapPlain = reshape([&
1_pInt,1_pInt, &
@ -71,6 +67,56 @@ module math
3_pInt,3_pInt &
],[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 :: &
math_init, &
math_qsort, &
@ -109,16 +155,14 @@ module math
math_equivStress33, &
math_trace33, &
math_det33, &
math_Plain33to9, &
math_Plain9to33, &
math_Mandel33to6, &
math_Mandel6to33, &
math_Plain3333to99, &
math_Plain99to3333, &
math_Mandel66toPlain66, &
math_Plain66toMandel66, &
math_Mandel3333to66, &
math_Mandel66to3333, &
math_33to9, &
math_9to33, &
math_sym33to6, &
math_6toSym33, &
math_3333to99, &
math_99to3333, &
math_sym3333to66, &
math_66toSym3333, &
math_Voigt66to3333, &
math_qRand, &
math_qMul, &
@ -416,7 +460,7 @@ pure function math_identity2nd(dimen)
real(pReal), dimension(dimen,dimen) :: math_identity2nd
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
@ -430,9 +474,11 @@ pure function math_identity4th(dimen)
integer(pInt), intent(in) :: dimen !< tensor dimension
integer(pInt) :: i,j,k,l
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) = &
0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k))
identity2nd = math_identity2nd(dimen)
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
@ -501,7 +547,7 @@ pure function math_tensorproduct(A,B)
real(pReal), dimension(size(A,1),size(B,1)) :: math_tensorproduct
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
@ -516,7 +562,7 @@ pure function math_tensorproduct33(A,B)
real(pReal), dimension(3), intent(in) :: A,B
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
@ -557,7 +603,7 @@ real(pReal) pure function math_mul33xx33(A,B)
integer(pInt) :: i,j
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)
end function math_mul33xx33
@ -574,8 +620,7 @@ pure function math_mul3333xx33(A,B)
real(pReal), dimension(3,3), intent(in) :: B
integer(pInt) :: i,j
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))
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))
end function math_mul3333xx33
@ -607,8 +652,7 @@ pure function math_mul33x33(A,B)
real(pReal), dimension(3,3), intent(in) :: A,B
integer(pInt) :: i,j
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)
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)
end function math_mul33x33
@ -623,9 +667,9 @@ pure function math_mul66x66(A,B)
real(pReal), dimension(6,6), intent(in) :: A,B
integer(pInt) :: i,j
forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt) 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)
forall(i=1_pInt:6_pInt,j=1_pInt:6_pInt) &
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)
end function math_mul66x66
@ -640,10 +684,10 @@ pure function math_mul99x99(A,B)
real(pReal), dimension(9,9), intent(in) :: A,B
integer(pInt) i,j
forall (i=1_pInt:9_pInt,j=1_pInt:9_pInt) 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,7)*B(7,j) + A(i,8)*B(8,j) + A(i,9)*B(9,j)
forall(i=1_pInt:9_pInt,j=1_pInt:9_pInt) &
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,7)*B(7,j) + A(i,8)*B(8,j) + A(i,9)*B(9,j)
end function math_mul99x99
@ -691,9 +735,8 @@ pure function math_mul66x6(A,B)
real(pReal), dimension(6), intent(in) :: B
integer(pInt) :: 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,4)*B(4) + A(i,5)*B(5) + A(i,6)*B(6)
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,4)*B(4) + A(i,5)*B(5) + A(i,6)*B(6)
end function math_mul66x6
@ -740,8 +783,8 @@ end function math_transpose33
!--------------------------------------------------------------------------------------------------
!> @brief Cramer inversion of 33 matrix (function)
! direct Cramer inversion of matrix A.
! returns all zeroes if not possible, i.e. if det close to zero
!> @details Direct Cramer inversion of matrix A. Returns all zeroes if not possible, i.e.
! if determinant is close to zero
!--------------------------------------------------------------------------------------------------
pure function math_inv33(A)
use prec, only: &
@ -777,10 +820,9 @@ end function math_inv33
!--------------------------------------------------------------------------------------------------
!> @brief Cramer inversion of 33 matrix (subroutine)
! direct Cramer inversion of matrix A.
! also returns determinant
! returns error if not possible, i.e. if det close to zero
! ToDo: has wrong order of arguments (out should be first)
!> @details Direct Cramer inversion of matrix A. Also returns determinant
! Returns an error if not possible, i.e. if determinant is close to zero
! ToDo: Output arguments should be first
!--------------------------------------------------------------------------------------------------
pure subroutine math_invert33(A, InvA, DetA, error)
use prec, only: &
@ -837,11 +879,11 @@ function math_invSym3333(A)
dgetrf, &
dgetri
temp66_real = math_Mandel3333to66(A)
temp66_real = math_sym3333to66(A)
call dgetrf(6,6,temp66_real,6,ipiv6,ierr)
call dgetri(6,temp66_real,6,ipiv6,work6,6,ierr)
if (ierr == 0_pInt) then
math_invSym3333 = math_Mandel66to3333(temp66_real)
math_invSym3333 = math_66toSym3333(temp66_real)
else
call IO_error(400_pInt, ext_msg = 'math_invSym3333')
endif
@ -850,36 +892,26 @@ end function math_invSym3333
!--------------------------------------------------------------------------------------------------
!> @brief invert quare matrix of arbitrary dimension
!> @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,2)), intent(out) :: invA
real(pReal), dimension(size(A,1),size(A,1)), intent(out) :: invA
logical, intent(out) :: error
integer(pInt) :: ierr
integer(pInt), dimension(size(A,1)) :: ipiv
real(pReal), dimension(size(A,1)) :: work
external :: &
dgetrf, &
dgetri
invA = A
call dgetrf(size(A,1),size(A,2),invA,size(A,1),ipiv,ierr)
call dgetri(size(A,1),InvA,size(A,1),ipiv,work,size(A,1),ierr)
error = merge(.true.,.false., ierr /= 0_pInt)
call math_invert(size(A,1), A, InvA, error)
end subroutine math_invert2
!--------------------------------------------------------------------------------------------------
!> @brief invert matrix of arbitrary dimension
! Obsolete: has wrong order of arguments and superflouous argumen myDim
! use math_inver2 instead
! ToDo: Wrong order of arguments and superfluous myDim argument.
! Use math_invert2 instead
!--------------------------------------------------------------------------------------------------
subroutine math_invert(myDim,A, InvA, error)
@ -901,7 +933,7 @@ subroutine math_invert(myDim,A, InvA, error)
invA = A
call dgetrf(myDim,myDim,invA,myDim,ipiv,ierr)
call dgetri(myDim,InvA,myDim,ipiv,work,myDim,ierr)
error = merge(.true.,.false., ierr /= 0_pInt) ! http://fortraninacworld.blogspot.de/2012/12/ternary-operator.html
error = merge(.true.,.false., ierr /= 0_pInt)
end subroutine math_invert
@ -984,15 +1016,14 @@ pure function math_equivStrain33(m)
real(pReal), dimension(3,3), intent(in) :: m
real(pReal), dimension(3) :: e,s
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), &
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
s = [m(1,2),m(2,3),m(1,3)]*2.0_pReal
math_equivStrain33 = TWOTHIRD*(1.50_pReal*(sum(e**2.0_pReal)) + &
0.75_pReal*(sum(s**2.0_pReal)))**(0.5_pReal)
math_equivStrain33 = 2.0_pReal/3.0_pReal &
* (1.50_pReal*(sum(e**2.0_pReal))+ 0.75_pReal*(sum(s**2.0_pReal)))**(0.5_pReal)
end function math_equivStrain33
@ -1064,173 +1095,188 @@ end function math_detSym33
!--------------------------------------------------------------------------------------------------
!> @brief convert 33 matrix into vector 9
!--------------------------------------------------------------------------------------------------
pure function math_Plain33to9(m33)
pure function math_33to9(m33)
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))
forall(i=1_pInt:9_pInt) math_33to9(i) = m33(mapPlain(1,i),mapPlain(2,i))
end function math_Plain33to9
end function math_33to9
!--------------------------------------------------------------------------------------------------
!> @brief convert Plain 9 back to 33 matrix
!> @brief convert 9 vector into 33 matrix
!--------------------------------------------------------------------------------------------------
pure function math_Plain9to33(v9)
pure function math_9to33(v9)
implicit none
real(pReal), dimension(3,3) :: math_Plain9to33
real(pReal), dimension(3,3) :: math_9to33
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)
forall(i=1_pInt:9_pInt) math_9to33(mapPlain(1,i),mapPlain(2,i)) = v9(i)
end function math_Plain9to33
end function math_9to33
!--------------------------------------------------------------------------------------------------
!> @brief convert symmetric 33 matrix into Mandel vector 6
!> @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_Mandel33to6(m33)
pure function math_sym33to6(m33,weighted)
implicit none
real(pReal), dimension(6) :: math_Mandel33to6
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
forall (i=1_pInt:6_pInt) math_Mandel33to6(i) = nrmMandel(i)*m33(mapMandel(1,i),mapMandel(2,i))
if(present(weighted)) then
w = merge(nrmMandel,1.0_pReal,weighted)
else
w = nrmMandel
endif
end function math_Mandel33to6
forall(i=1_pInt:6_pInt) math_sym33to6(i) = w(i)*m33(mapNye(1,i),mapNye(2,i))
end function math_sym33to6
!--------------------------------------------------------------------------------------------------
!> @brief convert Mandel 6 back to symmetric 33 matrix
!> @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_Mandel6to33(v6)
pure function math_6toSym33(v6,weighted)
implicit none
real(pReal), dimension(3,3) :: math_6toSym33
real(pReal), dimension(6), intent(in) :: v6
real(pReal), dimension(3,3) :: math_Mandel6to33
logical, optional, intent(in) :: weighted
real(pReal), dimension(6) :: w
integer(pInt) :: i
forall (i=1_pInt:6_pInt)
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
if(present(weighted)) then
w = merge(invnrmMandel,1.0_pReal,weighted)
else
w = invnrmMandel
endif
end function math_Mandel6to33
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 tensor into plain matrix 99
!> @brief convert 3333 matrix into 99 matrix
!--------------------------------------------------------------------------------------------------
pure function math_Plain3333to99(m3333)
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(9,9) :: math_Plain3333to99
integer(pInt) :: i,j
forall (i=1_pInt:9_pInt,j=1_pInt:9_pInt) math_Plain3333to99(i,j) = &
m3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j))
forall(i=1_pInt:9_pInt,j=1_pInt:9_pInt) &
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
real(pReal), dimension(3,3,3,3) :: math_99to3333
real(pReal), dimension(9,9), intent(in) :: m99
real(pReal), dimension(3,3,3,3) :: math_Plain99to3333
integer(pInt) :: i,j
forall (i=1_pInt:9_pInt,j=1_pInt:9_pInt) math_Plain99to3333(mapPlain(1,i),mapPlain(2,i),&
mapPlain(1,j),mapPlain(2,j)) = m99(i,j)
forall(i=1_pInt:9_pInt,j=1_pInt:9_pInt) &
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
real(pReal), dimension(6,6) :: math_sym3333to66
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
forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt) math_Mandel3333to66(i,j) = &
nrmMandel(i)*nrmMandel(j)*m3333(mapMandel(1,i),mapMandel(2,i),mapMandel(1,j),mapMandel(2,j))
if(present(weighted)) then
w = merge(nrmMandel,1.0_pReal,weighted)
else
w = nrmMandel
endif
end function math_Mandel3333to66
forall(i=1_pInt:6_pInt,j=1_pInt:6_pInt) &
math_sym3333to66(i,j) = w(i)*w(j)*m3333(mapNye(1,i),mapNye(2,i),mapNye(1,j),mapNye(2,j))
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
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
logical, optional, intent(in) :: weighted
real(pReal), dimension(6) :: w
integer(pInt) :: i,j
forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt)
math_Mandel66to3333(mapMandel(1,i),mapMandel(2,i),mapMandel(1,j),mapMandel(2,j)) = &
invnrmMandel(i)*invnrmMandel(j)*m66(i,j)
math_Mandel66to3333(mapMandel(2,i),mapMandel(1,i),mapMandel(1,j),mapMandel(2,j)) = &
invnrmMandel(i)*invnrmMandel(j)*m66(i,j)
math_Mandel66to3333(mapMandel(1,i),mapMandel(2,i),mapMandel(2,j),mapMandel(1,j)) = &
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
if(present(weighted)) then
w = merge(invnrmMandel,1.0_pReal,weighted)
else
w = invnrmMandel
endif
end function math_Mandel66to3333
do i=1_pInt,6_pInt; do j=1_pInt, 6_pInt
math_66toSym3333(mapNye(1,i),mapNye(2,i),mapNye(1,j),mapNye(2,j)) = w(i)*w(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_66toSym3333(mapNye(1,i),mapNye(2,i),mapNye(2,j),mapNye(1,j)) = w(i)*w(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)
enddo; enddo
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)
@ -1239,16 +1285,12 @@ pure function math_Voigt66to3333(m66)
real(pReal), dimension(6,6), intent(in) :: m66
integer(pInt) :: i,j
forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt)
math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(1,j),mapVoigt(2,j)) = &
invnrmVoigt(i)*invnrmVoigt(j)*m66(i,j)
math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(1,j),mapVoigt(2,j)) = &
invnrmVoigt(i)*invnrmVoigt(j)*m66(i,j)
math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(2,j),mapVoigt(1,j)) = &
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
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)) = m66(i,j)
math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(1,j),mapVoigt(2,j)) = m66(i,j)
math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(2,j),mapVoigt(1,j)) = m66(i,j)
math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(2,j),mapVoigt(1,j)) = m66(i,j)
enddo; enddo
end function math_Voigt66to3333
@ -1656,8 +1698,7 @@ pure function math_qToR(q)
real(pReal), dimension(3,3) :: math_qToR, T,S
integer(pInt) :: i, j
forall (i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) &
T(i,j) = q(i+1_pInt) * q(j+1_pInt)
forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) T(i,j) = q(i+1_pInt) * q(j+1_pInt)
S = reshape( [0.0_pReal, -q(4), q(3), &
q(4), 0.0_pReal, -q(2), &
@ -1957,7 +1998,7 @@ end function math_symmetricEulers
!--------------------------------------------------------------------------------------------------
!> @brief eigenvalues and eigenvectors of symmetric matrix m
! ToDo: has wrong order of arguments
! ToDo: has wrong oder of arguments
!--------------------------------------------------------------------------------------------------
subroutine math_eigenValuesVectorsSym(m,values,vectors,error)
@ -1981,10 +2022,10 @@ end subroutine math_eigenValuesVectorsSym
!--------------------------------------------------------------------------------------------------
!> @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
!> @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
!> @details See http://arxiv.org/abs/physics/0610206 (DSYEVH3)
! ToDo: has wrong order of arguments
! ToDo: has wrong oder of arguments
!--------------------------------------------------------------------------------------------------
subroutine math_eigenValuesVectorsSym33(m,values,vectors)
@ -2064,7 +2105,7 @@ end function math_eigenvectorBasisSym
!--------------------------------------------------------------------------------------------------
!> @brief eigenvector basis of symmetric 33 matrix m
!--------------------------------------------------------------------------------------------------
function math_eigenvectorBasisSym33(m)
pure function math_eigenvectorBasisSym33(m)
implicit none
real(pReal), dimension(3,3) :: math_eigenvectorBasisSym33
@ -2129,7 +2170,7 @@ end function math_eigenvectorBasisSym33
!--------------------------------------------------------------------------------------------------
!> @brief logarithm eigenvector basis of symmetric 33 matrix m
!--------------------------------------------------------------------------------------------------
function math_eigenvectorBasisSym33_log(m)
pure function math_eigenvectorBasisSym33_log(m)
implicit none
real(pReal), dimension(3,3) :: math_eigenvectorBasisSym33_log
@ -2190,6 +2231,7 @@ function math_eigenvectorBasisSym33_log(m)
end function math_eigenvectorBasisSym33_log
!--------------------------------------------------------------------------------------------------
!> @brief rotational part from polar decomposition of 33 tensor m
!--------------------------------------------------------------------------------------------------
@ -2634,13 +2676,12 @@ pure function math_rotate_forward3333(tensor,rot_tensor)
real(pReal), dimension(3,3,3,3), intent(in) :: tensor
integer(pInt) :: i,j,k,l,m,n,o,p
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 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) &
+ rot_tensor(i,m) * rot_tensor(j,n) &
* rot_tensor(k,o) * rot_tensor(l,p) * tensor(m,n,o,p)
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 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) &
+ rot_tensor(i,m) * rot_tensor(j,n) * rot_tensor(k,o) * rot_tensor(l,p) * tensor(m,n,o,p)
enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo
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)
case ('usepingpong')
usepingpong = IO_intValue(line,chunkPos,2_pInt) > 0_pInt
case ('timesyncing')
numerics_timeSyncing = IO_intValue(line,chunkPos,2_pInt) > 0_pInt
case ('unitlength')
numerics_unitlength = IO_floatValue(line,chunkPos,2_pInt)
@ -454,8 +452,6 @@ subroutine numerics_init
end select
#endif
numerics_timeSyncing = numerics_timeSyncing .and. all(numerics_integrator==2_pInt) ! timeSyncing only allowed for explicit Euler integrator
!--------------------------------------------------------------------------------------------------
! writing parameters to output
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)') ' aTol_crystalliteStress: ',aTol_crystalliteStress
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,es8.1,/)')' unitlength: ',numerics_unitlength

50
src/plastic_disloUCLA.f90
View File

@ -81,7 +81,7 @@ module plastic_disloUCLA
rhoEdge, &
rhoEdgeDip, &
accshear
end type
end type tDisloUCLAState
type, private :: tDisloUCLAdependentState
real(pReal), allocatable, dimension(:,:) :: &
@ -90,13 +90,13 @@ module plastic_disloUCLA
threshold_stress
end type tDisloUCLAdependentState
type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance)
type(tDisloUCLAState ), allocatable, dimension(:), private :: &
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
type(tParameters), allocatable, dimension(:), private :: param
type(tDisloUCLAState), allocatable, dimension(:), private :: &
dotState, &
state
type(tDisloUCLAdependentState), allocatable, dimension(:), private :: &
dependentState
type(tDisloUCLAdependentState), allocatable, dimension(:), private :: dependentState
public :: &
plastic_disloUCLA_init, &
@ -164,7 +164,6 @@ subroutine plastic_disloUCLA_init()
outputID
character(len=pStringLen) :: &
structure = '',&
extmsg = ''
character(len=65536), dimension(:), allocatable :: &
outputs
@ -197,8 +196,6 @@ subroutine plastic_disloUCLA_init()
dst => dependentState(phase_plasticityInstance(p)), &
config => config_phase(p))
structure = config%getString('lattice_structure')
!--------------------------------------------------------------------------------------------------
! optional parameters that need to be defined
prm%mu = lattice_mu(p)
@ -213,9 +210,10 @@ subroutine plastic_disloUCLA_init()
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%totalNslip = sum(prm%Nslip)
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))
if(structure=='bcc') then
if(trim(config%getString('lattice_structure')) == 'bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt)
@ -226,23 +224,23 @@ subroutine plastic_disloUCLA_init()
endif
prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, &
config%getFloats('interaction_slipslip'), &
structure(1:3))
prm%rho0 = config%getFloats('rhoedge0', requiredShape=shape(prm%Nslip))
prm%rhoDip0 = config%getFloats('rhoedgedip0', requiredShape=shape(prm%Nslip))
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))
config%getString('lattice_structure'))
prm%rho0 = config%getFloats('rhoedge0', requiredSize=size(prm%Nslip))
prm%rhoDip0 = config%getFloats('rhoedgedip0', requiredSize=size(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', requiredShape=shape(prm%Nslip))
prm%tau_Peierls = config%getFloats('tau_peierls', requiredShape=shape(prm%Nslip)) ! ToDo: Deprecated
prm%p = config%getFloats('p_slip', requiredShape=shape(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))])
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))])
prm%kink_height = config%getFloats('kink_height', requiredShape=shape(prm%Nslip))
prm%w = config%getFloats('kink_width', requiredShape=shape(prm%Nslip))
prm%omega = config%getFloats('omega', requiredShape=shape(prm%Nslip))
prm%B = config%getFloats('friction_coeff', requiredShape=shape(prm%Nslip))
prm%kink_height = config%getFloats('kink_height', requiredSize=size(prm%Nslip))
prm%w = config%getFloats('kink_width', requiredSize=size(prm%Nslip))
prm%omega = config%getFloats('omega', requiredSize=size(prm%Nslip))
prm%B = config%getFloats('friction_coeff', requiredSize=size(prm%Nslip))
prm%SolidSolutionStrength = config%getFloat('solidsolutionstrength') ! ToDo: Deprecated
prm%grainSize = config%getFloat('grainsize')
@ -250,7 +248,7 @@ subroutine plastic_disloUCLA_init()
prm%Qsd = config%getFloat('qsd')
prm%atomicVolume = config%getFloat('catomicvolume') * prm%burgers**3.0_pReal
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
prm%rho0 = math_expand(prm%rho0, prm%Nslip)

351
src/plastic_dislotwin.f90
View File

@ -104,12 +104,11 @@ module plastic_dislotwin
interaction_TwinSlip, & !< coefficients for twin-slip interaction for each interaction type
interaction_TwinTwin, & !< coefficients for twin-twin interaction for each interaction type
interaction_SlipTrans, & !< coefficients for slip-trans interaction for each interaction type
interaction_TransSlip, & !< coefficients for trans-slip interaction for each interaction type
interaction_TransTrans !< coefficients for trans-trans interaction for each interaction type
integer(pInt), dimension(:,:), allocatable :: &
fcc_twinNucleationSlipPair
fcc_twinNucleationSlipPair ! ToDo: Better name? Is also use for trans
real(pReal), dimension(:,:), allocatable :: &
forestProjectionEdge, &
forestProjection, &
C66
real(pReal), dimension(:,:,:), allocatable :: &
Schmid_trans, &
@ -124,7 +123,7 @@ module plastic_dislotwin
outputID !< ID of each post result output
logical :: &
isFCC !< twinning and transformation models are for fcc
fccTwinTransNucleation !< twinning and transformation models are for fcc
integer(pInt) :: &
totalNslip, & !< number of active slip systems for each family and instance
totalNtwin, & !< number of active twin systems for each family and instance
@ -190,7 +189,7 @@ contains
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine plastic_dislotwin_init(fileUnit)
subroutine plastic_dislotwin_init
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
@ -230,8 +229,6 @@ subroutine plastic_dislotwin_init(fileUnit)
use lattice
implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt) :: Ninstance,&
f,j,i,k,o,p, &
offset_slip, index_myFamily, index_otherFamily, &
@ -239,20 +236,11 @@ subroutine plastic_dislotwin_init(fileUnit)
integer(pInt) :: sizeState, sizeDotState
integer(pInt) :: NipcMyPhase
real(pReal), allocatable, dimension(:,:) :: temp1,temp2
integer(pInt), dimension(1,200), parameter :: lattice_ntranssystem = 12 ! HACK!!
integer(pInt), dimension(0), parameter :: emptyIntArray = [integer(pInt)::]
real(pReal), dimension(0), parameter :: emptyRealArray = [real(pReal)::]
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
type(tParameters) :: &
prm
type(tDislotwinState) :: &
stt, &
dot
type(tDislotwinMicrostructure) :: &
mse
integer(kind(undefined_ID)) :: &
outputID !< ID of each post result output
@ -293,42 +281,51 @@ subroutine plastic_dislotwin_init(fileUnit)
associate(prm => param(phase_plasticityInstance(p)), &
dot => dotState(phase_plasticityInstance(p)), &
stt => state(phase_plasticityInstance(p)), &
mse => microstructure(phase_plasticityInstance(p)))
mse => microstructure(phase_plasticityInstance(p)), &
config => config_phase(p))
! This data is read in already in lattice
prm%isFCC = merge(.true., .false., lattice_structure(p) == LATTICE_FCC_ID)
prm%mu = lattice_mu(p)
prm%nu = lattice_nu(p)
prm%C66 = lattice_C66(1:6,1:6,p)
structure = config_phase(p)%getString('lattice_structure')
structure = config%getString('lattice_structure')
!--------------------------------------------------------------------------------------------------
! slip related parameters
prm%Nslip = config_phase(p)%getInts('nslip',defaultVal=emptyIntArray)
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%totalNslip = sum(prm%Nslip)
slipActive: if (prm%totalNslip > 0_pInt) then
prm%fccTwinTransNucleation = merge(.true., .false., lattice_structure(p) == LATTICE_FCC_ID) &
.and. (prm%Nslip(1) == 12_pInt)
if(prm%fccTwinTransNucleation) &
prm%fcc_twinNucleationSlipPair = lattice_fcc_twinNucleationSlipPair
prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,structure(1:3),&
config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal))
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%forestProjection = lattice_forestProjection (prm%Nslip,structure(1:3),&
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, &
config_phase(p)%getFloats('interaction_slipslip'), &
config%getFloats('interaction_slipslip'), &
structure(1:3))
prm%rho0 = config_phase(p)%getFloats('rhoedge0', requiredShape=shape(prm%Nslip)) !ToDo: rename to rho_0
prm%rhoDip0 = config_phase(p)%getFloats('rhoedgedip0',requiredShape=shape(prm%Nslip)) !ToDo: rename to rho_dip_0
prm%v0 = config_phase(p)%getFloats('v0', requiredShape=shape(prm%Nslip))
prm%burgers_slip = config_phase(p)%getFloats('slipburgers',requiredShape=shape(prm%Nslip))
prm%Qedge = config_phase(p)%getFloats('qedge', requiredShape=shape(prm%Nslip)) !ToDo: rename (ask Karo)
prm%CLambdaSlip = config_phase(p)%getFloats('clambdaslip',requiredShape=shape(prm%Nslip))
prm%p = config_phase(p)%getFloats('p_slip', requiredShape=shape(prm%Nslip))
prm%q = config_phase(p)%getFloats('q_slip', requiredShape=shape(prm%Nslip))
prm%B = config_phase(p)%getFloats('b', requiredShape=shape(prm%Nslip), &
defaultVal=[(0.0_pReal, i=1,size(prm%Nslip))])
prm%tau_peierls = config_phase(p)%getFloats('tau_peierls',requiredShape=shape(prm%Nslip), &
prm%rho0 = config%getFloats('rhoedge0', requiredShape=shape(prm%Nslip)) !ToDo: rename to rho_0
prm%rhoDip0 = config%getFloats('rhoedgedip0',requiredShape=shape(prm%Nslip)) !ToDo: rename to rho_dip_0
prm%v0 = config%getFloats('v0', requiredShape=shape(prm%Nslip))
prm%burgers_slip = config%getFloats('slipburgers',requiredShape=shape(prm%Nslip))
prm%Qedge = config%getFloats('qedge', requiredShape=shape(prm%Nslip)) !ToDo: rename (ask Karo)
prm%CLambdaSlip = config%getFloats('clambdaslip',requiredShape=shape(prm%Nslip))
prm%p = config%getFloats('p_slip', requiredShape=shape(prm%Nslip))
prm%q = config%getFloats('q_slip', requiredShape=shape(prm%Nslip))
prm%B = config%getFloats('b', requiredShape=shape(prm%Nslip), &
defaultVal=[(0.0_pReal, i=1,size(prm%Nslip))])
prm%tau_peierls = config%getFloats('tau_peierls',requiredShape=shape(prm%Nslip), &
defaultVal=[(0.0_pReal, i=1,size(prm%Nslip))]) ! Deprecated
prm%CEdgeDipMinDistance = config_phase(p)%getFloat('cedgedipmindistance')
prm%CEdgeDipMinDistance = config%getFloat('cedgedipmindistance')
! expand: family => system
prm%rho0 = math_expand(prm%rho0, prm%Nslip)
@ -360,28 +357,33 @@ subroutine plastic_dislotwin_init(fileUnit)
!--------------------------------------------------------------------------------------------------
! twin related parameters
prm%Ntwin = config_phase(p)%getInts('ntwin', defaultVal=emptyIntArray)
prm%Ntwin = config%getInts('ntwin', defaultVal=emptyIntArray)
prm%totalNtwin = sum(prm%Ntwin)
if (prm%totalNtwin > 0_pInt) then
prm%Schmid_twin = lattice_SchmidMatrix_twin(prm%Ntwin,structure(1:3),&
config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal))
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%interaction_TwinTwin = lattice_interaction_TwinTwin(prm%Ntwin,&
config_phase(p)%getFloats('interaction_twintwin'), &
config%getFloats('interaction_twintwin'), &
structure(1:3))
prm%burgers_twin = config_phase(p)%getFloats('twinburgers')
prm%twinsize = config_phase(p)%getFloats('twinsize')
prm%r = config_phase(p)%getFloats('r_twin')
prm%burgers_twin = config%getFloats('twinburgers')
prm%twinsize = config%getFloats('twinsize')
prm%r = config%getFloats('r_twin')
prm%xc_twin = config_phase(p)%getFloat('xc_twin')
prm%L0_twin = config_phase(p)%getFloat('l0_twin')
prm%MaxTwinFraction = config_phase(p)%getFloat('maxtwinfraction') ! ToDo: only used in postResults
prm%Cthresholdtwin = config_phase(p)%getFloat('cthresholdtwin', defaultVal=0.0_pReal)
prm%Cmfptwin = config_phase(p)%getFloat('cmfptwin', defaultVal=0.0_pReal) ! ToDo: How to handle that???
prm%xc_twin = config%getFloat('xc_twin')
prm%L0_twin = config%getFloat('l0_twin')
prm%MaxTwinFraction = config%getFloat('maxtwinfraction') ! ToDo: only used in postResults
prm%Cthresholdtwin = config%getFloat('cthresholdtwin', defaultVal=0.0_pReal)
prm%Cmfptwin = config%getFloat('cmfptwin', defaultVal=0.0_pReal) ! ToDo: How to handle that???
prm%shear_twin = lattice_characteristicShear_Twin(prm%Ntwin,structure(1:3),&
config%getFloat('c/a',defaultVal=0.0_pReal))
if (.not. prm%isFCC) then
prm%Ndot0_twin = config_phase(p)%getFloats('ndot0_twin')
prm%C66_twin = lattice_C66_twin(prm%Ntwin,prm%C66,structure(1:3),&
config%getFloat('c/a',defaultVal=0.0_pReal))
if (.not. prm%fccTwinTransNucleation) then
prm%Ndot0_twin = config%getFloats('ndot0_twin')
prm%Ndot0_twin = math_expand(prm%Ndot0_twin,prm%Ntwin)
endif
@ -398,27 +400,42 @@ subroutine plastic_dislotwin_init(fileUnit)
!--------------------------------------------------------------------------------------------------
! transformation related parameters
prm%Ntrans = config_phase(p)%getInts('ntrans', defaultVal=emptyIntArray)
prm%Ntrans = config%getInts('ntrans', defaultVal=emptyIntArray)
prm%totalNtrans = sum(prm%Ntrans)
if (prm%totalNtrans > 0_pInt) then
prm%burgers_trans = config_phase(p)%getFloats('transburgers')
prm%burgers_trans = config%getFloats('transburgers')
prm%burgers_trans = math_expand(prm%burgers_trans,prm%Ntrans)
prm%Cthresholdtrans = config_phase(p)%getFloat('cthresholdtrans', defaultVal=0.0_pReal) ! ToDo: How to handle that???
prm%transStackHeight = config_phase(p)%getFloat('transstackheight', defaultVal=0.0_pReal) ! ToDo: How to handle that???
prm%Cmfptrans = config_phase(p)%getFloat('cmfptrans', defaultVal=0.0_pReal) ! ToDo: How to handle that???
prm%deltaG = config_phase(p)%getFloat('deltag')
prm%xc_trans = config_phase(p)%getFloat('xc_trans', defaultVal=0.0_pReal) ! ToDo: How to handle that???
prm%L0_trans = config_phase(p)%getFloat('l0_trans')
prm%Cthresholdtrans = config%getFloat('cthresholdtrans', defaultVal=0.0_pReal) ! ToDo: How to handle that???
prm%transStackHeight = config%getFloat('transstackheight', defaultVal=0.0_pReal) ! ToDo: How to handle that???
prm%Cmfptrans = config%getFloat('cmfptrans', defaultVal=0.0_pReal) ! ToDo: How to handle that???
prm%deltaG = config%getFloat('deltag')
prm%xc_trans = config%getFloat('xc_trans', defaultVal=0.0_pReal) ! ToDo: How to handle that???
prm%L0_trans = config%getFloat('l0_trans')
prm%interaction_TransTrans = lattice_interaction_TransTrans(prm%Ntrans,&
config%getFloats('interaction_transtrans'), &
structure(1:3))
prm%C66_trans = lattice_C66_trans(prm%Ntrans,prm%C66, &
config%getString('trans_lattice_structure'), &
0.0_pReal, &
config%getFloat('a_bcc', defaultVal=0.0_pReal), &
config%getFloat('a_fcc', defaultVal=0.0_pReal))
prm%Schmid_trans = lattice_SchmidMatrix_trans(prm%Ntrans, &
config%getString('trans_lattice_structure'), &
0.0_pReal, &
config%getFloat('a_bcc', defaultVal=0.0_pReal), &
config%getFloat('a_fcc', defaultVal=0.0_pReal))
prm%interaction_TransTrans = spread(config_phase(p)%getFloats('interaction_transtrans'),2,1)
if (lattice_structure(p) /= LATTICE_fcc_ID) then
prm%Ndot0_trans = config_phase(p)%getFloats('ndot0_trans')
prm%Ndot0_trans = config%getFloats('ndot0_trans')
prm%Ndot0_trans = math_expand(prm%Ndot0_trans,prm%Ntrans)
endif
prm%lamellarsizePerTransSystem = config_phase(p)%getFloats('lamellarsize')
prm%lamellarsizePerTransSystem = config%getFloats('lamellarsize')
prm%lamellarsizePerTransSystem = math_expand(prm%lamellarsizePerTransSystem,prm%Ntrans)
prm%s = config_phase(p)%getFloats('s_trans',defaultVal=[0.0_pReal])
prm%s = config%getFloats('s_trans',defaultVal=[0.0_pReal])
prm%s = math_expand(prm%s,prm%Ntrans)
else
allocate(prm%lamellarsizePerTransSystem(0))
@ -426,45 +443,48 @@ subroutine plastic_dislotwin_init(fileUnit)
endif
if (sum(prm%Ntwin) > 0_pInt .or. prm%totalNtrans > 0_pInt) then
prm%SFE_0K = config_phase(p)%getFloat('sfe_0k')
prm%dSFE_dT = config_phase(p)%getFloat('dsfe_dt')
prm%VcrossSlip = config_phase(p)%getFloat('vcrossslip')
prm%SFE_0K = config%getFloat('sfe_0k')
prm%dSFE_dT = config%getFloat('dsfe_dt')
prm%VcrossSlip = config%getFloat('vcrossslip')
endif
if (prm%totalNslip > 0_pInt .and. prm%totalNtwin > 0_pInt) then
prm%interaction_SlipTwin = lattice_interaction_SlipTwin(prm%Nslip,prm%Ntwin,&
config_phase(p)%getFloats('interaction_sliptwin'), &
config%getFloats('interaction_sliptwin'), &
structure(1:3))
prm%interaction_TwinSlip = lattice_interaction_TwinSlip(prm%Ntwin,prm%Nslip,&
config_phase(p)%getFloats('interaction_twinslip'), &
config%getFloats('interaction_twinslip'), &
structure(1:3))
if (prm%fccTwinTransNucleation .and. prm%totalNtwin > 12_pInt) write(6,*) 'mist' ! ToDo: implement better test. The model will fail also if ntwin is [6,6]
endif
if (prm%totalNslip > 0_pInt .and. prm%totalNtrans > 0_pInt) then
prm%interaction_TransSlip = spread(config_phase(p)%getFloats('interaction_transslip'),2,1)
prm%interaction_SlipTrans = spread(config_phase(p)%getFloats('interaction_sliptrans'),2,1)
prm%interaction_SlipTrans = lattice_interaction_SlipTrans(prm%Nslip,prm%Ntrans,&
config%getFloats('interaction_sliptrans'), &
structure(1:3))
if (prm%fccTwinTransNucleation .and. prm%totalNtrans > 12_pInt) write(6,*) 'mist' ! ToDo: implement better test. The model will fail also if ntrans is [6,6]
endif
prm%aTolRho = config_phase(p)%getFloat('atol_rho', defaultVal=0.0_pReal)
prm%aTolTwinFrac = config_phase(p)%getFloat('atol_twinfrac', defaultVal=0.0_pReal)
prm%aTolTransFrac = config_phase(p)%getFloat('atol_transfrac', defaultVal=0.0_pReal)
prm%aTolRho = config%getFloat('atol_rho', defaultVal=0.0_pReal)
prm%aTolTwinFrac = config%getFloat('atol_twinfrac', defaultVal=0.0_pReal)
prm%aTolTransFrac = config%getFloat('atol_transfrac', defaultVal=0.0_pReal)
prm%CAtomicVolume = config_phase(p)%getFloat('catomicvolume')
prm%GrainSize = config_phase(p)%getFloat('grainsize')
prm%CAtomicVolume = config%getFloat('catomicvolume')
prm%GrainSize = config%getFloat('grainsize')
prm%D0 = config_phase(p)%getFloat('d0')
prm%Qsd = config_phase(p)%getFloat('qsd')
prm%SolidSolutionStrength = config_phase(p)%getFloat('solidsolutionstrength')
if (config_phase(p)%keyExists('dipoleformationfactor')) call IO_error(1,ext_msg='use /nodipoleformation/')
prm%dipoleformation = .not. config_phase(p)%keyExists('/nodipoleformation/')
prm%sbVelocity = config_phase(p)%getFloat('shearbandvelocity',defaultVal=0.0_pReal)
prm%D0 = config%getFloat('d0')
prm%Qsd = config%getFloat('qsd')
prm%SolidSolutionStrength = config%getFloat('solidsolutionstrength') ! Deprecated
if (config%keyExists('dipoleformationfactor')) call IO_error(1,ext_msg='use /nodipoleformation/')
prm%dipoleformation = .not. config%keyExists('/nodipoleformation/')
prm%sbVelocity = config%getFloat('shearbandvelocity',defaultVal=0.0_pReal)
if (prm%sbVelocity > 0.0_pReal) then
prm%sbResistance = config_phase(p)%getFloat('shearbandresistance')
prm%sbQedge = config_phase(p)%getFloat('qedgepersbsystem')
prm%pShearBand = config_phase(p)%getFloat('p_shearband')
prm%qShearBand = config_phase(p)%getFloat('q_shearband')
prm%sbResistance = config%getFloat('shearbandresistance')
prm%sbQedge = config%getFloat('qedgepersbsystem')
prm%pShearBand = config%getFloat('p_shearband')
prm%qShearBand = config%getFloat('q_shearband')
endif
!if (Ndot0PerTwinFamily(f,p) < 0.0_pReal) &
@ -505,7 +525,7 @@ subroutine plastic_dislotwin_init(fileUnit)
prm%qShearBand <= 0.0_pReal) &
call IO_error(211_pInt,el=p,ext_msg='qShearBand ('//PLASTICITY_DISLOTWIN_label//')')
outputs = config_phase(p)%getStrings('(output)', defaultVal=emptyStringArray)
outputs = config%getStrings('(output)', defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i= 1_pInt, size(outputs)
outputID = undefined_ID
@ -600,96 +620,6 @@ subroutine plastic_dislotwin_init(fileUnit)
plasticState(p)%accumulatedSlip => &
plasticState(p)%state (offset_slip+1:offset_slip+plasticState(p)%nslip,1:NipcMyPhase)
allocate(temp1(prm%totalNslip,prm%totalNtrans),source =0.0_pReal)
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
do o = 1_pInt,size(prm%Ntrans,1)
index_otherFamily = sum(prm%Ntrans(1:o-1_pInt))
do k = 1_pInt,prm%Ntrans(o) ! loop over (active) systems in other family (trans)
temp1(index_myFamily+j,index_otherFamily+k) = &
prm%interaction_SlipTrans(lattice_interactionSlipTrans( &
sum(lattice_NslipSystem(1:f-1_pInt,p))+j, &
sum(lattice_NtransSystem(1:o-1_pInt,p))+k, &
p),1 )
enddo; enddo
enddo slipSystemsLoop
enddo mySlipFamilies
prm%interaction_SlipTrans = temp1; deallocate(temp1)
allocate(prm%C66_twin(6,6,prm%totalNtwin), source=0.0_pReal)
if (lattice_structure(p) == LATTICE_fcc_ID) &
allocate(prm%fcc_twinNucleationSlipPair(2,prm%totalNtwin),source = 0_pInt)
allocate(prm%shear_twin(prm%totalNtwin),source = 0.0_pReal)
i = 0_pInt
twinFamiliesLoop: do f = 1_pInt, size(prm%Ntwin,1)
index_myFamily = sum(prm%Ntwin(1:f-1_pInt)) ! index in truncated twin system list
twinSystemsLoop: do j = 1_pInt,prm%Ntwin(f)
i = i + 1_pInt
prm%shear_twin(i) = lattice_shearTwin(sum(lattice_Ntwinsystem(1:f-1,p))+j,p)
if (lattice_structure(p) == LATTICE_fcc_ID) prm%fcc_twinNucleationSlipPair(1:2,i) = &
lattice_fcc_twinNucleationSlipPair(1:2,sum(lattice_Ntwinsystem(1:f-1,p))+j)
!* Rotate twin elasticity matrices
index_otherFamily = sum(lattice_NtwinSystem(1:f-1_pInt,p)) ! index in full lattice twin list
prm%C66_twin(1:6,1:6,index_myFamily+j) = &
math_Mandel3333to66(math_rotate_forward3333(lattice_C3333(1:3,1:3,1:3,1:3,p),&
lattice_Qtwin(1:3,1:3,index_otherFamily+j,p)))
enddo twinSystemsLoop
enddo twinFamiliesLoop
allocate(temp1(prm%totalNtrans,prm%totalNslip), source =0.0_pReal)
allocate(temp2(prm%totalNtrans,prm%totalNtrans), source =0.0_pReal)
allocate(prm%C66_trans(6,6,prm%totalNtrans) ,source=0.0_pReal)
allocate(prm%Schmid_trans(3,3,prm%totalNtrans),source = 0.0_pReal)
i = 0_pInt
transFamiliesLoop: do f = 1_pInt,size(prm%Ntrans,1)
index_myFamily = sum(prm%Ntrans(1:f-1_pInt)) ! index in truncated trans system list
transSystemsLoop: do j = 1_pInt,prm%Ntrans(f)
i = i + 1_pInt
prm%Schmid_trans(1:3,1:3,i) = lattice_Strans(1:3,1:3,sum(lattice_Ntranssystem(1:f-1,p))+j,p)
!* Rotate trans elasticity matrices
index_otherFamily = sum(lattice_NtransSystem(1:f-1_pInt,p)) ! index in full lattice trans list
prm%C66_trans(1:6,1:6,index_myFamily+j) = &
math_Mandel3333to66(math_rotate_forward3333(lattice_trans_C3333(1:3,1:3,1:3,1:3,p),&
lattice_Qtrans(1:3,1:3,index_otherFamily+j,p)))
!* Interaction matrices
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)
temp1(index_myFamily+j,index_otherFamily+k) = &
prm%interaction_TransSlip(lattice_interactionTransSlip( &
sum(lattice_NtransSystem(1:f-1_pInt,p))+j, &
sum(lattice_NslipSystem(1:o-1_pInt,p))+k, &
p) ,1 )
enddo; enddo
do o = 1_pInt,size(prm%Ntrans,1)
index_otherFamily = sum(prm%Ntrans(1:o-1_pInt))
do k = 1_pInt,prm%Ntrans(o) ! loop over (active) systems in other family (trans)
temp2(index_myFamily+j,index_otherFamily+k) = &
prm%interaction_TransTrans(lattice_interactionTransTrans( &
sum(lattice_NtransSystem(1:f-1_pInt,p))+j, &
sum(lattice_NtransSystem(1:o-1_pInt,p))+k, &
p),1 )
enddo; enddo
enddo transSystemsLoop
enddo transFamiliesLoop
prm%interaction_TransSlip = temp1; deallocate(temp1)
prm%interaction_TransTrans = temp2; deallocate(temp2)
startIndex=1_pInt
endIndex=prm%totalNslip
@ -738,25 +668,23 @@ subroutine plastic_dislotwin_init(fileUnit)
plasticState(p)%state0 = plasticState(p)%state
dot%whole => plasticState(p)%dotState
allocate(mse%invLambdaSlip(prm%totalNslip,NipcMyPhase),source=0.0_pReal)
allocate(mse%invLambdaSlipTwin(prm%totalNslip,NipcMyPhase),source=0.0_pReal)
allocate(mse%invLambdaTwin(prm%totalNtwin,NipcMyPhase),source=0.0_pReal)
allocate(mse%invLambdaSlipTrans(prm%totalNtrans,NipcMyPhase),source=0.0_pReal)
allocate(mse%invLambdaTrans(prm%totalNtrans,NipcMyPhase),source=0.0_pReal)
allocate(mse%invLambdaSlip (prm%totalNslip, NipcMyPhase),source=0.0_pReal)
allocate(mse%invLambdaSlipTwin (prm%totalNslip, NipcMyPhase),source=0.0_pReal)
allocate(mse%invLambdaSlipTrans (prm%totalNslip, NipcMyPhase),source=0.0_pReal)
allocate(mse%mfp_slip (prm%totalNslip, NipcMyPhase),source=0.0_pReal)
allocate(mse%threshold_stress_slip (prm%totalNslip, NipcMyPhase),source=0.0_pReal)
allocate(mse%mfp_slip(prm%totalNslip,NipcMyPhase), source=0.0_pReal)
allocate(mse%mfp_twin(prm%totalNtwin,NipcMyPhase), source=0.0_pReal)
allocate(mse%mfp_trans(prm%totalNtrans,NipcMyPhase),source=0.0_pReal)
allocate(mse%invLambdaTwin (prm%totalNtwin, NipcMyPhase),source=0.0_pReal)
allocate(mse%mfp_twin (prm%totalNtwin, NipcMyPhase),source=0.0_pReal)
allocate(mse%threshold_stress_twin (prm%totalNtwin, NipcMyPhase),source=0.0_pReal)
allocate(mse%tau_r_twin (prm%totalNtwin, NipcMyPhase),source=0.0_pReal)
allocate(mse%twinVolume (prm%totalNtwin, NipcMyPhase),source=0.0_pReal)
allocate(mse%threshold_stress_slip(prm%totalNslip,NipcMyPhase), source=0.0_pReal)
allocate(mse%threshold_stress_twin(prm%totalNtwin,NipcMyPhase), source=0.0_pReal)
allocate(mse%invLambdaTrans (prm%totalNtrans,NipcMyPhase),source=0.0_pReal)
allocate(mse%mfp_trans (prm%totalNtrans,NipcMyPhase),source=0.0_pReal)
allocate(mse%threshold_stress_trans(prm%totalNtrans,NipcMyPhase),source=0.0_pReal)
allocate(mse%tau_r_twin(prm%totalNtwin,NipcMyPhase), source=0.0_pReal)
allocate(mse%tau_r_trans(prm%totalNtrans,NipcMyPhase), source=0.0_pReal)
allocate(mse%twinVolume(prm%totalNtwin,NipcMyPhase), source=0.0_pReal)
allocate(mse%martensiteVolume(prm%totalNtrans,NipcMyPhase), source=0.0_pReal)
allocate(mse%tau_r_trans (prm%totalNtrans,NipcMyPhase),source=0.0_pReal)
allocate(mse%martensiteVolume (prm%totalNtrans,NipcMyPhase),source=0.0_pReal)
end associate
enddo
@ -779,8 +707,6 @@ function plastic_dislotwin_homogenizedC(ipc,ip,el)
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
type(tParameters) :: prm
type(tDislotwinState) :: stt
integer(pInt) :: i, &
of
@ -838,10 +764,6 @@ subroutine plastic_dislotwin_microstructure(temperature,ipc,ip,el)
fOverStacksize, &
ftransOverLamellarSize
type(tParameters) :: prm !< parameters of present instance
type(tDislotwinState) :: stt !< state of present instance
type(tDislotwinMicrostructure) :: mse
of = phasememberAt(ipc,ip,el)
associate(prm => param(phase_plasticityInstance(material_phase(ipc,ip,el))),&
@ -864,7 +786,7 @@ subroutine plastic_dislotwin_microstructure(temperature,ipc,ip,el)
forall (i = 1_pInt:prm%totalNslip) &
mse%invLambdaSlip(i,of) = &
sqrt(dot_product((stt%rhoEdge(1_pInt:prm%totalNslip,of)+stt%rhoEdgeDip(1_pInt:prm%totalNslip,of)),&
prm%forestProjectionEdge(1:prm%totalNslip,i)))/prm%CLambdaSlip(i)
prm%forestProjection(1:prm%totalNslip,i)))/prm%CLambdaSlip(i)
!* 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation
!$OMP CRITICAL (evilmatmul)
@ -881,7 +803,7 @@ subroutine plastic_dislotwin_microstructure(temperature,ipc,ip,el)
!* 1/mean free distance between 2 martensite lamellar from different systems seen by a moving dislocation
if (prm%totalNtrans > 0_pInt .and. prm%totalNslip > 0_pInt) &
mse%invLambdaSlipTrans(1_pInt:prm%totalNslip,of) = &
mse%invLambdaSlipTrans(1_pInt:prm%totalNslip,of) = & ! ToDo: does not work if Ntrans is not 12
matmul(prm%interaction_SlipTrans,ftransOverLamellarSize)/(1.0_pReal-sumf_trans)
!* 1/mean free distance between 2 martensite stacks from different systems seen by a growing martensite (1/lambda_trans)
@ -953,10 +875,6 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,Temperature,instance,
math_symmetric33, &
math_mul33xx33, &
math_mul33x3
use material, only: &
material_phase, &
phase_plasticityInstance, &
phasememberAt
implicit none
real(pReal), dimension(3,3), intent(out) :: Lp
@ -999,9 +917,6 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,Temperature,instance,
0, 1, 1 &
],pReal),[ 3,6])
type(tParameters) :: prm !< parameters of present instance
type(tDislotwinState) :: ste !< state of present instance
associate(prm => param(instance), stt => state(instance), mse => microstructure(instance))
f_unrotated = 1.0_pReal &
@ -1068,7 +983,7 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,Temperature,instance,
significantTransStress: if (tau > tol_math_check) then
StressRatio_s = (mse%threshold_stress_trans(i,of)/tau)**prm%s(i)
isFCCtrans: if (prm%isFCC) then
isFCCtrans: if (prm%fccTwinTransNucleation) then
s1=prm%fcc_twinNucleationSlipPair(1,i)
s2=prm%fcc_twinNucleationSlipPair(2,i)
if (tau < mse%tau_r_trans(i,of)) then
@ -1112,10 +1027,7 @@ subroutine plastic_dislotwin_dotState(Mp,Temperature,instance,of)
math_Mandel6to33, &
pi
use material, only: &
material_phase, &
phase_plasticityInstance, &
plasticState, &
phasememberAt
plasticState
implicit none
real(pReal), dimension(3,3), intent(in):: &
@ -1135,12 +1047,6 @@ subroutine plastic_dislotwin_dotState(Mp,Temperature,instance,of)
real(pReal), dimension(plasticState(instance)%Nslip) :: &
gdot_slip
type(tParameters) :: prm
type(tDislotwinState) :: stt, dot
type(tDislotwinMicrostructure) :: mse
associate(prm => param(instance), stt => state(instance), &
dot => dotstate(instance), mse => microstructure(instance))
@ -1206,7 +1112,7 @@ subroutine plastic_dislotwin_dotState(Mp,Temperature,instance,of)
significantTwinStress: if (tau > tol_math_check) then
StressRatio_r = (mse%threshold_stress_twin(i,of)/tau)**prm%r(i)
isFCCtwin: if (prm%isFCC) then
isFCCtwin: if (prm%fccTwinTransNucleation) then
s1=prm%fcc_twinNucleationSlipPair(1,i)
s2=prm%fcc_twinNucleationSlipPair(2,i)
if (tau < mse%tau_r_twin(i,of)) then
@ -1232,7 +1138,7 @@ subroutine plastic_dislotwin_dotState(Mp,Temperature,instance,of)
significantTransStress: if (tau > tol_math_check) then
StressRatio_s = (mse%threshold_stress_trans(i,of)/tau)**prm%s(i)
isFCCtrans: if (prm%isFCC) then
isFCCtrans: if (prm%fccTwinTransNucleation) then
s1=prm%fcc_twinNucleationSlipPair(1,i)
s2=prm%fcc_twinNucleationSlipPair(2,i)
if (tau < mse%tau_r_trans(i,of)) then
@ -1375,7 +1281,7 @@ pure subroutine kinetics_twin(prm,stt,mse,of,Mp,temperature,gdot_slip,gdot_twin,
do i = 1_pInt, prm%totalNtwin
tau(i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i))
isFCC: if (prm%isFCC) then
isFCC: if (prm%fccTwinTransNucleation) then
s1=prm%fcc_twinNucleationSlipPair(1,i)
s2=prm%fcc_twinNucleationSlipPair(2,i)
if (tau(i) < mse%tau_r_twin(i,of)) then
@ -1447,7 +1353,7 @@ pure subroutine kinetics_trans(prm,stt,mse,of,Mp,temperature,gdot_slip,gdot_tran
do i = 1_pInt, prm%totalNtrans
tau(i) = math_mul33xx33(Mp,prm%Schmid_trans(1:3,1:3,i))
isFCC: if (prm%isFCC) then
isFCC: if (prm%fccTwinTransNucleation) then
s1=prm%fcc_twinNucleationSlipPair(1,i)
s2=prm%fcc_twinNucleationSlipPair(2,i)
if (tau(i) < mse%tau_r_trans(i,of)) then
@ -1500,11 +1406,6 @@ function plastic_dislotwin_postResults(Mp,Temperature,instance,of) result(postRe
PI, &
math_mul33xx33, &
math_Mandel6to33
use material, only: &
material_phase, &
plasticState, &
phase_plasticityInstance,&
phasememberAt
implicit none
real(pReal), dimension(3,3),intent(in) :: &
@ -1622,7 +1523,7 @@ function plastic_dislotwin_postResults(Mp,Temperature,instance,of) result(postRe
tau = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,j))
if ( tau > 0.0_pReal ) then
isFCCtwin: if (prm%isFCC) then
isFCCtwin: if (prm%fccTwinTransNucleation) then
s1=prm%fcc_twinNucleationSlipPair(1,j)
s2=prm%fcc_twinNucleationSlipPair(2,j)
if (tau < mse%tau_r_twin(j,of)) then

9
src/plastic_isotropic.f90
View File

@ -48,16 +48,17 @@ module plastic_isotropic
outputID
logical :: &
dilatation
end type
end type tParameters
type, private :: tIsotropicState
real(pReal), pointer, dimension(:) :: &
flowstress, &
accumulatedShear
end type
end type tIsotropicState
type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance)
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
type(tParameters), allocatable, dimension(:), private :: param
type(tIsotropicState), allocatable, dimension(:), private :: &
dotState, &
state

34
src/plastic_kinematichardening.f90
View File

@ -58,7 +58,7 @@ module plastic_kinehardening
Nslip !< number of active slip systems for each family
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID !< ID of each post result output
end type
end type tParameters
type, private :: tKinehardeningState
real(pReal), pointer, dimension(:,:) :: & !< vectors along NipcMyInstance
@ -68,10 +68,11 @@ module plastic_kinehardening
chi0, & !< backstress at last switch of stress sense
gamma0, & !< accumulated shear at last switch of stress sense
accshear !< accumulated (absolute) shear
end type
end type tKinehardeningState
type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance)
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
type(tParameters), allocatable, dimension(:), private :: param
type(tKinehardeningState), allocatable, dimension(:), private :: &
dotState, &
deltaState, &
@ -150,7 +151,6 @@ subroutine plastic_kinehardening_init
outputID
character(len=pStringLen) :: &
structure = '',&
extmsg = ''
character(len=65536), dimension(:), allocatable :: &
outputs
@ -186,8 +186,6 @@ subroutine plastic_kinehardening_init
endif
#endif
structure = config%getString('lattice_structure')
!--------------------------------------------------------------------------------------------------
! optional parameters that need to be defined
prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal)
@ -202,9 +200,10 @@ subroutine plastic_kinehardening_init
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%totalNslip = sum(prm%Nslip)
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))
if(structure=='bcc') then
if(trim(config%getString('lattice_structure')) == 'bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt)
@ -215,15 +214,15 @@ subroutine plastic_kinehardening_init
endif
prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, &
config%getFloats('interaction_slipslip'), &
structure(1:3))
config%getString('lattice_structure'))
prm%crss0 = config%getFloats('crss0', requiredShape=shape(prm%Nslip))
prm%tau1 = config%getFloats('tau1', requiredShape=shape(prm%Nslip))
prm%tau1_b = config%getFloats('tau1_b', requiredShape=shape(prm%Nslip))
prm%theta0 = config%getFloats('theta0', requiredShape=shape(prm%Nslip))
prm%theta1 = config%getFloats('theta1', requiredShape=shape(prm%Nslip))
prm%theta0_b = config%getFloats('theta0_b', requiredShape=shape(prm%Nslip))
prm%theta1_b = config%getFloats('theta1_b', requiredShape=shape(prm%Nslip))
prm%crss0 = config%getFloats('crss0', requiredSize=size(prm%Nslip))
prm%tau1 = config%getFloats('tau1', requiredSize=size(prm%Nslip))
prm%tau1_b = config%getFloats('tau1_b', requiredSize=size(prm%Nslip))
prm%theta0 = config%getFloats('theta0', requiredSize=size(prm%Nslip))
prm%theta1 = config%getFloats('theta1', requiredSize=size(prm%Nslip))
prm%theta0_b = config%getFloats('theta0_b', requiredSize=size(prm%Nslip))
prm%theta1_b = config%getFloats('theta1_b', requiredSize=size(prm%Nslip))
prm%gdot0 = config%getFloat('gdot0')
prm%n = config%getFloat('n_slip')
@ -301,7 +300,6 @@ subroutine plastic_kinehardening_init
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState, &
prm%totalNslip,0_pInt,0_pInt)
plasticState(p)%sizePostResults = sum(plastic_kinehardening_sizePostResult(:,phase_plasticityInstance(p)))
plasticState(p)%offsetDeltaState = sizeDotState
!--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and aTolState

31
src/plastic_phenopowerlaw.f90
View File

@ -72,7 +72,7 @@ module plastic_phenopowerlaw
Ntwin !< number of active twin systems for each family
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID !< ID of each post result output
end type !< container type for internal constitutive parameters
end type tParameters
type, private :: tPhenopowerlawState
real(pReal), pointer, dimension(:,:) :: &
@ -80,10 +80,11 @@ module plastic_phenopowerlaw
xi_twin, &
gamma_slip, &
gamma_twin
end type
end type tPhenopowerlawState
type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance)
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
type(tParameters), allocatable, dimension(:), private :: param
type(tPhenopowerlawState), allocatable, dimension(:), private :: &
dotState, &
state
@ -152,7 +153,6 @@ subroutine plastic_phenopowerlaw_init
outputID
character(len=pStringLen) :: &
structure = '',&
extmsg = ''
character(len=65536), dimension(:), allocatable :: &
outputs
@ -180,8 +180,6 @@ subroutine plastic_phenopowerlaw_init
stt => state(phase_plasticityInstance(p)), &
config => config_phase(p))
structure = config%getString('lattice_structure')
!--------------------------------------------------------------------------------------------------
! optional parameters that need to be defined
prm%twinB = config%getFloat('twin_b',defaultVal=1.0_pReal)
@ -203,9 +201,10 @@ subroutine plastic_phenopowerlaw_init
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%totalNslip = sum(prm%Nslip)
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))
if(structure=='bcc') then
if(trim(config%getString('lattice_structure')) == 'bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt)
@ -216,7 +215,7 @@ subroutine plastic_phenopowerlaw_init
endif
prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, &
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_sat = config%getFloats('tausat_slip', requiredSize=size(prm%Nslip))
@ -238,7 +237,7 @@ subroutine plastic_phenopowerlaw_init
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 (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
allocate(prm%interaction_SlipSlip(0,0))
allocate(prm%xi_slip_0(0))
@ -249,12 +248,12 @@ subroutine plastic_phenopowerlaw_init
prm%Ntwin = config%getInts('ntwin', defaultVal=emptyIntArray)
prm%totalNtwin = sum(prm%Ntwin)
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))
prm%interaction_TwinTwin = lattice_interaction_TwinTwin(prm%Ntwin,&
config%getFloats('interaction_twintwin'), &
structure(1:3))
prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,structure(1:3),&
config%getString('lattice_structure'))
prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,config%getString('lattice_structure'),&
config%getFloat('c/a'))
prm%xi_twin_0 = config%getFloats('tau0_twin',requiredSize=size(prm%Ntwin))
@ -281,10 +280,10 @@ subroutine plastic_phenopowerlaw_init
slipAndTwinActive: if (prm%totalNslip > 0_pInt .and. prm%totalNtwin > 0_pInt) then
prm%interaction_SlipTwin = lattice_interaction_SlipTwin(prm%Nslip,prm%Ntwin,&
config%getFloats('interaction_sliptwin'), &
structure(1:3))
config%getString('lattice_structure'))
prm%interaction_TwinSlip = lattice_interaction_TwinSlip(prm%Ntwin,prm%Nslip,&
config%getFloats('interaction_twinslip'), &
structure(1:3))
config%getString('lattice_structure'))
else slipAndTwinActive
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