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

overhaul to allow for multiphase/homogenization setup.

This commit is contained in:
Philip Eisenlohr 2009-03-04 11:48:54 +00:00
parent 5a09badb71
commit 9351508ae1
13 changed files with 1688 additions and 133 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

51
trunk/FEsolving.f90
View File

@ -9,8 +9,9 @@
integer(pInt) cycleCounter integer(pInt) cycleCounter
integer(pInt) theInc,theCycle,theLovl integer(pInt) theInc,theCycle,theLovl
real(pReal) theTime real(pReal) theTime
logical :: lastIncConverged = .false.,outdatedByNewInc = .false., outdatedFFN1 = .false. logical :: lastIncConverged = .false.,outdatedByNewInc = .false.,outdatedFFN1 = .false.
logical :: symmetricSolver = .false. logical :: symmetricSolver = .false.
logical :: parallelExecution = .true.
CONTAINS CONTAINS
@ -18,36 +19,36 @@
!*********************************************************** !***********************************************************
! determine wether a symmetric solver is used ! determine wether a symmetric solver is used
!*********************************************************** !***********************************************************
subroutine FE_get_solverSymmetry(unit) subroutine FE_init()
use prec, only: pInt use prec, only: pInt
use IO use IO
implicit none implicit none
integer(pInt) unit integer(pInt), parameter :: fileunit = 222
integer(pInt), dimension (133) :: pos integer(pInt), dimension (1+2*2) :: pos
character*300 line character(len=1024) line
610 FORMAT(A300)
rewind(unit)
do
read (unit,610,END=630) line
pos = IO_stringPos(line,1)
if( IO_lc(IO_stringValue(line,pos,1)) == 'solver' ) then
read (unit,610,END=630) line ! Garbage line
pos = IO_stringPos(line,2) ! limit to 64 nodes max (plus ID, type)
if(IO_intValue(line,pos,2) /= 1_pInt) then
symmetricSolver = .true.
!$OMP CRITICAL (write2out)
write (6,*)
write (6,*) 'Symmetric solver detected. d-Matrix will be symmetrized!'
!$OMP END CRITICAL (write2out)
endif
endif
enddo
630 return if (IO_open_inputFile(fileunit)) then
rewind(fileunit)
do
read (fileunit,'(a1024)',END=100) line
pos = IO_stringPos(line,1)
if( IO_lc(IO_stringValue(line,pos,1)) == 'solver' ) then
read (fileunit,'(a1024)',END=100) line ! Garbage line
pos = IO_stringPos(line,2)
symmetricSolver = (IO_intValue(line,pos,2) /= 1_pInt)
exit
endif
enddo
else
call IO_error(100) ! cannot open input file
endif
100 close(fileunit)
return
end subroutine end subroutine

202
trunk/IO.f90
View File

@ -244,6 +244,132 @@
END FUNCTION END FUNCTION
!********************************************************************
! identifies lines without content
!********************************************************************
PURE FUNCTION IO_isBlank (line)
use prec, only: pInt
implicit none
character(len=*), intent(in) :: line
character(len=*), parameter :: blank = achar(32)//achar(9)//achar(10)//achar(13) ! whitespaces
character(len=*), parameter :: comment = achar(35) ! comment id '#'
integer(pInt) posNonBlank, posComment
logical IO_isBlank
posNonBlank = verify(line,blank)
posComment = scan(line,comment)
IO_isBlank = posNonBlank == 0 .or. posNonBlank == posComment
return
END FUNCTION
!********************************************************************
! get tagged content of line
!********************************************************************
PURE FUNCTION IO_getTag (line,openChar,closechar)
use prec, only: pInt
implicit none
character(len=*), intent(in) :: line,openChar,closeChar
character(len=*), parameter :: sep=achar(32)//achar(9)//achar(10)//achar(13) ! whitespaces
character(len=len_trim(line)) IO_getTag
integer(pInt) left,right
IO_getTag = ''
left = scan(line,openChar)
right = scan(line,closeChar)
if (left == verify(line,sep) .and. right > left) & ! openChar is first and closeChar occurs
IO_getTag = line(left+1:right-1)
return
END FUNCTION
!*********************************************************************
FUNCTION IO_countSections(file,part)
!*********************************************************************
use prec, only: pInt
implicit none
!* Definition of variables
integer(pInt), intent(in) :: file
character(len=*), intent(in) :: part
integer(pInt) IO_countSections
character(len=1024) line
IO_countSections = 0
line = ''
rewind(file)
do while (IO_getTag(line,'<','>') /= part) ! search for part
read(file,'(a1024)',END=100) line
enddo
do
read(file,'(a1024)',END=100) line
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') & ! found [section] identifier
IO_countSections = IO_countSections + 1
enddo
100 return
END FUNCTION
!*********************************************************************
! return array of myTag counts within <part> for at most N[sections]
!*********************************************************************
FUNCTION IO_countTagInPart(file,part,myTag,Nsections)
use prec, only: pInt
implicit none
!* Definition of variables
integer(pInt), intent(in) :: file, Nsections
character(len=*), intent(in) :: part, myTag
integer(pInt), dimension(Nsections) :: IO_countTagInPart, counter
integer(pInt), parameter :: maxNchunks = 1
integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) section
character(len=1024) line,tag
counter = 0_pInt
section = 0_pInt
line = ''
rewind(file)
do while (IO_getTag(line,'<','>') /= part) ! search for part
read(file,'(a1024)',END=100) line
enddo
do
read(file,'(a1024)',END=100) line
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') & ! found [section] identifier
section = section + 1
if (section > 0) then
positions = IO_stringPos(line,maxNchunks)
tag = IO_lc(IO_stringValue(line,positions,1)) ! extract key
if (tag == myTag) & ! match
counter(section) = counter(section) + 1
endif
enddo
100 IO_countTagInPart = counter
return
END FUNCTION
!******************************************************************** !********************************************************************
! locate at most N space-separated parts in line ! locate at most N space-separated parts in line
! return array containing number of parts found and ! return array containing number of parts found and
@ -559,48 +685,90 @@
! and terminate the Marc run with exit #9xxx ! and terminate the Marc run with exit #9xxx
! in ABAQUS either time step is reduced or execution terminated ! in ABAQUS either time step is reduced or execution terminated
!******************************************************************** !********************************************************************
SUBROUTINE IO_error(ID) SUBROUTINE IO_error(ID,e,i,g,ext_msg)
use prec, only: pInt use prec, only: pInt
use debug use debug
implicit none implicit none
integer(pInt) ID integer(pInt), intent(in) :: ID
integer(pInt), optional, intent(in) :: e,i,g
character(len=*), optional, intent(in) :: ext_msg
character(len=80) msg character(len=80) msg
select case (ID) select case (ID)
case (100) case (0)
msg='Unable to open input file.' msg='Unable to open input file.'
case (100)
msg='Error reading from configuration file.'
case (105)
msg='Error reading from ODF file.'
case (110) case (110)
msg='No materials specified via State Variable 2.' msg='No homogenization specified via State Variable 2.'
case (120) case (120)
msg='No textures specified via State Variable 3.' msg='No microstructure specified via State Variable 3.'
case (130)
msg='Homogenization index out of bounds.'
case (140)
msg='Microstructure index out of bounds.'
case (150)
msg='Phase index out of bounds.'
case (160)
msg='Texture index out of bounds.'
case (170)
msg='Sum of phase fractions differs from 1.'
case (200) case (200)
msg='Error reading from material+texture file' msg='Unknown constitution specified.'
case (201)
msg='Unknown lattice type specified.'
case (202)
msg='Number of slip systems too small.'
case (203)
msg='Negative initial slip resistance.'
case (204)
msg='Non-positive reference shear rate.'
case (205)
msg='Non-positive stress exponent.'
case (206)
msg='Non-positive initial hardening slope.'
case (207)
msg='Non-positive saturation stress.'
case (208)
msg='Non-positive w0.'
case (209)
msg='Negative latent hardening ratio.'
case (300) case (300)
msg='This material can only be used with elements with three direct stress components' msg='This material can only be used with elements with three direct stress components.'
case (400)
msg='Unknown alloy number specified'
case (500) case (500)
msg='Unknown lattice type specified' msg='Unknown lattice type specified.'
case (600) case (600)
msg='Convergence not reached' msg='Convergence not reached.'
case (610)
msg='Stress loop not converged.'
case (650) case (650)
msg='Polar decomposition failed' msg='Polar decomposition failed.'
case (700) case (700)
msg='Singular matrix in stress iteration' msg='Singular matrix in stress iteration.'
case (800) case (800)
msg='GIA requires 8 grains per IP (bonehead, you!)' msg='GIA requires 8 grains per IP (bonehead, you!)'
case default case default
msg='Unknown error number' msg='Unknown error number...'
end select end select
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) 'MPIE Material Routine Ver. 0.0 by the coding team' write(6,*) '+------------------------------+'
write(6,*)
write(6,*) msg write(6,*) msg
write(6,*) if (present(ext_msg)) write(6,*) ext_msg
write(6,*) '+------------------------------+'
if (present(e)) then
if (present(i) .and. present(g)) then
write(6,'(a10,x,i6,x,a2,x,i2,x,a5,x,i4)') 'at element',e,'IP',i,'grain',g
else
write(6,'(a2,x,i6)') 'at',e
endif
write(6,*)
endif
call debug_info() call debug_info()
call flush(6) call flush(6)

182
trunk/concom2007r1 Normal file
View File

