Merge branch 'development' into 20-NewStyleDislotwin

This commit is contained in:
Martin Diehl 2018-10-01 21:16:11 +02:00
commit ddfbe19f02
13 changed files with 630 additions and 778 deletions

@ -1 +1 @@
Subproject commit df3b567c90dbcb291a5c2563d3984fce346cb55b
Subproject commit 02a391e38d1a44a8cfd2f169687a9d9f314a8146

View File

@ -1 +1 @@
v2.0.2-558-gdfba3c9b
v2.0.2-591-ga00d15b8

View File

@ -5,15 +5,10 @@
(output) orientation # quaternion
(output) eulerangles # orientation as Bunge triple in degree
(output) grainrotation # deviation from initial orientation as axis (1-3) and angle in degree (4) in crystal reference coordinates
(output) grainrotationx # deviation from initial orientation as angle in degrees around sample reference x axis
(output) grainrotationy # deviation from initial orientation as angle in degrees around sample reference y axis
(output) grainrotationz # deviation from initial orientation as angle in degrees around sample reference z axis
(output) f # deformation gradient tensor; synonyms: "defgrad"
(output) f # deformation gradient tensor
(output) fe # elastic deformation gradient tensor
(output) fp # plastic deformation gradient tensor
(output) e # total strain as Green-Lagrange tensor
(output) ee # elastic strain as Green-Lagrange tensor
(output) p # first Piola-Kichhoff stress tensor; synonyms: "firstpiola", "1stpiola"
(output) s # second Piola-Kichhoff stress tensor; synonyms: "tstar", "secondpiola", "2ndpiola"
(output) p # first Piola-Kichhoff stress tensor
(output) s # second Piola-Kichhoff stress tensor
(output) lp # plastic velocity gradient tensor
(output) elasmatrix # elastic stiffness matrix

View File

@ -18,8 +18,6 @@ mech none
(output) f # deformation gradient tensor; synonyms: "defgrad"
(output) fe # elastic deformation gradient tensor
(output) fp # plastic deformation gradient tensor
(output) e # total strain as Green-Lagrange tensor
(output) ee # elastic strain as Green-Lagrange tensor
(output) p # first Piola-Kichhoff stress tensor; synonyms: "firstpiola", "1stpiola"
(output) lp # plastic velocity gradient tensor
@ -56,7 +54,6 @@ h0_slipslip 75e6
interaction_slipslip 1 1 1.4 1.4 1.4 1.4
atol_resistance 1
#-------------------#
<microstructure>
#-------------------#

View File

@ -1,8 +1,8 @@
#! /usr/bin/env bash
if [ $1x != 3to2x ]; then
echo 'python2.7 to python'
find . -name '*.py' | xargs sed -i 's/usr\/bin\/env python2.7/usr\/bin\/env python/g'
find . -name '*.py' -type f | xargs sed -i 's/usr\/bin\/env python2.7/usr\/bin\/env python/g'
else
echo 'python to python2.7'
find . -name '*.py' | xargs sed -i 's/usr\/bin\/env python/usr\/bin\/env python2.7/g'
find . -name '*.py' -type f | xargs sed -i 's/usr\/bin\/env python/usr\/bin\/env python2.7/g'
fi

View File

