From 21cad3213710f89fc4fe9a295f696a17ce85d570 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Mon, 12 Oct 2009 16:01:49 +0000 Subject: [PATCH] # new interface for Abaqus # IO has some additional functionality for Abaqus parsing # ping pong scheme in FE interface now similar (and more human understandable) in both versions # mesh has better splitting of different tasks, plus operation on database whenever possible # FEsolver as new global var to indicate FEM solver type # computation mode reshuffling: 6 is now Marc special case of recycling... --- code/CPFEM.f90 | 62 +- code/FEsolving.f90 | 7 +- code/IO.f90 | 249 +++- code/homogenization_RGC.f90 | 2260 +++++++++++++++++------------------ code/mesh.f90 | 1347 ++++++++++++++++----- code/mpie_cpfem_abaqus.f90 | 270 +++++ code/mpie_cpfem_marc.f90 | 159 ++- 7 files changed, 2775 insertions(+), 1579 deletions(-) create mode 100644 code/mpie_cpfem_abaqus.f90 diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index 37380abd5..b26ab28c0 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -79,8 +79,6 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt terminallyIll, & cycleCounter, & theInc, & - theCycle, & - theLovl, & theTime, & FEsolving_execElem, & FEsolving_execIP @@ -120,8 +118,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt materialpoint_Temperature, & materialpoint_stressAndItsTangent, & materialpoint_postResults - use IO, only: IO_init - use cpfem_interface + use IO, only: IO_init + use cpfem_interface implicit none @@ -133,12 +131,12 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt real(pReal), intent(in) :: dt ! time increment real(pReal), dimension (3,3), intent(in) :: ffn, & ! deformation gradient for t=t0 ffn1 ! deformation gradient for t=t1 - integer(pInt), intent(in) :: mode ! computation mode 1: regular computation with aged results + integer(pInt), intent(in) :: mode ! computation mode 1: regular computation plus aging of results ! 2: regular computation ! 3: collection of FEM data - ! 4: recycling of former results (MARC speciality) - ! 5: record tangent from former converged inc - ! 6: restore tangent from former converged inc + ! 4: backup tangent from former converged inc + ! 5: restore tangent from former converged inc + ! 6: recycling of former results (MARC speciality) !*** output variables ***! real(pReal), dimension(ngens), intent(out) :: cauchyStress ! stress vector in Mandel notation @@ -176,7 +174,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt call debug_init() call math_init() call FE_init() - call mesh_init() + call mesh_init(IP, element) ! pass on coordinates to alter calcMode of first ip call lattice_init() call material_init() call constitutive_init() @@ -190,12 +188,14 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt cp_en = mesh_FEasCP('elem',element) if (cp_en == 1 .and. IP == 1) then + !$OMP CRITICAL (write2out) write(6,*) write(6,*) '#####################################' - write(6,'(a10,1x,f8.4,1x,a10,1x,i4,1x,a10,1x,i3,1x,a10,1x,i2,x,a10,1x,i2)') & - 'theTime',theTime,'theInc',theInc,'theCycle',theCycle,'theLovl',theLovl,& - 'mode',mode + write(6,'(a10,x,f8.4,x,a10,x,i4,x,a10,x,i3,x,a16,x,i2,x,a16,x,i2)') & + 'theTime',theTime,'theInc',theInc,'cycleCounter',cycleCounter,'computationMode',mode write(6,*) '#####################################' + call flush (6) + !$OMP END CRITICAL (write2out) endif ! according to our "mode" we decide what to do @@ -213,7 +213,9 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt j = 1:mesh_maxNips, & k = 1:mesh_NcpElems ) & constitutive_state0(i,j,k)%p = constitutive_state(i,j,k)%p ! microstructure of crystallites + !$OMP CRITICAL (write2out) write(6,'(a10,/,4(3(e20.8,x),/))') 'aged state',constitutive_state(1,1,1)%p + !$OMP END CRITICAL (write2out) do k = 1,mesh_NcpElems do j = 1,mesh_maxNips if (homogenization_sizeState(j,k) > 0_pInt) & @@ -225,7 +227,9 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt ! deformation gradient outdated or any actual deformation gradient differs more than relevantStrain from the stored one if (terminallyIll .or. outdatedFFN1 .or. any(abs(ffn1 - materialpoint_F(:,:,IP,cp_en)) > relevantStrain)) then if (.not. terminallyIll .and. .not. outdatedFFN1) then + !$OMP CRITICAL (write2out) write(6,'(a11,x,i5,x,i2,x,a10,/,3(3(f10.6,x),/))') 'outdated at',cp_en,IP,'FFN1 now:',ffn1(:,1),ffn1(:,2),ffn1(:,3) + !$OMP END CRITICAL (write2out) outdatedFFN1 = .true. endif CPFEM_cs(1:ngens,IP,cp_en) = CPFEM_odd_stress @@ -234,7 +238,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt ! deformation gradient is not outdated else ! set flag for Jacobian update - updateJaco = (mod(cycleCounter-4,4_pInt*iJacoStiffness)==0) + updateJaco = mod(cycleCounter,iJacoStiffness) == 0 ! no parallel computation if (.not. parallelExecution) then @@ -278,8 +282,13 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt endif endif - ! --+>> COLLECTION OF FEM DATA AND RETURN OF ODD STRESS AND JACOBIAN <<+-- - case (3) + ! --+>> COLLECTION OF FEM INPUT WITH RETURNING OF ODD STRESS AND JACOBIAN <<+-- + case (3,4,5) + if (mode == 4) then + CPFEM_dcsde_knownGood = CPFEM_dcsde ! --+>> BACKUP JACOBIAN FROM FORMER CONVERGED INC + else if (mode == 5) then + CPFEM_dcsde = CPFEM_dcsde_knownGood ! --+>> RESTORE CONSISTENT JACOBIAN FROM FORMER CONVERGED INC + end if materialpoint_Temperature(IP,cp_en) = Temperature materialpoint_F0(:,:,IP,cp_en) = ffn materialpoint_F(:,:,IP,cp_en) = ffn1 @@ -288,29 +297,26 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt CPFEM_calc_done = .false. ! --+>> RECYCLING OF FORMER RESULTS (MARC SPECIALTY) <<+-- - case (4) + case (6) ! do nothing - ! --+>> RECORD JACOBIAN FROM FORMER CONVERGED INC <<+-- - case (5) - CPFEM_dcsde_knownGood = CPFEM_dcsde - - ! --+>> RESTORE CONSISTENT JACOBIAN FROM FORMER CONVERGED INC <<+-- - case (6) - CPFEM_dcsde = CPFEM_dcsde_knownGood - end select ! return the local stress and the jacobian from storage cauchyStress(1:ngens) = CPFEM_cs(1:ngens,IP,cp_en) jacobian(1:ngens,1:ngens) = CPFEM_dcsdE(1:ngens,1:ngens,IP,cp_en) - if (IP == 1 .and. cp_en == 1) write(6,'(a,/,6(6(f10.3,x)/))') 'jacobian/GPa at ip 1 el 1',jacobian/1e9 + if (IP == 1 .and. cp_en == 1) then + !$OMP CRITICAL (write2out) + write(6,'(a,/,6(6(f10.3,x)/))') 'jacobian/GPa at ip 1 el 1',jacobian/1e9 + call flush(6) + !$OMP END CRITICAL (write2out) + endif + ! return temperature if (theInc > 0_pInt) Temperature = materialpoint_Temperature(IP,cp_en) ! homogenized result except for potentially non-isothermal starting condition. return end subroutine - - END MODULE CPFEM - \ No newline at end of file +END MODULE CPFEM + diff --git a/code/FEsolving.f90 b/code/FEsolving.f90 index 32ce3a40b..86b995b37 100644 --- a/code/FEsolving.f90 +++ b/code/FEsolving.f90 @@ -6,12 +6,13 @@ use prec, only: pInt,pReal implicit none - integer(pInt) cycleCounter - integer(pInt) theInc,theCycle,theLovl - real(pReal) theTime + integer(pInt) cycleCounter, theInc + real(pReal) theTime, theDelta logical :: lastIncConverged = .false.,outdatedByNewInc = .false.,outdatedFFN1 = .false.,terminallyIll = .false. logical :: symmetricSolver = .false. logical :: parallelExecution = .true. + logical :: lastMode = .true., cutBack = .false. + logical, dimension(:,:), allocatable :: calcMode integer(pInt), dimension(:,:), allocatable :: FEsolving_execIP integer(pInt), dimension(2) :: FEsolving_execElem diff --git a/code/IO.f90 b/code/IO.f90 index c807a4795..df00c572e 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -27,11 +27,16 @@ !******************************************************************** ! output version number !******************************************************************** -subroutine IO_init - write(6,*) - write(6,*) '<<<+- IO init -+>>>' - write(6,*) '$Id$' - write(6,*) +subroutine IO_init () + + !$OMP CRITICAL (write2out) + write(6,*) + write(6,*) '<<<+- IO init -+>>>' + write(6,*) '$Id$' + write(6,*) + call flush(6) + !$OMP END CRITICAL (write2out) + return endsubroutine @@ -64,6 +69,7 @@ endsubroutine !******************************************************************** logical function IO_open_inputFile(unit) + use cpfem_interface use prec, only: pReal, pInt implicit none @@ -74,11 +80,18 @@ endsubroutine inquire(6, name=outName) ! determine outputfileName extPos = len_trim(outName)-2 - if(outName(extPos:extPos+2)=='out') then - ext='dat' ! MARC - else - ext='inp' ! ABAQUS - endif +! if(outName(extPos:extPos+2)=='out') then +! ext='dat' ! MARC +! else +! ext='inp' ! ABAQUS +! endif + select case (FEsolver) + case ('Marc') + ext='dat' + case ('Abaqus') + ext='inp' + end select + open(unit,status='old',err=100,file=outName(1:extPos-1)//ext) IO_open_inputFile = .true. return @@ -461,7 +474,7 @@ endfunction implicit none character(len=*), intent(in) :: line - character(len=*), parameter :: sep=achar(32)//achar(9)//achar(10)//achar(13) ! whitespaces + character(len=*), parameter :: sep=achar(44)//achar(32)//achar(9)//achar(10)//achar(13) ! comma and whitespaces integer(pInt), intent(in) :: N integer(pInt) part integer(pInt) IO_stringPos(1+N*2) @@ -703,91 +716,211 @@ endfunction endsubroutine + +!******************************************************************** +! extract value from key=value pair and check whether key matches +!******************************************************************** + pure function IO_extractValue (line,key) + use prec, only: pReal,pInt + implicit none + + character(len=*), intent(in) :: line,key + character(len=*), parameter :: sep = achar(61) ! '=' + integer(pInt) pos + character(len=300) IO_extractValue + + IO_extractValue = '' + + pos = scan(line,sep) + if (pos > 0 .and. line(:pos-1) == key(:pos-1)) & ! key matches expected key + IO_extractValue = line(pos+1:) ! extract value + + return + + endfunction + + !******************************************************************** -! count items in consecutive lines of ints concatenated by "c" -! as last char or range of values a "to" b +! count lines containig data up to next *keyword !******************************************************************** - function IO_countContinousIntValues (unit) + function IO_countDataLines (unit) use prec, only: pReal,pInt implicit none - integer(pInt) IO_countContinousIntValues,unit - integer(pInt), dimension(67) :: pos ! allow for 32 values excl "c" - character(len=300) line + integer(pInt) IO_countDataLines,unit + integer(pInt), parameter :: maxNchunks = 64 + integer(pInt), dimension(1+2*maxNchunks) :: pos + character(len=300) line,tmp - IO_countContinousIntValues = 0 + IO_countDataLines = 0 do read(unit,'(A300)',end=100) line - pos = IO_stringPos(line,33) - if (IO_lc(IO_stringValue(line,pos,2)) == 'to' ) then ! found range indicator - IO_countContinousIntValues = IO_countContinousIntValues+1+IO_intValue(line,pos,3)-IO_intValue(line,pos,1) - exit + pos = IO_stringPos(line,maxNchunks) + tmp = IO_lc(IO_stringValue(line,pos,1)) + if (tmp(1:1) == '*' ) then ! found keyword + exit else - IO_countContinousIntValues = IO_countContinousIntValues+pos(1)-1 - if ( IO_lc(IO_stringValue(line,pos,pos(1))) /= 'c' ) then ! line finished, read last value - IO_countContinousIntValues = IO_countContinousIntValues+1 - exit - endif + IO_countDataLines = IO_countDataLines + 1_pInt endif enddo +100 backspace(unit) + return + + endfunction + + +!******************************************************************** +! count items in consecutive lines +! Marc: ints concatenated by "c" as last char or range of values a "to" b +! Abaqus: triplet of start,stop,inc +!******************************************************************** + function IO_countContinousIntValues (unit) + + use cpfem_interface + use prec, only: pReal,pInt + implicit none + + integer(pInt) unit,i,j,l,count + integer(pInt) IO_countContinousIntValues + integer(pInt), parameter :: maxNchunks = 64 + integer(pInt), dimension(1+2*maxNchunks) :: pos + character(len=300) line + + IO_countContinousIntValues = 0_pInt + + select case (FEsolver) + case ('Marc') + + do + read(unit,'(A300)',end=100) line + pos = IO_stringPos(line,maxNchunks) + if (IO_lc(IO_stringValue(line,pos,2)) == 'to' ) then ! found range indicator + IO_countContinousIntValues = 1 + IO_intValue(line,pos,3) - IO_intValue(line,pos,1) + exit ! only one single range indicator allowed + else + IO_countContinousIntValues = IO_countContinousIntValues+pos(1)-1 ! add line's count when assuming 'c' + if ( IO_lc(IO_stringValue(line,pos,pos(1))) /= 'c' ) then ! line finished, read last value + IO_countContinousIntValues = IO_countContinousIntValues+1 + exit ! data ended + endif + endif + enddo + + case('Abaqus') + + count = IO_countDataLines(unit) + do l = 1,count + backspace(unit) + enddo + + do l = 1,count + read(unit,'(A300)',end=100) line + pos = IO_stringPos(line,maxNchunks) + IO_countContinousIntValues = IO_countContinousIntValues + 1 + & ! assuming range generation + (IO_intValue(line,pos,2)-IO_intValue(line,pos,1))/max(1,IO_intValue(line,pos,3)) + enddo + + endselect + 100 return endfunction -!********************************************************************* -! read consecutive lines of ints concatenated by "c" as last char -! or range of values a "to" b -!********************************************************************* + +!******************************************************************** +! return integer list corrsponding to items in consecutive lines +! Marc: ints concatenated by "c" as last char, range of a "to" b, or named set +! Abaqus: triplet of start,stop,inc or named set +!******************************************************************** function IO_continousIntValues (unit,maxN,lookupName,lookupMap,lookupMaxN) + use cpfem_interface use prec, only: pReal,pInt implicit none - integer(pInt) unit,maxN,i + integer(pInt) unit,maxN,i,j,l,count,first,last integer(pInt), dimension(1+maxN) :: IO_continousIntValues - integer(pInt), dimension(67) :: pos ! allow for 32 values excl "c" + integer(pInt), parameter :: maxNchunks = 64 + integer(pInt), dimension(1+2*maxNchunks) :: pos character(len=64), dimension(:) :: lookupName integer(pInt) :: lookupMaxN integer(pInt), dimension(:,:) :: lookupMap character(len=300) line - IO_continousIntValues = 0_pInt - do - read(unit,'(A300)',end=100) line - pos = IO_stringPos(line,33) - if (verify(IO_stringValue(line,pos,1),"0123456789") > 0) then ! a non-int, i.e. set name - do i = 1,lookupMaxN ! loop over known set names - if (IO_stringValue(line,pos,1) == lookupName(i)) then ! found matching name - IO_continousIntValues = lookupMap(:,i) ! return resp. entity list + IO_continousIntValues = 0 + + select case (FEsolver) + case ('Marc') + + do + read(unit,'(A300)',end=100) line + pos = IO_stringPos(line,maxNchunks) + if (verify(IO_stringValue(line,pos,1),"0123456789") > 0) then ! a non-int, i.e. set name + do i = 1,lookupMaxN ! loop over known set names + if (IO_stringValue(line,pos,1) == lookupName(i)) then ! found matching name + IO_continousIntValues = lookupMap(:,i) ! return resp. entity list + exit + endif + enddo exit + else if (IO_lc(IO_stringValue(line,pos,2)) == 'to' ) then ! found range indicator + do i = IO_intValue(line,pos,1),IO_intValue(line,pos,3) + IO_continousIntValues(1) = IO_continousIntValues(1) + 1 + IO_continousIntValues(1+IO_continousIntValues(1)) = i + enddo + exit + else + do i = 1,pos(1)-1 ! interpret up to second to last value + IO_continousIntValues(1) = IO_continousIntValues(1) + 1 + IO_continousIntValues(1+IO_continousIntValues(1)) = IO_intValue(line,pos,i) + enddo + if ( IO_lc(IO_stringValue(line,pos,pos(1))) /= 'c' ) then ! line finished, read last value + IO_continousIntValues(1) = IO_continousIntValues(1)+1 + IO_continousIntValues(1+IO_continousIntValues(1)) = IO_intValue(line,pos,pos(1)) + exit + endif endif enddo - exit - else if (IO_lc(IO_stringValue(line,pos,2)) == 'to' ) then ! found range indicator - do i = IO_intValue(line,pos,1),IO_intValue(line,pos,3) - IO_continousIntValues(1) = IO_continousIntValues(1)+1 - IO_continousIntValues(1+IO_continousIntValues(1)) = i + + case('Abaqus') + + count = IO_countDataLines(unit) + do l = 1,count + backspace(unit) enddo - exit - else - do i = 1,pos(1)-1 ! interpret up to second to last value - IO_continousIntValues(1) = IO_continousIntValues(1)+1 - IO_continousIntValues(1+IO_continousIntValues(1)) = IO_intValue(line,pos,i) + + do l = 1,count + read(unit,'(A300)',end=100) line + pos = IO_stringPos(line,maxNchunks) + if (verify(IO_stringValue(line,pos,1),"0123456789") > 0) then ! a non-int, i.e. set names follow on this line + do i = 1,pos(1) ! loop over set names in line + do j = 1,lookupMaxN ! look thru known set names + if (IO_stringValue(line,pos,i) == lookupName(j)) then ! found matching name + first = 2 + IO_continousIntValues(1) ! where to start appending data + last = first + lookupMap(1,j) - 1 ! up to where to append data + IO_continousIntValues(first:last) = lookupMap(2:1+lookupMap(1,j),j) ! add resp. entity list + IO_continousIntValues(1) = IO_continousIntValues(1) + lookupMap(1,j) ! count them + endif + enddo + enddo + else ! assuming range generation + do i = IO_intValue(line,pos,1),IO_intValue(line,pos,2),max(1,IO_intValue(line,pos,3)) + IO_continousIntValues(1) = IO_continousIntValues(1) + 1 + IO_continousIntValues(1+IO_continousIntValues(1)) = i + enddo + endif enddo - if ( IO_lc(IO_stringValue(line,pos,pos(1))) /= 'c' ) then ! line finished, read last value - IO_continousIntValues(1) = IO_continousIntValues(1)+1 - IO_continousIntValues(1+IO_continousIntValues(1)) = IO_intValue(line,pos,pos(1)) - exit - endif - endif - enddo + + endselect + 100 return endfunction + !******************************************************************** ! write error statements to standard out ! and terminate the Marc run with exit #9xxx diff --git a/code/homogenization_RGC.f90 b/code/homogenization_RGC.f90 index 8dbe339f8..d4c86e6b7 100644 --- a/code/homogenization_RGC.f90 +++ b/code/homogenization_RGC.f90 @@ -1,1130 +1,1130 @@ - -!***************************************************** -!* Module: HOMOGENIZATION_RGC * -!***************************************************** -!* contains: * -!***************************************************** - -! [rgc] -! type rgc -! Ngrains p x q x r (cluster) -! (output) Ngrains - -MODULE homogenization_RGC - -!*** Include other modules *** - use prec, only: pReal,pInt - implicit none - - character (len=*), parameter :: homogenization_RGC_label = 'rgc' - - integer(pInt), dimension(:), allocatable :: homogenization_RGC_sizeState, & - homogenization_RGC_sizePostResults - integer(pInt), dimension(:,:), allocatable :: homogenization_RGC_Ngrains - real(pReal), dimension(:,:), allocatable :: homogenization_RGC_xiAlpha, & - homogenization_RGC_ciAlpha - character(len=64), dimension(:,:), allocatable :: homogenization_RGC_output - -CONTAINS -!**************************************** -!* - homogenization_RGC_init -!* - homogenization_RGC_stateInit -!* - homogenization_RGC_deformationPartition -!* - homogenization_RGC_stateUpdate -!* - homogenization_RGC_averageStressAndItsTangent -!* - homogenization_RGC_postResults -!**************************************** - - -!************************************** -!* Module initialization * -!************************************** -subroutine homogenization_RGC_init(& - file & ! file pointer to material configuration - ) - - use prec, only: pInt, pReal - use math, only: math_Mandel3333to66, math_Voigt66to3333 - use mesh, only: mesh_maxNips,mesh_NcpElems - use IO - use material - integer(pInt), intent(in) :: file - integer(pInt), parameter :: maxNchunks = 4 - integer(pInt), dimension(1+2*maxNchunks) :: positions - integer(pInt) section, maxNinstance, i,j,k,l, output - character(len=64) tag - character(len=1024) line - - write(6,*) - write(6,'(a20,a20,a12)') '<<<+- homogenization',homogenization_RGC_label,' init -+>>>' - write(6,*) '$Id$' - write(6,*) - - maxNinstance = count(homogenization_type == homogenization_RGC_label) - if (maxNinstance == 0) return - - allocate(homogenization_RGC_sizeState(maxNinstance)); homogenization_RGC_sizeState = 0_pInt - allocate(homogenization_RGC_sizePostResults(maxNinstance)); homogenization_RGC_sizePostResults = 0_pInt - allocate(homogenization_RGC_Ngrains(3,maxNinstance)); homogenization_RGC_Ngrains = 0_pInt - allocate(homogenization_RGC_ciAlpha(3,maxNinstance)); homogenization_RGC_ciAlpha = 0.0_pReal - allocate(homogenization_RGC_xiAlpha(3,maxNinstance)); homogenization_RGC_xiAlpha = 0.0_pReal - allocate(homogenization_RGC_output(maxval(homogenization_Noutput),maxNinstance)); homogenization_RGC_output = '' - - rewind(file) - line = '' - section = 0 - - do while (IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization) ! 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. homogenization_type(section) == homogenization_RGC_label) then ! one of my sections - i = homogenization_typeInstance(section) ! which instance of my type is present homogenization - positions = IO_stringPos(line,maxNchunks) - tag = IO_lc(IO_stringValue(line,positions,1)) ! extract key - select case(tag) - case ('(output)') - output = output + 1 - homogenization_RGC_output(output,i) = IO_lc(IO_stringValue(line,positions,2)) - case ('clustersize') - homogenization_RGC_Ngrains(1,i) = IO_intValue(line,positions,2) - homogenization_RGC_Ngrains(2,i) = IO_intValue(line,positions,3) - homogenization_RGC_Ngrains(3,i) = IO_intValue(line,positions,4) - case ('grainsizeparameter') - homogenization_RGC_xiAlpha(1,i) = IO_floatValue(line,positions,2) - homogenization_RGC_xiAlpha(2,i) = IO_floatValue(line,positions,3) - homogenization_RGC_xiAlpha(3,i) = IO_floatValue(line,positions,4) - case ('overproportionality') - homogenization_RGC_ciAlpha(1,i) = IO_floatValue(line,positions,2) - homogenization_RGC_ciAlpha(2,i) = IO_floatValue(line,positions,3) - homogenization_RGC_ciAlpha(3,i) = IO_floatValue(line,positions,4) - end select - endif - enddo - -100 do i = 1,maxNinstance ! sanity checks - enddo - - do i = 1,maxNinstance - do j = 1,maxval(homogenization_Noutput) - select case(homogenization_RGC_output(j,i)) - case('constitutivework') - homogenization_RGC_sizePostResults(i) = & - homogenization_RGC_sizePostResults(i) + 1 - case('penaltyenergy') - homogenization_RGC_sizePostResults(i) = & - homogenization_RGC_sizePostResults(i) + 1 - case('magnitudemismatch') - homogenization_RGC_sizePostResults(i) = & - homogenization_RGC_sizePostResults(i) + 1 - end select - enddo - - homogenization_RGC_sizeState(i) & - = 3*(homogenization_RGC_Ngrains(1,i)-1)*homogenization_RGC_Ngrains(2,i)*homogenization_RGC_Ngrains(3,i) & - + 3*homogenization_RGC_Ngrains(1,i)*(homogenization_RGC_Ngrains(2,i)-1)*homogenization_RGC_Ngrains(3,i) & - + 3*homogenization_RGC_Ngrains(1,i)*homogenization_RGC_Ngrains(2,i)*(homogenization_RGC_Ngrains(3,i)-1) & - + 3_pInt ! (1) Average constitutive work, (2) Average penalty energy, (3) Overall mismatch - enddo - - return - -endsubroutine - - -!********************************************************************* -!* initial homogenization state * -!********************************************************************* -function homogenization_RGC_stateInit(myInstance) - use prec, only: pReal,pInt - implicit none - -!* Definition of variables - integer(pInt), intent(in) :: myInstance - real(pReal), dimension(homogenization_RGC_sizeState(myInstance)) :: homogenization_RGC_stateInit - -!* Open a debugging file -! open(1978,file='homogenization_RGC_debugging.out',status='unknown') - homogenization_RGC_stateInit = 0.0_pReal - - return - -endfunction - - -!******************************************************************** -! partition material point def grad onto constituents -!******************************************************************** -subroutine homogenization_RGC_partitionDeformation(& - F, & ! partioned def grad per grain -! - F0, & ! initial partioned def grad per grain - avgF, & ! my average def grad - state, & ! my state - ip, & ! my integration point - el & ! my element - ) - use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips - use material, only: homogenization_maxNgrains,homogenization_Ngrains,homogenization_typeInstance - use FEsolving, only: theInc,theCycle,theTime - - implicit none - -!* Definition of variables - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: F0 - real(pReal), dimension (3,3), intent(in) :: avgF - type(p_vec), intent(in) :: state - integer(pInt), intent(in) :: ip,el -! - real(pReal), dimension (3) :: aVect,nVect - integer(pInt), dimension (4) :: intFace - integer(pInt), dimension (3) :: iGrain3 - integer(pInt) homID, iGrain,iFace,i,j - logical RGCdebug -! - integer(pInt), parameter :: nFace = 6 - - RGCdebug = .false. - -!* Debugging the overall deformation gradient - if (RGCdebug) then - write(6,'(x,a,i3,a,i3,a)')'========== Increment: ',theInc,' Cycle: ',theCycle,' ==========' - write(6,'(x,a32)')'Overall deformation gradient: ' - do i = 1,3 - write(6,'(x,3(e14.8,x))')(avgF(i,j), j = 1,3) - enddo - write(6,*)' ' - call flush(6) - endif - -!* Compute the deformation gradient of individual grains due to relaxations - homID = homogenization_typeInstance(mesh_element(3,el)) - F = 0.0_pReal - do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) - call homogenization_RGC_grain1to3(iGrain3,iGrain,homID) - do iFace = 1,nFace - call homogenization_RGC_getInterface(intFace,iFace,iGrain3) - call homogenization_RGC_relaxationVector(aVect,intFace,state,homID) - call homogenization_RGC_interfaceNormal(nVect,intFace) - forall (i=1:3,j=1:3) & - F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations - enddo - F(:,:,iGrain) = F(:,:,iGrain) + avgF(:,:) ! relaxed deformation gradient -!* Debugging the grain deformation gradients - if (RGCdebug) then - write(6,'(x,a32,x,i3)')'Deformation gradient of grain: ',iGrain - do i = 1,3 - write(6,'(x,3(e14.8,x))')(F(i,j,iGrain), j = 1,3) - enddo - write(6,*)' ' - call flush(6) - endif - enddo - - return - -endsubroutine - -!******************************************************************** -! update the internal state of the homogenization scheme -! and tell whether "done" and "happy" with result -!******************************************************************** -function homogenization_RGC_updateState(& - state, & ! my state -! - P, & ! array of current grain stresses - F, & ! array of current grain deformation gradients - F0, & ! array of initial grain deformation gradients - avgF, & ! average deformation gradient - dPdF, & ! array of current grain stiffnesses - ip, & ! my integration point - el & ! my element - ) - - use prec, only: pReal,pInt,p_vec - use math, only: math_invert - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips - use material, only: homogenization_maxNgrains,homogenization_typeInstance,homogenization_Ngrains - use numerics, only: absTol_RGC,relTol_RGC,absMax_RGC,relMax_RGC,pPert_RGC - use FEsolving, only: theInc,theCycle,theTime - - implicit none - -!* Definition of variables - type(p_vec), intent(inout) :: state - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: P,F,F0 - real(pReal), dimension (3,3,3,3,homogenization_maxNgrains), intent(in) :: dPdF - real(pReal), dimension (3,3), intent(in) :: avgF - integer(pInt), intent(in) :: ip,el - logical, dimension(2) :: homogenization_RGC_updateState - integer(pInt), dimension (4) :: intFaceN,intFaceP,faceID - integer(pInt), dimension (3) :: nGDim,iGr3N,iGr3P,stresLoc - integer(pInt), dimension (2) :: residLoc - integer(pInt) homID,i1,i2,i3,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ival,ipert,iGrain - real(pReal), dimension (3,3,homogenization_maxNgrains) :: R,pF,pR - real(pReal), dimension (homogenization_maxNgrains) :: NN,pNN - real(pReal), dimension (3) :: normP,normN,mornP,mornN - real(pReal) residMax,stresMax,constitutiveWork,penaltyEnergy - logical error,RGCdebug,RGCdebugJacobi,RGCcheck -! - integer(pInt), parameter :: nFace = 6 - real(pInt), parameter :: max_drelax = 0.5_pReal -! - real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix - real(pReal), dimension(:), allocatable :: resid,relax,p_relax,p_resid,drelax - - RGCcheck = (el == 1 .and. ip == 1) - RGCdebug = .false. - RGCdebugJacobi = .false. - -!* Get the dimension of the cluster (grains and interfaces) - homID = homogenization_typeInstance(mesh_element(3,el)) - nGDim = homogenization_RGC_Ngrains(:,homID) - nIntFaceTot = (nGDim(1)-1)*nGDim(2)*nGDim(3) + nGDim(1)*(nGDim(2)-1)*nGDim(3) & - + nGDim(1)*nGDim(2)*(nGDim(3)-1) - -!* Allocate the size of the arrays/matrices depending on the size of the cluster - allocate(resid(3*nIntFaceTot)); resid = 0.0_pReal - allocate(tract(nIntFaceTot,3)); tract = 0.0_pReal - allocate(relax(3*nIntFaceTot)); relax = state%p(1:3*nIntFaceTot) - allocate(drelax(3*nIntFaceTot)); drelax = 0.0_pReal -!* Debugging the obtained state - if (RGCdebug) then - write(6,'(x,a30)')'Obtained state: ' - do i = 1,3*nIntFaceTot - write(6,'(x,2(e14.8,x))')state%p(i) - enddo - write(6,*)' ' - endif - -!* Stress-like penalty related to mismatch or incompatibility at interfaces - call homogenization_RGC_stressPenalty(R,NN,F,ip,el,homID) -!* Debugging the mismatch, stress and penalty of grains - if (RGCdebug) then - do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) - write(6,'(x,a30,x,i3,x,a4,x,e14.8)')'Mismatch magnitude of grain(',iGrain,') :',NN(iGrain) - write(6,*)' ' - write(6,'(x,a30,x,i3)')'Stress and penalty of grain: ',iGrain - do i = 1,3 - write(6,'(x,3(e14.8,x),x,3(e14.8,x))')(P(i,j,iGrain), j = 1,3),(R(i,j,iGrain), j = 1,3) - enddo - write(6,*)' ' - enddo - endif - -!* Compute the residual stress from the balance of traction at all (interior) interfaces - do iNum = 1,nIntFaceTot - call homogenization_RGC_interface1to4(faceID,iNum,homID) -!* Identify the left/bottom/back grain (-|N) - iGr3N = faceID(2:4) - call homogenization_RGC_grain3to1(iGrN,iGr3N,homID) - call homogenization_RGC_getInterface(intFaceN,2*faceID(1),iGr3N) - call homogenization_RGC_interfaceNormal(normN,intFaceN) ! get the interface normal -!* Identify the right/up/front grain (+|P) - iGr3P = iGr3N - iGr3P(faceID(1)) = iGr3N(faceID(1))+1 - call homogenization_RGC_grain3to1(iGrP,iGr3P,homID) - call homogenization_RGC_getInterface(intFaceP,2*faceID(1)-1,iGr3P) - call homogenization_RGC_interfaceNormal(normP,intFaceP) ! get the interface normal - do i = 1,3 ! compute the traction balance at the interface - do j = 1,3 - tract(iNum,i) = tract(iNum,i) + (P(i,j,iGrP) + R(i,j,iGrP))*normP(j) & - + (P(i,j,iGrN) + R(i,j,iGrN))*normN(j) - resid(i+3*(iNum-1)) = tract(iNum,i) ! map into 1D residual array - enddo - enddo -!* Debugging the residual stress - if (RGCdebug) then - write(6,'(x,a30,x,i3)')'Traction at interface: ',iNum - write(6,'(x,3(e14.8,x))')(tract(iNum,j), j = 1,3) - write(6,*)' ' - endif - enddo - -!* Convergence check for stress residual - stresMax = maxval(abs(P)) - stresLoc = maxloc(abs(P)) - residMax = maxval(abs(tract)) - residLoc = maxloc(abs(tract)) -!* Debugging the convergent criteria - if (RGCcheck) then - write(6,'(x,a)')' ' - write(6,'(x,a,x,i2,x,i4)')'RGC residual check ...',ip,el - write(6,'(x,a15,x,e14.8,x,a7,i3,x,a12,i2,i2)')'Max stress: ',stresMax, & - '@ grain',stresLoc(3),'in component',stresLoc(1),stresLoc(2) - write(6,'(x,a15,x,e14.8,x,a7,i3,x,a12,i2)')'Max residual: ',residMax, & - '@ iface',residLoc(1),'in direction',residLoc(2) - call flush(6) - endif - homogenization_RGC_updateState = .false. -!*** If convergence reached => done and happy - if (residMax < relTol_RGC*stresMax .or. residMax < absTol_RGC) then - homogenization_RGC_updateState = .true. - if (RGCcheck) then - write(6,'(x,a55)')'... done and happy' - write(6,*)' ' - call flush(6) - endif -! write(6,'(x,a,x,i3,x,a6,x,i3,x,a12)')'RGC_updateState: ip',ip,'| el',el,'converged :)' -!* Then compute/update the state for postResult, i.e., ... -!* ... the (bulk) constitutive work and penalty energy - constitutiveWork = state%p(3*nIntFaceTot+1) - penaltyEnergy = state%p(3*nIntFaceTot+2) - do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) - do i = 1,3 - do j = 1,3 - constitutiveWork = constitutiveWork + P(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain)) - penaltyEnergy = penaltyEnergy + R(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain)) - enddo - enddo - enddo - state%p(3*nIntFaceTot+1) = constitutiveWork - state%p(3*nIntFaceTot+2) = penaltyEnergy -!* ... the overall mismatch - state%p(3*nIntFaceTot+3) = sum(NN) - if (RGCdebug) then - write(6,'(x,a30,x,e14.8)')'Constitutive work: ',constitutiveWork - write(6,'(x,a30,x,e14.8)')'Penalty energy: ',penaltyEnergy - write(6,'(x,a30,x,e14.8)')'Magnitude mismatch: ',sum(NN) - write(6,*)'' - call flush(6) - endif - deallocate(tract,resid,relax,drelax) - return -!*** If residual blows-up => done but unhappy - elseif (residMax > relMax_RGC*stresMax .or. residMax > absMax_RGC) then -!* Try to restart when residual blows up exceeding maximum bound - homogenization_RGC_updateState = (/.true.,.false./) ! ... with direct cut-back - write(6,'(x,a,x,i3,x,a,x,i3,x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback' - write(6,'(x,a,x,e14.8,x,a,x,e14.8)')'due to high residual =',residMax,'for stress =',stresMax -! state%p(1:3*nIntFaceTot) = 0.0_pReal ! ... with full Taylor assumption -! write(6,'(x,a,x,i3,x,a,x,i3,x,a)')'RGC_updateState: ip',ip,'| el',el,'relaxation vectors reset to zero' - if (RGCcheck) then - write(6,'(x,a55)')'... broken' - write(6,*)' ' - call flush(6) - endif -! write(6,'(x,a,x,i3,x,a6,x,i3,x,a9)')'RGC_updateState: ip',ip,'| el',el,'broken :(' - deallocate(tract,resid,relax,drelax) - return -!*** Otherwise, proceed with computing the Jacobian and state update - else - if (RGCcheck) then - write(6,'(x,a55)')'... not yet done' - write(6,*)' ' - call flush(6) - endif - endif - -!* Construct the Jacobian matrix of the constitutive stress tangent from dPdF - allocate(smatrix(3*nIntFaceTot,3*nIntFaceTot)); smatrix = 0.0_pReal - do iNum = 1,nIntFaceTot - call homogenization_RGC_interface1to4(faceID,iNum,homID) -!* Identify the left/bottom/back grain (-|N) - iGr3N = faceID(2:4) - call homogenization_RGC_grain3to1(iGrN,iGr3N,homID) - call homogenization_RGC_getInterface(intFaceN,2*faceID(1),iGr3N) - call homogenization_RGC_interfaceNormal(normN,intFaceN) ! get the interface normal - do iFace = 1,nFace - call homogenization_RGC_getInterface(intFaceN,iFace,iGr3N) - call homogenization_RGC_interfaceNormal(mornN,intFaceN) ! get influencing interfaces normal - call homogenization_RGC_interface4to1(iMun,intFaceN,homID) - if (iMun .gt. 0) then ! get the tangent - forall(i=1:3,j=1:3,k=1:3,l=1:3) & - smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) + dPdF(i,k,j,l,iGrN)*normN(k)*mornN(l) - endif - enddo -!* Identify the right/up/front grain (+|P) - iGr3P = iGr3N - iGr3P(faceID(1)) = iGr3N(faceID(1))+1 - call homogenization_RGC_grain3to1(iGrP,iGr3P,homID) - call homogenization_RGC_getInterface(intFaceP,2*faceID(1)-1,iGr3P) - call homogenization_RGC_interfaceNormal(normP,intFaceP) ! get the interface normal - do iFace = 1,nFace - call homogenization_RGC_getInterface(intFaceP,iFace,iGr3P) - call homogenization_RGC_interfaceNormal(mornP,intFaceP) ! get influencing interfaces normal - call homogenization_RGC_interface4to1(iMun,intFaceP,homID) - if (iMun .gt. 0) then ! get the tangent - forall(i=1:3,j=1:3,k=1:3,l=1:3) & - smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) + dPdF(i,k,j,l,iGrP)*normP(k)*mornP(l) - endif - enddo - enddo -!* Debugging the global Jacobian matrix of stress tangent - if (RGCdebugJacobi) then - write(6,'(x,a30)')'Jacobian matrix of stress' - do i = 1,3*nIntFaceTot - write(6,'(x,100(e10.4,x))')(smatrix(i,j), j = 1,3*nIntFaceTot) - enddo - write(6,*)' ' - call flush(6) - endif - -!* Compute the Jacobian of the stress-like penalty (penalty tangent) using perturbation technique - allocate(pmatrix(3*nIntFaceTot,3*nIntFaceTot)); pmatrix = 0.0_pReal - allocate(p_relax(3*nIntFaceTot)); p_relax = 0.0_pReal - allocate(p_resid(3*nIntFaceTot)); p_resid = 0.0_pReal - do ipert = 1,3*nIntFaceTot - p_relax = relax - p_relax(ipert) = relax(ipert) + pPert_RGC ! perturb the relaxation vector - state%p(1:3*nIntFaceTot) = p_relax - call homogenization_RGC_grainDeformation(pF,F0,avgF,state,el) - call homogenization_RGC_stressPenalty(pR,pNN,pF,ip,el,homID) - p_resid = 0.0_pReal - do iNum = 1,nIntFaceTot - call homogenization_RGC_interface1to4(faceID,iNum,homID) -!* Identify the left/bottom/back grain (-|N) - iGr3N = faceID(2:4) - call homogenization_RGC_grain3to1(iGrN,iGr3N,homID) - call homogenization_RGC_getInterface(intFaceN,2*faceID(1),iGr3N) - call homogenization_RGC_interfaceNormal(normN,intFaceN) ! get the corresponding normal -!* Identify the right/up/front grain (+|P) - iGr3P = iGr3N - iGr3P(faceID(1)) = iGr3N(faceID(1))+1 - call homogenization_RGC_grain3to1(iGrP,iGr3P,homID) - call homogenization_RGC_getInterface(intFaceP,2*faceID(1)-1,iGr3P) - call homogenization_RGC_interfaceNormal(normP,intFaceP) ! get the corresponding normal -!* Compute the perturbed traction at interface - do i = 1,3 - do j = 1,3 - p_resid(i+3*(iNum-1)) = p_resid(i+3*(iNum-1)) + (pR(i,j,iGrP) - R(i,j,iGrP))*normP(j) & - + (pR(i,j,iGrN) - R(i,j,iGrN))*normN(j) - enddo - enddo - enddo - pmatrix(:,ipert) = p_resid/pPert_RGC - enddo -!* Debugging the global Jacobian matrix of penalty tangent - if (RGCdebugJacobi) then - write(6,'(x,a30)')'Jacobian matrix of penalty' - do i = 1,3*nIntFaceTot - write(6,'(x,100(e10.4,x))')(pmatrix(i,j), j = 1,3*nIntFaceTot) - enddo - write(6,*)' ' - call flush(6) - endif - -!* The overall Jacobian matrix (due to constitutive and penalty tangents) - allocate(jmatrix(3*nIntFaceTot,3*nIntFaceTot)); jmatrix = smatrix + pmatrix - if (RGCdebugJacobi) then - write(6,'(x,a30)')'Jacobian matrix (total)' - do i = 1,3*nIntFaceTot - write(6,'(x,100(e10.4,x))')(jmatrix(i,j), j = 1,3*nIntFaceTot) - enddo - write(6,*)' ' - call flush(6) - endif - allocate(jnverse(3*nIntFaceTot,3*nIntFaceTot)); jnverse = 0.0_pReal - call math_invert(3*nIntFaceTot,jmatrix,jnverse,ival,error) -!* Debugging the inverse Jacobian matrix - if (RGCdebugJacobi) then - write(6,'(x,a30)')'Jacobian inverse' - do i = 1,3*nIntFaceTot - write(6,'(x,100(e10.4,x))')(jnverse(i,j), j = 1,3*nIntFaceTot) - enddo - write(6,*)' ' - call flush(6) - endif - -!* Calculate the state update (i.e., new relaxation vectors) - do i = 1,3*nIntFaceTot - do j = 1,3*nIntFaceTot - drelax(i) = drelax(i) - jnverse(i,j)*resid(j) - enddo - enddo - relax = relax + drelax - state%p(1:3*nIntFaceTot) = relax - if (any(abs(drelax(:)) > max_drelax)) then -! state%p(1:3*nIntFaceTot) = 0.0_pReal -! write(6,'(x,a,x,i3,x,a,x,i3,x,a)')'RGC_updateState: ip',ip,'| el',el,'relaxation vectors reset to zero' - homogenization_RGC_updateState = (/.true.,.false./) - write(6,'(x,a,x,i3,x,a,x,i3,x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback' - write(6,'(x,a,x,e14.8)')'due to large relaxation change =',maxval(abs(drelax)) - call flush(6) - endif -!* Debugging the return state - if (RGCdebugJacobi) then - write(6,'(x,a30)')'Returned state: ' - do i = 1,3*nIntFaceTot - write(6,'(x,2(e14.8,x))')state%p(i) - enddo - write(6,*)' ' - call flush(6) - endif - - deallocate(tract,resid,jmatrix,jnverse,relax,drelax,pmatrix,smatrix,p_relax,p_resid) - return - -endfunction - -!******************************************************************** -! derive average stress and stiffness from constituent quantities -!******************************************************************** -subroutine homogenization_RGC_averageStressAndItsTangent(& - avgP, & ! average stress at material point - dAvgPdAvgF, & ! average stiffness at material point -! - P, & ! array of current grain stresses - dPdF, & ! array of current grain stiffnesses - ip, & ! my integration point - el & ! my element - ) - - use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips - use material, only: homogenization_maxNgrains, homogenization_Ngrains,homogenization_typeInstance - use math, only: math_Plain3333to99 - implicit none - -!* Definition of variables - real(pReal), dimension (3,3), intent(out) :: avgP - real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: P - real(pReal), dimension (3,3,3,3,homogenization_maxNgrains), intent(in) :: dPdF - real(pReal), dimension (9,9) :: dPdF99 - integer(pInt), intent(in) :: ip,el -! - logical homogenization_RGC_stateUpdate,RGCdebug - integer(pInt) homID, i, j, Ngrains, iGrain - - RGCdebug = .false. !(ip == 1 .and. el == 1) - - homID = homogenization_typeInstance(mesh_element(3,el)) -!* Debugging the grain tangent - if (RGCdebug) then - do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) - dPdF99 = math_Plain3333to99(dPdF(:,:,:,:,iGrain)) - write(6,'(x,a30,x,i3)')'Stress tangent of grain: ',iGrain - do i = 1,9 - write(6,'(x,(e14.8,x))') (dPdF99(i,j), j = 1,9) - enddo - write(6,*)' ' - enddo - call flush(6) - endif - Ngrains = homogenization_Ngrains(mesh_element(3,el)) - avgP = sum(P,3)/dble(Ngrains) - dAvgPdAvgF = sum(dPdF,5)/dble(Ngrains) - - return - -endsubroutine - -!******************************************************************** -! derive average stress and stiffness from constituent quantities -!******************************************************************** -function homogenization_RGC_averageTemperature(& - Temperature, & ! temperature - ip, & ! my integration point - el & ! my element - ) - - use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips - use material, only: homogenization_maxNgrains, homogenization_Ngrains - implicit none - -!* Definition of variables - real(pReal), dimension (homogenization_maxNgrains), intent(in) :: Temperature - integer(pInt), intent(in) :: ip,el - real(pReal) homogenization_RGC_averageTemperature - integer(pInt) homID, i, Ngrains - -! homID = homogenization_typeInstance(mesh_element(3,el)) ! <> - Ngrains = homogenization_Ngrains(mesh_element(3,el)) - homogenization_RGC_averageTemperature = sum(Temperature(1:Ngrains))/dble(Ngrains) - - return - -endfunction - -!******************************************************************** -! return array of homogenization results for post file inclusion -!******************************************************************** -pure function homogenization_RGC_postResults(& - state, & ! my state - ip, & ! my integration point - el & ! my element - ) - - use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_element - use material, only: homogenization_typeInstance,homogenization_Noutput,homogenization_Ngrains - implicit none - -!* Definition of variables - type(p_vec), intent(in) :: state - integer(pInt), intent(in) :: ip,el -! - integer(pInt) homID,o,c,nIntFaceTot - real(pReal), dimension(homogenization_RGC_sizePostResults(homogenization_typeInstance(mesh_element(3,el)))) :: & - homogenization_RGC_postResults - - homID = homogenization_typeInstance(mesh_element(3,el)) - nIntFaceTot = (homogenization_RGC_Ngrains(1,homID)-1)*homogenization_RGC_Ngrains(2,homID)*homogenization_RGC_Ngrains(3,homID) + & - homogenization_RGC_Ngrains(1,homID)*(homogenization_RGC_Ngrains(2,homID)-1)*homogenization_RGC_Ngrains(3,homID) + & - homogenization_RGC_Ngrains(1,homID)*homogenization_RGC_Ngrains(2,homID)*(homogenization_RGC_Ngrains(3,homID)-1) - - c = 0_pInt - homogenization_RGC_postResults = 0.0_pReal - do o = 1,homogenization_Noutput(mesh_element(3,el)) - select case(homogenization_RGC_output(o,homID)) - case('constitutivework') - homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+1) - c = c + 1 - case('penaltyenergy') - homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+2) - c = c + 1 - case('magnitudemismatch') - homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+3) - c = c + 1 - end select - enddo - - return - -endfunction - -!******************************************************************** -! subroutine to calculate stress-like penalty due to mismatch -!******************************************************************** -subroutine homogenization_RGC_stressPenalty(& - rPen, & ! stress-like penalty - nMis, & ! total amount of mismatch -! - fDef, & ! relaxation vectors - ip, & ! integration point - el, & ! element - homID & ! homogenization ID - ) - use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_element - use constitutive, only: constitutive_homogenizedC - use math, only: math_civita - use material, only: homogenization_maxNgrains,homogenization_Ngrains - use numerics, only: xSmoo_RGC - implicit none - -!* Definition of variables - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: rPen - real(pReal), dimension (homogenization_maxNgrains), intent(out) :: nMis - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: fDef - integer(pInt), intent(in) :: ip,el - integer(pInt), dimension (4) :: intFace - integer(pInt), dimension (3) :: iGrain3,iGNghb3,nGDim - real(pReal), dimension (3,3) :: gDef,nDef - real(pReal), dimension (3) :: nVect - integer(pInt) homID,iGrain,iGNghb,iFace,i,j,k,l - real(pReal) muGrain,muGNghb,nDefNorm -! - integer(pInt), parameter :: nFace = 6 - real(pReal), parameter :: nDefToler = 1.0e-10 - - nGDim = homogenization_RGC_Ngrains(:,homID) - - rPen = 0.0_pReal - nMis = 0.0_pReal - do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) - call homogenization_RGC_equivalentShearMod(muGrain,constitutive_homogenizedC(iGrain,ip,el)) - call homogenization_RGC_grain1to3(iGrain3,iGrain,homID) -!* Compute the mismatch tensor at all six interfaces - do iFace = 1,nFace - call homogenization_RGC_getInterface(intFace,iFace,iGrain3) - call homogenization_RGC_interfaceNormal(nVect,intFace) ! get the interface normal - iGNghb3 = iGrain3 ! identify the grain neighbor - iGNghb3(abs(intFace(1))) = iGNghb3(abs(intFace(1))) + int(dble(intFace(1))/dble(abs(intFace(1)))) -!* The grain periodicity along e1 - if (iGNghb3(1) < 1) iGNghb3(1) = nGDim(1) - if (iGNghb3(1) > nGDim(1)) iGNghb3(1) = 1 -!* The grain periodicity along e2 - if (iGNghb3(2) < 1) iGNghb3(2) = nGDim(2) - if (iGNghb3(2) > nGDim(2)) iGNghb3(2) = 1 -!* The grain periodicity along e3 - if (iGNghb3(3) < 1) iGNghb3(3) = nGDim(3) - if (iGNghb3(3) > nGDim(3)) iGNghb3(3) = 1 - call homogenization_RGC_grain3to1(iGNghb,iGNghb3,homID) ! get the grain neighbor - call homogenization_RGC_equivalentShearMod(muGNghb,constitutive_homogenizedC(iGNghb,ip,el)) - gDef = 0.5_pReal*(fDef(:,:,iGNghb) - fDef(:,:,iGrain)) ! difference in F with the neighbor - nDefNorm = 0.0_pReal - nDef = 0.0_pReal - do i = 1,3 - do j = 1,3 - do k = 1,3 - do l = 1,3 - nDef(i,j) = nDef(i,j) - nVect(k)*gDef(i,l)*math_civita(j,k,l) ! compute the interface mismatch tensor - enddo - enddo - nDefNorm = nDefNorm + nDef(i,j)*nDef(i,j) - enddo - enddo - nDefNorm = max(nDefToler,sqrt(nDefNorm)) ! zero mismatch approximation if too small -!* Debugging the mismatch tensor -! if (ip == 1 .and. el == 1) then -! write(6,'(x,a20,i2,x,a20,x,i3)')'Mismatch to face: ',intFace(1),'neighbor grain: ',iGNghb -! do i = 1,3 -! write(6,'(x,3(e10.4,x))')(nDef(i,j), j = 1,3) -! enddo -! write(6,'(x,a20,e10.4))')'with magnitude: ',nDefNorm -! endif -!* Compute the stress-like penalty from all six interfaces - do i = 1,3 - do j = 1,3 - do k = 1,3 - do l = 1,3 - rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain + muGNghb)/homogenization_RGC_xiAlpha(abs(intFace(1)),homID) & - *cosh(homogenization_RGC_ciAlpha(abs(intFace(1)),homID)*nDefNorm) & - *0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_civita(k,l,j) & - *tanh(nDefNorm/xSmoo_RGC) - enddo - enddo - enddo - enddo -!* Total amount of mismatch experienced by the grain (at all six interfaces) - nMis(iGrain) = nMis(iGrain) + nDefNorm - enddo -!* Debugging the stress-like penalty -! if (ip == 1 .and. el == 1) then -! write(6,'(x,a20,i2)')'Penalty of grain: ',iGrain -! do i = 1,3 -! write(6,'(x,3(e10.4,x))')(rPen(i,j,iGrain), j = 1,3) -! enddo -! endif - enddo - - return - -endsubroutine - -!******************************************************************** -! subroutine to compute the equivalent shear modulus from the elasticity tensor -!******************************************************************** -subroutine homogenization_RGC_equivalentShearMod(& - shearMod, & ! equivalent (isotropic) shear modulus -! - elasTens & ! elasticity tensor in Mandel notation - ) - - use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips - use material, only: homogenization_typeInstance - implicit none - -!* Definition of variables - real(pReal), dimension (6,6), intent(in) :: elasTens - real(pReal), intent(out) :: shearMod - real(pReal) cEquiv_11,cEquiv_12,cEquiv_44 - -!* Compute the equivalent shear modulus using Turterltaub and Suiker, JMPS (2005) - cEquiv_11 = (elasTens(1,1) + elasTens(2,2) + elasTens(3,3))/3.0_pReal - cEquiv_12 = (elasTens(1,2) + elasTens(2,3) + elasTens(3,1) + & - elasTens(1,3) + elasTens(2,1) + elasTens(3,2))/6.0_pReal - cEquiv_44 = (elasTens(4,4) + elasTens(5,5) + elasTens(6,6))/3.0_pReal - shearMod = 0.2_pReal*(cEquiv_11 - cEquiv_12) + 0.6_pReal*cEquiv_44 - - return - -endsubroutine - -!******************************************************************** -! subroutine to collect relaxation vectors of an interface -!******************************************************************** -subroutine homogenization_RGC_relaxationVector(& - aVect, & ! relaxation vector of the interface -! - intFace, & ! set of interface ID in 4D array (normal and position) - state, & ! set of global relaxation vectors - homID & ! homogenization ID - ) - - use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips - use material, only: homogenization_typeInstance - implicit none - -!* Definition of variables - real(pReal), dimension (3), intent(out) :: aVect - integer(pInt), dimension (4), intent(in) :: intFace - type(p_vec), intent(in) :: state - integer(pInt), dimension (3) :: nGDim - integer(pInt) iNum,homID - -!* Collect the interface relaxation vector from the global state array - aVect = 0.0_pReal - nGDim = homogenization_RGC_Ngrains(:,homID) - call homogenization_RGC_interface4to1(iNum,intFace,homID) ! Get the position in global state array - if (iNum .gt. 0_pInt) aVect = state%p((3*iNum-2):(3*iNum)) ! Collect the corresponding entries - - return - -endsubroutine - -!******************************************************************** -! subroutine to identify the normal of an interface -!******************************************************************** -subroutine homogenization_RGC_interfaceNormal(& - nVect, & ! interface normal -! - intFace & ! interface ID in 4D array (normal and position) - ) - - use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips - implicit none - -!* Definition of variables - real(pReal), dimension (3), intent(out) :: nVect - integer(pInt), dimension (4), intent(in) :: intFace - integer(pInt) nPos - -!* Get the normal of the interface, identified from the value of intFace(1) - nVect = 0.0_pReal - nPos = abs(intFace(1)) - nVect(nPos) = intFace(1)/abs(intFace(1)) - - return - -endsubroutine - -!******************************************************************** -! subroutine to collect six faces of a grain in 4D (normal and position) -!******************************************************************** -subroutine homogenization_RGC_getInterface(& - intFace, & ! interface ID in 4D (normal and position) -! - iFace, & ! number of faces of grain - iGrain3 & ! grain ID in 3D array - ) - use prec, only: pReal,pInt,p_vec - implicit none - -!* Definition of variables - integer(pInt), dimension (4), intent(out) :: intFace - integer(pInt), dimension (3), intent(in) :: iGrain3 - integer(pInt), intent(in) :: iFace - integer(pInt) iDir - -!* Direction of interface normal - iDir = (int(dble(iFace-1)/2.0_pReal)+1)*(-1_pInt)**iFace - intFace(1) = iDir - -!* Identify the interface position by the direction of its normal - intFace(2:4) = iGrain3(:) - if (iDir .eq. -1_pInt) intFace(2) = intFace(2)-1 - if (iDir .eq. -2_pInt) intFace(3) = intFace(3)-1 - if (iDir .eq. -3_pInt) intFace(4) = intFace(4)-1 - - return - -endsubroutine - -!******************************************************************** -! subroutine to map grain ID from in 1D (array) to in 3D (position) -!******************************************************************** -subroutine homogenization_RGC_grain1to3(& - grain3, & ! grain ID in 3D array (pos.x,pos.y,pos.z) -! - grain1, & ! grain ID in 1D array - homID & ! homogenization ID - ) - - use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips - implicit none - -!* Definition of variables - integer(pInt), dimension (3), intent(out) :: grain3 - integer(pInt), intent(in) :: grain1,homID - integer(pInt), dimension (3) :: nGDim - -!* Get the grain position - nGDim = homogenization_RGC_Ngrains(:,homID) - grain3(3) = int(dble(grain1-1)/dble(nGDim(1))/dble(nGDim(2)))+1 - grain3(2) = mod(int(dble(grain1-1)/dble(nGDim(1))),nGDim(2))+1 - grain3(1) = mod((grain1-1),nGDim(1))+1 - - return - -endsubroutine - -!******************************************************************** -! subroutine to map grain ID from in 3D (position) to in 1D (array) -!******************************************************************** -subroutine homogenization_RGC_grain3to1(& - grain1, & ! grain ID in 1D array -! - grain3, & ! grain ID in 3D array (pos.x,pos.y,pos.z) - homID & ! homogenization ID - ) - - use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips - implicit none - -!* Definition of variables - integer(pInt), dimension (3), intent(in) :: grain3 - integer(pInt), intent(out) :: grain1 - integer(pInt), dimension (3) :: nGDim - integer(pInt) homID - -!* Get the grain ID - nGDim = homogenization_RGC_Ngrains(:,homID) - grain1 = grain3(1) + nGDim(1)*(grain3(2)-1) + nGDim(1)*nGDim(2)*(grain3(3)-1) - - return - -endsubroutine - -!******************************************************************** -! subroutine to map interface ID from 4D (normal and position) into 1D (array) -!******************************************************************** -subroutine homogenization_RGC_interface4to1(& - iFace1D, & ! interface ID in 1D array -! - iFace4D, & ! interface ID in 4D array (n.dir,pos.x,pos.y,pos.z) - homID & ! homogenization ID - ) - - use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips - implicit none - -!* Definition of variables - integer(pInt), dimension (4), intent(in) :: iFace4D - integer(pInt), intent(out) :: iFace1D - integer(pInt), dimension (3) :: nGDim,nIntFace - integer(pInt) homID - - nGDim = homogenization_RGC_Ngrains(:,homID) -!* Get the number of interfaces, which ... - nIntFace(1) = (nGDim(1)-1)*nGDim(2)*nGDim(3) ! ... normal //e1 - nIntFace(2) = nGDim(1)*(nGDim(2)-1)*nGDim(3) ! ... normal //e2 - nIntFace(3) = nGDim(1)*nGDim(2)*(nGDim(3)-1) ! ... normal //e3 - -!* For interface with normal //e1 - if (abs(iFace4D(1)) == 1_pInt) then - iFace1D = iFace4D(3) + nGDim(2)*(iFace4D(4)-1) + nGDim(2)*nGDim(3)*(iFace4D(2)-1) - if ((iFace4D(2) == 0_pInt) .or. (iFace4D(2) == nGDim(1))) iFace1D = 0_pInt -!* For interface with normal //e2 - elseif (abs(iFace4D(1)) == 2_pInt) then - iFace1D = iFace4D(4) + nGDim(3)*(iFace4D(2)-1) + nGDim(3)*nGDim(1)*(iFace4D(3)-1) + nIntFace(1) - if ((iFace4D(3) == 0_pInt) .or. (iFace4D(3) == nGDim(2))) iFace1D = 0_pInt -!* For interface with normal //e3 - elseif (abs(iFace4D(1)) == 3_pInt) then - iFace1D = iFace4D(2) + nGDim(1)*(iFace4D(3)-1) + nGDim(1)*nGDim(2)*(iFace4D(4)-1) + nIntFace(1) + nIntFace(2) - if ((iFace4D(4) == 0_pInt) .or. (iFace4D(4) == nGDim(3))) iFace1D = 0_pInt - endif - - return - -endsubroutine - -!******************************************************************** -! subroutine to map interface ID from 1D (array) into 4D (normal and position) -!******************************************************************** -subroutine homogenization_RGC_interface1to4(& - iFace4D, & ! interface ID in 4D array (n.dir,pos.x,pos.y,pos.z) -! - iFace1D, & ! interface ID in 1D array - homID & ! homogenization ID - ) - - use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips - implicit none - -!* Definition of variables - integer(pInt), dimension (4), intent(out) :: iFace4D - integer(pInt), intent(in) :: iFace1D - integer(pInt), dimension (3) :: nGDim,nIntFace - integer(pInt) homID - - nGDim = homogenization_RGC_Ngrains(:,homID) -!* Get the number of interfaces, which ... - nIntFace(1) = (nGDim(1)-1)*nGDim(2)*nGDim(3) ! ... normal //e1 - nIntFace(2) = nGDim(1)*(nGDim(2)-1)*nGDim(3) ! ... normal //e2 - nIntFace(3) = nGDim(1)*nGDim(2)*(nGDim(3)-1) ! ... normal //e3 - -!* For interface ID between 1 and nIntFace(1) - if (iFace1D > 0 .and. iFace1D <= nIntFace(1)) then - iFace4D(1) = 1 - iFace4D(3) = mod((iFace1D-1),nGDim(2))+1 - iFace4D(4) = mod(int(dble(iFace1D-1)/dble(nGDim(2))),nGDim(3))+1 - iFace4D(2) = int(dble(iFace1D-1)/dble(nGDim(2))/dble(nGDim(3)))+1 -!* For interface ID between nIntFace(1) and nIntFace(1) + nIntFace(2) - elseif (iFace1D > nIntFace(1) .and. iFace1D <= (nIntFace(2) + nIntFace(1))) then - iFace4D(1) = 2 - iFace4D(4) = mod((iFace1D-nIntFace(1)-1),nGDim(3))+1 - iFace4D(2) = mod(int(dble(iFace1D-nIntFace(1)-1)/dble(nGDim(3))),nGDim(1))+1 - iFace4D(3) = int(dble(iFace1D-nIntFace(1)-1)/dble(nGDim(3))/dble(nGDim(1)))+1 -!* For interface ID between nIntFace(1) + nIntFace(2) and nIntFace(1) + nIntFace(2) + nIntFace(3) - elseif (iFace1D > nIntFace(2) + nIntFace(1) .and. iFace1D <= (nIntFace(3) + nIntFace(2) + nIntFace(1))) then - iFace4D(1) = 3 - iFace4D(2) = mod((iFace1D-nIntFace(2)-nIntFace(1)-1),nGDim(1))+1 - iFace4D(3) = mod(int(dble(iFace1D-nIntFace(2)-nIntFace(1)-1)/dble(nGDim(1))),nGDim(2))+1 - iFace4D(4) = int(dble(iFace1D-nIntFace(2)-nIntFace(1)-1)/dble(nGDim(1))/dble(nGDim(2)))+1 - endif - - return - -endsubroutine - -!******************************************************************** -! calculating the grain deformation gradient -!******************************************************************** -subroutine homogenization_RGC_grainDeformation(& - F, & ! partioned def grad per grain -! - F0, & ! initial partioned def grad per grain - avgF, & ! my average def grad - state, & ! my state - el & ! my element - ) - use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_element - use material, only: homogenization_maxNgrains,homogenization_Ngrains,homogenization_typeInstance - implicit none - -!* Definition of variables - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: F0 - real(pReal), dimension (3,3), intent(in) :: avgF - type(p_vec), intent(in) :: state - integer(pInt), intent(in) :: el -! - real(pReal), dimension (3) :: aVect,nVect - integer(pInt), dimension (4) :: intFace - integer(pInt), dimension (3) :: iGrain3 - integer(pInt) homID, iGrain,iFace,i,j -! - integer(pInt), parameter :: nFace = 6 - -!* Compute the deformation gradient of individual grains due to relaxations - homID = homogenization_typeInstance(mesh_element(3,el)) - F = 0.0_pReal - do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) - call homogenization_RGC_grain1to3(iGrain3,iGrain,homID) - do iFace = 1,nFace - call homogenization_RGC_getInterface(intFace,iFace,iGrain3) - call homogenization_RGC_relaxationVector(aVect,intFace,state,homID) - call homogenization_RGC_interfaceNormal(nVect,intFace) - forall (i=1:3,j=1:3) & - F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations - enddo - F(:,:,iGrain) = F(:,:,iGrain) + avgF(:,:) ! relaxed deformation gradient - enddo - - return - -endsubroutine - -END MODULE + +!***************************************************** +!* Module: HOMOGENIZATION_RGC * +!***************************************************** +!* contains: * +!***************************************************** + +! [rgc] +! type rgc +! Ngrains p x q x r (cluster) +! (output) Ngrains + +MODULE homogenization_RGC + +!*** Include other modules *** + use prec, only: pReal,pInt + implicit none + + character (len=*), parameter :: homogenization_RGC_label = 'rgc' + + integer(pInt), dimension(:), allocatable :: homogenization_RGC_sizeState, & + homogenization_RGC_sizePostResults + integer(pInt), dimension(:,:), allocatable :: homogenization_RGC_Ngrains + real(pReal), dimension(:,:), allocatable :: homogenization_RGC_xiAlpha, & + homogenization_RGC_ciAlpha + character(len=64), dimension(:,:), allocatable :: homogenization_RGC_output + +CONTAINS +!**************************************** +!* - homogenization_RGC_init +!* - homogenization_RGC_stateInit +!* - homogenization_RGC_deformationPartition +!* - homogenization_RGC_stateUpdate +!* - homogenization_RGC_averageStressAndItsTangent +!* - homogenization_RGC_postResults +!**************************************** + + +!************************************** +!* Module initialization * +!************************************** +subroutine homogenization_RGC_init(& + file & ! file pointer to material configuration + ) + + use prec, only: pInt, pReal + use math, only: math_Mandel3333to66, math_Voigt66to3333 + use mesh, only: mesh_maxNips,mesh_NcpElems + use IO + use material + integer(pInt), intent(in) :: file + integer(pInt), parameter :: maxNchunks = 4 + integer(pInt), dimension(1+2*maxNchunks) :: positions + integer(pInt) section, maxNinstance, i,j,k,l, output + character(len=64) tag + character(len=1024) line + + write(6,*) + write(6,'(a20,a20,a12)') '<<<+- homogenization',homogenization_RGC_label,' init -+>>>' + write(6,*) '$Id$' + write(6,*) + + maxNinstance = count(homogenization_type == homogenization_RGC_label) + if (maxNinstance == 0) return + + allocate(homogenization_RGC_sizeState(maxNinstance)); homogenization_RGC_sizeState = 0_pInt + allocate(homogenization_RGC_sizePostResults(maxNinstance)); homogenization_RGC_sizePostResults = 0_pInt + allocate(homogenization_RGC_Ngrains(3,maxNinstance)); homogenization_RGC_Ngrains = 0_pInt + allocate(homogenization_RGC_ciAlpha(3,maxNinstance)); homogenization_RGC_ciAlpha = 0.0_pReal + allocate(homogenization_RGC_xiAlpha(3,maxNinstance)); homogenization_RGC_xiAlpha = 0.0_pReal + allocate(homogenization_RGC_output(maxval(homogenization_Noutput),maxNinstance)); homogenization_RGC_output = '' + + rewind(file) + line = '' + section = 0 + + do while (IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization) ! 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. homogenization_type(section) == homogenization_RGC_label) then ! one of my sections + i = homogenization_typeInstance(section) ! which instance of my type is present homogenization + positions = IO_stringPos(line,maxNchunks) + tag = IO_lc(IO_stringValue(line,positions,1)) ! extract key + select case(tag) + case ('(output)') + output = output + 1 + homogenization_RGC_output(output,i) = IO_lc(IO_stringValue(line,positions,2)) + case ('clustersize') + homogenization_RGC_Ngrains(1,i) = IO_intValue(line,positions,2) + homogenization_RGC_Ngrains(2,i) = IO_intValue(line,positions,3) + homogenization_RGC_Ngrains(3,i) = IO_intValue(line,positions,4) + case ('grainsizeparameter') + homogenization_RGC_xiAlpha(1,i) = IO_floatValue(line,positions,2) + homogenization_RGC_xiAlpha(2,i) = IO_floatValue(line,positions,3) + homogenization_RGC_xiAlpha(3,i) = IO_floatValue(line,positions,4) + case ('overproportionality') + homogenization_RGC_ciAlpha(1,i) = IO_floatValue(line,positions,2) + homogenization_RGC_ciAlpha(2,i) = IO_floatValue(line,positions,3) + homogenization_RGC_ciAlpha(3,i) = IO_floatValue(line,positions,4) + end select + endif + enddo + +100 do i = 1,maxNinstance ! sanity checks + enddo + + do i = 1,maxNinstance + do j = 1,maxval(homogenization_Noutput) + select case(homogenization_RGC_output(j,i)) + case('constitutivework') + homogenization_RGC_sizePostResults(i) = & + homogenization_RGC_sizePostResults(i) + 1 + case('penaltyenergy') + homogenization_RGC_sizePostResults(i) = & + homogenization_RGC_sizePostResults(i) + 1 + case('magnitudemismatch') + homogenization_RGC_sizePostResults(i) = & + homogenization_RGC_sizePostResults(i) + 1 + end select + enddo + + homogenization_RGC_sizeState(i) & + = 3*(homogenization_RGC_Ngrains(1,i)-1)*homogenization_RGC_Ngrains(2,i)*homogenization_RGC_Ngrains(3,i) & + + 3*homogenization_RGC_Ngrains(1,i)*(homogenization_RGC_Ngrains(2,i)-1)*homogenization_RGC_Ngrains(3,i) & + + 3*homogenization_RGC_Ngrains(1,i)*homogenization_RGC_Ngrains(2,i)*(homogenization_RGC_Ngrains(3,i)-1) & + + 3_pInt ! (1) Average constitutive work, (2) Average penalty energy, (3) Overall mismatch + enddo + + return + +endsubroutine + + +!********************************************************************* +!* initial homogenization state * +!********************************************************************* +function homogenization_RGC_stateInit(myInstance) + use prec, only: pReal,pInt + implicit none + +!* Definition of variables + integer(pInt), intent(in) :: myInstance + real(pReal), dimension(homogenization_RGC_sizeState(myInstance)) :: homogenization_RGC_stateInit + +!* Open a debugging file +! open(1978,file='homogenization_RGC_debugging.out',status='unknown') + homogenization_RGC_stateInit = 0.0_pReal + + return + +endfunction + + +!******************************************************************** +! partition material point def grad onto constituents +!******************************************************************** +subroutine homogenization_RGC_partitionDeformation(& + F, & ! partioned def grad per grain +! + F0, & ! initial partioned def grad per grain + avgF, & ! my average def grad + state, & ! my state + ip, & ! my integration point + el & ! my element + ) + use prec, only: pReal,pInt,p_vec + use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips + use material, only: homogenization_maxNgrains,homogenization_Ngrains,homogenization_typeInstance + use FEsolving, only: theInc,cycleCounter,theTime + + implicit none + +!* Definition of variables + real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F + real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: F0 + real(pReal), dimension (3,3), intent(in) :: avgF + type(p_vec), intent(in) :: state + integer(pInt), intent(in) :: ip,el +! + real(pReal), dimension (3) :: aVect,nVect + integer(pInt), dimension (4) :: intFace + integer(pInt), dimension (3) :: iGrain3 + integer(pInt) homID, iGrain,iFace,i,j + logical RGCdebug +! + integer(pInt), parameter :: nFace = 6 + + RGCdebug = .false. + +!* Debugging the overall deformation gradient + if (RGCdebug) then + write(6,'(x,a,i3,a,i3,a)')'========== Increment: ',theInc,' Cycle: ',cycleCounter,' ==========' + write(6,'(x,a32)')'Overall deformation gradient: ' + do i = 1,3 + write(6,'(x,3(e14.8,x))')(avgF(i,j), j = 1,3) + enddo + write(6,*)' ' + call flush(6) + endif + +!* Compute the deformation gradient of individual grains due to relaxations + homID = homogenization_typeInstance(mesh_element(3,el)) + F = 0.0_pReal + do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) + call homogenization_RGC_grain1to3(iGrain3,iGrain,homID) + do iFace = 1,nFace + call homogenization_RGC_getInterface(intFace,iFace,iGrain3) + call homogenization_RGC_relaxationVector(aVect,intFace,state,homID) + call homogenization_RGC_interfaceNormal(nVect,intFace) + forall (i=1:3,j=1:3) & + F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations + enddo + F(:,:,iGrain) = F(:,:,iGrain) + avgF(:,:) ! relaxed deformation gradient +!* Debugging the grain deformation gradients + if (RGCdebug) then + write(6,'(x,a32,x,i3)')'Deformation gradient of grain: ',iGrain + do i = 1,3 + write(6,'(x,3(e14.8,x))')(F(i,j,iGrain), j = 1,3) + enddo + write(6,*)' ' + call flush(6) + endif + enddo + + return + +endsubroutine + +!******************************************************************** +! update the internal state of the homogenization scheme +! and tell whether "done" and "happy" with result +!******************************************************************** +function homogenization_RGC_updateState(& + state, & ! my state +! + P, & ! array of current grain stresses + F, & ! array of current grain deformation gradients + F0, & ! array of initial grain deformation gradients + avgF, & ! average deformation gradient + dPdF, & ! array of current grain stiffnesses + ip, & ! my integration point + el & ! my element + ) + + use prec, only: pReal,pInt,p_vec + use math, only: math_invert + use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips + use material, only: homogenization_maxNgrains,homogenization_typeInstance,homogenization_Ngrains + use numerics, only: absTol_RGC,relTol_RGC,absMax_RGC,relMax_RGC,pPert_RGC + use FEsolving, only: theInc,cycleCounter,theTime + + implicit none + +!* Definition of variables + type(p_vec), intent(inout) :: state + real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: P,F,F0 + real(pReal), dimension (3,3,3,3,homogenization_maxNgrains), intent(in) :: dPdF + real(pReal), dimension (3,3), intent(in) :: avgF + integer(pInt), intent(in) :: ip,el + logical, dimension(2) :: homogenization_RGC_updateState + integer(pInt), dimension (4) :: intFaceN,intFaceP,faceID + integer(pInt), dimension (3) :: nGDim,iGr3N,iGr3P,stresLoc + integer(pInt), dimension (2) :: residLoc + integer(pInt) homID,i1,i2,i3,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ival,ipert,iGrain + real(pReal), dimension (3,3,homogenization_maxNgrains) :: R,pF,pR + real(pReal), dimension (homogenization_maxNgrains) :: NN,pNN + real(pReal), dimension (3) :: normP,normN,mornP,mornN + real(pReal) residMax,stresMax,constitutiveWork,penaltyEnergy + logical error,RGCdebug,RGCdebugJacobi,RGCcheck +! + integer(pInt), parameter :: nFace = 6 + real(pInt), parameter :: max_drelax = 0.5_pReal +! + real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix + real(pReal), dimension(:), allocatable :: resid,relax,p_relax,p_resid,drelax + + RGCcheck = (el == 1 .and. ip == 1) + RGCdebug = .false. + RGCdebugJacobi = .false. + +!* Get the dimension of the cluster (grains and interfaces) + homID = homogenization_typeInstance(mesh_element(3,el)) + nGDim = homogenization_RGC_Ngrains(:,homID) + nIntFaceTot = (nGDim(1)-1)*nGDim(2)*nGDim(3) + nGDim(1)*(nGDim(2)-1)*nGDim(3) & + + nGDim(1)*nGDim(2)*(nGDim(3)-1) + +!* Allocate the size of the arrays/matrices depending on the size of the cluster + allocate(resid(3*nIntFaceTot)); resid = 0.0_pReal + allocate(tract(nIntFaceTot,3)); tract = 0.0_pReal + allocate(relax(3*nIntFaceTot)); relax = state%p(1:3*nIntFaceTot) + allocate(drelax(3*nIntFaceTot)); drelax = 0.0_pReal +!* Debugging the obtained state + if (RGCdebug) then + write(6,'(x,a30)')'Obtained state: ' + do i = 1,3*nIntFaceTot + write(6,'(x,2(e14.8,x))')state%p(i) + enddo + write(6,*)' ' + endif + +!* Stress-like penalty related to mismatch or incompatibility at interfaces + call homogenization_RGC_stressPenalty(R,NN,F,ip,el,homID) +!* Debugging the mismatch, stress and penalty of grains + if (RGCdebug) then + do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) + write(6,'(x,a30,x,i3,x,a4,x,e14.8)')'Mismatch magnitude of grain(',iGrain,') :',NN(iGrain) + write(6,*)' ' + write(6,'(x,a30,x,i3)')'Stress and penalty of grain: ',iGrain + do i = 1,3 + write(6,'(x,3(e14.8,x),x,3(e14.8,x))')(P(i,j,iGrain), j = 1,3),(R(i,j,iGrain), j = 1,3) + enddo + write(6,*)' ' + enddo + endif + +!* Compute the residual stress from the balance of traction at all (interior) interfaces + do iNum = 1,nIntFaceTot + call homogenization_RGC_interface1to4(faceID,iNum,homID) +!* Identify the left/bottom/back grain (-|N) + iGr3N = faceID(2:4) + call homogenization_RGC_grain3to1(iGrN,iGr3N,homID) + call homogenization_RGC_getInterface(intFaceN,2*faceID(1),iGr3N) + call homogenization_RGC_interfaceNormal(normN,intFaceN) ! get the interface normal +!* Identify the right/up/front grain (+|P) + iGr3P = iGr3N + iGr3P(faceID(1)) = iGr3N(faceID(1))+1 + call homogenization_RGC_grain3to1(iGrP,iGr3P,homID) + call homogenization_RGC_getInterface(intFaceP,2*faceID(1)-1,iGr3P) + call homogenization_RGC_interfaceNormal(normP,intFaceP) ! get the interface normal + do i = 1,3 ! compute the traction balance at the interface + do j = 1,3 + tract(iNum,i) = tract(iNum,i) + (P(i,j,iGrP) + R(i,j,iGrP))*normP(j) & + + (P(i,j,iGrN) + R(i,j,iGrN))*normN(j) + resid(i+3*(iNum-1)) = tract(iNum,i) ! map into 1D residual array + enddo + enddo +!* Debugging the residual stress + if (RGCdebug) then + write(6,'(x,a30,x,i3)')'Traction at interface: ',iNum + write(6,'(x,3(e14.8,x))')(tract(iNum,j), j = 1,3) + write(6,*)' ' + endif + enddo + +!* Convergence check for stress residual + stresMax = maxval(abs(P)) + stresLoc = maxloc(abs(P)) + residMax = maxval(abs(tract)) + residLoc = maxloc(abs(tract)) +!* Debugging the convergent criteria + if (RGCcheck) then + write(6,'(x,a)')' ' + write(6,'(x,a,x,i2,x,i4)')'RGC residual check ...',ip,el + write(6,'(x,a15,x,e14.8,x,a7,i3,x,a12,i2,i2)')'Max stress: ',stresMax, & + '@ grain',stresLoc(3),'in component',stresLoc(1),stresLoc(2) + write(6,'(x,a15,x,e14.8,x,a7,i3,x,a12,i2)')'Max residual: ',residMax, & + '@ iface',residLoc(1),'in direction',residLoc(2) + call flush(6) + endif + homogenization_RGC_updateState = .false. +!*** If convergence reached => done and happy + if (residMax < relTol_RGC*stresMax .or. residMax < absTol_RGC) then + homogenization_RGC_updateState = .true. + if (RGCcheck) then + write(6,'(x,a55)')'... done and happy' + write(6,*)' ' + call flush(6) + endif +! write(6,'(x,a,x,i3,x,a6,x,i3,x,a12)')'RGC_updateState: ip',ip,'| el',el,'converged :)' +!* Then compute/update the state for postResult, i.e., ... +!* ... the (bulk) constitutive work and penalty energy + constitutiveWork = state%p(3*nIntFaceTot+1) + penaltyEnergy = state%p(3*nIntFaceTot+2) + do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) + do i = 1,3 + do j = 1,3 + constitutiveWork = constitutiveWork + P(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain)) + penaltyEnergy = penaltyEnergy + R(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain)) + enddo + enddo + enddo + state%p(3*nIntFaceTot+1) = constitutiveWork + state%p(3*nIntFaceTot+2) = penaltyEnergy +!* ... the overall mismatch + state%p(3*nIntFaceTot+3) = sum(NN) + if (RGCdebug) then + write(6,'(x,a30,x,e14.8)')'Constitutive work: ',constitutiveWork + write(6,'(x,a30,x,e14.8)')'Penalty energy: ',penaltyEnergy + write(6,'(x,a30,x,e14.8)')'Magnitude mismatch: ',sum(NN) + write(6,*)'' + call flush(6) + endif + deallocate(tract,resid,relax,drelax) + return +!*** If residual blows-up => done but unhappy + elseif (residMax > relMax_RGC*stresMax .or. residMax > absMax_RGC) then +!* Try to restart when residual blows up exceeding maximum bound + homogenization_RGC_updateState = (/.true.,.false./) ! ... with direct cut-back + write(6,'(x,a,x,i3,x,a,x,i3,x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback' + write(6,'(x,a,x,e14.8,x,a,x,e14.8)')'due to high residual =',residMax,'for stress =',stresMax +! state%p(1:3*nIntFaceTot) = 0.0_pReal ! ... with full Taylor assumption +! write(6,'(x,a,x,i3,x,a,x,i3,x,a)')'RGC_updateState: ip',ip,'| el',el,'relaxation vectors reset to zero' + if (RGCcheck) then + write(6,'(x,a55)')'... broken' + write(6,*)' ' + call flush(6) + endif +! write(6,'(x,a,x,i3,x,a6,x,i3,x,a9)')'RGC_updateState: ip',ip,'| el',el,'broken :(' + deallocate(tract,resid,relax,drelax) + return +!*** Otherwise, proceed with computing the Jacobian and state update + else + if (RGCcheck) then + write(6,'(x,a55)')'... not yet done' + write(6,*)' ' + call flush(6) + endif + endif + +!* Construct the Jacobian matrix of the constitutive stress tangent from dPdF + allocate(smatrix(3*nIntFaceTot,3*nIntFaceTot)); smatrix = 0.0_pReal + do iNum = 1,nIntFaceTot + call homogenization_RGC_interface1to4(faceID,iNum,homID) +!* Identify the left/bottom/back grain (-|N) + iGr3N = faceID(2:4) + call homogenization_RGC_grain3to1(iGrN,iGr3N,homID) + call homogenization_RGC_getInterface(intFaceN,2*faceID(1),iGr3N) + call homogenization_RGC_interfaceNormal(normN,intFaceN) ! get the interface normal + do iFace = 1,nFace + call homogenization_RGC_getInterface(intFaceN,iFace,iGr3N) + call homogenization_RGC_interfaceNormal(mornN,intFaceN) ! get influencing interfaces normal + call homogenization_RGC_interface4to1(iMun,intFaceN,homID) + if (iMun .gt. 0) then ! get the tangent + forall(i=1:3,j=1:3,k=1:3,l=1:3) & + smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) + dPdF(i,k,j,l,iGrN)*normN(k)*mornN(l) + endif + enddo +!* Identify the right/up/front grain (+|P) + iGr3P = iGr3N + iGr3P(faceID(1)) = iGr3N(faceID(1))+1 + call homogenization_RGC_grain3to1(iGrP,iGr3P,homID) + call homogenization_RGC_getInterface(intFaceP,2*faceID(1)-1,iGr3P) + call homogenization_RGC_interfaceNormal(normP,intFaceP) ! get the interface normal + do iFace = 1,nFace + call homogenization_RGC_getInterface(intFaceP,iFace,iGr3P) + call homogenization_RGC_interfaceNormal(mornP,intFaceP) ! get influencing interfaces normal + call homogenization_RGC_interface4to1(iMun,intFaceP,homID) + if (iMun .gt. 0) then ! get the tangent + forall(i=1:3,j=1:3,k=1:3,l=1:3) & + smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) + dPdF(i,k,j,l,iGrP)*normP(k)*mornP(l) + endif + enddo + enddo +!* Debugging the global Jacobian matrix of stress tangent + if (RGCdebugJacobi) then + write(6,'(x,a30)')'Jacobian matrix of stress' + do i = 1,3*nIntFaceTot + write(6,'(x,100(e10.4,x))')(smatrix(i,j), j = 1,3*nIntFaceTot) + enddo + write(6,*)' ' + call flush(6) + endif + +!* Compute the Jacobian of the stress-like penalty (penalty tangent) using perturbation technique + allocate(pmatrix(3*nIntFaceTot,3*nIntFaceTot)); pmatrix = 0.0_pReal + allocate(p_relax(3*nIntFaceTot)); p_relax = 0.0_pReal + allocate(p_resid(3*nIntFaceTot)); p_resid = 0.0_pReal + do ipert = 1,3*nIntFaceTot + p_relax = relax + p_relax(ipert) = relax(ipert) + pPert_RGC ! perturb the relaxation vector + state%p(1:3*nIntFaceTot) = p_relax + call homogenization_RGC_grainDeformation(pF,F0,avgF,state,el) + call homogenization_RGC_stressPenalty(pR,pNN,pF,ip,el,homID) + p_resid = 0.0_pReal + do iNum = 1,nIntFaceTot + call homogenization_RGC_interface1to4(faceID,iNum,homID) +!* Identify the left/bottom/back grain (-|N) + iGr3N = faceID(2:4) + call homogenization_RGC_grain3to1(iGrN,iGr3N,homID) + call homogenization_RGC_getInterface(intFaceN,2*faceID(1),iGr3N) + call homogenization_RGC_interfaceNormal(normN,intFaceN) ! get the corresponding normal +!* Identify the right/up/front grain (+|P) + iGr3P = iGr3N + iGr3P(faceID(1)) = iGr3N(faceID(1))+1 + call homogenization_RGC_grain3to1(iGrP,iGr3P,homID) + call homogenization_RGC_getInterface(intFaceP,2*faceID(1)-1,iGr3P) + call homogenization_RGC_interfaceNormal(normP,intFaceP) ! get the corresponding normal +!* Compute the perturbed traction at interface + do i = 1,3 + do j = 1,3 + p_resid(i+3*(iNum-1)) = p_resid(i+3*(iNum-1)) + (pR(i,j,iGrP) - R(i,j,iGrP))*normP(j) & + + (pR(i,j,iGrN) - R(i,j,iGrN))*normN(j) + enddo + enddo + enddo + pmatrix(:,ipert) = p_resid/pPert_RGC + enddo +!* Debugging the global Jacobian matrix of penalty tangent + if (RGCdebugJacobi) then + write(6,'(x,a30)')'Jacobian matrix of penalty' + do i = 1,3*nIntFaceTot + write(6,'(x,100(e10.4,x))')(pmatrix(i,j), j = 1,3*nIntFaceTot) + enddo + write(6,*)' ' + call flush(6) + endif + +!* The overall Jacobian matrix (due to constitutive and penalty tangents) + allocate(jmatrix(3*nIntFaceTot,3*nIntFaceTot)); jmatrix = smatrix + pmatrix + if (RGCdebugJacobi) then + write(6,'(x,a30)')'Jacobian matrix (total)' + do i = 1,3*nIntFaceTot + write(6,'(x,100(e10.4,x))')(jmatrix(i,j), j = 1,3*nIntFaceTot) + enddo + write(6,*)' ' + call flush(6) + endif + allocate(jnverse(3*nIntFaceTot,3*nIntFaceTot)); jnverse = 0.0_pReal + call math_invert(3*nIntFaceTot,jmatrix,jnverse,ival,error) +!* Debugging the inverse Jacobian matrix + if (RGCdebugJacobi) then + write(6,'(x,a30)')'Jacobian inverse' + do i = 1,3*nIntFaceTot + write(6,'(x,100(e10.4,x))')(jnverse(i,j), j = 1,3*nIntFaceTot) + enddo + write(6,*)' ' + call flush(6) + endif + +!* Calculate the state update (i.e., new relaxation vectors) + do i = 1,3*nIntFaceTot + do j = 1,3*nIntFaceTot + drelax(i) = drelax(i) - jnverse(i,j)*resid(j) + enddo + enddo + relax = relax + drelax + state%p(1:3*nIntFaceTot) = relax + if (any(abs(drelax(:)) > max_drelax)) then +! state%p(1:3*nIntFaceTot) = 0.0_pReal +! write(6,'(x,a,x,i3,x,a,x,i3,x,a)')'RGC_updateState: ip',ip,'| el',el,'relaxation vectors reset to zero' + homogenization_RGC_updateState = (/.true.,.false./) + write(6,'(x,a,x,i3,x,a,x,i3,x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback' + write(6,'(x,a,x,e14.8)')'due to large relaxation change =',maxval(abs(drelax)) + call flush(6) + endif +!* Debugging the return state + if (RGCdebugJacobi) then + write(6,'(x,a30)')'Returned state: ' + do i = 1,3*nIntFaceTot + write(6,'(x,2(e14.8,x))')state%p(i) + enddo + write(6,*)' ' + call flush(6) + endif + + deallocate(tract,resid,jmatrix,jnverse,relax,drelax,pmatrix,smatrix,p_relax,p_resid) + return + +endfunction + +!******************************************************************** +! derive average stress and stiffness from constituent quantities +!******************************************************************** +subroutine homogenization_RGC_averageStressAndItsTangent(& + avgP, & ! average stress at material point + dAvgPdAvgF, & ! average stiffness at material point +! + P, & ! array of current grain stresses + dPdF, & ! array of current grain stiffnesses + ip, & ! my integration point + el & ! my element + ) + + use prec, only: pReal,pInt,p_vec + use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips + use material, only: homogenization_maxNgrains, homogenization_Ngrains,homogenization_typeInstance + use math, only: math_Plain3333to99 + implicit none + +!* Definition of variables + real(pReal), dimension (3,3), intent(out) :: avgP + real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF + real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: P + real(pReal), dimension (3,3,3,3,homogenization_maxNgrains), intent(in) :: dPdF + real(pReal), dimension (9,9) :: dPdF99 + integer(pInt), intent(in) :: ip,el +! + logical homogenization_RGC_stateUpdate,RGCdebug + integer(pInt) homID, i, j, Ngrains, iGrain + + RGCdebug = .false. !(ip == 1 .and. el == 1) + + homID = homogenization_typeInstance(mesh_element(3,el)) +!* Debugging the grain tangent + if (RGCdebug) then + do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) + dPdF99 = math_Plain3333to99(dPdF(:,:,:,:,iGrain)) + write(6,'(x,a30,x,i3)')'Stress tangent of grain: ',iGrain + do i = 1,9 + write(6,'(x,(e14.8,x))') (dPdF99(i,j), j = 1,9) + enddo + write(6,*)' ' + enddo + call flush(6) + endif + Ngrains = homogenization_Ngrains(mesh_element(3,el)) + avgP = sum(P,3)/dble(Ngrains) + dAvgPdAvgF = sum(dPdF,5)/dble(Ngrains) + + return + +endsubroutine + +!******************************************************************** +! derive average stress and stiffness from constituent quantities +!******************************************************************** +function homogenization_RGC_averageTemperature(& + Temperature, & ! temperature + ip, & ! my integration point + el & ! my element + ) + + use prec, only: pReal,pInt,p_vec + use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips + use material, only: homogenization_maxNgrains, homogenization_Ngrains + implicit none + +!* Definition of variables + real(pReal), dimension (homogenization_maxNgrains), intent(in) :: Temperature + integer(pInt), intent(in) :: ip,el + real(pReal) homogenization_RGC_averageTemperature + integer(pInt) homID, i, Ngrains + +! homID = homogenization_typeInstance(mesh_element(3,el)) ! <> + Ngrains = homogenization_Ngrains(mesh_element(3,el)) + homogenization_RGC_averageTemperature = sum(Temperature(1:Ngrains))/dble(Ngrains) + + return + +endfunction + +!******************************************************************** +! return array of homogenization results for post file inclusion +!******************************************************************** +pure function homogenization_RGC_postResults(& + state, & ! my state + ip, & ! my integration point + el & ! my element + ) + + use prec, only: pReal,pInt,p_vec + use mesh, only: mesh_element + use material, only: homogenization_typeInstance,homogenization_Noutput,homogenization_Ngrains + implicit none + +!* Definition of variables + type(p_vec), intent(in) :: state + integer(pInt), intent(in) :: ip,el +! + integer(pInt) homID,o,c,nIntFaceTot + real(pReal), dimension(homogenization_RGC_sizePostResults(homogenization_typeInstance(mesh_element(3,el)))) :: & + homogenization_RGC_postResults + + homID = homogenization_typeInstance(mesh_element(3,el)) + nIntFaceTot = (homogenization_RGC_Ngrains(1,homID)-1)*homogenization_RGC_Ngrains(2,homID)*homogenization_RGC_Ngrains(3,homID) + & + homogenization_RGC_Ngrains(1,homID)*(homogenization_RGC_Ngrains(2,homID)-1)*homogenization_RGC_Ngrains(3,homID) + & + homogenization_RGC_Ngrains(1,homID)*homogenization_RGC_Ngrains(2,homID)*(homogenization_RGC_Ngrains(3,homID)-1) + + c = 0_pInt + homogenization_RGC_postResults = 0.0_pReal + do o = 1,homogenization_Noutput(mesh_element(3,el)) + select case(homogenization_RGC_output(o,homID)) + case('constitutivework') + homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+1) + c = c + 1 + case('penaltyenergy') + homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+2) + c = c + 1 + case('magnitudemismatch') + homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+3) + c = c + 1 + end select + enddo + + return + +endfunction + +!******************************************************************** +! subroutine to calculate stress-like penalty due to mismatch +!******************************************************************** +subroutine homogenization_RGC_stressPenalty(& + rPen, & ! stress-like penalty + nMis, & ! total amount of mismatch +! + fDef, & ! relaxation vectors + ip, & ! integration point + el, & ! element + homID & ! homogenization ID + ) + use prec, only: pReal,pInt,p_vec + use mesh, only: mesh_element + use constitutive, only: constitutive_homogenizedC + use math, only: math_civita + use material, only: homogenization_maxNgrains,homogenization_Ngrains + use numerics, only: xSmoo_RGC + implicit none + +!* Definition of variables + real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: rPen + real(pReal), dimension (homogenization_maxNgrains), intent(out) :: nMis + real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: fDef + integer(pInt), intent(in) :: ip,el + integer(pInt), dimension (4) :: intFace + integer(pInt), dimension (3) :: iGrain3,iGNghb3,nGDim + real(pReal), dimension (3,3) :: gDef,nDef + real(pReal), dimension (3) :: nVect + integer(pInt) homID,iGrain,iGNghb,iFace,i,j,k,l + real(pReal) muGrain,muGNghb,nDefNorm +! + integer(pInt), parameter :: nFace = 6 + real(pReal), parameter :: nDefToler = 1.0e-10 + + nGDim = homogenization_RGC_Ngrains(:,homID) + + rPen = 0.0_pReal + nMis = 0.0_pReal + do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) + call homogenization_RGC_equivalentShearMod(muGrain,constitutive_homogenizedC(iGrain,ip,el)) + call homogenization_RGC_grain1to3(iGrain3,iGrain,homID) +!* Compute the mismatch tensor at all six interfaces + do iFace = 1,nFace + call homogenization_RGC_getInterface(intFace,iFace,iGrain3) + call homogenization_RGC_interfaceNormal(nVect,intFace) ! get the interface normal + iGNghb3 = iGrain3 ! identify the grain neighbor + iGNghb3(abs(intFace(1))) = iGNghb3(abs(intFace(1))) + int(dble(intFace(1))/dble(abs(intFace(1)))) +!* The grain periodicity along e1 + if (iGNghb3(1) < 1) iGNghb3(1) = nGDim(1) + if (iGNghb3(1) > nGDim(1)) iGNghb3(1) = 1 +!* The grain periodicity along e2 + if (iGNghb3(2) < 1) iGNghb3(2) = nGDim(2) + if (iGNghb3(2) > nGDim(2)) iGNghb3(2) = 1 +!* The grain periodicity along e3 + if (iGNghb3(3) < 1) iGNghb3(3) = nGDim(3) + if (iGNghb3(3) > nGDim(3)) iGNghb3(3) = 1 + call homogenization_RGC_grain3to1(iGNghb,iGNghb3,homID) ! get the grain neighbor + call homogenization_RGC_equivalentShearMod(muGNghb,constitutive_homogenizedC(iGNghb,ip,el)) + gDef = 0.5_pReal*(fDef(:,:,iGNghb) - fDef(:,:,iGrain)) ! difference in F with the neighbor + nDefNorm = 0.0_pReal + nDef = 0.0_pReal + do i = 1,3 + do j = 1,3 + do k = 1,3 + do l = 1,3 + nDef(i,j) = nDef(i,j) - nVect(k)*gDef(i,l)*math_civita(j,k,l) ! compute the interface mismatch tensor + enddo + enddo + nDefNorm = nDefNorm + nDef(i,j)*nDef(i,j) + enddo + enddo + nDefNorm = max(nDefToler,sqrt(nDefNorm)) ! zero mismatch approximation if too small +!* Debugging the mismatch tensor +! if (ip == 1 .and. el == 1) then +! write(6,'(x,a20,i2,x,a20,x,i3)')'Mismatch to face: ',intFace(1),'neighbor grain: ',iGNghb +! do i = 1,3 +! write(6,'(x,3(e10.4,x))')(nDef(i,j), j = 1,3) +! enddo +! write(6,'(x,a20,e10.4))')'with magnitude: ',nDefNorm +! endif +!* Compute the stress-like penalty from all six interfaces + do i = 1,3 + do j = 1,3 + do k = 1,3 + do l = 1,3 + rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain + muGNghb)/homogenization_RGC_xiAlpha(abs(intFace(1)),homID) & + *cosh(homogenization_RGC_ciAlpha(abs(intFace(1)),homID)*nDefNorm) & + *0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_civita(k,l,j) & + *tanh(nDefNorm/xSmoo_RGC) + enddo + enddo + enddo + enddo +!* Total amount of mismatch experienced by the grain (at all six interfaces) + nMis(iGrain) = nMis(iGrain) + nDefNorm + enddo +!* Debugging the stress-like penalty +! if (ip == 1 .and. el == 1) then +! write(6,'(x,a20,i2)')'Penalty of grain: ',iGrain +! do i = 1,3 +! write(6,'(x,3(e10.4,x))')(rPen(i,j,iGrain), j = 1,3) +! enddo +! endif + enddo + + return + +endsubroutine + +!******************************************************************** +! subroutine to compute the equivalent shear modulus from the elasticity tensor +!******************************************************************** +subroutine homogenization_RGC_equivalentShearMod(& + shearMod, & ! equivalent (isotropic) shear modulus +! + elasTens & ! elasticity tensor in Mandel notation + ) + + use prec, only: pReal,pInt,p_vec + use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips + use material, only: homogenization_typeInstance + implicit none + +!* Definition of variables + real(pReal), dimension (6,6), intent(in) :: elasTens + real(pReal), intent(out) :: shearMod + real(pReal) cEquiv_11,cEquiv_12,cEquiv_44 + +!* Compute the equivalent shear modulus using Turterltaub and Suiker, JMPS (2005) + cEquiv_11 = (elasTens(1,1) + elasTens(2,2) + elasTens(3,3))/3.0_pReal + cEquiv_12 = (elasTens(1,2) + elasTens(2,3) + elasTens(3,1) + & + elasTens(1,3) + elasTens(2,1) + elasTens(3,2))/6.0_pReal + cEquiv_44 = (elasTens(4,4) + elasTens(5,5) + elasTens(6,6))/3.0_pReal + shearMod = 0.2_pReal*(cEquiv_11 - cEquiv_12) + 0.6_pReal*cEquiv_44 + + return + +endsubroutine + +!******************************************************************** +! subroutine to collect relaxation vectors of an interface +!******************************************************************** +subroutine homogenization_RGC_relaxationVector(& + aVect, & ! relaxation vector of the interface +! + intFace, & ! set of interface ID in 4D array (normal and position) + state, & ! set of global relaxation vectors + homID & ! homogenization ID + ) + + use prec, only: pReal,pInt,p_vec + use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips + use material, only: homogenization_typeInstance + implicit none + +!* Definition of variables + real(pReal), dimension (3), intent(out) :: aVect + integer(pInt), dimension (4), intent(in) :: intFace + type(p_vec), intent(in) :: state + integer(pInt), dimension (3) :: nGDim + integer(pInt) iNum,homID + +!* Collect the interface relaxation vector from the global state array + aVect = 0.0_pReal + nGDim = homogenization_RGC_Ngrains(:,homID) + call homogenization_RGC_interface4to1(iNum,intFace,homID) ! Get the position in global state array + if (iNum .gt. 0_pInt) aVect = state%p((3*iNum-2):(3*iNum)) ! Collect the corresponding entries + + return + +endsubroutine + +!******************************************************************** +! subroutine to identify the normal of an interface +!******************************************************************** +subroutine homogenization_RGC_interfaceNormal(& + nVect, & ! interface normal +! + intFace & ! interface ID in 4D array (normal and position) + ) + + use prec, only: pReal,pInt,p_vec + use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips + implicit none + +!* Definition of variables + real(pReal), dimension (3), intent(out) :: nVect + integer(pInt), dimension (4), intent(in) :: intFace + integer(pInt) nPos + +!* Get the normal of the interface, identified from the value of intFace(1) + nVect = 0.0_pReal + nPos = abs(intFace(1)) + nVect(nPos) = intFace(1)/abs(intFace(1)) + + return + +endsubroutine + +!******************************************************************** +! subroutine to collect six faces of a grain in 4D (normal and position) +!******************************************************************** +subroutine homogenization_RGC_getInterface(& + intFace, & ! interface ID in 4D (normal and position) +! + iFace, & ! number of faces of grain + iGrain3 & ! grain ID in 3D array + ) + use prec, only: pReal,pInt,p_vec + implicit none + +!* Definition of variables + integer(pInt), dimension (4), intent(out) :: intFace + integer(pInt), dimension (3), intent(in) :: iGrain3 + integer(pInt), intent(in) :: iFace + integer(pInt) iDir + +!* Direction of interface normal + iDir = (int(dble(iFace-1)/2.0_pReal)+1)*(-1_pInt)**iFace + intFace(1) = iDir + +!* Identify the interface position by the direction of its normal + intFace(2:4) = iGrain3(:) + if (iDir .eq. -1_pInt) intFace(2) = intFace(2)-1 + if (iDir .eq. -2_pInt) intFace(3) = intFace(3)-1 + if (iDir .eq. -3_pInt) intFace(4) = intFace(4)-1 + + return + +endsubroutine + +!******************************************************************** +! subroutine to map grain ID from in 1D (array) to in 3D (position) +!******************************************************************** +subroutine homogenization_RGC_grain1to3(& + grain3, & ! grain ID in 3D array (pos.x,pos.y,pos.z) +! + grain1, & ! grain ID in 1D array + homID & ! homogenization ID + ) + + use prec, only: pReal,pInt,p_vec + use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips + implicit none + +!* Definition of variables + integer(pInt), dimension (3), intent(out) :: grain3 + integer(pInt), intent(in) :: grain1,homID + integer(pInt), dimension (3) :: nGDim + +!* Get the grain position + nGDim = homogenization_RGC_Ngrains(:,homID) + grain3(3) = int(dble(grain1-1)/dble(nGDim(1))/dble(nGDim(2)))+1 + grain3(2) = mod(int(dble(grain1-1)/dble(nGDim(1))),nGDim(2))+1 + grain3(1) = mod((grain1-1),nGDim(1))+1 + + return + +endsubroutine + +!******************************************************************** +! subroutine to map grain ID from in 3D (position) to in 1D (array) +!******************************************************************** +subroutine homogenization_RGC_grain3to1(& + grain1, & ! grain ID in 1D array +! + grain3, & ! grain ID in 3D array (pos.x,pos.y,pos.z) + homID & ! homogenization ID + ) + + use prec, only: pReal,pInt,p_vec + use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips + implicit none + +!* Definition of variables + integer(pInt), dimension (3), intent(in) :: grain3 + integer(pInt), intent(out) :: grain1 + integer(pInt), dimension (3) :: nGDim + integer(pInt) homID + +!* Get the grain ID + nGDim = homogenization_RGC_Ngrains(:,homID) + grain1 = grain3(1) + nGDim(1)*(grain3(2)-1) + nGDim(1)*nGDim(2)*(grain3(3)-1) + + return + +endsubroutine + +!******************************************************************** +! subroutine to map interface ID from 4D (normal and position) into 1D (array) +!******************************************************************** +subroutine homogenization_RGC_interface4to1(& + iFace1D, & ! interface ID in 1D array +! + iFace4D, & ! interface ID in 4D array (n.dir,pos.x,pos.y,pos.z) + homID & ! homogenization ID + ) + + use prec, only: pReal,pInt,p_vec + use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips + implicit none + +!* Definition of variables + integer(pInt), dimension (4), intent(in) :: iFace4D + integer(pInt), intent(out) :: iFace1D + integer(pInt), dimension (3) :: nGDim,nIntFace + integer(pInt) homID + + nGDim = homogenization_RGC_Ngrains(:,homID) +!* Get the number of interfaces, which ... + nIntFace(1) = (nGDim(1)-1)*nGDim(2)*nGDim(3) ! ... normal //e1 + nIntFace(2) = nGDim(1)*(nGDim(2)-1)*nGDim(3) ! ... normal //e2 + nIntFace(3) = nGDim(1)*nGDim(2)*(nGDim(3)-1) ! ... normal //e3 + +!* For interface with normal //e1 + if (abs(iFace4D(1)) == 1_pInt) then + iFace1D = iFace4D(3) + nGDim(2)*(iFace4D(4)-1) + nGDim(2)*nGDim(3)*(iFace4D(2)-1) + if ((iFace4D(2) == 0_pInt) .or. (iFace4D(2) == nGDim(1))) iFace1D = 0_pInt +!* For interface with normal //e2 + elseif (abs(iFace4D(1)) == 2_pInt) then + iFace1D = iFace4D(4) + nGDim(3)*(iFace4D(2)-1) + nGDim(3)*nGDim(1)*(iFace4D(3)-1) + nIntFace(1) + if ((iFace4D(3) == 0_pInt) .or. (iFace4D(3) == nGDim(2))) iFace1D = 0_pInt +!* For interface with normal //e3 + elseif (abs(iFace4D(1)) == 3_pInt) then + iFace1D = iFace4D(2) + nGDim(1)*(iFace4D(3)-1) + nGDim(1)*nGDim(2)*(iFace4D(4)-1) + nIntFace(1) + nIntFace(2) + if ((iFace4D(4) == 0_pInt) .or. (iFace4D(4) == nGDim(3))) iFace1D = 0_pInt + endif + + return + +endsubroutine + +!******************************************************************** +! subroutine to map interface ID from 1D (array) into 4D (normal and position) +!******************************************************************** +subroutine homogenization_RGC_interface1to4(& + iFace4D, & ! interface ID in 4D array (n.dir,pos.x,pos.y,pos.z) +! + iFace1D, & ! interface ID in 1D array + homID & ! homogenization ID + ) + + use prec, only: pReal,pInt,p_vec + use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips + implicit none + +!* Definition of variables + integer(pInt), dimension (4), intent(out) :: iFace4D + integer(pInt), intent(in) :: iFace1D + integer(pInt), dimension (3) :: nGDim,nIntFace + integer(pInt) homID + + nGDim = homogenization_RGC_Ngrains(:,homID) +!* Get the number of interfaces, which ... + nIntFace(1) = (nGDim(1)-1)*nGDim(2)*nGDim(3) ! ... normal //e1 + nIntFace(2) = nGDim(1)*(nGDim(2)-1)*nGDim(3) ! ... normal //e2 + nIntFace(3) = nGDim(1)*nGDim(2)*(nGDim(3)-1) ! ... normal //e3 + +!* For interface ID between 1 and nIntFace(1) + if (iFace1D > 0 .and. iFace1D <= nIntFace(1)) then + iFace4D(1) = 1 + iFace4D(3) = mod((iFace1D-1),nGDim(2))+1 + iFace4D(4) = mod(int(dble(iFace1D-1)/dble(nGDim(2))),nGDim(3))+1 + iFace4D(2) = int(dble(iFace1D-1)/dble(nGDim(2))/dble(nGDim(3)))+1 +!* For interface ID between nIntFace(1) and nIntFace(1) + nIntFace(2) + elseif (iFace1D > nIntFace(1) .and. iFace1D <= (nIntFace(2) + nIntFace(1))) then + iFace4D(1) = 2 + iFace4D(4) = mod((iFace1D-nIntFace(1)-1),nGDim(3))+1 + iFace4D(2) = mod(int(dble(iFace1D-nIntFace(1)-1)/dble(nGDim(3))),nGDim(1))+1 + iFace4D(3) = int(dble(iFace1D-nIntFace(1)-1)/dble(nGDim(3))/dble(nGDim(1)))+1 +!* For interface ID between nIntFace(1) + nIntFace(2) and nIntFace(1) + nIntFace(2) + nIntFace(3) + elseif (iFace1D > nIntFace(2) + nIntFace(1) .and. iFace1D <= (nIntFace(3) + nIntFace(2) + nIntFace(1))) then + iFace4D(1) = 3 + iFace4D(2) = mod((iFace1D-nIntFace(2)-nIntFace(1)-1),nGDim(1))+1 + iFace4D(3) = mod(int(dble(iFace1D-nIntFace(2)-nIntFace(1)-1)/dble(nGDim(1))),nGDim(2))+1 + iFace4D(4) = int(dble(iFace1D-nIntFace(2)-nIntFace(1)-1)/dble(nGDim(1))/dble(nGDim(2)))+1 + endif + + return + +endsubroutine + +!******************************************************************** +! calculating the grain deformation gradient +!******************************************************************** +subroutine homogenization_RGC_grainDeformation(& + F, & ! partioned def grad per grain +! + F0, & ! initial partioned def grad per grain + avgF, & ! my average def grad + state, & ! my state + el & ! my element + ) + use prec, only: pReal,pInt,p_vec + use mesh, only: mesh_element + use material, only: homogenization_maxNgrains,homogenization_Ngrains,homogenization_typeInstance + implicit none + +!* Definition of variables + real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F + real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: F0 + real(pReal), dimension (3,3), intent(in) :: avgF + type(p_vec), intent(in) :: state + integer(pInt), intent(in) :: el +! + real(pReal), dimension (3) :: aVect,nVect + integer(pInt), dimension (4) :: intFace + integer(pInt), dimension (3) :: iGrain3 + integer(pInt) homID, iGrain,iFace,i,j +! + integer(pInt), parameter :: nFace = 6 + +!* Compute the deformation gradient of individual grains due to relaxations + homID = homogenization_typeInstance(mesh_element(3,el)) + F = 0.0_pReal + do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) + call homogenization_RGC_grain1to3(iGrain3,iGrain,homID) + do iFace = 1,nFace + call homogenization_RGC_getInterface(intFace,iFace,iGrain3) + call homogenization_RGC_relaxationVector(aVect,intFace,state,homID) + call homogenization_RGC_interfaceNormal(nVect,intFace) + forall (i=1:3,j=1:3) & + F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations + enddo + F(:,:,iGrain) = F(:,:,iGrain) + avgF(:,:) ! relaxed deformation gradient + enddo + + return + +endsubroutine + +END MODULE diff --git a/code/mesh.f90 b/code/mesh.f90 index 17485e52b..04899202d 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -42,10 +42,13 @@ ! order is +x,-x,+y,-y,+z,-z but meaning strongly depends on Elemtype ! --------------------------- integer(pInt) mesh_Nelems,mesh_NcpElems,mesh_NelemSets,mesh_maxNelemInSet + integer(pInt) mesh_Nmaterials integer(pInt) mesh_Nnodes,mesh_maxNnodes,mesh_maxNips,mesh_maxNipNeighbors,mesh_maxNsharedElems,mesh_maxNsubNodes integer(pInt), dimension(2) :: mesh_maxValStateVar = 0_pInt - character(len=64), dimension(:), allocatable :: mesh_nameElemSet - integer(pInt), dimension(:,:), allocatable :: mesh_mapElemSet + character(len=64), dimension(:), allocatable :: mesh_nameElemSet, & ! names of elementSet + mesh_nameMaterial, & ! names of material in solid section + mesh_mapMaterial ! name of elementSet for material + integer(pInt), dimension(:,:), allocatable :: mesh_mapElemSet ! list of elements in elementSet integer(pInt), dimension(:,:), allocatable, target :: mesh_mapFEtoCPelem, mesh_mapFEtoCPnode integer(pInt), dimension(:,:), allocatable :: mesh_element, mesh_sharedElem integer(pInt), dimension(:,:,:,:), allocatable :: mesh_ipNeighborhood @@ -62,8 +65,8 @@ integer(pInt), dimension(:,:,:), allocatable :: FE_subNodeParent integer(pInt), dimension(:,:,:,:), allocatable :: FE_subNodeOnIPFace - integer(pInt) :: hypoelasticTableStyle = 0 - integer(pInt) :: initialcondTableStyle = 0 + integer(pInt) :: hypoelasticTableStyle + integer(pInt) :: initialcondTableStyle integer(pInt), parameter :: FE_Nelemtypes = 7 integer(pInt), parameter :: FE_maxNnodes = 8 integer(pInt), parameter :: FE_maxNsubNodes = 56 @@ -182,53 +185,75 @@ !*********************************************************** ! initialization !*********************************************************** - subroutine mesh_init () + subroutine mesh_init (ip,element) + use cpfem_interface use prec, only: pInt use IO, only: IO_error,IO_open_InputFile - use FEsolving, only: parallelExecution, FEsolving_execElem, FEsolving_execIP + use FEsolving, only: parallelExecution, FEsolving_execElem, FEsolving_execIP, calcMode, lastMode implicit none integer(pInt), parameter :: fileUnit = 222 - integer(pInt) e + integer(pInt) e,element,ip write(6,*) write(6,*) '<<<+- mesh init -+>>>' write(6,*) '$Id$' write(6,*) - mesh_Nelems = 0_pInt - mesh_NcpElems = 0_pInt - mesh_Nnodes = 0_pInt - mesh_maxNips = 0_pInt - mesh_maxNnodes = 0_pInt - mesh_maxNipNeighbors = 0_pInt - mesh_maxNsharedElems = 0_pInt - mesh_maxNsubNodes = 0_pInt - mesh_NelemSets = 0_pInt - mesh_maxNelemInSet = 0_pInt - - ! get properties of the different types of elements - call mesh_get_FEdata() + call mesh_build_FEdata() ! --- get properties of the different types of elements -! call to various subroutines to parse the stuff from the input file... - if (IO_open_inputFile(fileUnit)) then + if (IO_open_inputFile(fileUnit)) then ! --- parse info from input file... - call mesh_get_meshDimensions(fileUnit) - call mesh_build_nodeMapping(fileUnit) - call mesh_build_elemMapping(fileUnit) - call mesh_build_elemSetMapping(fileUnit) - call mesh_get_nodeElemDimensions(fileUnit) - call mesh_build_nodes(fileUnit) - call mesh_build_elements(fileUnit) + select case (FEsolver) + case ('Marc') + write(6,*) '<<---mesh_marc_get_tableStyles--->>'; call flush(6) + call mesh_marc_get_tableStyles(fileUnit) + write(6,*) '<<---mesh_marc_count_nodesAndElements--->>'; call flush(6) + call mesh_marc_count_nodesAndElements(fileUnit) + write(6,*) '<<---mesh_marc_count_elementSets--->>'; call flush(6) + call mesh_marc_count_elementSets(fileUnit) + write(6,*) '<<---mesh_marc_map_elementSets--->>'; call flush(6) + call mesh_marc_map_elementSets(fileUnit) + write(6,*) '<<---mesh_marc_count_cpElements--->>'; call flush(6) + call mesh_marc_count_cpElements(fileUnit) + write(6,*) '<<---mesh_marc_map_elements--->>'; call flush(6) + call mesh_marc_map_elements(fileUnit) + write(6,*) '<<---mesh_marc_map_nodes--->>'; call flush(6) + call mesh_marc_map_nodes(fileUnit) + write(6,*) '<<---mesh_marc_build_nodes--->>'; call flush(6) + call mesh_marc_build_nodes(fileUnit) + write(6,*) '<<---mesh_marc_count_cpSizes--->>'; call flush(6) + call mesh_marc_count_cpSizes(fileunit) + write(6,*) '<<---mesh_marc_build_elements--->>'; call flush(6) + call mesh_marc_build_elements(fileUnit) + case ('Abaqus') + call mesh_abaqus_count_nodesAndElements(fileUnit) + call mesh_abaqus_count_elementSets(fileUnit) + call mesh_abaqus_count_materials(fileUnit) + call mesh_abaqus_map_elementSets(fileUnit) + call mesh_abaqus_map_materials(fileUnit) + call mesh_abaqus_count_cpElements(fileUnit) + call mesh_abaqus_map_elements(fileUnit) + call mesh_abaqus_map_nodes(fileUnit) + call mesh_abaqus_build_nodes(fileUnit) + call mesh_abaqus_count_cpSizes(fileunit) + call mesh_abaqus_build_elements(fileUnit) + end select + close (fileUnit) + + write(6,*) '<<---mesh_build_sharedElems--->>'; call flush(6) call mesh_build_sharedElems() + write(6,*) '<<---mesh_build_ipNeighborhood--->>'; call flush(6) call mesh_build_ipNeighborhood() + write(6,*) '<<---mesh_build_subNodeCoords--->>'; call flush(6) call mesh_build_subNodeCoords() + write(6,*) '<<---mesh_build_ipVolumes--->>'; call flush(6) call mesh_build_ipVolumes() + write(6,*) '<<---mesh_build_ipAreas--->>'; call flush(6) call mesh_build_ipAreas() call mesh_tell_statistics() - close (fileUnit) parallelExecution = (mesh_Nelems == mesh_NcpElems) ! plus potential killer from non-local constitutive else @@ -239,6 +264,10 @@ allocate(FEsolving_execIP(2,mesh_NcpElems)); FEsolving_execIP = 1_pInt forall (e = 1:mesh_NcpElems) FEsolving_execIP(2,e) = FE_Nips(mesh_element(2,e)) + allocate(calcMode(mesh_maxNips,mesh_NcpElems)) + calcMode = .false. ! pretend to have collected what first call is asking (F = I) + calcMode(ip,mesh_FEasCP('elem',element)) = .true. ! first ip,el needs to be already pingponged to "calc" + lastMode = .true. ! and its mode is already known... endsubroutine @@ -254,22 +283,22 @@ integer(pInt) FE_mapElemtype select case (what) - case ('7', 'C3D8') - FE_mapElemtype = 1 ! Three-dimensional Arbitrarily Distorted Brick - case ('134') - FE_mapElemtype = 2 ! Three-dimensional Four-node Tetrahedron - case ('11') - FE_mapElemtype = 3 ! Arbitrary Quadrilateral Plane-strain - case ('27') - FE_mapElemtype = 4 ! Plane Strain, Eight-node Distorted Quadrilateral + case ( '7', 'C3D8') + FE_mapElemtype = 1 ! Three-dimensional Arbitrarily Distorted Brick + case ('134', 'C3D4') + FE_mapElemtype = 2 ! Three-dimensional Four-node Tetrahedron + case ( '11', 'CPE4') + FE_mapElemtype = 3 ! Arbitrary Quadrilateral Plane-strain + case ( '27', 'CPE8') + FE_mapElemtype = 4 ! Plane Strain, Eight-node Distorted Quadrilateral case ('157') - FE_mapElemtype = 5 ! Three-dimensional, Low-order, Tetrahedron, Herrmann Formulations - case ('136') - FE_mapElemtype = 6 ! Three-dimensional Arbitrarily Distorted Pentahedral - case ('21') + FE_mapElemtype = 5 ! Three-dimensional, Low-order, Tetrahedron, Herrmann Formulations + case ('136', 'C3D6') + FE_mapElemtype = 6 ! Three-dimensional Arbitrarily Distorted Pentahedral + case ( '21', 'C3D20') FE_mapElemtype = 7 ! Three-dimensional Arbitrarily Distorted qudratic hexahedral - case default - FE_mapElemtype = 0 ! unknown element --> should raise an error upstream..! + case default + FE_mapElemtype = 0 ! unknown element --> should raise an error upstream..! endselect endfunction @@ -333,7 +362,7 @@ !*********************************************************** ! find face-matching element of same type -!!*********************************************************** +!*********************************************************** function mesh_faceMatch(face,elem) use prec, only: pInt @@ -344,22 +373,22 @@ integer(pInt), dimension(FE_NfaceNodes(face,mesh_element(2,elem))) :: nodeMap integer(pInt) minN,NsharedElems,lonelyNode,faceNode,i,n,t - minN = mesh_maxNsharedElems+1 ! init to worst case - mesh_faceMatch = 0_pInt ! intialize to "no match found" - t = mesh_element(2,elem) ! figure elemType + minN = mesh_maxNsharedElems+1 ! init to worst case + mesh_faceMatch = 0_pInt ! intialize to "no match found" + t = mesh_element(2,elem) ! figure elemType - do faceNode=1,FE_NfaceNodes(face,t) ! loop over nodes on face + do faceNode=1,FE_NfaceNodes(face,t) ! loop over nodes on face nodeMap(faceNode) = mesh_FEasCP('node',mesh_element(4+FE_nodeOnFace(faceNode,face,t),elem)) ! CP id of face node NsharedElems = mesh_sharedElem(1,nodeMap(faceNode)) ! figure # shared elements for this node if (NsharedElems < minN) then - minN = NsharedElems ! remember min # shared elems - lonelyNode = faceNode ! remember most lonely node + minN = NsharedElems ! remember min # shared elems + lonelyNode = faceNode ! remember most lonely node endif enddo candidate: do i=1,minN ! iterate over lonelyNode's shared elements mesh_faceMatch = mesh_sharedElem(1+i,nodeMap(lonelyNode)) ! present candidate elem - if (mesh_faceMatch == elem) then ! my own element ? - mesh_faceMatch = 0_pInt ! disregard + if (mesh_faceMatch == elem) then ! my own element ? + mesh_faceMatch = 0_pInt ! disregard cycle candidate endif do faceNode=1,FE_NfaceNodes(face,t) ! check remaining face nodes to match @@ -370,7 +399,7 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements cycle candidate ! next candidate elem endif enddo - exit ! surviving candidate + exit ! surviving candidate enddo candidate return @@ -384,7 +413,7 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements ! assign globals: ! FE_nodesAtIP, FE_ipNeighbor, FE_subNodeParent, FE_subNodeOnIPFace !******************************************************************** - subroutine mesh_get_FEdata () + subroutine mesh_build_FEdata () use prec, only: pInt implicit none @@ -1020,51 +1049,39 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements endsubroutine + !******************************************************************** -! get count of elements, nodes, and cp elements in mesh -! for subsequent array allocations +! figure out table styles (Marc only) ! -! assign globals: -! _Nelems, _Nnodes, _NcpElems +! initialcondTableStyle, hypoelasticTableStyle !******************************************************************** - subroutine mesh_get_meshDimensions (unit) + subroutine mesh_marc_get_tableStyles (unit) use prec, only: pInt use IO implicit none - integer(pInt) unit,i,pos(41) - character*300 line + integer(pInt), parameter :: maxNchunks = 6 + integer(pInt), dimension (1+2*maxNchunks) :: pos + integer(pInt) unit + character(len=300) line + + initialcondTableStyle = 0_pInt + hypoelasticTableStyle = 0_pInt + 610 FORMAT(A300) rewind(unit) do read (unit,610,END=620) line - pos = IO_stringPos(line,20) - - select case ( IO_lc(IO_StringValue(line,pos,1))) - case('table') - if (pos(1) == 6) then - initialcondTableStyle = IO_IntValue (line,pos,4) - hypoelasticTableStyle = IO_IntValue (line,pos,5) - endif - case('sizing') - mesh_Nelems = IO_IntValue (line,pos,3) - mesh_Nnodes = IO_IntValue (line,pos,4) - case('define') - select case (IO_lc(IO_StringValue(line,pos,2))) - case('element') ! Count the number of encountered element sets - mesh_NelemSets=mesh_NelemSets+1 - mesh_maxNelemInSet = max(mesh_maxNelemInSet,IO_countContinousIntValues(unit)) - endselect - case('hypoelastic') - do i=1,3+hypoelasticTableStyle ! Skip 3 or 4 lines - read (unit,610,END=620) line - enddo - mesh_NcpElems = mesh_NcpElems + IO_countContinousIntValues(unit) - endselect + pos = IO_stringPos(line,maxNchunks) + if ( IO_lc(IO_stringValue(line,pos,1)) == 'table' .and. pos(1) == 6) then + initialcondTableStyle = IO_intValue(line,pos,4) + hypoelasticTableStyle = IO_intValue(line,pos,5) + exit + endif enddo 620 return @@ -1072,133 +1089,501 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements endsubroutine -!!******************************************************************** -! get maximum count of nodes, IPs, IP neighbors, and shared elements -! for subsequent array allocations -! -! assign globals: -! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsharedElems !******************************************************************** -subroutine mesh_get_nodeElemDimensions (unit) - +! count overall number of nodes and elements in mesh +! +! mesh_Nelems, mesh_Nnodes +!******************************************************************** + subroutine mesh_marc_count_nodesAndElements (unit) + use prec, only: pInt use IO implicit none - - integer(pInt), parameter :: maxNchunks = 66 - integer(pInt), dimension (mesh_Nnodes) :: node_count - integer(pInt), dimension (:), allocatable :: node_seen - integer(pInt) unit,i,j,n,t,e + + integer(pInt), parameter :: maxNchunks = 4 integer(pInt), dimension (1+2*maxNchunks) :: pos - character*300 line - + + integer(pInt) unit + character(len=300) line + + mesh_Nnodes = 0_pInt + mesh_Nelems = 0_pInt + 610 FORMAT(A300) - - node_count = 0_pInt - allocate(node_seen(maxval(FE_Nnodes))) - + rewind(unit) - do - read (unit,610,END=630) line - pos = IO_stringPos(line,1) - if( IO_lc(IO_stringValue(line,pos,1)) == 'connectivity' ) then - read (unit,610,END=630) line ! Garbage line - do i=1,mesh_Nelems ! read all elements - read (unit,610,END=630) line - pos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max (plus ID, type) - e = mesh_FEasCP('elem',IO_intValue(line,pos,1)) - if (e /= 0) then - t = FE_mapElemtype(IO_StringValue(line,pos,2)) - mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(t)) - mesh_maxNips = max(mesh_maxNips,FE_Nips(t)) - mesh_maxNipNeighbors = max(mesh_maxNipNeighbors,FE_NipNeighbors(t)) - mesh_maxNsubNodes = max(mesh_maxNsubNodes,FE_NsubNodes(t)) - node_seen = 0_pInt - do j=1,FE_Nnodes(t) - n = mesh_FEasCP('node',IO_IntValue (line,pos,j+2)) - if (all(node_seen /= n)) node_count(n) = node_count(n)+1 - node_seen(j) = n - enddo - call IO_skipChunks(unit,FE_NoriginalNodes(t)-(pos(1)-2)) ! read on if FE_Nnodes exceeds node count present on current line - endif - enddo + do + read (unit,610,END=620) line + pos = IO_stringPos(line,maxNchunks) + + if ( IO_lc(IO_StringValue(line,pos,1)) == 'sizing') then + mesh_Nelems = IO_IntValue (line,pos,3) + mesh_Nnodes = IO_IntValue (line,pos,4) exit endif enddo + +620 return -630 mesh_maxNsharedElems = maxval(node_count) - - return endsubroutine +!******************************************************************** +! count overall number of nodes and elements in mesh +! +! mesh_Nelems, mesh_Nnodes +!******************************************************************** + subroutine mesh_abaqus_count_nodesAndElements (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 2 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt) unit + logical inPart + + mesh_Nnodes = 0 + mesh_Nelems = 0 + +610 FORMAT(A300) + + inPart = .false. + rewind(unit) + do + read (unit,610,END=620) line + pos = IO_stringPos(line,maxNchunks) + if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & + IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. + + if (inPart) then + select case ( IO_lc(IO_stringValue(line,pos,1))) + case('*node') + if( & + IO_lc(IO_stringValue(line,pos,2)) /= 'output' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'print' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'file' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'response' & + ) & + mesh_Nnodes = mesh_Nnodes + IO_countDataLines(unit) + case('*element') + if( & + IO_lc(IO_stringValue(line,pos,2)) /= 'output' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'matrix' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'response' & + ) then + mesh_Nelems = mesh_Nelems + IO_countDataLines(unit) + endif + endselect + endif + enddo + +620 return + + endsubroutine + + +!******************************************************************** +! count overall number of element sets in mesh +! +! mesh_NelemSets, mesh_maxNelemInSet +!******************************************************************** + subroutine mesh_marc_count_elementSets (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 2 + integer(pInt), dimension (1+2*maxNchunks) :: pos + + integer(pInt) unit + character(len=300) line + + mesh_NelemSets = 0_pInt + mesh_maxNelemInSet = 0_pInt + +610 FORMAT(A300) + + rewind(unit) + do + read (unit,610,END=620) line + pos = IO_stringPos(line,maxNchunks) + + if ( IO_lc(IO_StringValue(line,pos,1)) == 'define' .and. & + IO_lc(IO_StringValue(line,pos,2)) == 'element' ) then + mesh_NelemSets = mesh_NelemSets + 1_pInt + mesh_maxNelemInSet = max(mesh_maxNelemInSet, & + IO_countContinousIntValues(unit)) + endif + enddo + +620 return + + endsubroutine + + +!******************************************************************** +! count overall number of element sets in mesh +! +! mesh_NelemSets, mesh_maxNelemInSet +!******************************************************************** + subroutine mesh_abaqus_count_elementSets (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 2 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt) unit + logical inPart + + mesh_NelemSets = 0_pInt + mesh_maxNelemInSet = mesh_Nelems ! have to be conservative, since Abaqus allows for recursive definitons + +610 FORMAT(A300) + + inPart = .false. + rewind(unit) + do + read (unit,610,END=620) line + pos = IO_stringPos(line,maxNchunks) + if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & + IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. + + if ( inPart .and. IO_lc(IO_stringValue(line,pos,1)) == '*elset' ) & + mesh_NelemSets = mesh_NelemSets + 1_pInt + enddo + +620 return + + endsubroutine + + +!******************************************************************** +! count overall number of solid sections sets in mesh (Abaqus only) +! +! mesh_Nmaterials +!******************************************************************** + subroutine mesh_abaqus_count_materials (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 2 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt) unit + logical inPart + + mesh_Nmaterials = 0_pInt + +610 FORMAT(A300) + + inPart = .false. + rewind(unit) + do + read (unit,610,END=620) line + pos = IO_stringPos(line,maxNchunks) + if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & + IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. + + if ( inPart .and. & + IO_lc(IO_StringValue(line,pos,1)) == '*solid' .and. & + IO_lc(IO_StringValue(line,pos,2)) == 'section' ) & + mesh_Nmaterials = mesh_Nmaterials + 1_pInt + enddo + +620 return + + endsubroutine + + +!******************************************************************** +! count overall number of cpElements in mesh +! +! mesh_NcpElems +!******************************************************************** + subroutine mesh_marc_count_cpElements (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 1 + integer(pInt), dimension (1+2*maxNchunks) :: pos + + integer(pInt) unit,i + character(len=300) line + + mesh_NcpElems = 0_pInt + +610 FORMAT(A300) + + rewind(unit) + do + read (unit,610,END=620) line + pos = IO_stringPos(line,maxNchunks) + + if ( IO_lc(IO_stringValue(line,pos,1)) == 'hypoelastic') then + do i=1,3+hypoelasticTableStyle ! Skip 3 or 4 lines + read (unit,610,END=620) line + enddo + mesh_NcpElems = mesh_NcpElems + IO_countContinousIntValues(unit) + exit + endif + enddo + +620 return + + endsubroutine + + +!******************************************************************** +! count overall number of cpElements in mesh +! +! mesh_NcpElems +!******************************************************************** + subroutine mesh_abaqus_count_cpElements (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 2 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt) unit,i + logical materialFound + character (len=64) materialName,elemSetName + + mesh_NcpElems = 0 + +610 FORMAT(A300) + + rewind(unit) + do + read (unit,610,END=620) line + pos = IO_stringPos(line,maxNchunks) + select case ( IO_lc(IO_stringValue(line,pos,1)) ) + case('*material') + materialName = IO_extractValue(IO_lc(IO_stringValue(line,pos,2)),'name') ! extract name=value + materialFound = materialName /= '' ! valid name? + case('*user') + if (IO_lc(IO_StringValue(line,pos,2)) == 'material' .and. materialFound) then + do i = 1,mesh_Nmaterials ! look thru material names + if (materialName == mesh_nameMaterial(i)) exit ! found one + enddo + if (i <= mesh_Nmaterials) then ! found one? + elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet + do i = 1,mesh_NelemSets ! look thru all elemSet definitions + if (elemSetName == mesh_nameElemSet(i)) & ! matched? + mesh_NcpElems = mesh_NcpElems + mesh_mapElemSet(1,i) ! add those elem count + enddo + endif + materialFound = .false. + endif + endselect + enddo + +620 return + + endsubroutine + + +!******************************************************************** +! map element sets +! +! allocate globals: mesh_nameElemSet, mesh_mapElemSet +!******************************************************************** + subroutine mesh_marc_map_elementSets (unit) + + use prec, only: pInt + use IO + + implicit none + + integer(pInt), parameter :: maxNchunks = 4 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt) unit,elemSet + + allocate (mesh_nameElemSet(mesh_NelemSets)) ; mesh_nameElemSet = '' + allocate (mesh_mapElemSet(1+mesh_maxNelemInSet,mesh_NelemSets)) ; mesh_mapElemSet = 0_pInt + +610 FORMAT(A300) + + elemSet = 0_pInt + rewind(unit) + do + read (unit,610,END=640) line + pos = IO_stringPos(line,maxNchunks) + if( (IO_lc(IO_stringValue(line,pos,1)) == 'define' ) .and. & + (IO_lc(IO_stringValue(line,pos,2)) == 'element' ) ) then + elemSet = elemSet+1 + mesh_nameElemSet(elemSet) = IO_stringValue(line,pos,4) + mesh_mapElemSet(:,elemSet) = IO_continousIntValues(unit,mesh_maxNelemInSet,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) + endif + enddo + +640 return + + endsubroutine + + !******************************************************************** ! Build element set mapping ! ! allocate globals: mesh_nameElemSet, mesh_mapElemSet !******************************************************************** - subroutine mesh_build_elemSetMapping (unit) + subroutine mesh_abaqus_map_elementSets (unit) use prec, only: pInt use IO implicit none - integer unit, elem_set - character*300 line - integer(pInt), dimension (9) :: pos ! count plus 4 entities on a line + integer(pInt), parameter :: maxNchunks = 4 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt) unit,elemSet + logical inPart + + allocate (mesh_nameElemSet(mesh_NelemSets)) ; mesh_nameElemSet = '' + allocate (mesh_mapElemSet(1+mesh_maxNelemInSet,mesh_NelemSets)) ; mesh_mapElemSet = 0_pInt 610 FORMAT(A300) - allocate (mesh_nameElemSet(mesh_NelemSets)) - allocate (mesh_mapElemSet(1+mesh_maxNelemInSet,mesh_NelemSets)) ; mesh_mapElemSet = 0_pInt - elem_set = 0_pInt - + elemSet = 0_pInt + inPart = .false. rewind(unit) do read (unit,610,END=640) line - pos = IO_stringPos(line,4) - if( (IO_lc(IO_stringValue(line,pos,1)) == 'define' ).and. & - (IO_lc(IO_stringValue(line,pos,2)) == 'element' ) )then - elem_set = elem_set+1 - mesh_nameElemSet(elem_set) = IO_stringValue(line,pos,4) - mesh_mapElemSet(:,elem_set) = IO_continousIntValues(unit,mesh_maxNelemInSet,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) + pos = IO_stringPos(line,maxNchunks) + if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & + IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. + + if ( inPart .and. IO_lc(IO_stringValue(line,pos,1)) == '*elset' ) then + elemSet = elemSet + 1_pInt + mesh_nameElemSet(elemSet) = IO_extractValue(IO_lc(IO_stringValue(line,pos,2)),'elset') + mesh_mapElemSet(:,elemSet) = IO_continousIntValues(unit,mesh_Nelems,mesh_nameElemSet,mesh_mapElemSet,elemSet-1) endif enddo + 640 return + endsubroutine - + !******************************************************************** -! Build node mapping from FEM to CP +! map solid section (Abaqus only) ! -! allocate globals: -! _mapFEtoCPnode +! allocate globals: mesh_nameMaterial, mesh_mapMaterial !******************************************************************** - subroutine mesh_build_nodeMapping (unit) + subroutine mesh_abaqus_map_materials (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 20 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt) unit,i,count + logical inPart,materialFound + character(len=64) elemSetName,materialName + + allocate (mesh_nameMaterial(mesh_Nmaterials)) ; mesh_nameMaterial = '' + allocate (mesh_mapMaterial(mesh_Nmaterials)) ; mesh_mapMaterial = '' + +610 FORMAT(A300) + + count = 0 + inPart = .false. + rewind(unit) + do + read (unit,610,END=620) line + pos = IO_stringPos(line,maxNchunks) + if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & + IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. + + if ( inPart .and. & + IO_lc(IO_StringValue(line,pos,1)) == '*solid' .and. & + IO_lc(IO_StringValue(line,pos,2)) == 'section' ) then + + elemSetName = '' + materialName = '' + + do i = 3,pos(1) + if (IO_extractValue(IO_lc(IO_stringValue(line,pos,i)),'elset') /= '') & + elemSetName = IO_extractValue(IO_lc(IO_stringValue(line,pos,i)),'elset') + if (IO_extractValue(IO_lc(IO_stringValue(line,pos,i)),'material') /= '') & + materialName = IO_extractValue(IO_lc(IO_stringValue(line,pos,i)),'material') + enddo + + if (elemSetName /= '' .and. materialName /= '') then + count = count + 1_pInt + mesh_nameMaterial(count) = materialName ! name of material used for this section + mesh_mapMaterial(count) = elemSetName ! mapped to respective element set + endif + endif + enddo + +620 return + + endsubroutine + + + + !******************************************************************** +! map nodes from FE id to internal (consecutive) representation +! +! allocate globals: mesh_mapFEtoCPnode +!******************************************************************** + subroutine mesh_marc_map_nodes (unit) use prec, only: pInt use math, only: qsort use IO + implicit none integer(pInt), parameter :: maxNchunks = 1 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + integer(pInt), dimension (mesh_Nnodes) :: node_count integer(pInt) unit,i - integer(pInt), dimension (1+2*maxNchunks) :: pos - character*300 line + + allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt 610 FORMAT(A300) - allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt - node_count(:) = 0_pInt + node_count = 0_pInt rewind(unit) do read (unit,610,END=650) line pos = IO_stringPos(line,maxNchunks) if( IO_lc(IO_stringValue(line,pos,1)) == 'coordinates' ) then - read (unit,610,END=650) line ! skip crap line - do i=1,mesh_Nnodes + read (unit,610,END=650) line ! skip crap line + do i = 1,mesh_Nnodes read (unit,610,END=650) line mesh_mapFEtoCPnode(1,i) = IO_fixedIntValue (line,(/0,10/),1) mesh_mapFEtoCPnode(2,i) = i @@ -1210,16 +1595,17 @@ subroutine mesh_get_nodeElemDimensions (unit) 650 call qsort(mesh_mapFEtoCPnode,1,size(mesh_mapFEtoCPnode,2)) return + endsubroutine + !******************************************************************** -! Build element mapping from FEM to CP +! map nodes from FE id to internal (consecutive) representation ! -! allocate globals: -! _mapFEtoCPelem +! allocate globals: mesh_mapFEtoCPnode !******************************************************************** - subroutine mesh_build_elemMapping (unit) + subroutine mesh_abaqus_map_nodes (unit) use prec, only: pInt use math, only: qsort @@ -1227,37 +1613,277 @@ subroutine mesh_get_nodeElemDimensions (unit) implicit none - integer unit, i,CP_elem - character*300 line - integer(pInt), dimension (3) :: pos - integer(pInt), dimension (1+mesh_NcpElems) :: contInts + integer(pInt), parameter :: maxNchunks = 2 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + integer(pInt) unit,i,count,cpNode + logical inPart + + allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt 610 FORMAT(A300) - allocate (mesh_mapFEtoCPelem(2,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt - CP_elem = 0_pInt - + cpNode = 0_pInt + inPart = .false. rewind(unit) do - read (unit,610,END=660) line - pos = IO_stringPos(line,1) - if( IO_lc(IO_stringValue(line,pos,1)) == 'hypoelastic' ) then - do i=1,3+hypoelasticTableStyle ! skip three (or four if new table style!) lines - read (unit,610,END=660) line + read (unit,610,END=650) line + pos = IO_stringPos(line,maxNchunks) + if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & + IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. + + if( inPart .and. & + IO_lc(IO_stringValue(line,pos,1)) == '*node' .and. & + ( IO_lc(IO_stringValue(line,pos,2)) /= 'output' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'print' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'file' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'response' ) & + ) then + count = IO_countDataLines(unit) + do i = 1,count + backspace(unit) enddo - contInts = IO_continousIntValues(unit,mesh_NcpElems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) - do i = 1,contInts(1) - CP_elem = CP_elem+1 - mesh_mapFEtoCPelem(1,CP_elem) = contInts(1+i) - mesh_mapFEtoCPelem(2,CP_elem) = CP_elem + do i = 1,count + read (unit,610,END=650) line + pos = IO_stringPos(line,maxNchunks) + cpNode = cpNode + 1_pInt + mesh_mapFEtoCPnode(1,cpNode) = IO_intValue(line,pos,1) + mesh_mapFEtoCPnode(2,cpNode) = cpNode enddo endif enddo -660 call qsort(mesh_mapFEtoCPelem,1,size(mesh_mapFEtoCPelem,2)) ! should be mesh_NcpElems +650 call qsort(mesh_mapFEtoCPnode,1,size(mesh_mapFEtoCPnode,2)) return + + endsubroutine + + +!******************************************************************** +! map elements from FE id to internal (consecutive) representation +! +! allocate globals: mesh_mapFEtoCPelem +!******************************************************************** + subroutine mesh_marc_map_elements (unit) + + use prec, only: pInt + use math, only: qsort + use IO + + implicit none + + integer(pInt), parameter :: maxNchunks = 1 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt), dimension (1+mesh_NcpElems) :: contInts + integer(pInt) unit,i,cpElem + + allocate (mesh_mapFEtoCPelem(2,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt + +610 FORMAT(A300) + + cpElem = 0_pInt + rewind(unit) + do + read (unit,610,END=660) line + pos = IO_stringPos(line,maxNchunks) + if( IO_lc(IO_stringValue(line,pos,1)) == 'hypoelastic' ) then + do i=1,3+hypoelasticTableStyle ! skip three (or four if new table style!) lines + read (unit,610,END=660) line + enddo + contInts = IO_continousIntValues(unit,mesh_NcpElems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) + do i = 1,contInts(1) + cpElem = cpElem+1 + mesh_mapFEtoCPelem(1,cpElem) = contInts(1+i) + mesh_mapFEtoCPelem(2,cpElem) = cpElem + enddo + endif + enddo + +660 call qsort(mesh_mapFEtoCPelem,1,size(mesh_mapFEtoCPelem,2)) ! should be mesh_NcpElems + + return + endsubroutine + + +!******************************************************************** +! map elements from FE id to internal (consecutive) representation +! +! allocate globals: mesh_mapFEtoCPelem +!******************************************************************** + subroutine mesh_abaqus_map_elements (unit) + + use prec, only: pInt + use math, only: qsort + use IO + + implicit none + + integer(pInt), parameter :: maxNchunks = 2 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt) unit,i,j,cpElem + logical materialFound + character (len=64) materialName,elemSetName + + allocate (mesh_mapFEtoCPelem(2,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt + +610 FORMAT(A300) + + cpElem = 0_pInt + rewind(unit) + do + read (unit,610,END=660) line + pos = IO_stringPos(line,maxNchunks) + select case ( IO_lc(IO_stringValue(line,pos,1)) ) + case('*material') + materialName = IO_extractValue(IO_lc(IO_stringValue(line,pos,2)),'name') ! extract name=value + materialFound = materialName /= '' ! valid name? + case('*user') + if (IO_lc(IO_stringValue(line,pos,2)) == 'material' .and. materialFound) then + do i = 1,mesh_Nmaterials ! look thru material names + if (materialName == mesh_nameMaterial(i)) exit ! found one + enddo + if (i <= mesh_Nmaterials) then ! found one? + elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet + do i = 1,mesh_NelemSets ! look thru all elemSet definitions + if (elemSetName == mesh_nameElemSet(i)) then ! matched? + do j = 1,mesh_mapElemSet(1,i) + cpElem = cpElem + 1_pInt + mesh_mapFEtoCPelem(1,cpElem) = mesh_mapElemSet(1+j,i) ! store FE id + mesh_mapFEtoCPelem(2,cpElem) = cpElem ! store our id + enddo + endif + enddo + endif + materialFound = .false. + endif + endselect + enddo + +660 call qsort(mesh_mapFEtoCPelem,1,size(mesh_mapFEtoCPelem,2)) ! should be mesh_NcpElems + + return + endsubroutine + + +!******************************************************************** +! get maximum count of nodes, IPs, IP neighbors, and subNodes +! among cpElements +! +! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsubNodes +!******************************************************************** +subroutine mesh_marc_count_cpSizes (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 2 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt) unit,i,t,e + + mesh_maxNnodes = 0_pInt + mesh_maxNips = 0_pInt + mesh_maxNipNeighbors = 0_pInt + mesh_maxNsubNodes = 0_pInt + +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)) == 'connectivity' ) then + read (unit,610,END=630) line ! Garbage line + do i=1,mesh_Nelems ! read all elements + read (unit,610,END=630) line + pos = IO_stringPos(line,maxNchunks) ! limit to id and type + e = mesh_FEasCP('elem',IO_intValue(line,pos,1)) + if (e /= 0) then + t = FE_mapElemtype(IO_stringValue(line,pos,2)) + mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(t)) + mesh_maxNips = max(mesh_maxNips,FE_Nips(t)) + mesh_maxNipNeighbors = max(mesh_maxNipNeighbors,FE_NipNeighbors(t)) + mesh_maxNsubNodes = max(mesh_maxNsubNodes,FE_NsubNodes(t)) + call IO_skipChunks(unit,FE_NoriginalNodes(t)-(pos(1)-2)) ! read on if FE_Nnodes exceeds node count present on current line + endif + enddo + exit + endif + enddo + +630 return + endsubroutine + + +!******************************************************************** +! get maximum count of nodes, IPs, IP neighbors, and subNodes +! among cpElements +! +! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsubNodes +!******************************************************************** + subroutine mesh_abaqus_count_cpSizes (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 2 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt) unit,i,count,t + logical inPart + + mesh_maxNnodes = 0_pInt + mesh_maxNips = 0_pInt + mesh_maxNipNeighbors = 0_pInt + mesh_maxNsubNodes = 0_pInt + +610 FORMAT(A300) + + inPart = .false. + rewind(unit) + do + read (unit,610,END=620) line + pos = IO_stringPos(line,2) + if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & + IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. + + if( inPart .and. & + IO_lc(IO_stringValue(line,pos,1)) == '*element' .and. & + ( IO_lc(IO_stringValue(line,pos,2)) /= 'output' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'matrix' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'response' ) & + ) then + t = FE_mapElemtype(IO_extractValue(IO_stringValue(line,pos,2),'type')) ! remember elem type + count = IO_countDataLines(unit) + do i = 1,count + backspace(unit) + enddo + do i = 1,count + read (unit,610,END=620) line + pos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max + if (mesh_FEasCP('elem',IO_intValue(line,pos,1)) /= 0) then ! disregard non CP elems + mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(t)) + mesh_maxNips = max(mesh_maxNips,FE_Nips(t)) + mesh_maxNipNeighbors = max(mesh_maxNipNeighbors,FE_NipNeighbors(t)) + mesh_maxNsubNodes = max(mesh_maxNsubNodes,FE_NsubNodes(t)) + endif + enddo + endif + enddo + +620 return + endsubroutine @@ -1267,34 +1893,33 @@ subroutine mesh_get_nodeElemDimensions (unit) ! allocate globals: ! _node !******************************************************************** - subroutine mesh_build_nodes (unit) + subroutine mesh_marc_build_nodes (unit) use prec, only: pInt use IO implicit none - integer unit,i,j,m - integer(pInt), dimension(3) :: pos integer(pInt), dimension(5), parameter :: node_ends = (/0,10,30,50,70/) - character*300 line + integer(pInt), parameter :: maxNchunks = 1 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line - allocate ( mesh_node (3,mesh_Nnodes) ) - mesh_node(:,:) = 0_pInt + integer(pInt) unit,i,j,m + + allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0_pInt 610 FORMAT(A300) rewind(unit) do read (unit,610,END=670) line - pos = IO_stringPos(line,1) + pos = IO_stringPos(line,maxNchunks) if( IO_lc(IO_stringValue(line,pos,1)) == 'coordinates' ) then - read (unit,610,END=670) line ! skip crap line + read (unit,610,END=670) line ! skip crap line do i=1,mesh_Nnodes read (unit,610,END=670) line - m = mesh_FEasCP('node',IO_fixedIntValue (line,node_ends,1)) - do j=1,3 - mesh_node(j,m) = IO_fixedNoEFloatValue (line,node_ends,j+1) - enddo + m = mesh_FEasCP('node',IO_fixedIntValue(line,node_ends,1)) + forall (j = 1:3) mesh_node(j,m) = IO_fixedNoEFloatValue (line,node_ends,j+1) enddo exit endif @@ -1305,97 +1930,258 @@ subroutine mesh_get_nodeElemDimensions (unit) endsubroutine +!******************************************************************** +! store x,y,z coordinates of all nodes in mesh +! +! allocate globals: +! _node +!******************************************************************** + subroutine mesh_abaqus_build_nodes (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 4 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line + + integer(pInt) unit,i,j,m,count + logical inPart + + allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0_pInt + +610 FORMAT(A300) + + inPart = .false. + rewind(unit) + do + read (unit,610,END=670) line + pos = IO_stringPos(line,maxNchunks) + if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & + IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. + + if( inPart .and. & + IO_lc(IO_stringValue(line,pos,1)) == '*node' .and. & + ( IO_lc(IO_stringValue(line,pos,2)) /= 'output' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'print' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'file' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'response' ) & + ) then + count = IO_countDataLines(unit) ! how many nodes are defined here? + do i = 1,count + backspace(unit) ! rewind to first entry + enddo + do i = 1,count + read (unit,610,END=670) line + pos = IO_stringPos(line,maxNchunks) + m = mesh_FEasCP('node',IO_intValue(line,pos,1)) + forall (j=1:3) mesh_node(j,m) = IO_floatValue(line,pos,j+1) + enddo + endif + enddo + +670 return + + endsubroutine + + !******************************************************************** ! store FEid, type, mat, tex, and node list per element ! ! allocate globals: ! _element !******************************************************************** - subroutine mesh_build_elements (unit) + subroutine mesh_marc_build_elements (unit) use prec, only: pInt use IO implicit none - integer unit,i,j,sv,val,e + integer(pInt), parameter :: maxNchunks = 66 + integer(pInt), dimension (1+2*maxNchunks) :: pos + character(len=300) line - integer(pInt), parameter :: maxNchunks = 66 - integer(pInt), dimension(1+2*maxNchunks) :: pos integer(pInt), dimension(1+mesh_NcpElems) :: contInts - character*300 line + integer(pInt) unit,i,j,sv,val,e allocate (mesh_element (4+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt 610 FORMAT(A300) - rewind(unit) do read (unit,610,END=620) line - pos = IO_stringPos(line,2) + pos = IO_stringPos(line,1) if( IO_lc(IO_stringValue(line,pos,1)) == 'connectivity' ) then read (unit,610,END=620) line ! Garbage line - do i=1,mesh_Nelems - read (unit,610,END=620) line - pos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max (plus ID, type) - e = mesh_FEasCP('elem',IO_intValue(line,pos,1)) - if (e /= 0) then ! disregard non CP elems - mesh_element (1,e) = IO_IntValue (line,pos,1) ! FE id - mesh_element (2,e) = FE_mapElemtype(IO_StringValue (line,pos,2)) ! elem type - forall (j=1:FE_Nnodes(mesh_element(2,e))) & - mesh_element(j+4,e) = IO_IntValue (line,pos,j+2) ! copy FE ids of nodes - call IO_skipChunks(unit,FE_NoriginalNodes(mesh_element(2,e))-(pos(1)-2)) ! read on if FE_Nnodes exceeds node count present on current line - endif + do i = 1,mesh_Nelems + read (unit,610,END=620) line + pos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max (plus ID, type) + e = mesh_FEasCP('elem',IO_intValue(line,pos,1)) + if (e /= 0) then ! disregard non CP elems + mesh_element (1,e) = IO_IntValue (line,pos,1) ! FE id + mesh_element (2,e) = FE_mapElemtype(IO_StringValue(line,pos,2)) ! elem type + forall (j = 1:FE_Nnodes(mesh_element(2,e))) & + mesh_element(j+4,e) = IO_IntValue(line,pos,j+2) ! copy FE ids of nodes + call IO_skipChunks(unit,FE_NoriginalNodes(mesh_element(2,e))-(pos(1)-2)) ! read on if FE_Nnodes exceeds node count present on current line + endif enddo - exit endif enddo - rewind(unit) ! just in case "initial state" apears before "connectivity" +620 rewind(unit) ! just in case "initial state" apears before "connectivity" read (unit,610,END=620) line do pos = IO_stringPos(line,2) - if( (IO_lc(IO_stringValue(line,pos,1)) == 'initial').and. & + if( (IO_lc(IO_stringValue(line,pos,1)) == 'initial') .and. & (IO_lc(IO_stringValue(line,pos,2)) == 'state') ) then - if (initialcondTableStyle == 2) read (unit,610,END=620) line ! read extra line for new style - read (unit,610,END=620) line ! read line with index of state var + if (initialcondTableStyle == 2) read (unit,610,END=620) line ! read extra line for new style + read (unit,610,END=630) line ! read line with index of state var pos = IO_stringPos(line,1) - sv = IO_IntValue (line,pos,1) ! figure state variable index - if( (sv == 2).or.(sv == 3) ) then ! only state vars 2 and 3 of interest - read (unit,610,END=620) line ! read line with value of state var + sv = IO_IntValue(line,pos,1) ! figure state variable index + if( (sv == 2).or.(sv == 3) ) then ! only state vars 2 and 3 of interest + read (unit,610,END=620) line ! read line with value of state var pos = IO_stringPos(line,1) - do while (scan(IO_stringValue(line,pos,1),'+-',back=.true.)>1) ! is noEfloat value? - val = NINT(IO_fixedNoEFloatValue (line,(/0,20/),1)) ! state var's value - mesh_maxValStateVar(sv-1) = max(val,mesh_maxValStateVar(sv-1)) ! remember max val of material and texture index + do while (scan(IO_stringValue(line,pos,1),'+-',back=.true.)>1) ! is noEfloat value? + val = NINT(IO_fixedNoEFloatValue(line,(/0,20/),1)) ! state var's value + mesh_maxValStateVar(sv-1) = max(val,mesh_maxValStateVar(sv-1)) ! remember max val of homogenization and microstructure index if (initialcondTableStyle == 2) then - read (unit,610,END=620) line ! read extra line - read (unit,610,END=620) line ! read extra line + read (unit,610,END=630) line ! read extra line + read (unit,610,END=630) line ! read extra line endif contInts = IO_continousIntValues(unit,mesh_Nelems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) ! get affected elements do i = 1,contInts(1) e = mesh_FEasCP('elem',contInts(1+i)) mesh_element(1+sv,e) = val enddo - if (initialcondTableStyle == 0) read (unit,610,END=620) line ! ignore IP range for old table style - read (unit,610,END=620) line + if (initialcondTableStyle == 0) read (unit,610,END=620) line ! ignore IP range for old table style + read (unit,610,END=630) line pos = IO_stringPos(line,1) enddo endif else - read (unit,610,END=620) line + read (unit,610,END=630) line endif enddo -620 return +630 return endsubroutine - + + !******************************************************************** -! build list of elements shared by each node in mesh +! store FEid, type, mat, tex, and node list per element ! ! allocate globals: +! _element +!******************************************************************** + subroutine mesh_abaqus_build_elements (unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt), parameter :: maxNchunks = 65 + integer(pInt), dimension (1+2*maxNchunks) :: pos + + integer(pInt) unit,i,j,count,e,t,homog,micro + logical inPart,materialFound + character (len=64) materialName,elemSetName + character(len=300) line + + allocate (mesh_element (4+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt + +610 FORMAT(A300) + + inPart = .false. + rewind(unit) + do + read (unit,610,END=620) line + pos = IO_stringPos(line,2) + if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & + IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. + + if( inPart .and. & + IO_lc(IO_stringValue(line,pos,1)) == '*element' .and. & + ( IO_lc(IO_stringValue(line,pos,2)) /= 'output' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'matrix' .and. & + IO_lc(IO_stringValue(line,pos,2)) /= 'response' ) & + ) then + t = FE_mapElemtype(IO_extractValue(IO_stringValue(line,pos,2),'type')) ! remember elem type + count = IO_countDataLines(unit) + do i = 1,count + backspace(unit) + enddo + do i = 1,count + read (unit,610,END=620) line + pos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max + e = mesh_FEasCP('elem',IO_intValue(line,pos,1)) + if (e /= 0) then ! disregard non CP elems + mesh_element(1,e) = IO_intValue(line,pos,1) ! FE id + mesh_element(2,e) = t ! elem type + forall (j=1:FE_Nnodes(t)) & + mesh_element(4+j,e) = IO_intValue(line,pos,1+j) ! copy FE ids of nodes to position 5: + call IO_skipChunks(unit,FE_NoriginalNodes(t)-(pos(1)-1)) ! read on (even multiple lines) if FE_NoriginalNodes exceeds required node count + endif + enddo + endif + enddo + + +620 rewind(unit) ! just in case "*material" definitions apear before "*element" + + materialFound = .false. + do + read (unit,610,END=630) line + pos = IO_stringPos(line,maxNchunks) + select case ( IO_lc(IO_StringValue(line,pos,1))) + case('*material') + materialName = IO_extractValue(IO_lc(IO_StringValue(line,pos,2)),'name') ! extract name=value + materialFound = materialName /= '' ! valid name? + case('*user') + if ( IO_lc(IO_StringValue(line,pos,2)) == 'material' .and. & + materialFound ) then + do i = 1,mesh_Nmaterials ! look thru material names + if (materialName == mesh_nameMaterial(i)) exit ! found one + enddo + if (i <= mesh_Nmaterials) then ! found one? + elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet + read (unit,610,END=630) line ! read homogenization and microstructure + pos = IO_stringPos(line,2) + homog = NINT(IO_floatValue(line,pos,1)) + micro = NINT(IO_floatValue(line,pos,2)) + do i = 1,mesh_NelemSets ! look thru all elemSet definitions + if (elemSetName == mesh_nameElemSet(i)) then ! matched? + do j = 1,mesh_mapElemSet(1,i) + e = mesh_FEasCP('elem',mesh_mapElemSet(1+j,i)) + mesh_element(3,e) = homog ! store homogenization + mesh_element(4,e) = micro ! store microstructure + mesh_maxValStateVar(1) = max(mesh_maxValStateVar(1),homog) + mesh_maxValStateVar(2) = max(mesh_maxValStateVar(2),micro) + enddo + endif + enddo + endif + materialFound = .false. + endif + endselect + enddo + +630 return + + endsubroutine + + +!******************************************************************** +! get maximum count of shared elements among cpElements and +! build list of elements shared by each node in mesh +! +! _maxNsharedElems ! _sharedElem !******************************************************************** subroutine mesh_build_sharedElems () @@ -1404,20 +2190,38 @@ subroutine mesh_get_nodeElemDimensions (unit) use IO implicit none - integer(pint) e,n,j + integer(pint) e,t,n,j + integer(pInt), dimension (mesh_Nnodes) :: node_count integer(pInt), dimension (:), allocatable :: node_seen allocate(node_seen(maxval(FE_Nnodes))) + node_count = 0_pInt + + do e = 1,mesh_NcpElems + + t = mesh_element(2,e) ! my type + write(6,*) e,t; call flush(6) + node_seen = 0_pInt ! reset node duplicates + do j = 1,FE_Nnodes(t) ! check each node of element + n = mesh_FEasCP('node',mesh_element(4+j,e)) + if (all(node_seen /= n)) node_count(n) = node_count(n) + 1_pInt ! if FE node not yet encountered -> count it + node_seen(j) = n ! remember this node to be counted already + enddo + + enddo + + mesh_maxNsharedElems = maxval(node_count) ! most shared node + allocate ( mesh_sharedElem( 1+mesh_maxNsharedElems,mesh_Nnodes) ) - mesh_sharedElem(:,:) = 0_pInt + mesh_sharedElem = 0_pInt do e = 1,mesh_NcpElems node_seen = 0_pInt do j = 1,FE_Nnodes(mesh_element(2,e)) n = mesh_FEasCP('node',mesh_element(4+j,e)) if (all(node_seen /= n)) then - mesh_sharedElem(1,n) = mesh_sharedElem(1,n) + 1 - mesh_sharedElem(1+mesh_sharedElem(1,n),n) = e + mesh_sharedElem(1,n) = mesh_sharedElem(1,n) + 1 ! count for each node the connected elements + mesh_sharedElem(1+mesh_sharedElem(1,n),n) = e ! store the respective element id endif node_seen(j) = n enddo @@ -1448,36 +2252,37 @@ subroutine mesh_get_nodeElemDimensions (unit) allocate(mesh_ipNeighborhood(2,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipNeighborhood = 0_pInt - do e = 1,mesh_NcpElems ! loop over cpElems - t = mesh_element(2,e) ! get elemType - do i = 1,FE_Nips(t) ! loop over IPs of elem - do n = 1,FE_NipNeighbors(t) ! loop over neighbors of IP + do e = 1,mesh_NcpElems ! loop over cpElems + t = mesh_element(2,e) ! get elemType + do i = 1,FE_Nips(t) ! loop over IPs of elem + do n = 1,FE_NipNeighbors(t) ! loop over neighbors of IP neighbor = FE_ipNeighbor(n,i,t) - if (neighbor > 0) then ! intra-element IP + if (neighbor > 0) then ! intra-element IP neighboringElem = e neighboringIP = neighbor - else ! neighboring element's IP + else ! neighboring element's IP neighboringElem = 0_pInt neighboringIP = 0_pInt - matchingElem = mesh_faceMatch(-neighbor,e) ! get CP elem id of face match + matchingElem = mesh_faceMatch(-neighbor,e) ! get CP elem id of face match if (matchingElem > 0 .and. mesh_element(2,matchingElem) == t) then ! found match of same type? - if (FE_nodesAtIP(2,1,i,t) == 0) then ! single linked node - matchNode1: do j = 1,FE_Nnodes(t) ! check against all neighbor's nodes + if (FE_nodesAtIP(2,1,i,t) == 0) then ! single linked node + matchNode1: do j = 1,FE_Nnodes(t) ! check against all neighbor's nodes if (mesh_element(4+FE_nodesAtIP(1,1,i,t),e)==mesh_element(4+j,matchingElem)) then - linkedNode(1) = j ! which neighboring node matches my first nodeAtIP (indexed globally) + linkedNode(1) = j ! which neighboring node matches my first nodeAtIP (indexed globally) linkedNode(2) = 0_pInt exit matchNode1 endif enddo matchNode1 matchFace1: do j = 1,FE_Nips(t) - if ((linkedNode(1) == FE_nodesAtIP(1,1,j,t)) .and. (FE_nodesAtIP(2,1,j,t) == 0))then + if ( (linkedNode(1) == FE_nodesAtIP(1,1,j,t)) .and. & + (FE_nodesAtIP(2,1,j,t) == 0) ) then neighboringElem = matchingElem neighboringIP = j exit matchFace1 endif enddo matchFace1 - else ! double linked node - matchNode2: do j = 1,FE_Nnodes(t) ! check against all neighbor's nodes + else ! double linked node + matchNode2: do j = 1,FE_Nnodes(t) ! check against all neighbor's nodes if (mesh_element(4+FE_nodesAtIP(1,1,i,t),e)==mesh_element(4+j,matchingElem)) linkedNode(1) = j ! which neighboring node matches my first nodeAtIP (indexed globally) if (mesh_element(4+FE_nodesAtIP(2,1,i,t),e)==mesh_element(4+j,matchingElem)) linkedNode(2) = j ! which neighboring node matches my second nodeAtIP (indexed globally) enddo matchNode2 @@ -1518,7 +2323,7 @@ subroutine mesh_get_nodeElemDimensions (unit) implicit none integer(pInt) e,t,n,p - + allocate(mesh_subNodeCoord(3,mesh_maxNnodes+mesh_maxNsubNodes,mesh_NcpElems)) ; mesh_subNodeCoord = 0.0_pReal do e = 1,mesh_NcpElems ! loop over cpElems @@ -1558,38 +2363,38 @@ subroutine mesh_get_nodeElemDimensions (unit) integer(pInt), parameter :: Ntriangles = FE_NipFaceNodes-2 ! each interface is made up of this many triangles integer(pInt), dimension(mesh_maxNnodes+mesh_maxNsubNodes) :: gravityNode ! flagList to find subnodes determining center of grav real(pReal), dimension(3,mesh_maxNnodes+mesh_maxNsubNodes) :: gravityNodePos ! coordinates of subnodes determining center of grav - real(pReal), dimension (3,FE_NipFaceNodes) :: nPos ! coordinates of nodes on IP face + real(pReal), dimension(3,FE_NipFaceNodes) :: nPos ! coordinates of nodes on IP face real(pReal), dimension(Ntriangles,FE_NipFaceNodes) :: volume ! volumes of possible tetrahedra real(pReal), dimension(3) :: centerOfGravity allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems)) ; mesh_ipVolume = 0.0_pReal allocate(mesh_ipCenterOfGravity(3,mesh_maxNips,mesh_NcpElems)) ; mesh_ipCenterOfGravity = 0.0_pReal - do e = 1,mesh_NcpElems ! loop over cpElems - t = mesh_element(2,e) ! get elemType - do i = 1,FE_Nips(t) ! loop over IPs of elem - gravityNode = 0_pInt ! reset flagList - gravityNodePos = 0.0_pReal ! reset coordinates - do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP - do n = 1,FE_NipFaceNodes ! loop over nodes on interface + do e = 1,mesh_NcpElems ! loop over cpElems + t = mesh_element(2,e) ! get elemType + do i = 1,FE_Nips(t) ! loop over IPs of elem + gravityNode = 0_pInt ! reset flagList + gravityNodePos = 0.0_pReal ! reset coordinates + do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP + do n = 1,FE_NipFaceNodes ! loop over nodes on interface gravityNode(FE_subNodeOnIPFace(n,f,i,t)) = 1 gravityNodePos(:,FE_subNodeOnIPFace(n,f,i,t)) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) enddo enddo - - do j = 1,mesh_maxNnodes+mesh_maxNsubNodes-1 ! walk through entire flagList except last - if (gravityNode(j) > 0_pInt) then ! valid node index - do k = j+1,mesh_maxNnodes+mesh_maxNsubNodes ! walk through remainder of list + + do j = 1,mesh_maxNnodes+mesh_maxNsubNodes-1 ! walk through entire flagList except last + if (gravityNode(j) > 0_pInt) then ! valid node index + do k = j+1,mesh_maxNnodes+mesh_maxNsubNodes ! walk through remainder of list if (all(abs(gravityNodePos(:,j) - gravityNodePos(:,k)) < 1.0e-100_pReal)) then ! found duplicate - gravityNode(j) = 0_pInt ! delete first instance + gravityNode(j) = 0_pInt ! delete first instance gravityNodePos(:,j) = 0.0_pReal - exit ! continue with next suspect + exit ! continue with next suspect endif enddo endif enddo centerOfGravity = sum(gravityNodePos,2)/count(gravityNode > 0) - + do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP and add tetrahedra which connect to CoG forall (n = 1:FE_NipFaceNodes) nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles) & ! start at each interface node and build valid triangles to cover interface @@ -1621,28 +2426,28 @@ subroutine mesh_get_nodeElemDimensions (unit) implicit none integer(pInt) e,f,t,i,j,n - integer(pInt), parameter :: Ntriangles = FE_NipFaceNodes-2 ! each interface is made up of this many triangles - real(pReal), dimension (3,FE_NipFaceNodes) :: nPos ! coordinates of nodes on IP face + integer(pInt), parameter :: Ntriangles = FE_NipFaceNodes-2 ! each interface is made up of this many triangles + real(pReal), dimension (3,FE_NipFaceNodes) :: nPos ! coordinates of nodes on IP face real(pReal), dimension(3,Ntriangles,FE_NipFaceNodes) :: normal real(pReal), dimension(Ntriangles,FE_NipFaceNodes) :: area allocate(mesh_ipArea(mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipArea = 0.0_pReal allocate(mesh_ipAreaNormal(3,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipAreaNormal = 0.0_pReal - do e = 1,mesh_NcpElems ! loop over cpElems - t = mesh_element(2,e) ! get elemType - do i = 1,FE_Nips(t) ! loop over IPs of elem - do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP + do e = 1,mesh_NcpElems ! loop over cpElems + t = mesh_element(2,e) ! get elemType + do i = 1,FE_Nips(t) ! loop over IPs of elem + do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP forall (n = 1:FE_NipFaceNodes) nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) - forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles) ! start at each interface node and build valid triangles to cover interface + forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles) ! start at each interface node and build valid triangles to cover interface normal(:,j,n) = math_vectorproduct(nPos(:,1+mod(n+j-1,FE_NipFaceNodes)) - nPos(:,n), & ! calc their normal vectors nPos(:,1+mod(n+j-0,FE_NipFaceNodes)) - nPos(:,n)) area(j,n) = dsqrt(sum(normal(:,j,n)*normal(:,j,n))) ! and area end forall forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles, area(j,n) > 0.0_pReal) & - normal(:,j,n) = normal(:,j,n) / area(j,n) ! make unit normal + normal(:,j,n) = normal(:,j,n) / area(j,n) ! make unit normal - mesh_ipArea(f,i,e) = sum(area) / (FE_NipFaceNodes*2.0_pReal) ! area of parallelograms instead of triangles + mesh_ipArea(f,i,e) = sum(area) / (FE_NipFaceNodes*2.0_pReal) ! area of parallelograms instead of triangles mesh_ipAreaNormal(:,f,i,e) = sum(sum(normal,3),2) / count(area > 0.0_pReal) ! average of all valid normals enddo enddo @@ -1748,13 +2553,13 @@ subroutine mesh_get_nodeElemDimensions (unit) write (6,fmt) i,"| ",mesh_HomogMicro(i,:) ! loop over all (possibly assigned) microstrcutures enddo write (6,*) + call flush(6) !$OMP END CRITICAL (write2out) - return endsubroutine - + END MODULE mesh diff --git a/code/mpie_cpfem_abaqus.f90 b/code/mpie_cpfem_abaqus.f90 new file mode 100644 index 000000000..ab6837805 --- /dev/null +++ b/code/mpie_cpfem_abaqus.f90 @@ -0,0 +1,270 @@ +!* $Id$ +!******************************************************************** +! Material subroutine for Abaqus +! +! written by P. Eisenlohr, +! F. Roters, +! K. Janssens +! +! MPI fuer Eisenforschung, Duesseldorf +! PSI, Switzerland +!******************************************************************** + +MODULE cpfem_interface + +character(len=64), parameter :: FEsolver = 'Abaqus' + +CONTAINS + +subroutine mpie_cpfem_init () + +!$OMP CRITICAL (write2out) + write(6,*) + write(6,*) '<<<+- mpie_cpfem_abaqus init -+>>>' + write(6,*) + call flush(6) +!$OMP END CRITICAL (write2out) + return +end subroutine + +END MODULE + + include "prec.f90" ! uses nothing else + include "IO.f90" ! uses prec + include "numerics.f90" ! uses prec, IO + include "math.f90" ! uses prec, numerics + include "debug.f90" ! uses prec, numerics + include "FEsolving.f90" ! uses prec, IO + include "mesh.f90" ! uses prec, math, IO, FEsolving + include "material.f90" ! uses prec, math, IO, mesh + include "lattice.f90" ! uses prec, math, IO, material + include "constitutive_phenopowerlaw.f90" ! uses prec, math, IO, latt ice, material, debug + include "constitutive_j2.f90" ! uses prec, math, IO, latt ice, material, debug + include "constitutive_dislobased.f90" ! uses prec, math, IO, latt ice, material, debug + include "constitutive_nonlocal.f90" ! uses prec, math, IO, latt ice, material, debug + include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug + include "crystallite.f90" ! uses prec, math, IO, numerics + include "homogenization_isostrain.f90" ! uses prec, math, IO, + include "homogenization_RGC.f90" ! uses prec, math, IO, numerics, mesh: added <<>> + include "homogenization.f90" ! uses prec, math, IO, numerics + include "CPFEM.f90" ! uses prec, math, IO, numerics, debug, FEsolving, mesh, lattice, constitutive, crystallite + +subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,& + RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,& + TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,NDI,NSHR,NTENS,& + NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,& + DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC) + + use prec, only: pReal, & + pInt + use FEsolving, only: cycleCounter, & + theInc, & + cutBack, & + calcMode, & + lastMode, & + theTime, & + theDelta, & + lastIncConverged, & + outdatedByNewInc, & + outdatedFFN1, & + terminallyIll, & + symmetricSolver + use math, only: invnrmMandel + use debug, only: debug_info, & + debug_reset + use mesh, only: mesh_FEasCP + use CPFEM, only: CPFEM_general,CPFEM_init_done + use homogenization, only: materialpoint_sizeResults, materialpoint_results + + + implicit none + + CHARACTER*80 CMNAME + integer(pInt) ndi, nshr, ntens, nstatv, nprops, noel, npt,& + kslay, kspt, kstep, kinc + real(pReal) STRESS(NTENS),STATEV(NSTATV),& + DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),& + STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),& + PROPS(NPROPS),COORDS(3),DROT(3,3),& + DFGRD0(3,3),DFGRD1(3,3) + real(pReal) SSE, SPD, SCD, RPL, DRPLDT, DTIME, TEMP,& + DTEMP, PNEWDT, CELENT + +! local variables + real(pReal), dimension(6) :: stress_h + real(pReal), dimension(6,6) :: ddsdde_h + integer(pInt) computationMode, i, cp_en + + if (noel == 1 .and. npt == 1) then +!$OMP CRITICAL (write2out) + write(6,*) 'el',noel,'ip',npt + write(6,*) 'got kinc as',kinc + write(6,*) 'got dStran',dstran + call flush(6) +!$OMP END CRITICAL (write2out) + endif + + if ( .not. CPFEM_init_done ) then + + computationMode = 2 ! calc + init +!$OMP CRITICAL (write2out) + write(6,'(i6,x,i2,x,a)') noel,npt,'first call special case..!'; call flush(6) +!$OMP END CRITICAL (write2out) + + else + cp_en = mesh_FEasCP('elem',noel) + if (theTime < time(2) .or. theInc /= kinc) then ! reached convergence + lastIncConverged = .true. + outdatedByNewInc = .true. + terminallyIll = .false. + cycleCounter = 0 +!$OMP CRITICAL (write2out) + write (6,'(i6,x,i2,x,a)') noel,npt,'lastIncConverged + outdated'; call flush(6) +!$OMP END CRITICAL (write2out) + endif + + if ( dtime < theDelta ) then ! cutBack + calcMode = .true. ! pretend last step was calculation + cutBack = .true. + terminallyIll = .false. + cycleCounter = 0 +!$OMP CRITICAL (write2out) + write(6,'(i6,x,i2,x,a)') noel,npt,'cutback detected..!'; call flush(6) +!$OMP END CRITICAL (write2out) + endif + + calcMode(npt,cp_en) = .not. calcMode(npt,cp_en) ! ping pong (calc <--> collect) + + if ( calcMode(npt,cp_en) ) then ! now calc + if ( lastMode .ne. calcMode(npt,cp_en) ) then ! first after ping pong + call debug_reset() ! resets debugging + outdatedFFN1 = .false. + cycleCounter = cycleCounter + 1 + endif + if ( outdatedByNewInc ) then + outdatedByNewInc = .false. + computationMode = 1 ! calc and age results + else + computationMode = 2 ! plain calc + endif + else ! now collect + if ( lastMode .ne. calcMode(npt,cp_en) ) call debug_info() ! first after ping pong reports debugging + if ( lastIncConverged ) then + lastIncConverged = .false. + computationMode = 4 ! collect and backup Jacobian after convergence + elseif ( cutBack ) then + cutBack = .false. + computationMode = 5 ! collect and restore Jacobian after cutback + else + computationMode = 3 ! plain collect + endif + endif + + endif + + theTime = time(2) ! record current starting time + theDelta = dtime ! record current time increment + theInc = kinc ! record current increment number + if (CPFEM_init_done) lastMode = calcMode(npt,cp_en) ! record calculationMode + +!$OMP CRITICAL (write2out) + write(6,'(a16,x,i2,x,a,i5,a,i5,x,i5,a)') 'computationMode',computationMode,'(',cp_en,':',noel,npt,')'; call flush(6) +!$OMP END CRITICAL (write2out) + + call CPFEM_general(computationMode,dfgrd0,dfgrd1,temp,dtime,noel,npt,stress,ddsdde,ntens) + +! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 +! straight: 11, 22, 33, 12, 23, 13 + forall(i=1:ntens) ddsdde(1:ntens,i) = invnrmMandel(i)*ddsdde(1:ntens,i)*invnrmMandel(1:ntens) + stress(1:ntens) = stress(1:ntens)*invnrmMandel(1:ntens) + if(symmetricSolver) ddsdde(1:ntens,1:ntens) = 0.5_pReal*(ddsdde(1:ntens,1:ntens) + transpose(ddsdde(1:ntens,1:ntens))) +! ABAQUS: 11, 22, 33, 12, 13, 23 + if(ntens == 6) then + stress_h=stress + stress(5)=stress_h(6) + stress(6)=stress_h(5) + ddsdde_h=ddsdde + ddsdde(:,5)=ddsdde_h(:,6) + ddsdde(:,6)=ddsdde_h(:,5) + ddsdde_h=ddsdde + ddsdde(5,:)=ddsdde_h(6,:) + ddsdde(6,:)=ddsdde_h(5,:) + end if + + statev = materialpoint_results(1:min(nstatv,materialpoint_sizeResults),npt,mesh_FEasCP('elem', noel)) + + if ( terminallyIll ) pnewdt = 0.5_pReal ! force cutback directly ? + + return + end subroutine + +!******************************************************************** +! This subroutine replaces the corresponding Marc subroutine +!******************************************************************** + subroutine quit(mpie_error) + + use prec, only: pReal, & + pInt + implicit none + integer(pInt) mpie_error + + call xit + end subroutine + +!############################################################################ +! + +! include "KJ_Disp.f" + subroutine disp(u,kstep,kinc,time,node,noel,jdof,coords) + +! hardwired aba_param.inc + implicit real*8(a-h,o-z) + parameter (nprecd=2) +! + dimension u(3),time(2), coords(3) + + real ktime,ktcl,ktmax,ktmin,ktdeltaup,ktdeltadown + real klmin,klmax,kldelta,klbegincycle + real kdeltau, kru, kret +! When using stdb_abqerr for debugging +! dimension intv(2),realv(4) +! character*8 charv(1) + +! ratchet displacement per cycle + kru = 0.1 +! ratcheting ends at time kret + kret = 400. +! displacement amplitude + kdeltau = 0.5 +! time cycle length: + ktcl = 4. + ktmax = ktcl/4 + ktmin = 3.*ktmax + ktdeltadown = ktmin - ktmax + ktdeltaup = ktcl - ktdeltadown +! load minimum & maximum: + klmin = -kdeltau + klmax = kdeltau + kldelta = klmax - klmin + klbegincycle = klmin + kldelta * (ktcl-ktmin) / ktdeltaup +! load as a function of (total time); trianglar loading cycle + ktime = time(2) + kru = kru * MIN(kret, ktime) / ktcl + if ( ktime .lt. ktmax ) then +! special case for path to first maximum + u(1) = kru + (ktime/ktmax) * klmax + else + do while ( ktime .ge. ktcl ) + ktime = ktime - ktcl + end do + if ( ktime .le. ktmax ) then + u(1) = kru + klbegincycle + ktime * (klmax-klbegincycle) / ktmax + else if ( ktime .lt. ktmin ) then + u(1) = kru + klmax - (ktime-ktmax) * kldelta / (ktmin-ktmax) + else + u(1) = kru + klmin + (ktime-ktmin) * (klbegincycle-klmin)/(ktcl-ktmin) + end if + endif + return + end +! diff --git a/code/mpie_cpfem_marc.f90 b/code/mpie_cpfem_marc.f90 index 2b3b59510..9c9a23485 100644 --- a/code/mpie_cpfem_marc.f90 +++ b/code/mpie_cpfem_marc.f90 @@ -38,6 +38,8 @@ ! MODULE cpfem_interface +character(len=64), parameter :: FEsolver = 'Marc' + CONTAINS subroutine mpie_cpfem_init @@ -61,7 +63,7 @@ END MODULE include "lattice.f90" ! uses prec, math, IO, material include "constitutive_phenopowerlaw.f90" ! uses prec, math, IO, lattice, material, debug include "constitutive_j2.f90" ! uses prec, math, IO, lattice, material, debug - include "constitutive_dislotwin.f90" ! uses prec, math, IO, lattice, material, debug + include "constitutive_dislobased.f90" ! uses prec, math, IO, lattice, material, debug include "constitutive_nonlocal.f90" ! uses prec, math, IO, lattice, material, debug include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug include "crystallite.f90" ! uses prec, math, IO, numerics @@ -148,18 +150,22 @@ subroutine hypela2(& pInt use FEsolving, only: cycleCounter, & theInc, & - theCycle, & - theLovl, & + cutBack, & + calcMode, & + lastMode, & theTime, & + theDelta, & lastIncConverged, & outdatedByNewInc, & outdatedFFN1, & terminallyIll, & symmetricSolver - use CPFEM, only: CPFEM_general use math, only: invnrmMandel use debug, only: debug_info, & debug_reset + use mesh, only: mesh_FEasCP + use CPFEM, only: CPFEM_general,CPFEM_init_done + use homogenization, only: materialpoint_sizeResults, materialpoint_results implicit none ! ** Start of generated type statements ** @@ -179,47 +185,74 @@ subroutine hypela2(& include "concom%%MARCVERSION%%" ! concom is needed for inc, subinc, ncycle, lovl include "creeps%%MARCVERSION%%" ! creeps is needed for timinc (time increment) - integer(pInt) computationMode, i + integer(pInt) computationMode, i, cp_en - if (inc == 0) then - cycleCounter = 4 + if ( .not. CPFEM_init_done ) then + + computationMode = 2 ! calc + init +!$OMP CRITICAL (write2out) + write (6,'(i6,x,i2,x,a)') n(1),nn,'first call special case..!'; call flush(6) +!$OMP END CRITICAL (write2out) + + else if (lovl == 4) then ! Marc requires stiffness in separate call + computationMode = 6 ! --> just return known value else - if (theCycle > ncycle .or. theInc /= inc) then - cycleCounter = 0 ! reset counter for each cutback or new inc + cp_en = mesh_FEasCP('elem',n(1)) + if (theTime < cptim .or. theInc /= inc) then ! reached convergence + lastIncConverged = .true. + outdatedByNewInc = .true. terminallyIll = .false. + cycleCounter = 0 +!$OMP CRITICAL (write2out) + write (6,'(i6,x,i2,x,a)') n(1),nn,'lastIncConverged + outdated'; call flush(6) +!$OMP END CRITICAL (write2out) endif - 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 - call debug_reset() + + if ( timinc < theDelta ) then ! cutBack + calcMode = .true. ! pretend last step was calculation + cutBack = .true. + terminallyIll = .false. + cycleCounter = 0 +!$OMP CRITICAL (write2out) + write(6,'(i6,x,i2,x,a)') n(1),nn,'cutback detected..!'; call flush(6) +!$OMP END CRITICAL (write2out) endif - endif - if (cptim > theTime .or. theInc /= inc) then ! reached convergence - lastIncConverged = .true. - outdatedByNewInc = .true. - write (6,*) n(1),nn,'lastIncConverged + outdated' + + calcMode(nn,cp_en) = .not. calcMode(nn,cp_en) ! ping pong (calc <--> collect) + + if ( calcMode(nn,cp_en) ) then ! now calc + if ( lastMode .ne. calcMode(nn,cp_en) ) then ! first after ping pong + call debug_reset() ! resets debugging + outdatedFFN1 = .false. + cycleCounter = cycleCounter + 1 + endif + if ( outdatedByNewInc ) then + outdatedByNewInc = .false. + computationMode = 1 ! calc and age results + else + computationMode = 2 ! plain calc + endif + else ! now collect + if ( lastMode .ne. calcMode(nn,cp_en) ) call debug_info() ! first after ping pong reports debugging + if ( lastIncConverged ) then + lastIncConverged = .false. + computationMode = 4 ! collect and backup Jacobian after convergence + elseif ( cutBack ) then + cutBack = .false. + computationMode = 5 ! collect and restore Jacobian after cutback + else + computationMode = 3 ! plain collect + endif + endif + endif - 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 - computationMode = 5 ! recycle and record former consistent tangent - lastIncConverged = .false. + if (lovl /= 4) then ! pass-by for Marc recycle + theTime = cptim ! record current starting time + theDelta = timinc ! record current time increment + theInc = inc ! record current increment number + if (CPFEM_init_done) lastMode = calcMode(nn,cp_en) ! record calculationMode endif - if (computationMode == 2 .and. outdatedByNewInc) then - computationMode = 1 ! compute and age former results - outdatedByNewInc = .false. - endif - - theTime = cptim ! record current starting time - theInc = inc ! record current increment number - theCycle = ncycle ! record current cycle count - theLovl = lovl ! record current lovl call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,s,d,ngens) @@ -269,55 +302,3 @@ subroutine plotv(& return end subroutine - - - -! subroutine utimestep(timestep,timestepold,icall,time,timeloadcase) -!******************************************************************** -! This routine modifies the addaptive time step of Marc -!******************************************************************** -! use prec, only: pReal,pInt -! use CPFEM, only : CPFEM_timefactor_max -! implicit none -! -! real(pReal) timestep, timestepold, time,timeloadcase -! integer(pInt) icall -! -! user subroutine for modifying the time step in auto step -! -! timestep : the current time step as suggested by marc -! to be modified in this routine -! timestepold : the current time step before it was modified by marc -! icall : =1 for setting the initial time step -! =2 if this routine is called during an increment -! =3 if this routine is called at the beginning -! of the increment -! time : time at the start of the current increment -! timeloadcase: time period of the current load case -! -! it is in general not recommended to increase the time step -! during the increment. -! this routine is called right after the time step has (possibly) -! been updated by marc. -! -! user coding -! reduce timestep during increment in case mpie_timefactor is too large -! if(icall==2_pInt) then -! if(mpie_timefactor_max>1.25_pReal) then -! timestep=min(timestep,timestepold*0.8_pReal) -! end if -! return -! modify timestep at beginning of new increment -! else if(icall==3_pInt) then -! if(mpie_timefactor_max<=0.8_pReal) then -! timestep=min(timestep,timestepold*1.25_pReal) -! else if (mpie_timefactor_max<=1.0_pReal) then -! timestep=min(timestep,timestepold/mpie_timefactor_max) -! else if (mpie_timefactor_max<=1.25_pReal) then -! timestep=min(timestep,timestepold*1.01_pReal) -! else -! timestep=min(timestep,timestepold*0.8_pReal) -! end if -! end if -! return -! end \ No newline at end of file