@ -0,0 +1,182 @@
! reformated to free format
!***********************************************************************
!
! File: concom.cmn
!
! MSC.Marc include file
!
integer(pInt) &
iacous, iasmbl, iautth, ibear, icompl, iconj, icreep, ideva, idyn, idynt,&
ielas, ielcma, ielect, iform, ifour, iharm, ihcps, iheat, iheatt, ihresp,&
ijoule, ilem, ilnmom, iloren, inc, incext, incsub, ipass, iplres, ipois,&
ipoist, irpflo, ismall, ismalt, isoil, ispect, ispnow, istore, iswep, ithcrp,&
itherm, iupblg, iupdat, jacflg, jel, jparks, largst, lfond, loadup, loaduq,&
lodcor, lovl, lsub, magnet, ncycle, newtnt, newton, noshr, linear, ivscpl,&
icrpim, iradrt, ipshft, itshr, iangin, iupmdr, iconjf, jincfl, jpermg, jhour,&
isolvr, jritz, jtable, jshell, jdoubl, jform, jcentr, imini, kautth, iautof,&
ibukty, iassum, icnstd, icnstt, kmakmas, imethvp,iradrte,iradrtp, iupdate,iupdatp,&
ncycnt, marmen ,idynme, ihavca, ispf, kmini, imixed, largtt, kdoela, iautofg,&
ipshftp,idntrc, ipore, jtablm, jtablc, isnecma,itrnspo,imsdif, jtrnspo,mcnear,&
imech, imecht, ielcmat, ielectt,magnett, imsdift,noplas, jtabls, jactch, jtablth,&
kgmsto ,jpzo, ifricsh, iremkin,iremfor, ishearp,jspf, machining, jlshell,icompsol,&
iupblgfo,jcondir,nstcrp, nactive,ipassref, nstspnt,ibeart,icheckmpc, noline, icuring,&
ishrink,ioffsflg,isetoff, ioffsetm,iharmt, inc_incdat,iautspc,ibrake, icbush ,istream_input,&
iprsinp,ivlsinp,ifirst_time,ipin_m,jgnstr_glb,imarc_return,iqvcinp,nqvceid,istpnx,imicro1
common/marc_concom/&
iacous, iasmbl, iautth, ibear, icompl, iconj, icreep, ideva(50), idyn, idynt,&
ielas, ielcma, ielect, iform, ifour, iharm, ihcps, iheat, iheatt, ihresp,&
ijoule, ilem, ilnmom, iloren, inc, incext, incsub, ipass, iplres, ipois,&
ipoist, irpflo, ismall, ismalt, isoil, ispect, ispnow, istore, iswep, ithcrp,&
itherm, iupblg, iupdat, jacflg, jel, jparks, largst, lfond, loadup, loaduq,&
lodcor, lovl, lsub, magnet, ncycle, newtnt, newton, noshr, linear, ivscpl,&
icrpim, iradrt, ipshft, itshr, iangin, iupmdr, iconjf, jincfl, jpermg, jhour,&
isolvr, jritz, jtable, jshell, jdoubl, jform, jcentr, imini, kautth, iautof,&
ibukty, iassum, icnstd, icnstt, kmakmas, imethvp,iradrte,iradrtp, iupdate,iupdatp,&
ncycnt, marmen ,idynme, ihavca, ispf, kmini, imixed, largtt, kdoela, iautofg,&
ipshftp,idntrc, ipore, jtablm, jtablc, isnecma,itrnspo,imsdif, jtrnspo,mcnear,&
imech, imecht, ielcmat, ielectt,magnett, imsdift,noplas, jtabls, jactch, jtablth,&
kgmsto ,jpzo, ifricsh, iremkin,iremfor, ishearp,jspf, machining, jlshell,icompsol,&
iupblgfo,jcondir,nstcrp, nactive,ipassref, nstspnt,ibeart,icheckmpc, noline, icuring,&
ishrink,ioffsflg,isetoff, ioffsetm,iharmt, inc_incdat,iautspc,ibrake, icbush ,istream_input,&
iprsinp,ivlsinp,ifirst_time,ipin_m,jgnstr_glb,imarc_return,iqvcinp,nqvceid,istpnx,imicro1
!
! comments of variables:
!
! ideva(50) - debug print out flag
! 1 print element stiffness matrices, mass matrix
! 2 output matrices used in tying
! 3 force the solution of a nonpositive definite matrix
! 4 print info of connections to each node
! 5 info of gap convergence, internal heat generated, contact
! touching and separation
! 6 nodal value array during rezoning
! 7 tying info in CONRAD GAP option, fluid element numbers in
! CHANNEL option
! 8 output incremental displacements in local coord. system
! 9 latent heat output
! 10 stress-strain in local coord. system
! 11 additional info on interlaminar stress
! 12 output right hand side and solution vector
! 13 info of CPU resources used and memory available on NT
! 14 info of mesh adaption process, 2D outline information
! info of penetration checking for remeshing
! save .fem files after afmesh3d meshing
! 15 surface energy balance flag
! 16 print info regarding pyrolysis
! 17 print info of "streamline topology"
! 18 print mesh data changes after remeshing
! 19 print material flow stress data read in from *.mat file
! if unit flag is on, print out flow stress after conversion
! 20 print information on table input
! 21 print out information regarding kinematic boundary conditions
! 22 print out information regarding dist loads, point loads, film
! and foundations
! 23 print out information about automatic domain decomposition
! 24 print out iteration information in SuperForm status report file
! 25 print out information for ablation
! 26 print out information for films - Table input
! 27 print out the tying forces
! 28 print out for CASI solver, convection,
! 29 DDM single file debug printout
! 30 print out cavity debug info
! 31 print out welding related info
! 32 prints categorized DDM memory usage
! 33 print out the cutting info regarding machining feature
! 34 print out the list of quantities which can be defined via a table
! and for each quantity the supported independent variables
! 35 print out detailed coupling region info
! 36 print out solver debug info level 1 (Least Detailed)
! 37 print out solver debug info level 1 (Medium Detailed)
! 38 print out solver debug info level 1 (Very Detailed)
! 39 print detailed memory allocation info
! 40 print out marc-adams debug info
! 41 output rezone mapping post file for debugging
! 42 output post file after calling oprofos() for debugging
! 43 debug printout for vcct
! 44 debug printout for progressive failure
! 45 print out automatically generated midside node coordinates (arecrd)
! 46 print out message about routine and location, where the ibort is raised (ibort_inc)
! 47 print out summary message of element variables on a
! group-basis after all the automatic changes have been
! made (em_ellibp)
! 48 Automatically generate check results based on max and min vals.
! These vals are stored in the checkr file, which is inserted
! into the *dat file by the generate_check_results script from /marc/tools
! 49 Automatically generate check results based on the real calculated values
! at the sppecified check result locations.
! These vals are stored in the checkr file, which is inserted
! into the *dat file by the update_check_results script from /marc/tools
!
!
! jactch = 1 or 2 if elements are activated or deactivated
! = 3 if elements are adaptively remeshed or rezoned
! = 0 normally / reset to 0 when assembly is done
! ifricsh = 0 call to fricsh in otest not needed
! = 1 call to fricsh (nodal friction) in otest needed
! iremkin = 0 remove deactivated kinematic boundary conditions
! immediately - only in new input format (this is default)
! = 1 remove deactivated kinematic boundary conditions
! gradually - only in new input format
! iremfor = 0 remove force boundary conditions immediately -
! only in new input format (this is default)
! = 1 remove force boundary conditions gradually -
! only in new input format (this is default)
! ishearp set to 1 if shear panel elements are present in the model
!
! jspf = 0 not in spf loadcase
! > 0 in spf loadcase (jspf=1 during first increment)
! machining = 1 if the metal cutting feature is used, for memory allocation purpose
! = 0 (default) if no metal cutting feature required
!
! jlshell = 1 if there is a shell element in the mesh
! icompsol = 1 if there is a composite solid element in the mesh
! iupblgfo = 1 if follower force for point loads
! jcondir = 1 if contact priority option is used
! nstcrp = 0 (default) steady state creep flag (undocumented feature.
! if not 0, turns off special ncycle = 0 code in radial.f)
! nactive = number of active passes, if =1 then it's not a coupled analysis
! ipassref = reference ipass, if not in a multiphysics pass ipass=ipassref
! icheckmpc = value of mpc-check parameter option
! noline = set to 1 in osolty if no line seacrh should be done in ogetst
! icuring = set to 1 if the curing is included for the heat transfer analysis.
! ishrink = set to 1 if shrinkage strain is included for mechancial analysis.
! ioffsflg = 1 for small displacement beam/shell offsets
! = 2 for large displacement beam/shell offsets
! isetoff = 0 - do not apply beam/shell offsets
! = 1 - apply beam/shell offsets
! ioffsetm = min. value of offset flag
! inc_incdat = flag to record increment number of a new loadcase in incdat.f
! iautspc = flag for AutoSPC option
! ibrake = brake squeal in this increment
! icbush = set to 1 if cbush elements present in model
! istream_input = set to 1 for streaming input calling Marc as library
! iprsinp = set to 1 if pressure input, introduced so other variables
! such as h could be a function of pressure
! ivlsinp = set to 1 if velocity input, introduced so other variables
! such as h could be a function of velocity
! ipin_m = # of beam element with PIN flag
! jgnstr_glb = global control over pre or fast integrated composite shells
! imarc_return = Marc return flag for streaming input control
! iqvcimp = if non-zero, then the number of QVECT boundary conditions
! nqvceid = number of QVECT boundary conditions, where emisivity/absorbtion id entered
! istpnx = 1 if to stop at end of increment
! imicro1 = 1 if micro1 interface is used
! iaxisymm = set to 1 if axisymmetric analysis
! jbreakglue = set to 1 if breaking glued option is used
! iglstif = 1 if ddm and global stiffness matrix formed (sgi solver 6 or solver9)
! jfastasm = 1 do fast assembly using SuperForm code
! iwear = set to 1 if wear model, set to 2 if wear model and coordinates updated
! iwearcf = set to 1 to store nodal coefficient of friction for wear calculation
! imixmeth = set=1 then use nonlinear mixture material - allocate memory
! ielcmadyn = flag for magnetodynamics
! 0 - electromagnetics using newmark beta
! 1 - transient magnetics using backward euler
! idinout = flag to control if inside out elements should be deactivated
! igena_meth = 0 - generalized alpha parameters depend on whether or not contact
! is flagged (dynamic,7)
! 10 - generalized alpha parameters are optimized for a contact
! analysis (dynamic,8)
! 11 - generalized alpha parameters are optimized for an analysis
! without contact (dynamic,8)
!
!***********************************************************************

447
trunk/constitutive_phenomenological.f90 Normal file
View File

