diff --git a/trunk/FEsolving.f90 b/trunk/FEsolving.f90 index 21f519a3f..577c795b0 100644 --- a/trunk/FEsolving.f90 +++ b/trunk/FEsolving.f90 @@ -9,8 +9,9 @@ integer(pInt) cycleCounter integer(pInt) theInc,theCycle,theLovl real(pReal) theTime - logical :: lastIncConverged = .false.,outdatedByNewInc = .false., outdatedFFN1 = .false. + logical :: lastIncConverged = .false.,outdatedByNewInc = .false.,outdatedFFN1 = .false. logical :: symmetricSolver = .false. + logical :: parallelExecution = .true. CONTAINS @@ -18,36 +19,36 @@ !*********************************************************** ! determine wether a symmetric solver is used !*********************************************************** - subroutine FE_get_solverSymmetry(unit) + subroutine FE_init() use prec, only: pInt use IO implicit none - integer(pInt) unit - integer(pInt), dimension (133) :: pos - character*300 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 + integer(pInt), parameter :: fileunit = 222 + integer(pInt), dimension (1+2*2) :: pos + character(len=1024) line -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 diff --git a/trunk/IO.f90 b/trunk/IO.f90 index 702908630..4ef768faa 100644 --- a/trunk/IO.f90 +++ b/trunk/IO.f90 @@ -244,6 +244,132 @@ 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 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 ! return array containing number of parts found and @@ -559,48 +685,90 @@ ! and terminate the Marc run with exit #9xxx ! 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 debug 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 select case (ID) - case (100) + case (0) msg='Unable to open input file.' + case (100) + msg='Error reading from configuration file.' + case (105) + msg='Error reading from ODF file.' case (110) - msg='No materials specified via State Variable 2.' + msg='No homogenization specified via State Variable 2.' 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) - 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) - msg='This material can only be used with elements with three direct stress components' - case (400) - msg='Unknown alloy number specified' + msg='This material can only be used with elements with three direct stress components.' case (500) - msg='Unknown lattice type specified' + msg='Unknown lattice type specified.' case (600) - msg='Convergence not reached' + msg='Convergence not reached.' + case (610) + msg='Stress loop not converged.' case (650) - msg='Polar decomposition failed' + msg='Polar decomposition failed.' case (700) - msg='Singular matrix in stress iteration' + msg='Singular matrix in stress iteration.' case (800) msg='GIA requires 8 grains per IP (bonehead, you!)' case default - msg='Unknown error number' + msg='Unknown error number...' end select !$OMP CRITICAL (write2out) - write(6,*) 'MPIE Material Routine Ver. 0.0 by the coding team' - write(6,*) + write(6,*) '+------------------------------+' 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 flush(6) diff --git a/trunk/concom2007r1 b/trunk/concom2007r1 new file mode 100644 index 000000000..dd10f58d5 --- /dev/null +++ b/trunk/concom2007r1 @@ -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) +! +!*********************************************************************** diff --git a/trunk/constitutive_phenomenological.f90 b/trunk/constitutive_phenomenological.f90 new file mode 100644 index 000000000..79d1c3918 --- /dev/null +++ b/trunk/constitutive_phenomenological.f90 @@ -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 + 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 diff --git a/trunk/creeps2007r1 b/trunk/creeps2007r1 new file mode 100644 index 000000000..31dff6093 --- /dev/null +++ b/trunk/creeps2007r1 @@ -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 +! +!*********************************************************************** diff --git a/trunk/debug.f90 b/trunk/debug.f90 index 4de13b164..724108623 100644 --- a/trunk/debug.f90 +++ b/trunk/debug.f90 @@ -6,10 +6,11 @@ implicit none - integer(pInt), dimension(nCutback+1) :: debug_cutbackDistribution - integer(pInt), dimension(nInner) :: debug_InnerLoopDistribution - integer(pInt), dimension(nOuter) :: debug_OuterLoopDistribution + integer(pInt), dimension(nCutback+1) :: debug_cutbackDistribution = 0_pInt + integer(pInt), dimension(nInner) :: debug_InnerLoopDistribution = 0_pInt + integer(pInt), dimension(nOuter) :: debug_OuterLoopDistribution = 0_pInt logical :: debugger = .false. + logical :: distribution_init = .false. CONTAINS diff --git a/trunk/material.config b/trunk/material.config new file mode 100644 index 000000000..ec732cce4 --- /dev/null +++ b/trunk/material.config @@ -0,0 +1,103 @@ +##################### + +##################### + +[SX] +type Taylor +Ngrains 1 + +[RGC] +type RGC +Ngrains 8 + +[Taylor4] +type Taylor +Ngrains 4 + +##################### + +##################### + +[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 + + +##################### + +##################### + +[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 + +##################### + +##################### + +[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 + + + diff --git a/trunk/material.f90 b/trunk/material.f90 new file mode 100644 index 000000000..00bf09405 --- /dev/null +++ b/trunk/material.f90 @@ -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 diff --git a/trunk/math.f90 b/trunk/math.f90 index c91286f9f..8a4eab5a7 100644 --- a/trunk/math.f90 +++ b/trunk/math.f90 @@ -415,6 +415,23 @@ end subroutine 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 !************************************************************************** @@ -1366,7 +1383,7 @@ math_sampleFiberOri = math_RtoEuler(math_mul33x33(pRot,math_mul33x33(fRot,oRot)) use prec, only: pReal, pInt implicit none - character(len=80), intent(in) :: sym + integer(pInt), intent(in) :: sym real(pReal), dimension(3), intent(in) :: Euler real(pReal), dimension(3,3) :: math_symmetricEulers 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) 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 case default ! return blank diff --git a/trunk/mesh.f90 b/trunk/mesh.f90 index 6b138df2b..d1d66092b 100644 --- a/trunk/mesh.f90 +++ b/trunk/mesh.f90 @@ -681,8 +681,8 @@ use prec, only: pInt use IO, only: IO_error,IO_open_InputFile - - use FEsolving, only: FE_get_solverSymmetry + use FEsolving, only: parallelExecution + implicit none integer(pInt), parameter :: fileUnit = 222 @@ -703,7 +703,6 @@ ! call to various subroutines to parse the stuff from the input file... if (IO_open_inputFile(fileUnit)) then - call FE_get_solverSymmetry(fileUnit) call mesh_get_meshDimensions(fileUnit) call mesh_build_nodeMapping(fileUnit) call mesh_build_elemMapping(fileUnit) @@ -718,6 +717,8 @@ call mesh_build_ipAreas() call mesh_tell_statistics() close (fileUnit) + + parallelExecution = (mesh_Nelems == mesh_NcpElems) ! plus potential killer from non-local constitutive else call IO_error(100) ! cannot open input file 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 -! _ipVolume +! _ipArea, _ipAreaNormal !*********************************************************** 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() - use prec, only: pInt + use prec, only: pInt + use math, only: math_range use IO, only: IO_error 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_maxNsharedElems, " : max number of CP elements sharing a node" write (6,*) - write (6,*) "Input Parser: MATERIAL/TEXTURE" + write (6,*) "Input Parser: HOMOGENIZATION/MICROSTRUCTURE" write (6,*) - write (6,*) mesh_maxValStateVar(1), " : maximum material index" - write (6,*) mesh_maxValStateVar(2), " : maximum texture index" + write (6,*) mesh_maxValStateVar(1), " : maximum homogenization index" + write (6,*) mesh_maxValStateVar(2), " : maximum microstructure index" 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))" 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 diff --git a/trunk/mpie_cpfem_marc2007r1.f90 b/trunk/mpie_cpfem_marc2007r1.f90 index 6f4e6ea0b..ca9f8f2f8 100644 --- a/trunk/mpie_cpfem_marc2007r1.f90 +++ b/trunk/mpie_cpfem_marc2007r1.f90 @@ -4,18 +4,18 @@ ! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts ! MPI fuer Eisenforschung, Duesseldorf ! -! last modified: 27.11.2008 +! last modified: 22.11.2008 !******************************************************************** ! Usage: ! - choose material as hypela2 -! - set statevariable 2 to index of material -! - set statevariable 3 to index of texture -! - choose output of user variables if desired -! - make sure the file "mattex.mpie" exists in the working +! - set statevariable 2 to index of homogenization +! - set statevariable 3 to index of microstructure +! - make sure the file "material.config" exists in the working ! directory ! - use nonsymmetric option for solver (e.g. direct ! profile or multifrontal sparse, the latter seems ! to be faster!) +! - in case of ddm a symmetric solver has to be used !******************************************************************** ! Marc subroutines used: ! - hypela2 @@ -34,15 +34,16 @@ include "FEsolving.f90" ! uses prec, IO include "mesh.f90" ! uses prec, IO, math, FEsolving 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 "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO 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,& nshear,disp,dispt,coord,ffn,frotn,strechn,eigvn,ffn1,& frotn1,strechn1,eigvn1,ncrd,itel,ndeg,ndm,& nnode,jtype,lclass,ifr,ifu) + !******************************************************************** ! This is the Marc material routine !******************************************************************** @@ -123,45 +124,34 @@ !2 continue !3 continue ! - use prec, only: pReal,pInt, ijaco use FEsolving use CPFEM, only: CPFEM_general use math, only: invnrmMandel - - implicit real(pReal) (a-h,o-z) - integer(pInt) computationMode - + use debug +! + 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,*),& 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 -! Beware of changes in newer Marc versions -- these are from 2005r3 -! 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 +! Marc common blocks are in fixed format so they have to be reformated to free format (f90) +! Beware of changes in newer Marc versions -! creeps is needed for timinc (time increment) -! include 'creeps' - 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 + include "concom2007r1" ! concom is needed for inc, subinc, ncycle, lovl + include "creeps2007r1" ! creeps is needed for timinc (time increment) + + integer(pInt) computationMode,i + +! write(6,'(3(3(f10.3,x),/))') ffn1(:,1),ffn1(:,2),ffn1(:,3) if (inc == 0) then cycleCounter = 4 @@ -170,16 +160,22 @@ if (theCycle /= ncycle .or. theLovl /= lovl) then cycleCounter = cycleCounter+1 ! ping pong 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 if (cptim > theTime .or. theInc /= inc) then ! reached convergence lastIncConverged = .true. outdatedByNewInc = .true. + write (6,*) n(1),nn,'lastIncConverged + outdated' endif - if (mod(cycleCounter,2) /= 0) computationMode = 4 ! recycle - if (mod(cycleCounter,4) == 2) computationMode = 3 ! collect - if (mod(cycleCounter,4) == 0) computationMode = 2 ! compute + if (mod(cycleCounter,2) /= 0) computationMode = 4 ! recycle in odd cycles + if (mod(cycleCounter,4) == 2) computationMode = 3 ! collect in 2,6,10,... + if (mod(cycleCounter,4) == 0) computationMode = 2 ! compute in 0,4,8,... if (computationMode == 4 .and. ncycle == 0 .and. .not. lastIncConverged) & computationMode = 6 ! recycle but restore known good consistent tangent if (computationMode == 4 .and. lastIncConverged) then @@ -191,10 +187,6 @@ outdatedByNewInc = .false. endif - if (computationMode == 2 .and. outdatedFFN1) then - computationMode = 4 ! return odd results to force new vyvle - endif - theTime = cptim ! record current starting time theInc = inc ! record current increment number 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) 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))) + return END SUBROUTINE @@ -220,22 +213,22 @@ ! select a variable contour plotting (user subroutine). ! ! v variable -! s (idss) stress array +! s (idss) stress array ! sp stresses in preferred direction -! etot total strain (generalized) -! eplas total plastic strain -! ecreep total creep strain -! t current temperature +! etot total strain (generalized) +! eplas total plastic strain +! ecreep total creep strain +! t current temperature ! m element number ! nn integration point number ! layer layer number -! ndi (3) number of direct stress components -! nshear (3) number of shear stress components +! ndi (3) number of direct stress components +! nshear (3) number of shear stress components ! !******************************************************************** use prec, only: pReal,pInt use CPFEM, only: CPFEM_results, CPFEM_Nresults - use constitutive, only: constitutive_maxNresults + use constitutive, only: constitutive_maxSizePostResults use mesh, only: mesh_FEasCP implicit none ! @@ -244,10 +237,9 @@ integer(pInt) m, nn, layer, ndi, nshear, jpltcd ! ! assign result variable - v=CPFEM_results(mod(jpltcd-1_pInt, CPFEM_Nresults+constitutive_maxNresults)+1_pInt,& - (jpltcd-1_pInt)/(CPFEM_Nresults+constitutive_maxNresults)+1_pInt,& - nn, mesh_FEasCP('elem', m)) - + v = CPFEM_results(mod(jpltcd-1_pInt, CPFEM_Nresults+constitutive_maxSizePostResults)+1_pInt,& + (jpltcd-1_pInt)/(CPFEM_Nresults+constitutive_maxSizePostResults)+1_pInt,& + nn, mesh_FEasCP('elem', m)) return END SUBROUTINE ! @@ -300,4 +292,4 @@ ! end if ! end if ! return -! end +! end \ No newline at end of file diff --git a/trunk/mpie_cpfem_marc2008r1.f90 b/trunk/mpie_cpfem_marc2008r1.f90 index ffbe3dd91..d0c67025d 100644 --- a/trunk/mpie_cpfem_marc2008r1.f90 +++ b/trunk/mpie_cpfem_marc2008r1.f90 @@ -36,10 +36,10 @@ include "mesh.f90" ! uses prec, IO, math, FEsolving include "lattice.f90" ! uses prec, math 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 -! -! + + logical, parameter :: parallelExecution = .true. + SUBROUTINE hypela2(d,g,e,de,s,t,dt,ngens,n,nn,kcus,matus,ndi,& nshear,disp,dispt,coord,ffn,frotn,strechn,eigvn,ffn1,& frotn1,strechn1,eigvn1,ncrd,itel,ndeg,ndm,& @@ -128,6 +128,7 @@ use FEsolving use CPFEM, only: CPFEM_general use math, only: invnrmMandel + use debug ! implicit none @@ -146,9 +147,9 @@ ! 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 ! concom is needed for inc, subinc, ncycle, lovl - include "concom_f90" + include "concom2008r1" ! creeps is needed for timinc (time increment) - include "creeps_f90" + include "creeps2008r1" ! integer(pInt) computationMode,i ! @@ -159,6 +160,10 @@ if (theCycle /= ncycle .or. theLovl /= lovl) then cycleCounter = cycleCounter+1 ! ping pong 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 if (cptim > theTime .or. theInc /= inc) then ! reached convergence @@ -179,9 +184,6 @@ computationMode = 1 ! compute and age former results outdatedByNewInc = .false. endif - if (computationMode == 2 .and. outdatedFFN1) then - computationMode = 4 ! return odd results to force new vyvle - endif theTime = cptim ! record current starting time theInc = inc ! record current increment number @@ -208,17 +210,17 @@ ! select a variable contour plotting (user subroutine). ! ! v variable -! s (idss) stress array +! s (idss) stress array ! sp stresses in preferred direction -! etot total strain (generalized) -! eplas total plastic strain -! ecreep total creep strain -! t current temperature +! etot total strain (generalized) +! eplas total plastic strain +! ecreep total creep strain +! t current temperature ! m element number ! nn integration point number ! layer layer number -! ndi (3) number of direct stress components -! nshear (3) number of shear stress components +! ndi (3) number of direct stress components +! nshear (3) number of shear stress components ! !******************************************************************** use prec, only: pReal,pInt @@ -232,8 +234,8 @@ integer(pInt) m, nn, layer, ndi, nshear, jpltcd ! ! assign result variable - v=CPFEM_results(mod(jpltcd-1_pInt, CPFEM_Nresults+constitutive_maxNresults)+1_pInt,& - (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_maxSizePostResults)+1_pInt,& nn, mesh_FEasCP('elem', m)) return END SUBROUTINE diff --git a/trunk/prec.f90 b/trunk/prec.f90 index 91f0e9319..7babb3be0 100644 --- a/trunk/prec.f90 +++ b/trunk/prec.f90 @@ -9,6 +9,10 @@ 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 + type :: p_vec + real(pReal), dimension(:), pointer :: p + end type p_vec + ! *** Strain increment considered significant *** 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 integer(pInt), parameter :: nOuter = 20_pInt ! outer loop limit 20 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 :: abstol_Inner = 1.0e-8_pReal ! absolute tolerance in inner loop (Lp) !