@ -1,14 +1,14 @@
# -*- coding: UTF-8 no BOM -*-
import os,sys
import os, sys
import numpy as np
# ------------------------------------------------------------------
# python 3 has no unicode object, this ensures that the code works on Python 2&3
try:
test=isinstance('test', unicode)
except(NameError):
unicode=str
test = isinstance('test', unicode)
except NameError:
unicode = str
# ------------------------------------------------------------------
class ASCIItable():
@ -63,16 +63,17 @@ class ASCIItable():
x):
try:
return float(x)
except:
except ValueError:
return 0.0
# ------------------------------------------------------------------
def _removeCRLF(self,
string):
string):
"""Delete any carriage return and line feed from string"""
try:
return string.replace('\n','').replace('\r','')
except:
return string
except AttributeError:
return str(string)
# ------------------------------------------------------------------
@ -93,22 +94,22 @@ class ASCIItable():
# ------------------------------------------------------------------
def input_close(self):
try:
# try:
if self.__IO__['in'] != sys.stdin: self.__IO__['in'].close()
except:
pass
# except:
# pass
# ------------------------------------------------------------------
def output_write(self,
what):
"""Aggregate a single row (string) or list of (possibly containing further lists of) rows into output"""
if not isinstance(what, (str, unicode)):
if isinstance(what, (str, unicode)):
self.__IO__['output'] += [what]
else:
try:
for item in what: self.output_write(item)
except:
except TypeError:
self.__IO__['output'] += [str(what)]
else:
self.__IO__['output'] += [what]
return self.__IO__['buffered'] or self.output_flush()
@ -129,10 +130,10 @@ class ASCIItable():
# ------------------------------------------------------------------
def output_close(self,
dismiss = False):
try:
if self.__IO__['out'] != sys.stdout: self.__IO__['out'].close()
except:
pass
# try:
if self.__IO__['out'] != sys.stdout: self.__IO__['out'].close()
# except:
# pass
if dismiss and os.path.isfile(self.__IO__['out'].name):
os.remove(self.__IO__['out'].name)
elif self.__IO__['inPlace']:
@ -150,7 +151,7 @@ class ASCIItable():
try:
self.__IO__['in'].seek(0)
except:
except IOError:
pass
firstline = self.__IO__['in'].readline().strip()
@ -170,7 +171,7 @@ class ASCIItable():
else: # other table format
try:
self.__IO__['in'].seek(0) # try to rewind
except:
except IOError:
self.__IO__['readBuffer'] = [firstline] # or at least save data in buffer
while self.data_read(advance = False, respectLabels = False):
@ -197,7 +198,9 @@ class ASCIItable():
"""Write current header information (info + labels)"""
head = ['{}\theader'.format(len(self.info)+self.__IO__['labeled'])] if header else []
head.append(self.info)
if self.__IO__['labeled']: head.append('\t'.join(map(self._quote,self.tags)))
if self.__IO__['labeled']:
head.append('\t'.join(map(self._quote,self.tags)))
if len(self.tags) == 0: raise ValueError('no labels present.')
return self.output_write(head)
@ -257,13 +260,13 @@ class ASCIItable():
what,
reset = False):
"""Add item or list to existing set of labels (and switch on labeling)"""
if not isinstance(what, (str, unicode)):
if isinstance(what, (str, unicode)):
self.tags += [self._removeCRLF(what)]
else:
try:
for item in what: self.labels_append(item)
except:
except TypeError:
self.tags += [self._removeCRLF(str(what))]
else:
self.tags += [self._removeCRLF(what)]
self.__IO__['labeled'] = True # switch on processing (in particular writing) of tags
if reset: self.__IO__['tags'] = list(self.tags) # subsequent data_read uses current tags as data size
@ -410,13 +413,13 @@ class ASCIItable():
def info_append(self,
what):
"""Add item or list to existing set of infos"""
if not isinstance(what, (str, unicode)):
if isinstance(what, (str, unicode)):
self.info += [self._removeCRLF(what)]
else:
try:
for item in what: self.info_append(item)
except:
except TypeError:
self.info += [self._removeCRLF(str(what))]
else:
self.info += [self._removeCRLF(what)]
# ------------------------------------------------------------------
def info_clear(self):
@ -468,10 +471,8 @@ class ASCIItable():
"""Read whole data of all (given) labels as numpy array"""
from collections import Iterable
try:
self.data_rewind() # try to wind back to start of data
except:
pass # assume/hope we are at data start already...
try: self.data_rewind() # try to wind back to start of data
except: pass # assume/hope we are at data start already...
if labels is None or labels == []:
use = None # use all columns (and keep labels intact)
@ -530,13 +531,13 @@ class ASCIItable():
# ------------------------------------------------------------------
def data_append(self,
what):
if not isinstance(what, (str, unicode)):
if isinstance(what, (str, unicode)):
self.data += [what]
else:
try:
for item in what: self.data_append(item)
except:
except TypeError:
self.data += [str(what)]
else:
self.data += [what]
# ------------------------------------------------------------------
def data_set(self,
@ -581,7 +582,7 @@ class ASCIItable():
if len(items) > 2:
if items[1].lower() == 'of':
items = np.ones(datatype(items[0]))*datatype(items[2])
elif items[1].lower() == 'to':
elif items[1].lower() == 'to':
items = np.linspace(datatype(items[0]),datatype(items[2]),
abs(datatype(items[2])-datatype(items[0]))+1,dtype=int)
else: items = list(map(datatype,items))

View File

@ -65,7 +65,7 @@ if np.any(np.array(options.size) <= 0.0):
if filename == []: filename = [None]
table = damask.ASCIItable(outname = filename[0],
buffered = False)
buffered = False, labeled=False)
damask.util.report(scriptName,filename[0])