@ -0,0 +1,447 @@
!*****************************************************
!* Module: CONSTITUTIVE_PHENOMENOLOGICAL *
!*****************************************************
!* contains: *
!* - constitutive equations *
!* - parameters definition *
!*****************************************************
MODULE constitutive_phenomenological
!*** Include other modules ***
use prec, only: pReal,pInt
implicit none
character (len=*), parameter :: constitutive_phenomenological_label = 'phenomenological'
integer(pInt), dimension(:), allocatable :: constitutive_phenomenological_sizeDotState, &
constitutive_phenomenological_sizeState, &
constitutive_phenomenological_sizePostResults
character(len=64), dimension(:,:), allocatable :: constitutive_phenomenological_output
integer(pInt), dimension(:), allocatable :: constitutive_phenomenological_structure
integer(pInt), dimension(:), allocatable :: constitutive_phenomenological_Nslip
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_C11
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_C12
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_C13
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_C33
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_C44
real(pReal), dimension(:,:,:), allocatable :: constitutive_phenomenological_Cslip_66
!* Visco-plastic constitutive_phenomenological parameters
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_s0_slip
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_gdot0_slip
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_n_slip
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_h0
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_s_sat
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_w0
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_latent
real(pReal), dimension(:,:,:), allocatable :: constitutive_phenomenological_HardeningMatrix
CONTAINS
!****************************************
!* - constitutive_init
!* - constitutive_homogenizedC
!* - constitutive_microstructure
!* - constitutive_LpAndItsTangent
!* - consistutive_dotState
!* - consistutive_postResults
!****************************************
subroutine constitutive_phenomenological_init(file)
!**************************************
!* Module initialization *
!**************************************
use prec, only: pInt, pReal
use math, only: math_Mandel3333to66, math_Voigt66to3333
use IO
use material
integer(pInt), intent(in) :: file
integer(pInt), parameter :: maxNchunks = 7
integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) section, maxNinstance, i,j,k,l, output
character(len=64) tag
character(len=1024) line
maxNinstance = count(phase_constitution == constitutive_phenomenological_label)
if (maxNinstance == 0) return
allocate(constitutive_phenomenological_sizeDotState(maxNinstance)) ; constitutive_phenomenological_sizeDotState = 0_pInt
allocate(constitutive_phenomenological_sizeState(maxNinstance)) ; constitutive_phenomenological_sizeState = 0_pInt
allocate(constitutive_phenomenological_sizePostResults(maxNinstance)); constitutive_phenomenological_sizePostResults = 0_pInt
allocate(constitutive_phenomenological_output(maxval(phase_Noutput), &
maxNinstance)) ; constitutive_phenomenological_output = ''
allocate(constitutive_phenomenological_structure(maxNinstance)) ; constitutive_phenomenological_structure = 0_pInt
allocate(constitutive_phenomenological_Nslip(maxNinstance)) ; constitutive_phenomenological_Nslip = 0_pInt
allocate(constitutive_phenomenological_C11(maxNinstance)) ; constitutive_phenomenological_C11 = 0.0_pReal
allocate(constitutive_phenomenological_C12(maxNinstance)) ; constitutive_phenomenological_C12 = 0.0_pReal
allocate(constitutive_phenomenological_C13(maxNinstance)) ; constitutive_phenomenological_C13 = 0.0_pReal
allocate(constitutive_phenomenological_C33(maxNinstance)) ; constitutive_phenomenological_C33 = 0.0_pReal
allocate(constitutive_phenomenological_C44(maxNinstance)) ; constitutive_phenomenological_C44 = 0.0_pReal
allocate(constitutive_phenomenological_Cslip_66(6,6,maxNinstance)) ; constitutive_phenomenological_Cslip_66 = 0.0_pReal
allocate(constitutive_phenomenological_s0_slip(maxNinstance)) ; constitutive_phenomenological_s0_slip = 0.0_pReal
allocate(constitutive_phenomenological_gdot0_slip(maxNinstance)) ; constitutive_phenomenological_gdot0_slip = 0.0_pReal
allocate(constitutive_phenomenological_n_slip(maxNinstance)) ; constitutive_phenomenological_n_slip = 0.0_pReal
allocate(constitutive_phenomenological_h0(maxNinstance)) ; constitutive_phenomenological_h0 = 0.0_pReal
allocate(constitutive_phenomenological_s_sat(maxNinstance)) ; constitutive_phenomenological_s_sat = 0.0_pReal
allocate(constitutive_phenomenological_w0(maxNinstance)) ; constitutive_phenomenological_w0 = 0.0_pReal
allocate(constitutive_phenomenological_latent(maxNinstance)) ; constitutive_phenomenological_latent = 1.0_pReal
rewind(file)
line = ''
section = 0
do while (IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to <phase>
read(file,'(a1024)',END=100) line
enddo
do ! read thru sections of phase part
read(file,'(a1024)',END=100) line
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1
output = 0 ! reset output counter
endif
if (section > 0 .and. phase_constitution(section) == constitutive_phenomenological_label) then ! one of my sections
i = phase_constitutionInstance(section) ! which instance of my constitution is present phase
positions = IO_stringPos(line,maxNchunks)
tag = IO_lc(IO_stringValue(line,positions,1)) ! extract key
select case(tag)
case ('(output)')
output = output + 1
constitutive_phenomenological_output(output,i) = IO_lc(IO_stringValue(line,positions,2))
case ('lattice_structure')
constitutive_phenomenological_structure(i) = IO_intValue(line,positions,2)
case ('nslip')
constitutive_phenomenological_Nslip(i) = IO_intValue(line,positions,2)
case ('c11')
constitutive_phenomenological_C11(i) = IO_floatValue(line,positions,2)
case ('c12')
constitutive_phenomenological_C12(i) = IO_floatValue(line,positions,2)
case ('c13')
constitutive_phenomenological_C13(i) = IO_floatValue(line,positions,2)
case ('c33')
constitutive_phenomenological_C33(i) = IO_floatValue(line,positions,2)
case ('c44')
constitutive_phenomenological_C44(i) = IO_floatValue(line,positions,2)
case ('s0_slip')
constitutive_phenomenological_s0_slip(i) = IO_floatValue(line,positions,2)
case ('gdot0_slip')
constitutive_phenomenological_gdot0_slip(i) = IO_floatValue(line,positions,2)
case ('n_slip')
constitutive_phenomenological_n_slip(i) = IO_floatValue(line,positions,2)
case ('h0')
constitutive_phenomenological_h0(i) = IO_floatValue(line,positions,2)
case ('s_sat')
constitutive_phenomenological_s_sat(i) = IO_floatValue(line,positions,2)
case ('w0')
constitutive_phenomenological_w0(i) = IO_floatValue(line,positions,2)
case ('latent_ratio')
constitutive_phenomenological_latent(i) = IO_floatValue(line,positions,2)
end select
endif
enddo
100 do i = 1,maxNinstance ! sanity checks
if (constitutive_phenomenological_structure(i) < 1 .or. &
constitutive_phenomenological_structure(i) > 3) call IO_error(201)
if (constitutive_phenomenological_Nslip(i) < 1) call IO_error(202)
if (constitutive_phenomenological_s0_slip(i) < 0.0_pReal) call IO_error(203)
if (constitutive_phenomenological_gdot0_slip(i) <= 0.0_pReal) call IO_error(204)
if (constitutive_phenomenological_n_slip(i) <= 0.0_pReal) call IO_error(205)
if (constitutive_phenomenological_h0(i) <= 0.0_pReal) call IO_error(206)
if (constitutive_phenomenological_s_sat(i) <= 0.0_pReal) call IO_error(207)
if (constitutive_phenomenological_w0(i) <= 0.0_pReal) call IO_error(208)
if (constitutive_phenomenological_latent(i) < 0.0_pReal) call IO_error(209)
enddo
allocate(constitutive_phenomenological_hardeningMatrix(maxval(constitutive_phenomenological_Nslip),&
maxval(constitutive_phenomenological_Nslip),&
maxNinstance))
do i = 1,maxNinstance
constitutive_phenomenological_sizeDotState(i) = constitutive_phenomenological_Nslip(i)
constitutive_phenomenological_sizeState(i) = constitutive_phenomenological_Nslip(i)
do j = 1,maxval(phase_Noutput)
select case(constitutive_phenomenological_output(j,i))
case('slipresistance')
constitutive_phenomenological_sizePostResults(i) = &
constitutive_phenomenological_sizePostResults(i) + constitutive_phenomenological_Nslip(i)
case('rateofshear')
constitutive_phenomenological_sizePostResults(i) = &
constitutive_phenomenological_sizePostResults(i) + constitutive_phenomenological_Nslip(i)
end select
enddo
select case (constitutive_phenomenological_structure(i))
case(1:2) ! cubic(s)
forall(k=1:3)
forall(j=1:3) &
constitutive_phenomenological_Cslip_66(k,j,i) = constitutive_phenomenological_C12(i)
constitutive_phenomenological_Cslip_66(k,k,i) = constitutive_phenomenological_C11(i)
constitutive_phenomenological_Cslip_66(k+3,k+3,i) = constitutive_phenomenological_C44(i)
end forall
case(3) ! hcp
constitutive_phenomenological_Cslip_66(1,1,i) = constitutive_phenomenological_C11(i)
constitutive_phenomenological_Cslip_66(2,2,i) = constitutive_phenomenological_C11(i)
constitutive_phenomenological_Cslip_66(3,3,i) = constitutive_phenomenological_C33(i)
constitutive_phenomenological_Cslip_66(1,2,i) = constitutive_phenomenological_C12(i)
constitutive_phenomenological_Cslip_66(2,1,i) = constitutive_phenomenological_C12(i)
constitutive_phenomenological_Cslip_66(1,3,i) = constitutive_phenomenological_C13(i)
constitutive_phenomenological_Cslip_66(3,1,i) = constitutive_phenomenological_C13(i)
constitutive_phenomenological_Cslip_66(2,3,i) = constitutive_phenomenological_C13(i)
constitutive_phenomenological_Cslip_66(3,2,i) = constitutive_phenomenological_C13(i)
constitutive_phenomenological_Cslip_66(4,4,i) = constitutive_phenomenological_C44(i)
constitutive_phenomenological_Cslip_66(5,5,i) = constitutive_phenomenological_C44(i)
constitutive_phenomenological_Cslip_66(6,6,i) = 0.5_pReal*(constitutive_phenomenological_C11(i)- &
constitutive_phenomenological_C12(i))
end select
constitutive_phenomenological_Cslip_66(:,:,i) = &
math_Mandel3333to66(math_Voigt66to3333(constitutive_phenomenological_Cslip_66(:,:,i)))
constitutive_phenomenological_hardeningMatrix(:,:,i) = constitutive_phenomenological_latent(i)
forall (j = 1:constitutive_phenomenological_Nslip(i)) &
constitutive_phenomenological_hardeningMatrix(j,j,i) = 1.0_pReal
enddo
return
end subroutine
function constitutive_phenomenological_stateInit(ipc,ip,el)
!*********************************************************************
!* initial microstructural state *
!*********************************************************************
use prec, only: pReal,pInt
use material, only: material_phase, phase_constitutionInstance
implicit none
!* Definition of variables
integer(pInt), intent(in) :: ipc,ip,el
integer(pInt) matID
real(pReal), dimension(constitutive_phenomenological_Nslip(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: &
constitutive_phenomenological_stateInit
matID = phase_constitutionInstance(material_phase(ipc,ip,el))
constitutive_phenomenological_stateInit = constitutive_phenomenological_s0_slip(matID)
return
end function
function constitutive_phenomenological_homogenizedC(state,ipc,ip,el)
!*********************************************************************
!* homogenized elacticity matrix *
!* INPUT: *
!* - state : state variables *
!* - ipc : component-ID of current integration point *
!* - ip : current integration point *
!* - el : current element *
!*********************************************************************
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
implicit none
!* Definition of variables
integer(pInt), intent(in) :: ipc,ip,el
integer(pInt) matID
real(pReal), dimension(6,6) :: constitutive_phenomenological_homogenizedC
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
matID = phase_constitutionInstance(material_phase(ipc,ip,el))
constitutive_phenomenological_homogenizedC = constitutive_phenomenological_Cslip_66(:,:,matID)
return
end function
subroutine constitutive_phenomenological_microstructure(Temperature,state,ipc,ip,el)
!*********************************************************************
!* calculate derived quantities from state (not used here) *
!* INPUT: *
!* - Tp : temperature *
!* - ipc : component-ID of current integration point *
!* - ip : current integration point *
!* - el : current element *
!*********************************************************************
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
implicit none
!* Definition of variables
integer(pInt) ipc,ip,el, matID
real(pReal) Temperature
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
matID = phase_constitutionInstance(material_phase(ipc,ip,el))
end subroutine
subroutine constitutive_phenomenological_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,state,ipc,ip,el)
!*********************************************************************
!* plastic velocity gradient and its tangent *
!* INPUT: *
!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) *
!* - ipc : component-ID at current integration point *
!* - ip : current integration point *
!* - el : current element *
!* OUTPUT: *
!* - Lp : plastic velocity gradient *
!* - dLp_dTstar : derivative of Lp (4th-rank tensor) *
!*********************************************************************
use prec, only: pReal,pInt,p_vec
use math, only: math_Plain3333to99
use lattice, only: lattice_Sslip,lattice_Sslip_v
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
implicit none
!* Definition of variables
integer(pInt) ipc,ip,el
integer(pInt) matID,i,k,l,m,n
real(pReal) Temperature
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
real(pReal), dimension(6) :: Tstar_v
real(pReal), dimension(3,3) :: Lp
real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333
real(pReal), dimension(9,9) :: dLp_dTstar
real(pReal), dimension(constitutive_phenomenological_Nslip(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: &
gdot_slip,dgdot_dtauslip,tau_slip
matID = phase_constitutionInstance(material_phase(ipc,ip,el))
!* Calculation of Lp
Lp = 0.0_pReal
do i = 1,constitutive_phenomenological_Nslip(matID)
tau_slip(i) = dot_product(Tstar_v,lattice_Sslip_v(:,i,constitutive_phenomenological_structure(matID)))
gdot_slip(i) = constitutive_phenomenological_gdot0_slip(matID)*(abs(tau_slip(i))/state(ipc,ip,el)%p(i))**&
constitutive_phenomenological_n_slip(matID)*sign(1.0_pReal,tau_slip(i))
Lp = Lp + gdot_slip(i)*lattice_Sslip(:,:,i,constitutive_phenomenological_structure(matID))
enddo
!* Calculation of the tangent of Lp
dLp_dTstar3333 = 0.0_pReal
dLp_dTstar = 0.0_pReal
do i = 1,constitutive_phenomenological_Nslip(matID)
dgdot_dtauslip(i) = constitutive_phenomenological_gdot0_slip(matID)*(abs(tau_slip(i))/state(ipc,ip,el)%p(i))**&
(constitutive_phenomenological_n_slip(matID)-1.0_pReal)*&
constitutive_phenomenological_n_slip(matID)/state(ipc,ip,el)%p(i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
dgdot_dtauslip(i)*lattice_Sslip(k,l,i,constitutive_phenomenological_structure(matID))* &
lattice_Sslip(m,n,i,constitutive_phenomenological_structure(matID))
enddo
dLp_dTstar = math_Plain3333to99(dLp_dTstar3333)
return
end subroutine
function constitutive_phenomenological_dotState(Tstar_v,Temperature,state,ipc,ip,el)
!*********************************************************************
!* rate of change of microstructure *
!* INPUT: *
!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) *
!* - ipc : component-ID at current integration point *
!* - ip : current integration point *
!* - el : current element *
!* OUTPUT: *
!* - constitutive_dotState : evolution of state variable *
!*********************************************************************
use prec, only: pReal,pInt,p_vec
use lattice, only: lattice_Sslip_v
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
implicit none
!* Definition of variables
integer(pInt) ipc,ip,el
integer(pInt) matID,i,n
real(pReal) Temperature,tau_slip,gdot_slip
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
real(pReal), dimension(6) :: Tstar_v
real(pReal), dimension(constitutive_phenomenological_Nslip(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: &
constitutive_phenomenological_dotState,self_hardening
matID = phase_constitutionInstance(material_phase(ipc,ip,el))
n = constitutive_phenomenological_Nslip(matID)
!* Self-Hardening of each system
do i = 1,n
tau_slip = dot_product(Tstar_v,lattice_Sslip_v(:,i,constitutive_phenomenological_structure(matID)))
gdot_slip = constitutive_phenomenological_gdot0_slip(matID)*(abs(tau_slip)/state(ipc,ip,el)%p(i))**&
constitutive_phenomenological_n_slip(matID)*sign(1.0_pReal,tau_slip)
self_hardening(i) = constitutive_phenomenological_h0(matID)*(1.0_pReal-state(ipc,ip,el)%p(i)/&
constitutive_phenomenological_s_sat(matID))**constitutive_phenomenological_w0(matID)*abs(gdot_slip)
enddo
!$OMP CRITICAL (evilmatmul)
constitutive_phenomenological_dotState = matmul(constitutive_phenomenological_hardeningMatrix(1:n,1:n,matID),self_hardening)
!$OMP END CRITICAL (evilmatmul)
return
end function
pure function constitutive_phenomenological_postResults(Tstar_v,Temperature,dt,state,ipc,ip,el)
!*********************************************************************
!* return array of constitutive results *
!* INPUT: *
!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) *
!* - dt : current time increment *
!* - ipc : component-ID at current integration point *
!* - ip : current integration point *
!* - el : current element *
!*********************************************************************
use prec, only: pReal,pInt,p_vec
use lattice, only: lattice_Sslip_v
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance,phase_Noutput
implicit none
!* Definition of variables
integer(pInt), intent(in) :: ipc,ip,el
real(pReal), intent(in) :: dt,Temperature
real(pReal), dimension(6), intent(in) :: Tstar_v
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state
integer(pInt) matID,o,i,c,n
real(pReal) tau_slip, active_rate
real(pReal), dimension(constitutive_phenomenological_sizePostResults(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: &
constitutive_phenomenological_postResults
matID = phase_constitutionInstance(material_phase(ipc,ip,el))
n = constitutive_phenomenological_Nslip(matID)
c = 0_pInt
constitutive_phenomenological_postResults = 0.0_pReal
do o = 1,phase_Noutput(material_phase(ipc,ip,el))
select case(constitutive_phenomenological_output(o,matID))
case ('slipresistance')
constitutive_phenomenological_postResults(c+1:c+n) = state(ipc,ip,el)%p(1:n)
c = c + n
case ('rateofshear')
do i = 1,n
tau_slip = dot_product(Tstar_v,lattice_Sslip_v(:,i,constitutive_phenomenological_structure(matID)))
constitutive_phenomenological_postResults(c+i) = sign(1.0_pReal,tau_slip)*constitutive_phenomenological_gdot0_slip(matID)*&
(abs(tau_slip)/state(ipc,ip,el)%p(i))**&
constitutive_phenomenological_n_slip(matID)
enddo
c = c + n
end select
enddo
return
end function
END MODULE

24
trunk/creeps2007r1 Normal file
View File

@ -0,0 +1,24 @@
! reformated to free format
!***********************************************************************
!
! File: creeps.cmn
!
! MSC.Marc include file
!
real(pReal) cptim,timinc,timinc_p,timinc_s,timincm,timinc_a,timinc_b,creept
integer(pInt) icptim,icfte,icfst,icfeq,icftm,icetem,mcreep,jcreep,icpa,icftmp,icfstr,&
icfqcp,icfcpm,icrppr,icrcha,icpb,iicpmt,iicpa
real(pReal) time_beg_lcase,time_beg_inc,fractol,time_beg_pst
!
common/marc_creeps/cptim,timinc,timinc_p,timinc_s,timincm,timinc_a,timinc_b,creept(33),icptim,icfte,icfst,&
icfeq,icftm,icetem,mcreep,jcreep,icpa,icftmp,icfstr,icfqcp,icfcpm,icrppr,icrcha,icpb,iicpmt,iicpa
common/marc_creeps2/time_beg_lcase,time_beg_inc,fractol,time_beg_pst
!
! time_beg_lcase time at the beginning of the current load case
! time_beg_inc time at the beginning of the current increment
! fractol fraction of loadcase or increment time when we
! consider it to be finished
! time_beg_pst time corresponding to first increment to be
! read in from thermal post file for auto step
!
!***********************************************************************

7
trunk/debug.f90
View File

@ -6,10 +6,11 @@
implicit none implicit none
integer(pInt), dimension(nCutback+1) :: debug_cutbackDistribution integer(pInt), dimension(nCutback+1) :: debug_cutbackDistribution = 0_pInt
integer(pInt), dimension(nInner) :: debug_InnerLoopDistribution integer(pInt), dimension(nInner) :: debug_InnerLoopDistribution = 0_pInt
integer(pInt), dimension(nOuter) :: debug_OuterLoopDistribution integer(pInt), dimension(nOuter) :: debug_OuterLoopDistribution = 0_pInt
logical :: debugger = .false. logical :: debugger = .false.
logical :: distribution_init = .false.
CONTAINS CONTAINS

103
trunk/material.config Normal file
View File

@ -0,0 +1,103 @@
#####################
<homogenization>
#####################
[SX]
type Taylor
Ngrains 1
[RGC]
type RGC
Ngrains 8
[Taylor4]
type Taylor
Ngrains 4
#####################
<microstructure>
#####################
[Aluminum_CubeSX]
(constitutent) material 1 texture 2 fraction 1.0
[Copper_rCubeSX]
(constitutent) material 2 texture 3 fraction 1.0
[DPsteel]
(constitutent) material 3 texture 1 fraction 0.8
(constitutent) material 4 texture 1 fraction 0.2
#####################
<phase>
#####################
[Aluminum]
[Copper]
[Ferrite]
[Martensite]
[TWIP steel FeMnC]
constitution phenomenological
lattice_structure 1
Nslip 12
(output) slipResistance
(output) rateOfShear
C11 183.9e9 # elastic constants in Pa
C12 101.9e9
C44 115.4e9
### phenomenological constitutive parameters ###
s0_slip 85.0e6 # initial slip resistance
gdot0_slip 0.001 # reference shear rate
n_slip 100.0 # stress exponent
h0 355.0e6 # initial hardening slope
s_sat 265.0e6 # saturation stress
w0 1.0 # exponent
latent_ratio 1.4 # latent/self hardening ratio
### dislocation density-based constitutive parameters ###
burgers 2.56e-10 # Burgers vector [m]
Qedge 5.5e-19 # Activation energy for dislocation glide [J/K] (0.5*G*b^3)
Qsd 4.7e-19 # Activation energy for self diffusion [J/K] (gamma-iron)
diff0 4.0e-5 # prefactor vacancy diffusion coeffficent (gamma-iron)
grain_size 2.0e-5 # Average grain size [m]
interaction_coefficients 1.0 2.2 3.0 1.6 3.8 4.5 # Dislocation interaction coefficients
rho0 6.0e12 # Initial dislocation density [m/m^3]
c1 0.1 # Passing stress adjustment
c2 2.0 # Jump width adjustment
c3 1.0 # Activation volume adjustment
c4 50.0 # Average slip distance adjustment for lock formation
c5 1.0 # Average slip distance adjustment when grain boundaries
c7 8.0 # Athermal recovery adjustment
c8 1.0e10 # Thermal recovery adjustment (plays no role for me)
stack_size 5.0e-8 # Average twin thickness (stacks) [m]
f_sat 1.0 # Total twin volume fraction saturation
c6 # Average slip distance adjustment when twin boundaries [???]
site_scaling 1.0e-6 # Scaling potential nucleation sites
q1 1.0 # Scaling the P-K force on the twinning dislocation
q2 1.0 # Scaling the resolved shear stress
#####################
<texture>
#####################
[random]
[CubeSX]
(gauss) phi1 0.0 Phi 0.0 phi2 0.0 scatter 0.0 fraction 1.0
[rCubeSX]
(gauss) phi1 45.0 Phi 0.0 phi2 0.0 scatter 0.0 fraction 1.0

610
trunk/material.f90 Normal file
View File

@ -0,0 +1,610 @@
!************************************
!* Module: MATERIAL *
!************************************
!* contains: *
!* - parsing of material.config *
!************************************
MODULE material
!*** Include other modules ***
use prec, only: pReal,pInt
implicit none
character(len=64), parameter :: material_configFile = 'material.config'
!*************************************
!* Definition of material properties *
!*************************************
!* Number of materials
integer(pInt) material_Nhomogenization, & ! number of homogenizations
material_Nmicrostructure, & ! number of microstructures
material_Nphase, & ! number of phases
material_Ntexture, & ! number of textures
microstructure_maxNconstituents, & ! max number of constituents in any phase
homogenization_maxNgrains, & ! max number of grains in any homogenization
texture_maxNgauss, & ! max number of Gauss components in any texture
texture_maxNfiber ! max number of Fiber components in any texture
character(len=64), dimension(:), allocatable :: homogenization_name, & ! name of each homogenization
homogenization_type, & ! type of each homogenization
microstructure_name, & ! name of each microstructure
phase_name, & ! name of each phase
phase_constitution, & ! constitution of each phase
texture_name ! name of each texture
character(len=256),dimension(:), allocatable :: texture_ODFfile ! name of each ODF file
integer(pInt), dimension(:), allocatable :: homogenization_Ngrains, & ! number of grains in each homogenization
microstructure_Nconstituents, & ! number of constituents in each microstructure
phase_constitutionInstance, & ! instance of particular constitution of each phase
phase_Noutput, & ! number of '(output)' items per phase
texture_symmetry, & ! number of symmetric orientations per texture
texture_Ngauss, & ! number of Gauss components per texture
texture_Nfiber ! number of Fiber components per texture
integer(pInt), dimension(:,:), allocatable :: microstructure_phase, & ! phase IDs of each microstructure
microstructure_texture ! texture IDs of each microstructure
real(pReal), dimension(:,:), allocatable :: microstructure_fraction ! vol fraction of each constituent in microstructure
integer(pInt), dimension(:,:,:), allocatable :: material_volFrac, & ! vol fraction of grain within phase (?)
material_phase ! phase of each grain,IP,element
real(pReal), dimension(:,:,:,:), allocatable :: material_EulerAngles ! initial orientation of each grain,IP,element
real(pReal), dimension(:,:,:), allocatable :: texture_Gauss, & ! data of each Gauss component
texture_Fiber ! data of each Fiber component
CONTAINS
!*********************************************************************
subroutine material_init()
!*********************************************************************
!* Module initialization *
!**************************************
use prec, only: pReal,pInt
use IO, only: IO_error, IO_open_file
implicit none
!* Definition of variables
integer(pInt), parameter :: fileunit = 200
integer(pInt) i
if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file
write(6,*) 'parsing homogenization...'
call material_parseHomogenization(fileunit)
write(6,*) 'parsing microstrcuture...'
call material_parseMicrostructure(fileunit)
write(6,*) 'parsing texture...'
call material_parseTexture(fileunit)
write(6,*) 'parsing phase...'
call material_parsePhase(fileunit)
close(fileunit)
if (minval(microstructure_phase) < 1 .or. maxval(microstructure_phase) > material_Nphase) call IO_error(150)
if (minval(microstructure_texture) < 1 .or. maxval(microstructure_texture) > material_Ntexture) call IO_error(160)
do i = 1,material_Nmicrostructure
if (sum(microstructure_fraction(:,i)) /= 1.0_pReal) call IO_error(170,i)
enddo
write (6,*)
write (6,*) 'MATERIAL configuration'
write (6,*)
write (6,*) 'Homogenization'
do i = 1,material_Nhomogenization
write (6,'(a32,x,a16,x,i4)') homogenization_name(i),homogenization_type(i),homogenization_Ngrains(i)
enddo
write (6,*)
write (6,*) 'Microstructure'
do i = 1,material_Nmicrostructure
write (6,'(a32,x,i4)') microstructure_name(i),microstructure_Nconstituents(i)
enddo
write(6,*) 'populating grains...'
call material_populateGrains()
write(6,*) 'populating grains finished...'
end subroutine
!*********************************************************************
subroutine material_parseHomogenization(file)
!*********************************************************************
use prec, only: pInt
use IO
implicit none
character(len=32), parameter :: myPart = 'homogenization'
integer(pInt), intent(in) :: file
integer(pInt), parameter :: maxNchunks = 2
integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) Nsections, section
character(len=64) tag
character(len=1024) line
Nsections= IO_countSections(file,myPart)
material_Nhomogenization = Nsections
allocate(homogenization_name(Nsections)); homogenization_name = ''
allocate(homogenization_type(Nsections)); homogenization_type = ''
allocate(homogenization_Ngrains(Nsections)); homogenization_Ngrains = 0_pInt
rewind(file)
line = ''
section = 0
do while (IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart
read(file,'(a1024)',END=100) line
enddo
do
read(file,'(a1024)',END=100) line
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1
homogenization_name(section) = IO_getTag(line,'[',']')
endif
if (section > 0) then
positions = IO_stringPos(line,maxNchunks)
tag = IO_lc(IO_stringValue(line,positions,1)) ! extract key
select case(tag)
case ('type')
homogenization_type(section) = IO_stringValue(line,positions,2)
case ('ngrains')
homogenization_Ngrains(section) = IO_intValue(line,positions,2)
end select
endif
enddo
100 homogenization_maxNgrains = maxval(homogenization_Ngrains)
return
end subroutine
!*********************************************************************
subroutine material_parseMicrostructure(file)
!*********************************************************************
use prec, only: pInt
use IO
implicit none
character(len=32), parameter :: myPart = 'microstructure'
integer(pInt), intent(in) :: file
integer(pInt), parameter :: maxNchunks = 7
integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) Nsections, section, constituent, i
character(len=64) tag
character(len=1024) line
Nsections = IO_countSections(file,myPart)
material_Nmicrostructure = Nsections
allocate(microstructure_name(Nsections)); microstructure_name = ''
allocate(microstructure_Nconstituents(Nsections))
microstructure_Nconstituents = IO_countTagInPart(file,myPart,'(constituent)',Nsections)
microstructure_maxNconstituents = maxval(microstructure_Nconstituents)
allocate(microstructure_phase (microstructure_maxNconstituents,Nsections)); microstructure_phase = 0_pInt
allocate(microstructure_texture (microstructure_maxNconstituents,Nsections)); microstructure_texture = 0_pInt
allocate(microstructure_fraction(microstructure_maxNconstituents,Nsections)); microstructure_fraction = 0.0_pReal
rewind(file)
line = ''
section = 0
do while (IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart
read(file,'(a1024)',END=100) line
enddo
do
read(file,'(a1024)',END=100) line
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1
constituent = 0
microstructure_name(section) = IO_getTag(line,'[',']')
endif
if (section > 0) then
positions = IO_stringPos(line,maxNchunks)
tag = IO_lc(IO_stringValue(line,positions,1)) ! extract key
select case(tag)
case ('(constituent)')
constituent = constituent + 1
do i=2,6,2
tag = IO_lc(IO_stringValue(line,positions,i))
select case (tag)
case('phase')
microstructure_phase(constituent,section) = IO_intValue(line,positions,i+1)
case('texture')
microstructure_texture(constituent,section) = IO_intValue(line,positions,i+1)
case('fraction')
microstructure_fraction(constituent,section) = IO_floatValue(line,positions,i+1)
end select
enddo
end select
endif
enddo
100 return
end subroutine
!*********************************************************************
subroutine material_parsePhase(file)
!*********************************************************************
use prec, only: pInt
use IO
implicit none
character(len=32), parameter :: myPart = 'phase'
integer(pInt), intent(in) :: file
integer(pInt), parameter :: maxNchunks = 2
integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) Nsections, section, s
character(len=64) tag
character(len=1024) line
Nsections= IO_countSections(file,myPart)
material_Nphase = Nsections
allocate(phase_name(Nsections)); phase_name = ''
allocate(phase_constitution(Nsections)); phase_constitution = ''
allocate(phase_constitutionInstance(Nsections)); phase_constitutionInstance = 0_pInt
allocate(phase_Noutput(Nsections))
phase_Noutput = IO_countTagInPart(file,myPart,'(output)',Nsections)
rewind(file)
line = ''
section = 0
do while (IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart
read(file,'(a1024)',END=100) line
enddo
do
read(file,'(a1024)',END=100) line
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1
phase_name(section) = IO_getTag(line,'[',']')
endif
if (section > 0) then
positions = IO_stringPos(line,maxNchunks)
tag = IO_lc(IO_stringValue(line,positions,1)) ! extract key
select case(tag)
case ('constitution')
phase_constitution(section) = IO_lc(IO_stringValue(line,positions,2))
do s = 1,section
if (phase_constitution(s) == phase_constitution(section)) &
phase_constitutionInstance(section) = phase_constitutionInstance(section) + 1 ! count instances
enddo
end select
endif
enddo
100 return
end subroutine
!*********************************************************************
subroutine material_parseTexture(file)
!*********************************************************************
use prec, only: pInt, pReal
use IO
use math, only: inRad
implicit none
character(len=32), parameter :: myPart = 'texture'
integer(pInt), intent(in) :: file
integer(pInt), parameter :: maxNchunks = 13
integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) Nsections, section, gauss, fiber, i
character(len=64) tag
character(len=1024) line
Nsections = IO_countSections(file,myPart)
material_Ntexture = Nsections
allocate(texture_name(Nsections)); texture_name = ''
allocate(texture_ODFfile(Nsections)); texture_ODFfile = ''
allocate(texture_symmetry(Nsections)); texture_symmetry = 1_pInt
allocate(texture_Ngauss(Nsections)); texture_Ngauss = 0_pInt
allocate(texture_Nfiber(Nsections)); texture_Nfiber = 0_pInt
texture_Ngauss = IO_countTagInPart(file,myPart,'(gauss)',Nsections)
texture_Nfiber = IO_countTagInPart(file,myPart,'(fiber)',Nsections)
texture_maxNgauss = maxval(texture_Ngauss)
texture_maxNfiber = maxval(texture_Nfiber)
allocate(texture_Gauss (5,texture_maxNgauss,Nsections)); texture_Gauss = 0.0_pReal
allocate(texture_Fiber (6,texture_maxNfiber,Nsections)); texture_Fiber = 0.0_pReal
rewind(file)
line = ''
section = 0
do while (IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart
read(file,'(a1024)',END=100) line
enddo
do
read(file,'(a1024)',END=100) line
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1
gauss = 0
fiber = 0
texture_name(section) = IO_getTag(line,'[',']')
endif
if (section > 0) then
positions = IO_stringPos(line,maxNchunks)
tag = IO_lc(IO_stringValue(line,positions,1)) ! extract key
select case(tag)
case ('hybridia')
texture_ODFfile(section) = IO_stringValue(line,positions,2)
case ('symmetry')
tag = IO_lc(IO_stringValue(line,positions,2))
select case (tag)
case('orthotropic')
texture_symmetry(section) = 4
case('monoclinic')
texture_symmetry(section) = 2
case default
texture_symmetry(section) = 1
end select
case ('(gauss)')
gauss = gauss + 1
do i = 2,10,2
tag = IO_lc(IO_stringValue(line,positions,i))
select case (tag)
case('phi1')
texture_Gauss(1,gauss,section) = IO_floatValue(line,positions,i+1)*inRad
case('phi')
texture_Gauss(2,gauss,section) = IO_floatValue(line,positions,i+1)*inRad
case('phi2')
texture_Gauss(3,gauss,section) = IO_floatValue(line,positions,i+1)*inRad
case('scatter')
texture_Gauss(4,gauss,section) = IO_floatValue(line,positions,i+1)*inRad
case('fraction')
texture_Gauss(5,gauss,section) = IO_floatValue(line,positions,i+1)
end select
enddo
case ('(fiber)')
fiber = fiber + 1
do i = 2,12,2
tag = IO_lc(IO_stringValue(line,positions,i))
select case (tag)
case('aplha1')
texture_Fiber(1,fiber,section) = IO_floatValue(line,positions,i+1)*inRad
case('alpha2')
texture_Fiber(2,fiber,section) = IO_floatValue(line,positions,i+1)*inRad
case('beta1')
texture_Fiber(3,fiber,section) = IO_floatValue(line,positions,i+1)*inRad
case('beta2')
texture_Fiber(4,fiber,section) = IO_floatValue(line,positions,i+1)*inRad
case('scatter')
texture_Fiber(5,fiber,section) = IO_floatValue(line,positions,i+1)*inRad
case('fraction')
texture_Fiber(6,fiber,section) = IO_floatValue(line,positions,i+1)
end select
enddo
end select
endif
enddo
100 return
end subroutine
!*********************************************************************
subroutine material_populateGrains()
!*********************************************************************
use prec, only: pInt, pReal
use math, only: math_sampleRandomOri, math_sampleGaussOri, math_sampleFiberOri, math_symmetricEulers
use mesh, only: mesh_element, mesh_maxNips, mesh_NcpElems, mesh_ipVolume, FE_Nips
use IO, only: IO_error, IO_hybridIA
implicit none
integer(pInt), dimension (:,:), allocatable :: Ngrains
integer(pInt), dimension (microstructure_maxNconstituents) :: NgrainsOfConstituent
real(pReal), dimension (:,:), allocatable :: volume
real(pReal), dimension (:), allocatable :: volFracOfGrain, phaseOfGrain
real(pReal), dimension (:,:), allocatable :: orientationOfGrain
real(pReal), dimension (3) :: orientation
real(pReal), dimension (3,3) :: symOrientation
integer(pInt) t,n,e,i,g,j,m,homog,micro,sgn
integer(pInt) phaseID,textureID,dGrains,myNgrains,myNorientations,NorientationsOfConstituent, &
grain,constituentGrain,symExtension
real(pReal) extreme,rnd
allocate(material_volFrac(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; material_volFrac = 0.0_pReal
allocate(material_phase(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; material_phase = 0_pInt
allocate(material_EulerAngles(3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; material_EulerAngles = 0.0_pReal
allocate(Ngrains(material_Nhomogenization,material_Nmicrostructure)); Ngrains = 0_pInt
allocate(volume(material_Nhomogenization,material_Nmicrostructure)); volume = 0.0_pReal
! count grains and total volume per homog/micro pair
do e = 1,mesh_NcpElems
homog = mesh_element(3,e)
micro = mesh_element(4,e)
if (homog < 1 .or. homog > material_Nhomogenization) & ! out of bounds
call IO_error(130,e,0,0)
if (micro < 1 .or. micro > material_Nmicrostructure) & ! out of bounds
call IO_error(140,e,0,0)
Ngrains(homog,micro) = Ngrains(homog,micro) + homogenization_Ngrains(homog) * FE_Nips(mesh_element(2,e))
volume(homog,micro) = volume(homog,micro) + sum(mesh_ipVolume(:,e))
enddo
allocate(volFracOfGrain(maxval(Ngrains))) ! reserve memory for maximum case
allocate(phaseOfGrain(maxval(Ngrains))) ! reserve memory for maximum case
allocate(orientationOfGrain(3,maxval(Ngrains))) ! reserve memory for maximum case
write (6,*)
write (6,*) 'MATERIAL grain population'
do homog = 1,material_Nhomogenization ! loop over homogenizations
dGrains = homogenization_Ngrains(homog) ! grain number per material point
do micro = 1,material_Nmicrostructure ! all pairs of homog and micro
if (Ngrains(homog,micro) > 0) then ! an active pair of homog and micro
myNgrains = Ngrains(homog,micro) ! assign short name
write (6,*)
write (6,'(a32,x,a32,x,i6)') homogenization_name(homog),microstructure_name(micro),myNgrains
! ----------------------------------------------------------------------------
volFracOfGrain = 0.0_pReal
grain = 0_pInt ! microstructure grain index
do e = 1,mesh_NcpElems ! check each element
if (mesh_element(3,e) == homog .and. mesh_element(4,e) == micro) then ! my combination of homog and micro
forall (i = 1:FE_Nips(mesh_element(2,e))) & ! loop over IPs
volFracOfGrain(grain+(i-1)*dGrains+1:grain+i*dGrains) = &
mesh_ipVolume(i,e)/volume(homog,micro)/dGrains ! assign IPvolfrac/Ngrains to grains
grain = grain + FE_Nips(mesh_element(2,e)) * dGrains ! wind forward by Nips*NgrainsPerIP
endif
enddo
write (6,*) 'now at grain count',grain
! ----------------------------------------------------------------------------
NgrainsOfConstituent = 0_pInt
forall (i = 1:microstructure_Nconstituents(micro)) &
NgrainsOfConstituent(i) = nint(microstructure_fraction(i,micro) * myNgrains, pInt)
write (6,*) 'NgrainsOfConstituent',NgrainsOfConstituent
do while (sum(NgrainsOfConstituent) /= myNgrains) ! total grain count over constituents wrong?
sgn = sign(1_pInt, myNgrains - sum(NgrainsOfConstituent)) ! direction of required change
extreme = 0.0_pReal
t = 0_pInt
do i = 1,microstructure_Nconstituents(micro) ! find largest deviator
if (sgn*log(NgrainsOfConstituent(i)/myNgrains/microstructure_fraction(i,micro)) > extreme) then
extreme = sgn*log(NgrainsOfConstituent(i)/myNgrains/microstructure_fraction(i,micro))
t = i
endif
enddo
NgrainsOfConstituent(t) = NgrainsOfConstituent(t) + sgn ! change that by one
end do
write (6,*) 'fixed NgrainsOfConstituent',NgrainsOfConstituent
! ----------------------------------------------------------------------------
phaseOfGrain = 0_pInt
orientationOfGrain = 0.0_pReal
grain = 0_pInt ! reset microstructure grain index
do i = 1,microstructure_Nconstituents(micro) ! loop over constituents
phaseID = microstructure_phase(i,micro)
textureID = microstructure_texture(i,micro)
phaseOfGrain(grain+1:grain+NgrainsOfConstituent(i)) = phaseID ! assign resp. phase
myNorientations = ceiling(float(NgrainsOfConstituent(i))/texture_symmetry(textureID)) ! max number of unique orientations (excl. symmetry)
write (6,'(a32,x,i6,x,f5.3,x,i6)') &
phase_name(phaseID),NgrainsOfConstituent(i),real(NgrainsOfConstituent(i)/myNgrains),myNorientations
constituentGrain = 0_pInt ! constituent grain index
! ---------
if (texture_ODFfile(textureID) == '') then ! dealing with texture components
! ---------
do t = 1,texture_Ngauss(textureID) ! loop over Gauss components
write (6,*) 'gauss',t,int(myNorientations*texture_Gauss(5,t,textureID))
do g = 1,int(myNorientations*texture_Gauss(5,t,textureID)) ! loop over required grain count
orientationOfGrain(:,grain+constituentGrain+g) = &
math_sampleGaussOri(texture_Gauss(1:3,t,textureID),&
texture_Gauss( 4,t,textureID))
enddo
constituentGrain = constituentGrain + int(myNorientations*texture_Gauss(5,t,textureID))
write (6,*) 'now at constitutent grain',constituentGrain
enddo
do t = 1,texture_Nfiber(textureID) ! loop over fiber components
write (6,*) 'fiber',t,int(myNorientations*texture_Fiber(6,t,textureID))
do g = 1,int(myNorientations*texture_Fiber(6,t,textureID)) ! loop over required grain count
orientationOfGrain(:,grain+constituentGrain+g) = &
math_sampleFiberOri(texture_Fiber(1:2,t,textureID),&
texture_Fiber(3:4,t,textureID),&
texture_Fiber( 5,t,textureID))
enddo
constituentGrain = constituentGrain + int(myNorientations*texture_fiber(6,t,textureID))
write (6,*) 'now at constitutent grain',constituentGrain
enddo
write (6,*) 'looping',constituentGrain+1,myNorientations
do j = constituentGrain+1,myNorientations ! fill remainder with random
orientationOfGrain(:,grain+j) = math_sampleRandomOri()
enddo
write (6,*) 'done...'
! ---------
else ! hybrid IA
! ---------
orientationOfGrain(:,grain+1:grain+myNorientations) = IO_hybridIA(myNorientations,texture_ODFfile(textureID))
if (all(orientationOfGrain(:,grain+1) == -1.0_pReal)) call IO_error(105)
constituentGrain = constituentGrain + myNorientations
endif
! ----------------------------------------------------------------------------
symExtension = texture_symmetry(textureID) - 1_pInt
if (symExtension > 0_pInt) then ! sample symmetry
constituentGrain = NgrainsOfConstituent(i)-myNorientations ! calc remainder of array
do j = 1,myNorientations ! loop over each "real" orientation
symOrientation = math_symmetricEulers(texture_symmetry(textureID),orientationOfGrain(:,j)) ! get symmetric equivalents
e = min(symExtension,constituentGrain) ! are we at end of constituent grain array?
if (e > 0_pInt) then
orientationOfGrain(:,grain+myNorientations+1+(j-1)*symExtension:&
grain+myNorientations+e+(j-1)*symExtension) = &
symOrientation(:,1:e)
constituentGrain = constituentGrain - e ! remainder shrinks by e
endif
enddo
endif
grain = grain + NgrainsOfConstituent(i) ! advance microstructure grain index
end do ! constituent
! ----------------------------------------------------------------------------
do i=1,myNgrains-1 ! walk thru grains
call random_number(rnd)
t = nint(rnd*(myNgrains-i)+i+0.5_pReal,pInt) ! select a grain in remaining list
m = phaseOfGrain(t) ! exchange current with random
phaseOfGrain(t) = phaseOfGrain(i)
phaseOfGrain(i) = m
orientation = orientationOfGrain(:,t)
orientationOfGrain(:,t) = orientationOfGrain(:,i)
orientationOfGrain(:,i) = orientation
end do
!calc fraction after weighing with volumePerGrain
!exchange in MC steps to improve result...
! ----------------------------------------------------------------------------
grain = 0_pInt ! microstructure grain index
do e = 1,mesh_NcpElems ! check each element
if (mesh_element(3,e) == homog .and. mesh_element(4,e) == micro) then ! my combination of homog and micro
forall (i = 1:FE_Nips(mesh_element(2,e)), g = 1:dGrains) ! loop over IPs and grains
material_volFrac(g,i,e) = volFracOfGrain(grain+(i-1)*dGrains+g)
material_phase(g,i,e) = phaseOfGrain(grain+(i-1)*dGrains+g)
material_EulerAngles(:,g,i,e) = orientationOfGrain(:,grain+(i-1)*dGrains+g)
end forall
write (6,*) e
write (6,*) material_phase(:,:,e)
write (6,*) material_EulerAngles(:,:,:,e)
grain = grain + FE_Nips(mesh_element(2,e)) * dGrains ! wind forward by Nips*NgrainsPerIP
endif
enddo
endif ! active homog,micro pair
enddo
enddo
deallocate(volFracOfGrain)
deallocate(phaseOfGrain)
deallocate(orientationOfGrain)
return
end subroutine
END MODULE

23
trunk/math.f90
View File

@ -415,6 +415,23 @@ end subroutine
END FUNCTION END FUNCTION
!**************************************************************************
! range of integers starting at one
!**************************************************************************
PURE FUNCTION math_range(N)
use prec, only: pInt
implicit none
integer(pInt), intent(in) :: N
integer(pInt) i
integer(pInt), dimension(N) :: math_range
forall (i=1:N) math_range(i) = i
return
END FUNCTION
!************************************************************************** !**************************************************************************
! second rank identity tensor of specified dimension ! second rank identity tensor of specified dimension
!************************************************************************** !**************************************************************************
@ -1366,7 +1383,7 @@ math_sampleFiberOri = math_RtoEuler(math_mul33x33(pRot,math_mul33x33(fRot,oRot))
use prec, only: pReal, pInt use prec, only: pReal, pInt
implicit none implicit none
character(len=80), intent(in) :: sym integer(pInt), intent(in) :: sym
real(pReal), dimension(3), intent(in) :: Euler real(pReal), dimension(3), intent(in) :: Euler
real(pReal), dimension(3,3) :: math_symmetricEulers real(pReal), dimension(3,3) :: math_symmetricEulers
integer(pInt) i,j integer(pInt) i,j
@ -1386,9 +1403,9 @@ math_sampleFiberOri = math_RtoEuler(math_mul33x33(pRot,math_mul33x33(fRot,oRot))
forall (i=1:3,j=1:3) math_symmetricEulers(j,i) = modulo(math_symmetricEulers(j,i),2.0_pReal*pi) forall (i=1:3,j=1:3) math_symmetricEulers(j,i) = modulo(math_symmetricEulers(j,i),2.0_pReal*pi)
select case (sym) select case (sym)
case ('orthotropic') ! all done case (4) ! all done
case ('monoclinic') ! return only first case (2) ! return only first
math_symmetricEulers(:,2:3) = 0.0_pReal math_symmetricEulers(:,2:3) = 0.0_pReal
case default ! return blank case default ! return blank

22
trunk/mesh.f90
View File

@ -681,8 +681,8 @@
use prec, only: pInt use prec, only: pInt
use IO, only: IO_error,IO_open_InputFile use IO, only: IO_error,IO_open_InputFile
use FEsolving, only: parallelExecution
use FEsolving, only: FE_get_solverSymmetry
implicit none implicit none
integer(pInt), parameter :: fileUnit = 222 integer(pInt), parameter :: fileUnit = 222
@ -703,7 +703,6 @@
! call to various subroutines to parse the stuff from the input file... ! call to various subroutines to parse the stuff from the input file...
if (IO_open_inputFile(fileUnit)) then if (IO_open_inputFile(fileUnit)) then
call FE_get_solverSymmetry(fileUnit)
call mesh_get_meshDimensions(fileUnit) call mesh_get_meshDimensions(fileUnit)
call mesh_build_nodeMapping(fileUnit) call mesh_build_nodeMapping(fileUnit)
call mesh_build_elemMapping(fileUnit) call mesh_build_elemMapping(fileUnit)
@ -718,6 +717,8 @@
call mesh_build_ipAreas() call mesh_build_ipAreas()
call mesh_tell_statistics() call mesh_tell_statistics()
close (fileUnit) close (fileUnit)
parallelExecution = (mesh_Nelems == mesh_NcpElems) ! plus potential killer from non-local constitutive
else else
call IO_error(100) ! cannot open input file call IO_error(100) ! cannot open input file
endif endif
@ -1442,10 +1443,10 @@ matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matc
!*********************************************************** !***********************************************************
! calculation of IP volume ! calculation of IP interface areas
! !
! allocate globals ! allocate globals
! _ipVolume ! _ipArea, _ipAreaNormal
!*********************************************************** !***********************************************************
SUBROUTINE mesh_build_ipAreas() SUBROUTINE mesh_build_ipAreas()
@ -1492,7 +1493,8 @@ matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matc
!*********************************************************** !***********************************************************
SUBROUTINE mesh_tell_statistics() SUBROUTINE mesh_tell_statistics()
use prec, only: pInt use prec, only: pInt
use math, only: math_range
use IO, only: IO_error use IO, only: IO_error
implicit none implicit none
@ -1540,11 +1542,13 @@ matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matc
write (6,*) mesh_maxNsubNodes, " : max number of (additional) subnodes in any CP element" write (6,*) mesh_maxNsubNodes, " : max number of (additional) subnodes in any CP element"
write (6,*) mesh_maxNsharedElems, " : max number of CP elements sharing a node" write (6,*) mesh_maxNsharedElems, " : max number of CP elements sharing a node"
write (6,*) write (6,*)
write (6,*) "Input Parser: MATERIAL/TEXTURE" write (6,*) "Input Parser: HOMOGENIZATION/MICROSTRUCTURE"
write (6,*) write (6,*)
write (6,*) mesh_maxValStateVar(1), " : maximum material index" write (6,*) mesh_maxValStateVar(1), " : maximum homogenization index"
write (6,*) mesh_maxValStateVar(2), " : maximum texture index" write (6,*) mesh_maxValStateVar(2), " : maximum microstructure index"
write (6,*) write (6,*)
write (fmt,"(a,i3,a)") "(9(x),a1,x,",mesh_maxValStateVar(2),"(i8))"
write (6,fmt) "+",math_range(mesh_maxValStateVar(2))
write (fmt,"(a,i3,a)") "(i8,x,a1,x,",mesh_maxValStateVar(2),"(i8))" write (fmt,"(a,i3,a)") "(i8,x,a1,x,",mesh_maxValStateVar(2),"(i8))"
do i=1,mesh_maxValStateVar(1) ! loop over all (possibly assigned) materials do i=1,mesh_maxValStateVar(1) ! loop over all (possibly assigned) materials
write (6,fmt) i,"|",mesh_MatTex(i,:) ! loop over all (possibly assigned) textures write (6,fmt) i,"|",mesh_MatTex(i,:) ! loop over all (possibly assigned) textures

108
trunk/mpie_cpfem_marc2007r1.f90
View File

@ -4,18 +4,18 @@
! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts ! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts
! MPI fuer Eisenforschung, Duesseldorf ! MPI fuer Eisenforschung, Duesseldorf
! !
! last modified: 27.11.2008 ! last modified: 22.11.2008
!******************************************************************** !********************************************************************
! Usage: ! Usage:
! - choose material as hypela2 ! - choose material as hypela2
! - set statevariable 2 to index of material ! - set statevariable 2 to index of homogenization
! - set statevariable 3 to index of texture ! - set statevariable 3 to index of microstructure
! - choose output of user variables if desired ! - make sure the file "material.config" exists in the working
! - make sure the file "mattex.mpie" exists in the working
! directory ! directory
! - use nonsymmetric option for solver (e.g. direct ! - use nonsymmetric option for solver (e.g. direct
! profile or multifrontal sparse, the latter seems ! profile or multifrontal sparse, the latter seems
! to be faster!) ! to be faster!)
! - in case of ddm a symmetric solver has to be used
!******************************************************************** !********************************************************************
! Marc subroutines used: ! Marc subroutines used:
! - hypela2 ! - hypela2
@ -34,15 +34,16 @@
include "FEsolving.f90" ! uses prec, IO include "FEsolving.f90" ! uses prec, IO
include "mesh.f90" ! uses prec, IO, math, FEsolving include "mesh.f90" ! uses prec, IO, math, FEsolving
include "lattice.f90" ! uses prec, math include "lattice.f90" ! uses prec, math
include "material.f90" ! uses prec, math, IO, mesh
include "constitutive_phenomenological.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug
! include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO
include "CPFEM.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite include "CPFEM.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite
!
SUBROUTINE hypela2(d,g,e,de,s,t,dt,ngens,n,nn,kcus,matus,ndi,& SUBROUTINE hypela2(d,g,e,de,s,t,dt,ngens,n,nn,kcus,matus,ndi,&
nshear,disp,dispt,coord,ffn,frotn,strechn,eigvn,ffn1,& nshear,disp,dispt,coord,ffn,frotn,strechn,eigvn,ffn1,&
frotn1,strechn1,eigvn1,ncrd,itel,ndeg,ndm,& frotn1,strechn1,eigvn1,ncrd,itel,ndeg,ndm,&
nnode,jtype,lclass,ifr,ifu) nnode,jtype,lclass,ifr,ifu)
!******************************************************************** !********************************************************************
! This is the Marc material routine ! This is the Marc material routine
!******************************************************************** !********************************************************************
@ -123,45 +124,34 @@
!2 continue !2 continue
!3 continue !3 continue
! !
use prec, only: pReal,pInt, ijaco use prec, only: pReal,pInt, ijaco
use FEsolving use FEsolving
use CPFEM, only: CPFEM_general use CPFEM, only: CPFEM_general
use math, only: invnrmMandel use math, only: invnrmMandel
use debug
implicit real(pReal) (a-h,o-z) !
integer(pInt) computationMode implicit none
! ** Start of generated type statements **
real(pReal) coord, d, de, disp, dispt, dt, e, eigvn, eigvn1, ffn, ffn1
real(pReal) frotn, frotn1, g
integer(pInt) ifr, ifu, itel, jtype, kcus, lclass, matus, n, ncrd, ndeg
integer(pInt) ndi, ndm, ngens, nn, nnode, nshear
real(pReal) s, strechn, strechn1, t
! ** End of generated type statements **
!
dimension e(*),de(*),t(*),dt(*),g(*),d(ngens,*),s(*), n(2),coord(ncrd,*),disp(ndeg,*),matus(2),dispt(ndeg,*),ffn(itel,*),& dimension e(*),de(*),t(*),dt(*),g(*),d(ngens,*),s(*), n(2),coord(ncrd,*),disp(ndeg,*),matus(2),dispt(ndeg,*),ffn(itel,*),&
frotn(itel,*),strechn(itel),eigvn(itel,*),ffn1(itel,*),frotn1(itel,*),strechn1(itel),eigvn1(itel,*),kcus(2) frotn(itel,*),strechn(itel),eigvn(itel,*),ffn1(itel,*),frotn1(itel,*),strechn1(itel),eigvn1(itel,*),kcus(2)
! Marc common blocks are in fixed format so they have to be pasted in here ! Marc common blocks are in fixed format so they have to be reformated to free format (f90)
! Beware of changes in newer Marc versions -- these are from 2005r3 ! Beware of changes in newer Marc versions
! concom is needed for inc, subinc, ncycle, lovl
! include 'concom'
common/marc_concom/ &
iacous, iasmbl, iautth, ibear, icompl, iconj, icreep, ideva(50), idyn, idynt,&
ielas, ielcma, ielect, iform, ifour, iharm, ihcps, iheat, iheatt, ihresp,&
ijoule, ilem, ilnmom, iloren, inc, incext, incsub, ipass, iplres, ipois,&
ipoist, irpflo, ismall, ismalt, isoil, ispect, ispnow, istore, iswep, ithcrp,&
itherm, iupblg, iupdat, jacflg, jel, jparks, largst, lfond, loadup, loaduq,&
lodcor, lovl, lsub, magnet, ncycle, newtnt, newton, noshr, linear, ivscpl,&
icrpim, iradrt, ipshft, itshr, iangin, iupmdr, iconjf, jincfl, jpermg, jhour,&
isolvr, jritz, jtable, jshell, jdoubl, jform, jcentr, imini, kautth, iautof,&
ibukty, iassum, icnstd, icnstt, kmakmas, imethvp,iradrte,iradrtp, iupdate,iupdatp,&
ncycnt, marmen ,idynme, ihavca, ispf, kmini, imixed, largtt, kdoela, iautofg,&
ipshftp,idntrc, ipore, jtablm, jtablc, isnecma,itrnspo,imsdif, jtrnspo,mcnear,&
imech, imecht, ielcmat, ielectt,magnett, imsdift,noplas, jtabls, jactch, jtablth,&
kgmsto ,jpzo, ifricsh, iremkin,iremfor, ishearp,jspf, machining, jlshell,icompsol,&
iupblgfo,jcondir,nstcrp, nactive,ipassref, nstspnt,ibeart,icheckmpc, noline, icuring,&
ishrink,ioffsflg,isetoff, ioffsetm,iharmt, inc_incdat,iautspc,ibrake, icbush ,istream_input,&
iprsinp,ivlsinp,ifirst_time,ipin_m,jgnstr_glb, imarc_return,iqvcinp,nqvceid,istpnx,imicro1
! creeps is needed for timinc (time increment) include "concom2007r1" ! concom is needed for inc, subinc, ncycle, lovl
! include 'creeps' include "creeps2007r1" ! creeps is needed for timinc (time increment)
common/marc_creeps/ &
cptim,timinc,timinc_p,timinc_s,timincm,timinc_a,timinc_b,creept(33),icptim,icfte,icfst,& integer(pInt) computationMode,i
icfeq,icftm,icetem,mcreep,jcreep,icpa,icftmp,icfstr,icfqcp,icfcpm,icrppr,icrcha,icpb,iicpmt,iicpa
! write(6,'(3(3(f10.3,x),/))') ffn1(:,1),ffn1(:,2),ffn1(:,3)
if (inc == 0) then if (inc == 0) then
cycleCounter = 4 cycleCounter = 4
@ -170,16 +160,22 @@
if (theCycle /= ncycle .or. theLovl /= lovl) then if (theCycle /= ncycle .or. theLovl /= lovl) then
cycleCounter = cycleCounter+1 ! ping pong cycleCounter = cycleCounter+1 ! ping pong
outdatedFFN1 = .false. outdatedFFN1 = .false.
write (6,*) n(1),nn,'cycleCounter',cycleCounter
call debug_info() ! output of debugging/performance statistics of former
debug_cutbackDistribution = 0_pInt ! initialize debugging data
debug_InnerLoopDistribution = 0_pInt
debug_OuterLoopDistribution = 0_pInt
endif endif
endif endif
if (cptim > theTime .or. theInc /= inc) then ! reached convergence if (cptim > theTime .or. theInc /= inc) then ! reached convergence
lastIncConverged = .true. lastIncConverged = .true.
outdatedByNewInc = .true. outdatedByNewInc = .true.
write (6,*) n(1),nn,'lastIncConverged + outdated'
endif endif
if (mod(cycleCounter,2) /= 0) computationMode = 4 ! recycle if (mod(cycleCounter,2) /= 0) computationMode = 4 ! recycle in odd cycles
if (mod(cycleCounter,4) == 2) computationMode = 3 ! collect if (mod(cycleCounter,4) == 2) computationMode = 3 ! collect in 2,6,10,...
if (mod(cycleCounter,4) == 0) computationMode = 2 ! compute if (mod(cycleCounter,4) == 0) computationMode = 2 ! compute in 0,4,8,...
if (computationMode == 4 .and. ncycle == 0 .and. .not. lastIncConverged) & if (computationMode == 4 .and. ncycle == 0 .and. .not. lastIncConverged) &
computationMode = 6 ! recycle but restore known good consistent tangent computationMode = 6 ! recycle but restore known good consistent tangent
if (computationMode == 4 .and. lastIncConverged) then if (computationMode == 4 .and. lastIncConverged) then
@ -191,10 +187,6 @@
outdatedByNewInc = .false. outdatedByNewInc = .false.
endif endif
if (computationMode == 2 .and. outdatedFFN1) then
computationMode = 4 ! return odd results to force new vyvle
endif
theTime = cptim ! record current starting time theTime = cptim ! record current starting time
theInc = inc ! record current increment number theInc = inc ! record current increment number
theCycle = ncycle ! record current cycle count theCycle = ncycle ! record current cycle count
@ -207,6 +199,7 @@
forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*d(1:ngens,i)*invnrmMandel(1:ngens) forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*d(1:ngens,i)*invnrmMandel(1:ngens)
s(1:ngens) = s(1:ngens)*invnrmMandel(1:ngens) s(1:ngens) = s(1:ngens)*invnrmMandel(1:ngens)
if(symmetricSolver) d(1:ngens,1:ngens) = 0.5_pReal*(d(1:ngens,1:ngens)+transpose(d(1:ngens,1:ngens))) if(symmetricSolver) d(1:ngens,1:ngens) = 0.5_pReal*(d(1:ngens,1:ngens)+transpose(d(1:ngens,1:ngens)))
return return
END SUBROUTINE END SUBROUTINE
@ -220,22 +213,22 @@
! select a variable contour plotting (user subroutine). ! select a variable contour plotting (user subroutine).
! !
! v variable ! v variable
! s (idss) stress array ! s (idss) stress array
! sp stresses in preferred direction ! sp stresses in preferred direction
! etot total strain (generalized) ! etot total strain (generalized)
! eplas total plastic strain ! eplas total plastic strain
! ecreep total creep strain ! ecreep total creep strain
! t current temperature ! t current temperature
! m element number ! m element number
! nn integration point number ! nn integration point number
! layer layer number ! layer layer number
! ndi (3) number of direct stress components ! ndi (3) number of direct stress components
! nshear (3) number of shear stress components ! nshear (3) number of shear stress components
! !
!******************************************************************** !********************************************************************
use prec, only: pReal,pInt use prec, only: pReal,pInt
use CPFEM, only: CPFEM_results, CPFEM_Nresults use CPFEM, only: CPFEM_results, CPFEM_Nresults
use constitutive, only: constitutive_maxNresults use constitutive, only: constitutive_maxSizePostResults
use mesh, only: mesh_FEasCP use mesh, only: mesh_FEasCP
implicit none implicit none
! !
@ -244,10 +237,9 @@
integer(pInt) m, nn, layer, ndi, nshear, jpltcd integer(pInt) m, nn, layer, ndi, nshear, jpltcd
! !
! assign result variable ! assign result variable
v=CPFEM_results(mod(jpltcd-1_pInt, CPFEM_Nresults+constitutive_maxNresults)+1_pInt,& v = CPFEM_results(mod(jpltcd-1_pInt, CPFEM_Nresults+constitutive_maxSizePostResults)+1_pInt,&
(jpltcd-1_pInt)/(CPFEM_Nresults+constitutive_maxNresults)+1_pInt,& (jpltcd-1_pInt)/(CPFEM_Nresults+constitutive_maxSizePostResults)+1_pInt,&
nn, mesh_FEasCP('elem', m)) nn, mesh_FEasCP('elem', m))
return return
END SUBROUTINE END SUBROUTINE
! !
@ -300,4 +292,4 @@
! end if ! end if
! end if ! end if
! return ! return
! end ! end

36
trunk/mpie_cpfem_marc2008r1.f90
View File

@ -36,10 +36,10 @@
include "mesh.f90" ! uses prec, IO, math, FEsolving include "mesh.f90" ! uses prec, IO, math, FEsolving
include "lattice.f90" ! uses prec, math include "lattice.f90" ! uses prec, math
include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug
! include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO
include "CPFEM.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite include "CPFEM.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite
!
! logical, parameter :: parallelExecution = .true.
SUBROUTINE hypela2(d,g,e,de,s,t,dt,ngens,n,nn,kcus,matus,ndi,& SUBROUTINE hypela2(d,g,e,de,s,t,dt,ngens,n,nn,kcus,matus,ndi,&
nshear,disp,dispt,coord,ffn,frotn,strechn,eigvn,ffn1,& nshear,disp,dispt,coord,ffn,frotn,strechn,eigvn,ffn1,&
frotn1,strechn1,eigvn1,ncrd,itel,ndeg,ndm,& frotn1,strechn1,eigvn1,ncrd,itel,ndeg,ndm,&
@ -128,6 +128,7 @@
use FEsolving use FEsolving
use CPFEM, only: CPFEM_general use CPFEM, only: CPFEM_general
use math, only: invnrmMandel use math, only: invnrmMandel
use debug
! !
implicit none implicit none
@ -146,9 +147,9 @@
! Marc common blocks are in fixed format so they have to be reformated to free format (f90) ! Marc common blocks are in fixed format so they have to be reformated to free format (f90)
! Beware of changes in newer Marc versions -- these are from 2005r3 ! Beware of changes in newer Marc versions -- these are from 2005r3
! concom is needed for inc, subinc, ncycle, lovl ! concom is needed for inc, subinc, ncycle, lovl
include "concom_f90" include "concom2008r1"
! creeps is needed for timinc (time increment) ! creeps is needed for timinc (time increment)
include "creeps_f90" include "creeps2008r1"
! !
integer(pInt) computationMode,i integer(pInt) computationMode,i
! !
@ -159,6 +160,10 @@
if (theCycle /= ncycle .or. theLovl /= lovl) then if (theCycle /= ncycle .or. theLovl /= lovl) then
cycleCounter = cycleCounter+1 ! ping pong cycleCounter = cycleCounter+1 ! ping pong
outdatedFFN1 = .false. outdatedFFN1 = .false.
call debug_info() ! output of debugging/performance statistics of former
debug_cutbackDistribution = 0_pInt ! initialize debugging data
debug_InnerLoopDistribution = 0_pInt
debug_OuterLoopDistribution = 0_pInt
endif endif
endif endif
if (cptim > theTime .or. theInc /= inc) then ! reached convergence if (cptim > theTime .or. theInc /= inc) then ! reached convergence
@ -179,9 +184,6 @@
computationMode = 1 ! compute and age former results computationMode = 1 ! compute and age former results
outdatedByNewInc = .false. outdatedByNewInc = .false.
endif endif
if (computationMode == 2 .and. outdatedFFN1) then
computationMode = 4 ! return odd results to force new vyvle
endif
theTime = cptim ! record current starting time theTime = cptim ! record current starting time
theInc = inc ! record current increment number theInc = inc ! record current increment number
@ -208,17 +210,17 @@
! select a variable contour plotting (user subroutine). ! select a variable contour plotting (user subroutine).
! !
! v variable ! v variable
! s (idss) stress array ! s (idss) stress array
! sp stresses in preferred direction ! sp stresses in preferred direction
! etot total strain (generalized) ! etot total strain (generalized)
! eplas total plastic strain ! eplas total plastic strain
! ecreep total creep strain ! ecreep total creep strain
! t current temperature ! t current temperature
! m element number ! m element number
! nn integration point number ! nn integration point number
! layer layer number ! layer layer number
! ndi (3) number of direct stress components ! ndi (3) number of direct stress components
! nshear (3) number of shear stress components ! nshear (3) number of shear stress components
! !
!******************************************************************** !********************************************************************
use prec, only: pReal,pInt use prec, only: pReal,pInt
@ -232,8 +234,8 @@
integer(pInt) m, nn, layer, ndi, nshear, jpltcd integer(pInt) m, nn, layer, ndi, nshear, jpltcd
! !
! assign result variable ! assign result variable
v=CPFEM_results(mod(jpltcd-1_pInt, CPFEM_Nresults+constitutive_maxNresults)+1_pInt,& v=CPFEM_results(mod(jpltcd-1_pInt, CPFEM_Nresults+constitutive_maxSizePostResults)+1_pInt,&
(jpltcd-1_pInt)/(CPFEM_Nresults+constitutive_maxNresults)+1_pInt,& (jpltcd-1_pInt)/(CPFEM_Nresults+constitutive_maxSizePostResults)+1_pInt,&
nn, mesh_FEasCP('elem', m)) nn, mesh_FEasCP('elem', m))
return return
END SUBROUTINE END SUBROUTINE

6
trunk/prec.f90
View File

@ -9,6 +9,10 @@
integer, parameter :: pReal = selected_real_kind(15,300) ! 15 significant digits, up to 1e+-300 integer, parameter :: pReal = selected_real_kind(15,300) ! 15 significant digits, up to 1e+-300
integer, parameter :: pInt = selected_int_kind(9) ! up to +- 1e9 integer, parameter :: pInt = selected_int_kind(9) ! up to +- 1e9
type :: p_vec
real(pReal), dimension(:), pointer :: p
end type p_vec
! *** Strain increment considered significant *** ! *** Strain increment considered significant ***
real(pReal), parameter :: relevantStrain = 1.0e-7_pReal real(pReal), parameter :: relevantStrain = 1.0e-7_pReal
@ -19,7 +23,7 @@
real(pReal), parameter :: pert_Fg = 1.0e-6_pReal ! strain perturbation for FEM Jacobi real(pReal), parameter :: pert_Fg = 1.0e-6_pReal ! strain perturbation for FEM Jacobi
integer(pInt), parameter :: nOuter = 20_pInt ! outer loop limit 20 integer(pInt), parameter :: nOuter = 20_pInt ! outer loop limit 20
integer(pInt), parameter :: nInner = 200_pInt ! inner loop limit 200 integer(pInt), parameter :: nInner = 200_pInt ! inner loop limit 200
real(pReal), parameter :: reltol_Outer = 1.0e-6_pReal ! relative tolerance in outer loop (state) real(pReal), parameter :: reltol_Outer = 1.0e-5_pReal ! relative tolerance in outer loop (state)
real(pReal), parameter :: reltol_Inner = 1.0e-6_pReal ! relative tolerance in inner loop (Lp) real(pReal), parameter :: reltol_Inner = 1.0e-6_pReal ! relative tolerance in inner loop (Lp)
real(pReal), parameter :: abstol_Inner = 1.0e-8_pReal ! absolute tolerance in inner loop (Lp) real(pReal), parameter :: abstol_Inner = 1.0e-8_pReal ! absolute tolerance in inner loop (Lp)
! !