View File

@ -440,7 +440,9 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
math_Mandel33to6, &
math_Plain99to3333
use material, only: &
phasememberAt, &
phase_plasticity, &
phase_plasticityInstance, &
material_phase, &
material_homog, &
temperature, &
@ -490,7 +492,7 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
ho, & !< homogenization
tme !< thermal member position
integer(pInt) :: &
i, j
i, j, instance, of
ho = material_homog(ip,el)
tme = thermalMapping(ho)%p(ip,el)
@ -509,8 +511,9 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp),ipc,ip,el)
dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
of = phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMp, Mp,instance,of)
case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp),ipc,ip,el)
@ -817,6 +820,7 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac
debug_constitutive, &
debug_levelBasic
use math, only: &
math_mul33x33, &
math_Mandel6to33, &
math_Mandel33to6, &
math_mul33x33
@ -824,6 +828,8 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac
mesh_NcpElems, &
mesh_maxNips
use material, only: &
phasememberAt, &
phase_plasticityInstance, &
phase_plasticity, &
phase_source, &
phase_Nsources, &
@ -881,38 +887,41 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac
real(pReal), intent(in), dimension(6) :: &
S6 !< 2nd Piola Kirchhoff stress (vector notation)
real(pReal), dimension(3,3) :: &
Mstar
Mp
integer(pInt) :: &
ho, & !< homogenization
tme, & !< thermal member position
s !< counter in source loop
s, & !< counter in source loop
instance, of
ho = material_homog( ip,el)
tme = thermalMapping(ho)%p(ip,el)
Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6))
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6))
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_ISOTROPIC_ID) plasticityType
call plastic_isotropic_dotState (math_Mandel33to6(Mstar),ipc,ip,el)
call plastic_isotropic_dotState (math_Mandel33to6(Mp),ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
call plastic_phenopowerlaw_dotState(math_Mandel33to6(Mstar),ipc,ip,el)
of = phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
call plastic_phenopowerlaw_dotState(Mp,instance,of)
case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_dotState(math_Mandel33to6(Mstar),ipc,ip,el)
call plastic_kinehardening_dotState(math_Mandel33to6(Mp),ipc,ip,el)
case (PLASTICITY_DISLOTWIN_ID) plasticityType
call plastic_dislotwin_dotState (Mstar,temperature(ho)%p(tme), &
call plastic_dislotwin_dotState (math_Mandel33to6(Mp),temperature(ho)%p(tme), &
ipc,ip,el)
case (PLASTICITY_DISLOUCLA_ID) plasticityType
call plastic_disloucla_dotState (math_Mandel33to6(Mstar),temperature(ho)%p(tme), &
call plastic_disloucla_dotState (math_Mandel33to6(Mp),temperature(ho)%p(tme), &
ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_dotState (math_Mandel33to6(Mstar),FeArray,FpArray,temperature(ho)%p(tme), &
call plastic_nonlocal_dotState (math_Mandel33to6(Mp),FeArray,FpArray,temperature(ho)%p(tme), &
subdt,subfracArray,ip,el)
end select plasticityType
@ -1026,13 +1035,18 @@ end subroutine constitutive_collectDeltaState
!--------------------------------------------------------------------------------------------------
!> @brief returns array of constitutive results
!--------------------------------------------------------------------------------------------------
function constitutive_postResults(S6, FeArray, ipc, ip, el)
function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el)
use prec, only: &
pReal
use math, only: &
math_Mandel6to33, &
math_mul33x33
use mesh, only: &
mesh_NcpElems, &
mesh_maxNips
use material, only: &
phasememberAt, &
phase_plasticityInstance, &
plasticState, &
sourceState, &
phase_plasticity, &
@ -1083,19 +1097,25 @@ function constitutive_postResults(S6, FeArray, ipc, ip, el)
real(pReal), dimension(plasticState(material_phase(ipc,ip,el))%sizePostResults + &
sum(sourceState(material_phase(ipc,ip,el))%p(:)%sizePostResults)) :: &
constitutive_postResults
real(pReal), intent(in), dimension(3,3) :: &
Fi !< intermediate deformation gradient
real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
FeArray !< elastic deformation gradient
real(pReal), intent(in), dimension(6) :: &
S6 !< 2nd Piola Kirchhoff stress (vector notation)
real(pReal), dimension(3,3) :: &
Mp !< Mandel stress
integer(pInt) :: &
startPos, endPos
integer(pInt) :: &
ho, & !< homogenization
tme, & !< thermal member position
s !< counter in source loop
s, of, instance !< counter in source loop
constitutive_postResults = 0.0_pReal
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6))
ho = material_homog( ip,el)
tme = thermalMapping(ho)%p(ip,el)
@ -1104,10 +1124,13 @@ function constitutive_postResults(S6, FeArray, ipc, ip, el)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_ISOTROPIC_ID) plasticityType
constitutive_postResults(startPos:endPos) = plastic_isotropic_postResults(S6,ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
constitutive_postResults(startPos:endPos) = &
plastic_phenopowerlaw_postResults(S6,ipc,ip,el)
plastic_isotropic_postResults(S6,ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
of = phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
constitutive_postResults(startPos:endPos) = &
plastic_phenopowerlaw_postResults(Mp,instance,of)
case (PLASTICITY_KINEHARDENING_ID) plasticityType
constitutive_postResults(startPos:endPos) = &
plastic_kinehardening_postResults(S6,ipc,ip,el)

File diff suppressed because it is too large Load Diff

View File

@ -448,8 +448,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
subStepSizeHomog, &
stepIncreaseHomog, &
nMPstate
use math, only: &
math_transpose33
use FEsolving, only: &
FEsolving_execElem, &
FEsolving_execIP, &
@ -496,6 +494,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
crystallite_converged, &
crystallite_stressAndItsTangent, &
crystallite_orientations
#ifdef DEBUG
use debug, only: &
debug_level, &
debug_homogenization, &
@ -504,6 +503,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
debug_levelSelective, &
debug_e, &
debug_i
#endif
implicit none
real(pReal), intent(in) :: dt !< time increment
@ -517,18 +517,16 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
mySource, &
myNgrains
!--------------------------------------------------------------------------------------------------
! initialize to starting condition
#ifdef DEBUG
if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(/a,i5,1x,i2)') '<< HOMOG >> Material Point start at el ip ', debug_e, debug_i
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< HOMOG >> F0', &
math_transpose33(materialpoint_F0(1:3,1:3,debug_i,debug_e))
transpose(materialpoint_F0(1:3,1:3,debug_i,debug_e))
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< HOMOG >> F', &
math_transpose33(materialpoint_F(1:3,1:3,debug_i,debug_e))
!$OMP END CRITICAL (write2out)
transpose(materialpoint_F(1:3,1:3,debug_i,debug_e))
endif
#endif
!--------------------------------------------------------------------------------------------------
! initialize restoration points of ...
@ -608,10 +606,8 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
!---------------------------------------------------------------------------------------------------
! calculate new subStep and new subFrac
materialpoint_subFrac(i,e) = materialpoint_subFrac(i,e) + materialpoint_subStep(i,e)
!$OMP FLUSH(materialpoint_subFrac)
materialpoint_subStep(i,e) = min(1.0_pReal-materialpoint_subFrac(i,e), &
stepIncreaseHomog*materialpoint_subStep(i,e)) ! introduce flexibility for step increase/acceleration
!$OMP FLUSH(materialpoint_subStep)
steppingNeeded: if (materialpoint_subStep(i,e) > subStepMinHomog) then
@ -671,7 +667,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
hydrogenfluxState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = &
hydrogenfluxState(mappingHomogenization(2,i,e))%State( :,mappingHomogenization(1,i,e))! ...internal hydrogen transport state
materialpoint_subF0(1:3,1:3,i,e) = materialpoint_subF(1:3,1:3,i,e) ! ...def grad
!$OMP FLUSH(materialpoint_subF0)
endif steppingNeeded
else converged
@ -689,7 +684,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
!$OMP END CRITICAL (setTerminallyIll)
else ! cutback makes sense
materialpoint_subStep(i,e) = subStepSizeHomog * materialpoint_subStep(i,e) ! crystallite had severe trouble, so do a significant cutback
!$OMP FLUSH(materialpoint_subStep)
#ifdef DEBUG
if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt &
@ -752,8 +746,9 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
if (materialpoint_subStep(i,e) > subStepMinHomog) then
materialpoint_requested(i,e) = .true.
materialpoint_subF(1:3,1:3,i,e) = materialpoint_subF0(1:3,1:3,i,e) + &
materialpoint_subStep(i,e) * (materialpoint_F(1:3,1:3,i,e) - materialpoint_F0(1:3,1:3,i,e))
materialpoint_subF(1:3,1:3,i,e) = materialpoint_subF0(1:3,1:3,i,e) &
+ materialpoint_subStep(i,e) * (materialpoint_F(1:3,1:3,i,e) &
- materialpoint_F0(1:3,1:3,i,e))
materialpoint_subdt(i,e) = materialpoint_subStep(i,e) * dt
materialpoint_doneAndHappy(1:2,i,e) = [.false.,.true.]
endif
@ -810,13 +805,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
materialpoint_doneAndHappy(1:2,i,e) = homogenization_updateState(i,e)
materialpoint_converged(i,e) = all(materialpoint_doneAndHappy(1:2,i,e)) ! converged if done and happy
endif
!$OMP FLUSH(materialpoint_converged)
if (materialpoint_converged(i,e)) then
if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (distributionMPState)
!$OMP END CRITICAL (distributionMPState)
endif
endif
endif
enddo IpLooping3
enddo elementLooping3
@ -838,9 +826,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
enddo elementLooping4
!$OMP END PARALLEL DO
else
!$OMP CRITICAL (write2out)
write(6,'(/,a,/)') '<< HOMOG >> Material Point terminally ill'
!$OMP END CRITICAL (write2out)
endif
end subroutine materialpoint_stressAndItsTangent

View File

@ -942,10 +942,10 @@ pure function homogenization_RGC_postResults(ip,el,avgP,avgF)
do o = 1_pInt,homogenization_Noutput(mesh_element(3,el))
select case(homogenization_RGC_outputID(o,homID))
case (avgdefgrad_ID)
homogenization_RGC_postResults(c+1_pInt:c+9_pInt) = reshape(avgF,[9])
homogenization_RGC_postResults(c+1_pInt:c+9_pInt) = reshape(transpose(avgF),[9])
c = c + 9_pInt
case (avgfirstpiola_ID)
homogenization_RGC_postResults(c+1_pInt:c+9_pInt) = reshape(avgP,[9])
homogenization_RGC_postResults(c+1_pInt:c+9_pInt) = reshape(transpose(avgP),[9])
c = c + 9_pInt
case (ipcoords_ID)
homogenization_RGC_postResults(c+1_pInt:c+3_pInt) = mesh_ipCoordinates(1:3,ip,el) ! current ip coordinates

View File

@ -300,10 +300,10 @@ pure function homogenization_isostrain_postResults(ip,el,avgP,avgF)
homogenization_isostrain_postResults(c+1_pInt) = real(homogenization_isostrain_Ngrains(homID),pReal)
c = c + 1_pInt
case (avgdefgrad_ID)
homogenization_isostrain_postResults(c+1_pInt:c+9_pInt) = reshape(avgF,[9])
homogenization_isostrain_postResults(c+1_pInt:c+9_pInt) = reshape(transpose(avgF),[9])
c = c + 9_pInt
case (avgfirstpiola_ID)
homogenization_isostrain_postResults(c+1_pInt:c+9_pInt) = reshape(avgP,[9])
homogenization_isostrain_postResults(c+1_pInt:c+9_pInt) = reshape(transpose(avgP),[9])
c = c + 9_pInt
case (ipcoords_ID)
homogenization_isostrain_postResults(c+1_pInt:c+3_pInt) = mesh_ipCoordinates(1:3,ip,el) ! current ip coordinates

File diff suppressed because it is too large Load Diff