diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 06c8bd44c..919680d6e 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -6,7 +6,6 @@ !-------------------------------------------------------------------------------------------------- module homogenization use prec, only: & - pInt, & pReal !-------------------------------------------------------------------------------------------------- @@ -21,7 +20,7 @@ module homogenization materialpoint_dPdF !< tangent of first P--K stress at IP real(pReal), dimension(:,:,:), allocatable, public :: & materialpoint_results !< results array of material point - integer(pInt), public, protected :: & + integer, public, protected :: & materialpoint_sizeResults, & homogenization_maxSizePostResults, & thermal_maxSizePostResults, & @@ -92,10 +91,10 @@ subroutine homogenization_init worldrank implicit none - integer(pInt), parameter :: FILEUNIT = 200_pInt - integer(pInt) :: e,i,p - integer(pInt), dimension(:,:), pointer :: thisSize - integer(pInt), dimension(:) , pointer :: thisNoutput + integer, parameter :: FILEUNIT = 200 + integer :: e,i,p + integer, dimension(:,:), pointer :: thisSize + integer, dimension(:) , pointer :: thisNoutput character(len=64), dimension(:,:), pointer :: thisOutput character(len=32) :: outputName !< name of output, intermediate fix until HDF5 output is ready logical :: valid @@ -232,9 +231,9 @@ subroutine homogenization_init !-------------------------------------------------------------------------------------------------- ! allocate and initialize global state and postresutls variables - homogenization_maxSizePostResults = 0_pInt - thermal_maxSizePostResults = 0_pInt - damage_maxSizePostResults = 0_pInt + homogenization_maxSizePostResults = 0 + thermal_maxSizePostResults = 0 + damage_maxSizePostResults = 0 do p = 1,size(config_homogenization) homogenization_maxSizePostResults = max(homogenization_maxSizePostResults,homogState (p)%sizePostResults) thermal_maxSizePostResults = max(thermal_maxSizePostResults, thermalState (p)%sizePostResults) @@ -252,7 +251,7 @@ subroutine homogenization_init write(6,'(/,a)') ' <<<+- homogenization init -+>>>' - if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt) then + if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0) then #ifdef TODO write(6,'(a32,1x,7(i8,1x))') 'homogenization_state0: ', shape(homogenization_state0) write(6,'(a32,1x,7(i8,1x))') 'homogenization_subState0: ', shape(homogenization_subState0) @@ -275,7 +274,7 @@ subroutine homogenization_init flush(6) if (debug_g < 1 .or. debug_g > homogenization_Ngrains(mesh_element(3,debug_e))) & - call IO_error(602_pInt,ext_msg='constituent', el=debug_e, g=debug_g) + call IO_error(602,ext_msg='constituent', el=debug_e, g=debug_g) end subroutine homogenization_init @@ -344,7 +343,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) implicit none real(pReal), intent(in) :: dt !< time increment logical, intent(in) :: updateJaco !< initiating Jacobian update - integer(pInt) :: & + integer :: & NiterationHomog, & NiterationMPstate, & g, & !< grain number @@ -354,7 +353,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) myNgrains #ifdef DEBUG - if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt) then + if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0) then write(6,'(/a,i5,1x,i2)') '<< HOMOG >> Material Point start at el ip ', debug_e, debug_i write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< HOMOG >> F0', & @@ -372,7 +371,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) plasticState (phaseAt(g,i,e))%partionedState0(:,phasememberAt(g,i,e)) = & plasticState (phaseAt(g,i,e))%state0( :,phasememberAt(g,i,e)) - do mySource = 1_pInt, phase_Nsources(phaseAt(g,i,e)) + do mySource = 1, phase_Nsources(phaseAt(g,i,e)) sourceState(phaseAt(g,i,e))%p(mySource)%partionedState0(:,phasememberAt(g,i,e)) = & sourceState(phaseAt(g,i,e))%p(mySource)%state0( :,phasememberAt(g,i,e)) enddo @@ -393,19 +392,19 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) materialpoint_requested(i,e) = .true. ! everybody requires calculation endforall forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - homogState(material_homogenizationAt(e))%sizeState > 0_pInt) & + homogState(material_homogenizationAt(e))%sizeState > 0) & homogState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) = & homogState(material_homogenizationAt(e))%State0( :,mappingHomogenization(1,i,e)) ! ...internal homogenization state forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - thermalState(material_homogenizationAt(e))%sizeState > 0_pInt) & + thermalState(material_homogenizationAt(e))%sizeState > 0) & thermalState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) = & thermalState(material_homogenizationAt(e))%State0( :,mappingHomogenization(1,i,e)) ! ...internal thermal state forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - damageState(material_homogenizationAt(e))%sizeState > 0_pInt) & + damageState(material_homogenizationAt(e))%sizeState > 0) & damageState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) = & damageState(material_homogenizationAt(e))%State0( :,mappingHomogenization(1,i,e)) ! ...internal damage state enddo - NiterationHomog = 0_pInt + NiterationHomog = 0 cutBackLooping: do while (.not. terminallyIll .and. & any(materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMinHomog)) @@ -417,9 +416,9 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) converged: if ( materialpoint_converged(i,e) ) then #ifdef DEBUG - if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt & + if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0 & .and. ((e == debug_e .and. i == debug_i) & - .or. .not. iand(debug_level(debug_homogenization),debug_levelSelective) /= 0_pInt)) then + .or. .not. iand(debug_level(debug_homogenization),debug_levelSelective) /= 0)) then write(6,'(a,1x,f12.8,1x,a,1x,f12.8,1x,a,i8,1x,i2/)') '<< HOMOG >> winding forward from', & materialpoint_subFrac(i,e), 'to current materialpoint_subFrac', & materialpoint_subFrac(i,e)+materialpoint_subStep(i,e),'in materialpoint_stressAndItsTangent at el ip',e,i @@ -456,29 +455,29 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) do g = 1,myNgrains plasticState (phaseAt(g,i,e))%partionedState0(:,phasememberAt(g,i,e)) = & plasticState (phaseAt(g,i,e))%state( :,phasememberAt(g,i,e)) - do mySource = 1_pInt, phase_Nsources(phaseAt(g,i,e)) + do mySource = 1, phase_Nsources(phaseAt(g,i,e)) sourceState(phaseAt(g,i,e))%p(mySource)%partionedState0(:,phasememberAt(g,i,e)) = & sourceState(phaseAt(g,i,e))%p(mySource)%state( :,phasememberAt(g,i,e)) enddo enddo forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - homogState(material_homogenizationAt(e))%sizeState > 0_pInt) & + homogState(material_homogenizationAt(e))%sizeState > 0) & homogState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) = & homogState(material_homogenizationAt(e))%State( :,mappingHomogenization(1,i,e)) ! ...internal homogenization state forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - thermalState(material_homogenizationAt(e))%sizeState > 0_pInt) & + thermalState(material_homogenizationAt(e))%sizeState > 0) & thermalState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) = & thermalState(material_homogenizationAt(e))%State( :,mappingHomogenization(1,i,e)) ! ...internal thermal state forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - damageState(material_homogenizationAt(e))%sizeState > 0_pInt) & + damageState(material_homogenizationAt(e))%sizeState > 0) & damageState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) = & damageState(material_homogenizationAt(e))%State( :,mappingHomogenization(1,i,e)) ! ...internal damage state materialpoint_subF0(1:3,1:3,i,e) = materialpoint_subF(1:3,1:3,i,e) ! ...def grad endif steppingNeeded else converged - if ( (myNgrains == 1_pInt .and. materialpoint_subStep(i,e) <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite + if ( (myNgrains == 1 .and. materialpoint_subStep(i,e) <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite subStepSizeHomog * materialpoint_subStep(i,e) <= subStepMinHomog ) then ! would require too small subStep ! cutback makes no sense !$OMP FLUSH(terminallyIll) @@ -494,9 +493,9 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) materialpoint_subStep(i,e) = subStepSizeHomog * materialpoint_subStep(i,e) ! crystallite had severe trouble, so do a significant cutback #ifdef DEBUG - if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt & + if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0 & .and. ((e == debug_e .and. i == debug_i) & - .or. .not. iand(debug_level(debug_homogenization), debug_levelSelective) /= 0_pInt)) then + .or. .not. iand(debug_level(debug_homogenization), debug_levelSelective) /= 0)) then write(6,'(a,1x,f12.8,a,i8,1x,i2/)') & '<< HOMOG >> cutback step in materialpoint_stressAndItsTangent with new materialpoint_subStep:',& materialpoint_subStep(i,e),' at el ip',e,i @@ -518,21 +517,21 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) do g = 1, myNgrains plasticState (phaseAt(g,i,e))%state( :,phasememberAt(g,i,e)) = & plasticState (phaseAt(g,i,e))%partionedState0(:,phasememberAt(g,i,e)) - do mySource = 1_pInt, phase_Nsources(phaseAt(g,i,e)) + do mySource = 1, phase_Nsources(phaseAt(g,i,e)) sourceState(phaseAt(g,i,e))%p(mySource)%state( :,phasememberAt(g,i,e)) = & sourceState(phaseAt(g,i,e))%p(mySource)%partionedState0(:,phasememberAt(g,i,e)) enddo enddo forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - homogState(material_homogenizationAt(e))%sizeState > 0_pInt) & + homogState(material_homogenizationAt(e))%sizeState > 0) & homogState(material_homogenizationAt(e))%State( :,mappingHomogenization(1,i,e)) = & homogState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) ! ...internal homogenization state forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - thermalState(material_homogenizationAt(e))%sizeState > 0_pInt) & + thermalState(material_homogenizationAt(e))%sizeState > 0) & thermalState(material_homogenizationAt(e))%State( :,mappingHomogenization(1,i,e)) = & thermalState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) ! ...internal thermal state forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - damageState(material_homogenizationAt(e))%sizeState > 0_pInt) & + damageState(material_homogenizationAt(e))%sizeState > 0) & damageState(material_homogenizationAt(e))%State( :,mappingHomogenization(1,i,e)) = & damageState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) ! ...internal damage state endif @@ -550,7 +549,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) enddo elementLooping1 !$OMP END PARALLEL DO - NiterationMPstate = 0_pInt + NiterationMPstate = 0 convergenceLooping: do while (.not. terminallyIll .and. & any( materialpoint_requested(:,FEsolving_execELem(1):FEsolving_execElem(2)) & @@ -606,7 +605,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) enddo convergenceLooping - NiterationHomog = NiterationHomog + 1_pInt + NiterationHomog = NiterationHomog + 1 enddo cutBackLooping @@ -652,7 +651,7 @@ subroutine materialpoint_postResults crystallite_postResults implicit none - integer(pInt) :: & + integer :: & thePos, & theSize, & myNgrains, & @@ -666,21 +665,21 @@ subroutine materialpoint_postResults myNgrains = homogenization_Ngrains(mesh_element(3,e)) myCrystallite = microstructure_crystallite(mesh_element(4,e)) IpLooping: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - thePos = 0_pInt + thePos = 0 theSize = homogState (material_homogenizationAt(e))%sizePostResults & + thermalState (material_homogenizationAt(e))%sizePostResults & + damageState (material_homogenizationAt(e))%sizePostResults materialpoint_results(thePos+1,i,e) = real(theSize,pReal) ! tell size of homogenization results - thePos = thePos + 1_pInt + thePos = thePos + 1 - if (theSize > 0_pInt) then ! any homogenization results to mention? + if (theSize > 0) then ! any homogenization results to mention? materialpoint_results(thePos+1:thePos+theSize,i,e) = postResults(i,e) ! tell homogenization results thePos = thePos + theSize endif materialpoint_results(thePos+1,i,e) = real(myNgrains,pReal) ! tell number of grains at materialpoint - thePos = thePos + 1_pInt + thePos = thePos + 1 grainLooping :do g = 1,myNgrains theSize = 1 + crystallite_sizePostResults(myCrystallite) + & @@ -716,7 +715,7 @@ subroutine partitionDeformation(ip,el) homogenization_RGC_partitionDeformation implicit none - integer(pInt), intent(in) :: & + integer, intent(in) :: & ip, & !< integration point el !< element number @@ -769,7 +768,7 @@ function updateState(ip,el) damage_local_updateState implicit none - integer(pInt), intent(in) :: & + integer, intent(in) :: & ip, & !< integration point el !< element number logical, dimension(2) :: updateState @@ -831,7 +830,7 @@ subroutine averageStressAndItsTangent(ip,el) homogenization_RGC_averageStressAndItsTangent implicit none - integer(pInt), intent(in) :: & + integer, intent(in) :: & ip, & !< integration point el !< element number @@ -900,20 +899,20 @@ function postResults(ip,el) damage_nonlocal_postResults implicit none - integer(pInt), intent(in) :: & + integer, intent(in) :: & ip, & !< integration point el !< element number real(pReal), dimension( homogState (material_homogenizationAt(el))%sizePostResults & + thermalState (material_homogenizationAt(el))%sizePostResults & + damageState (material_homogenizationAt(el))%sizePostResults) :: & postResults - integer(pInt) :: & + integer :: & startPos, endPos ,& of, instance, homog postResults = 0.0_pReal - startPos = 1_pInt + startPos = 1 endPos = homogState(material_homogenizationAt(el))%sizePostResults chosenHomogenization: select case (homogenization_type(mesh_element(3,el))) @@ -924,7 +923,7 @@ function postResults(ip,el) end select chosenHomogenization - startPos = endPos + 1_pInt + startPos = endPos + 1 endPos = endPos + thermalState(material_homogenizationAt(el))%sizePostResults chosenThermal: select case (thermal_type(mesh_element(3,el))) @@ -939,7 +938,7 @@ function postResults(ip,el) end select chosenThermal - startPos = endPos + 1_pInt + startPos = endPos + 1 endPos = endPos + damageState(material_homogenizationAt(el))%sizePostResults chosenDamage: select case (damage_type(mesh_element(3,el))) diff --git a/src/homogenization_RGC.f90 b/src/homogenization_RGC.f90 index 6a513193b..26eadf4f9 100644 --- a/src/homogenization_RGC.f90 +++ b/src/homogenization_RGC.f90 @@ -13,9 +13,9 @@ module homogenization_RGC implicit none private - integer(pInt), dimension(:,:), allocatable,target, public :: & + integer, dimension(:,:), allocatable,target, public :: & homogenization_RGC_sizePostResult - character(len=64), dimension(:,:), allocatable,target, public :: & + character(len=64), dimension(:,:), allocatable,target, public :: & homogenization_RGC_output ! name of each post result output enum, bind(c) @@ -30,7 +30,7 @@ module homogenization_RGC end enum type, private :: tParameters - integer(pInt), dimension(:), allocatable :: & + integer, dimension(:), allocatable :: & Nconstituents real(pReal) :: & xiAlpha, & @@ -38,8 +38,8 @@ module homogenization_RGC real(pReal), dimension(:), allocatable :: & dAlpha, & angles - integer(pInt) :: & - of_debug = 0_pInt + integer :: & + of_debug = 0 integer(kind(undefined_ID)), dimension(:), allocatable :: & outputID end type tParameters @@ -121,7 +121,7 @@ subroutine homogenization_RGC_init() config_homogenization implicit none - integer(pInt) :: & + integer :: & Ninstance, & h, i, & NofMyHomog, outputSize, & @@ -152,11 +152,11 @@ subroutine homogenization_RGC_init() allocate(state0(Ninstance)) allocate(dependentState(Ninstance)) - allocate(homogenization_RGC_sizePostResult(maxval(homogenization_Noutput),Ninstance),source=0_pInt) + allocate(homogenization_RGC_sizePostResult(maxval(homogenization_Noutput),Ninstance),source=0) allocate(homogenization_RGC_output(maxval(homogenization_Noutput),Ninstance)) homogenization_RGC_output='' - do h = 1_pInt, size(homogenization_type) + do h = 1, size(homogenization_type) if (homogenization_type(h) /= HOMOGENIZATION_RGC_ID) cycle associate(prm => param(homogenization_typeInstance(h)), & stt => state(homogenization_typeInstance(h)), & @@ -172,7 +172,7 @@ subroutine homogenization_RGC_init() prm%Nconstituents = config%getInts('clustersize',requiredSize=3) if (homogenization_Ngrains(h) /= product(prm%Nconstituents)) & - call IO_error(211_pInt,ext_msg='clustersize ('//HOMOGENIZATION_RGC_label//')') + call IO_error(211,ext_msg='clustersize ('//HOMOGENIZATION_RGC_label//')') prm%xiAlpha = config%getFloat('scalingparameter') prm%ciAlpha = config%getFloat('overproportionality') @@ -183,28 +183,28 @@ subroutine homogenization_RGC_init() outputs = config%getStrings('(output)',defaultVal=emptyStringArray) allocate(prm%outputID(0)) - do i=1_pInt, size(outputs) + do i=1, size(outputs) outputID = undefined_ID select case(outputs(i)) case('constitutivework') outputID = constitutivework_ID - outputSize = 1_pInt + outputSize = 1 case('penaltyenergy') outputID = penaltyenergy_ID - outputSize = 1_pInt + outputSize = 1 case('volumediscrepancy') outputID = volumediscrepancy_ID - outputSize = 1_pInt + outputSize = 1 case('averagerelaxrate') outputID = averagerelaxrate_ID - outputSize = 1_pInt + outputSize = 1 case('maximumrelaxrate') outputID = maximumrelaxrate_ID - outputSize = 1_pInt + outputSize = 1 case('magnitudemismatch') outputID = magnitudemismatch_ID - outputSize = 3_pInt + outputSize = 3 end select @@ -217,9 +217,9 @@ subroutine homogenization_RGC_init() enddo NofMyHomog = count(material_homogenizationAt == h) - nIntFaceTot = 3_pInt*( (prm%Nconstituents(1)-1_pInt)*prm%Nconstituents(2)*prm%Nconstituents(3) & - + prm%Nconstituents(1)*(prm%Nconstituents(2)-1_pInt)*prm%Nconstituents(3) & - + prm%Nconstituents(1)*prm%Nconstituents(2)*(prm%Nconstituents(3)-1_pInt)) + nIntFaceTot = 3*( (prm%Nconstituents(1)-1)*prm%Nconstituents(2)*prm%Nconstituents(3) & + + prm%Nconstituents(1)*(prm%Nconstituents(2)-1)*prm%Nconstituents(3) & + + prm%Nconstituents(1)*prm%Nconstituents(2)*(prm%Nconstituents(3)-1)) sizeState = nIntFaceTot & + size(['avg constitutive work ','average penalty energy']) @@ -266,36 +266,36 @@ subroutine homogenization_RGC_partitionDeformation(F,avgF,instance,of) real(pReal), dimension (:,:,:), intent(out) :: F !< partioned F per grain real(pReal), dimension (:,:), intent(in) :: avgF !< averaged F - integer(pInt), intent(in) :: & + integer, intent(in) :: & instance, & of - real(pReal), dimension (3) :: aVect,nVect - integer(pInt), dimension (4) :: intFace - integer(pInt), dimension (3) :: iGrain3 - integer(pInt) :: iGrain,iFace,i,j + real(pReal), dimension(3) :: aVect,nVect + integer, dimension(4) :: intFace + integer, dimension(3) :: iGrain3 + integer :: iGrain,iFace,i,j associate(prm => param(instance)) !-------------------------------------------------------------------------------------------------- ! compute the deformation gradient of individual grains due to relaxations F = 0.0_pReal - do iGrain = 1_pInt,product(prm%Nconstituents) + do iGrain = 1,product(prm%Nconstituents) iGrain3 = grain1to3(iGrain,prm%Nconstituents) - do iFace = 1_pInt,6_pInt + do iFace = 1,6 intFace = getInterface(iFace,iGrain3) ! identifying 6 interfaces of each grain aVect = relaxationVector(intFace,instance,of) ! get the relaxation vectors for each interface from global relaxation vector array nVect = interfaceNormal(intFace,instance,of) - forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) & + forall (i=1:3,j=1:3) & F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! calculating deformation relaxations due to interface relaxation enddo F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF ! resulting relaxed deformation gradient #ifdef DEBUG - if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then + if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0) then write(6,'(1x,a32,1x,i3)')'Deformation gradient of grain: ',iGrain - do i = 1_pInt,3_pInt - write(6,'(1x,3(e15.8,1x))')(F(i,j,iGrain), j = 1_pInt,3_pInt) + do i = 1,3 + write(6,'(1x,3(e15.8,1x))')(F(i,j,iGrain), j = 1,3) enddo write(6,*)' ' flush(6) @@ -340,32 +340,32 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) implicit none - real(pReal), dimension (:,:,:), intent(in) :: & + real(pReal), dimension(:,:,:), intent(in) :: & P,& !< array of P F,& !< array of F F0 !< array of initial F - real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< array of current grain stiffness - real(pReal), dimension (3,3), intent(in) :: avgF !< average F - real(pReal), intent(in) :: dt !< time increment - integer(pInt), intent(in) :: & + real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< array of current grain stiffness + real(pReal), dimension(3,3), intent(in) :: avgF !< average F + real(pReal), intent(in) :: dt !< time increment + integer, intent(in) :: & ip, & !< integration point number el !< element number - logical, dimension(2) :: homogenization_RGC_updateState + logical, dimension(2) :: homogenization_RGC_updateState - integer(pInt), dimension (4) :: intFaceN,intFaceP,faceID - integer(pInt), dimension (3) :: nGDim,iGr3N,iGr3P - integer(pInt) :: instance,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ipert,iGrain,nGrain, of - real(pReal), dimension (3,3,size(P,3)) :: R,pF,pR,D,pD - real(pReal), dimension (3,size(P,3)) :: NN,devNull - real(pReal), dimension (3) :: normP,normN,mornP,mornN + integer, dimension(4) :: intFaceN,intFaceP,faceID + integer, dimension(3) :: nGDim,iGr3N,iGr3P + integer :: instance,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ipert,iGrain,nGrain, of + real(pReal), dimension(3,3,size(P,3)) :: R,pF,pR,D,pD + real(pReal), dimension(3,size(P,3)) :: NN,devNull + real(pReal), dimension(3) :: normP,normN,mornP,mornN real(pReal) :: residMax,stresMax logical :: error real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix,rmatrix real(pReal), dimension(:), allocatable :: resid,relax,p_relax,p_resid,drelax #ifdef DEBUG - integer(pInt), dimension (3) :: stresLoc - integer(pInt), dimension (2) :: residLoc + integer, dimension(3) :: stresLoc + integer, dimension(2) :: residLoc #endif zeroTimeStep: if(dEq0(dt)) then @@ -382,21 +382,21 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) ! get the dimension of the cluster (grains and interfaces) nGDim = prm%Nconstituents nGrain = product(nGDim) - nIntFaceTot = (nGDim(1)-1_pInt)*nGDim(2)*nGDim(3) & - + nGDim(1)*(nGDim(2)-1_pInt)*nGDim(3) & - + nGDim(1)*nGDim(2)*(nGDim(3)-1_pInt) + 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 global relaxation arrays/jacobian matrices depending on the size of the cluster - allocate(resid(3_pInt*nIntFaceTot), source=0.0_pReal) + allocate(resid(3*nIntFaceTot), source=0.0_pReal) allocate(tract(nIntFaceTot,3), source=0.0_pReal) relax = stt%relaxationVector(:,of) drelax = stt%relaxationVector(:,of) - st0%relaxationVector(:,of) #ifdef DEBUG - if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then + if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0) then write(6,'(1x,a30)')'Obtained state: ' - do i = 1_pInt,size(stt%relaxationVector(:,of)) + do i = 1,size(stt%relaxationVector(:,of)) write(6,'(1x,2(e15.8,1x))') stt%relaxationVector(i,of) enddo write(6,*)' ' @@ -412,15 +412,15 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) call volumePenalty(D,dst%volumeDiscrepancy(of),avgF,F,nGrain,instance,of) #ifdef DEBUG - if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then - do iGrain = 1_pInt,nGrain + if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0) then + do iGrain = 1,nGrain write(6,'(1x,a30,1x,i3,1x,a4,3(1x,e15.8))')'Mismatch magnitude of grain(',iGrain,') :',& NN(1,iGrain),NN(2,iGrain),NN(3,iGrain) write(6,'(/,1x,a30,1x,i3)')'Stress and penalties of grain: ',iGrain - do i = 1_pInt,3_pInt - write(6,'(1x,3(e15.8,1x),1x,3(e15.8,1x),1x,3(e15.8,1x))')(P(i,j,iGrain), j = 1_pInt,3_pInt), & - (R(i,j,iGrain), j = 1_pInt,3_pInt), & - (D(i,j,iGrain), j = 1_pInt,3_pInt) + do i = 1,3 + write(6,'(1x,3(e15.8,1x),1x,3(e15.8,1x),1x,3(e15.8,1x))')(P(i,j,iGrain), j = 1,3), & + (R(i,j,iGrain), j = 1,3), & + (D(i,j,iGrain), j = 1,3) enddo write(6,*)' ' enddo @@ -429,40 +429,40 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) !------------------------------------------------------------------------------------------------ ! computing the residual stress from the balance of traction at all (interior) interfaces - do iNum = 1_pInt,nIntFaceTot - faceID = interface1to4(iNum,param(instance)%Nconstituents) ! identifying the interface ID in local coordinate system (4-dimensional index) + do iNum = 1,nIntFaceTot + faceID = interface1to4(iNum,param(instance)%Nconstituents) ! identifying the interface ID in local coordinate system (4-dimensional index) !-------------------------------------------------------------------------------------------------- ! identify the left/bottom/back grain (-|N) iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate system (3-dimensional index) - iGrN = grain3to1(iGr3N,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index) - intFaceN = getInterface(2_pInt*faceID(1),iGr3N) + iGrN = grain3to1(iGr3N,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index) + intFaceN = getInterface(2*faceID(1),iGr3N) normN = interfaceNormal(intFaceN,instance,of) !-------------------------------------------------------------------------------------------------- ! identify the right/up/front grain (+|P) iGr3P = iGr3N - iGr3P(faceID(1)) = iGr3N(faceID(1))+1_pInt ! identifying the grain ID in local coordinate system (3-dimensional index) - iGrP = grain3to1(iGr3P,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index) - intFaceP = getInterface(2_pInt*faceID(1)-1_pInt,iGr3P) + iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate system (3-dimensional index) + iGrP = grain3to1(iGr3P,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index) + intFaceP = getInterface(2*faceID(1)-1,iGr3P) normP = interfaceNormal(intFaceP,instance,of) !-------------------------------------------------------------------------------------------------- ! compute the residual of traction at the interface (in local system, 4-dimensional index) - do i = 1_pInt,3_pInt - tract(iNum,i) = sign(viscModus_RGC*(abs(drelax(i+3*(iNum-1_pInt)))/(refRelaxRate_RGC*dt))**viscPower_RGC, & - drelax(i+3*(iNum-1_pInt))) ! contribution from the relaxation viscosity - do j = 1_pInt,3_pInt + do i = 1,3 + tract(iNum,i) = sign(viscModus_RGC*(abs(drelax(i+3*(iNum-1)))/(refRelaxRate_RGC*dt))**viscPower_RGC, & + drelax(i+3*(iNum-1))) ! contribution from the relaxation viscosity + do j = 1,3 tract(iNum,i) = tract(iNum,i) + (P(i,j,iGrP) + R(i,j,iGrP) + D(i,j,iGrP))*normP(j) & ! contribution from material stress P, mismatch penalty R, and volume penalty D projected into the interface + (P(i,j,iGrN) + R(i,j,iGrN) + D(i,j,iGrN))*normN(j) - resid(i+3_pInt*(iNum-1_pInt)) = tract(iNum,i) ! translate the local residual into global 1-dimensional residual array + resid(i+3*(iNum-1)) = tract(iNum,i) ! translate the local residual into global 1-dimensional residual array enddo enddo #ifdef DEBUG - if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then + if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0) then write(6,'(1x,a30,1x,i3)')'Traction at interface: ',iNum - write(6,'(1x,3(e15.8,1x))')(tract(iNum,j), j = 1_pInt,3_pInt) + write(6,'(1x,3(e15.8,1x))')(tract(iNum,j), j = 1,3) write(6,*)' ' endif #endif @@ -474,10 +474,10 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) residMax = maxval(abs(tract)) ! get the maximum of the residual #ifdef DEBUG - if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & + if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0 & .and. prm%of_debug == of) then - stresLoc = int(maxloc(abs(P)),pInt) ! get the location of the maximum stress - residLoc = int(maxloc(abs(tract)),pInt) ! get the position of the maximum residual + stresLoc = maxloc(abs(P)) + residLoc = maxloc(abs(tract)) write(6,'(1x,a)')' ' write(6,'(1x,a,1x,i2,1x,i4)')'RGC residual check ...',ip,el write(6,'(1x,a15,1x,e15.8,1x,a7,i3,1x,a12,i2,i2)')'Max stress: ',stresMax, & @@ -495,15 +495,15 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) if (residMax < relTol_RGC*stresMax .or. residMax < absTol_RGC) then homogenization_RGC_updateState = .true. #ifdef DEBUG - if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & + if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0 & .and. prm%of_debug == of) write(6,'(1x,a55,/)')'... done and happy' flush(6) #endif !-------------------------------------------------------------------------------------------------- ! compute/update the state for postResult, i.e., all energy densities computed by time-integration - do iGrain = 1_pInt,product(prm%Nconstituents) - do i = 1_pInt,3_pInt;do j = 1_pInt,3_pInt + do iGrain = 1,product(prm%Nconstituents) + do i = 1,3;do j = 1,3 stt%work(of) = stt%work(of) & + P(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal) stt%penaltyEnergy(of) = stt%penaltyEnergy(of) & @@ -512,11 +512,11 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) enddo dst%mismatch(1:3,of) = sum(NN,2)/real(nGrain,pReal) - dst%relaxationRate_avg(of) = sum(abs(drelax))/dt/real(3_pInt*nIntFaceTot,pReal) + dst%relaxationRate_avg(of) = sum(abs(drelax))/dt/real(3*nIntFaceTot,pReal) dst%relaxationRate_max(of) = maxval(abs(drelax))/dt #ifdef DEBUG - if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & + if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0 & .and. prm%of_debug == of) then write(6,'(1x,a30,1x,e15.8)') 'Constitutive work: ',stt%work(of) write(6,'(1x,a30,3(1x,e15.8))')'Magnitude mismatch: ',dst%mismatch(1,of), & @@ -538,7 +538,7 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) homogenization_RGC_updateState = [.true.,.false.] ! with direct cut-back #ifdef DEBUG - if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & + if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0 & .and. prm%of_debug == of) write(6,'(1x,a,/)') '... broken' flush(6) #endif @@ -547,7 +547,7 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) else ! proceed with computing the Jacobian and state update #ifdef DEBUG - if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & + if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0 & .and. prm%of_debug == of) write(6,'(1x,a,/)') '... not yet done' flush(6) #endif @@ -561,21 +561,21 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) !-------------------------------------------------------------------------------------------------- ! ... of the constitutive stress tangent, assembled from dPdF or material constitutive model "smatrix" allocate(smatrix(3*nIntFaceTot,3*nIntFaceTot), source=0.0_pReal) - do iNum = 1_pInt,nIntFaceTot + do iNum = 1,nIntFaceTot faceID = interface1to4(iNum,param(instance)%Nconstituents) ! assembling of local dPdF into global Jacobian matrix !-------------------------------------------------------------------------------------------------- ! identify the left/bottom/back grain (-|N) iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate sytem iGrN = grain3to1(iGr3N,param(instance)%Nconstituents) ! translate into global grain ID - intFaceN = getInterface(2_pInt*faceID(1),iGr3N) ! identifying the connecting interface in local coordinate system + intFaceN = getInterface(2*faceID(1),iGr3N) ! identifying the connecting interface in local coordinate system normN = interfaceNormal(intFaceN,instance,of) - do iFace = 1_pInt,6_pInt + do iFace = 1,6 intFaceN = getInterface(iFace,iGr3N) ! identifying all interfaces that influence relaxation of the above interface mornN = interfaceNormal(intFaceN,instance,of) iMun = interface4to1(intFaceN,param(instance)%Nconstituents) ! translate the interfaces ID into local 4-dimensional index if (iMun > 0) then ! get the corresponding tangent - do i=1_pInt,3_pInt; do j=1_pInt,3_pInt; do k=1_pInt,3_pInt; do l=1_pInt,3_pInt + do i=1,3; do j=1,3; do k=1,3; do 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) enddo;enddo;enddo;enddo @@ -587,16 +587,16 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) !-------------------------------------------------------------------------------------------------- ! identify the right/up/front grain (+|P) iGr3P = iGr3N - iGr3P(faceID(1)) = iGr3N(faceID(1))+1_pInt ! identifying the grain ID in local coordinate sytem + iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate sytem iGrP = grain3to1(iGr3P,param(instance)%Nconstituents) ! translate into global grain ID - intFaceP = getInterface(2_pInt*faceID(1)-1_pInt,iGr3P) ! identifying the connecting interface in local coordinate system + intFaceP = getInterface(2*faceID(1)-1,iGr3P) ! identifying the connecting interface in local coordinate system normP = interfaceNormal(intFaceP,instance,of) - do iFace = 1_pInt,6_pInt + do iFace = 1,6 intFaceP = getInterface(iFace,iGr3P) ! identifying all interfaces that influence relaxation of the above interface mornP = interfaceNormal(intFaceP,instance,of) iMun = interface4to1(intFaceP,param(instance)%Nconstituents) ! translate the interfaces ID into local 4-dimensional index - if (iMun > 0_pInt) then ! get the corresponding tangent - do i=1_pInt,3_pInt; do j=1_pInt,3_pInt; do k=1_pInt,3_pInt; do l=1_pInt,3_pInt + if (iMun > 0) then ! get the corresponding tangent + do i=1,3; do j=1,3; do k=1,3; do 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) enddo;enddo;enddo;enddo @@ -605,10 +605,10 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) enddo #ifdef DEBUG - if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then + if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0) then write(6,'(1x,a30)')'Jacobian matrix of stress' - do i = 1_pInt,3_pInt*nIntFaceTot - write(6,'(1x,100(e11.4,1x))')(smatrix(i,j), j = 1_pInt,3_pInt*nIntFaceTot) + do i = 1,3*nIntFaceTot + write(6,'(1x,100(e11.4,1x))')(smatrix(i,j), j = 1,3*nIntFaceTot) enddo write(6,*)' ' flush(6) @@ -622,7 +622,7 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) allocate(p_relax(3*nIntFaceTot), source=0.0_pReal) allocate(p_resid(3*nIntFaceTot), source=0.0_pReal) - do ipert = 1_pInt,3_pInt*nIntFaceTot + do ipert = 1,3*nIntFaceTot p_relax = relax p_relax(ipert) = relax(ipert) + pPert_RGC ! perturb the relaxation vector stt%relaxationVector(:,of) = p_relax @@ -633,28 +633,28 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) !-------------------------------------------------------------------------------------------------- ! computing the global stress residual array from the perturbed state p_resid = 0.0_pReal - do iNum = 1_pInt,nIntFaceTot + do iNum = 1,nIntFaceTot faceID = interface1to4(iNum,param(instance)%Nconstituents) ! identifying the interface ID in local coordinate system (4-dimensional index) !-------------------------------------------------------------------------------------------------- ! identify the left/bottom/back grain (-|N) iGr3N = faceID(2:4) ! identify the grain ID in local coordinate system (3-dimensional index) iGrN = grain3to1(iGr3N,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index) - intFaceN = getInterface(2_pInt*faceID(1),iGr3N) ! identify the interface ID of the grain + intFaceN = getInterface(2*faceID(1),iGr3N) ! identify the interface ID of the grain normN = interfaceNormal(intFaceN,instance,of) !-------------------------------------------------------------------------------------------------- ! identify the right/up/front grain (+|P) iGr3P = iGr3N - iGr3P(faceID(1)) = iGr3N(faceID(1))+1_pInt ! identify the grain ID in local coordinate system (3-dimensional index) + iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identify the grain ID in local coordinate system (3-dimensional index) iGrP = grain3to1(iGr3P,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index) - intFaceP = getInterface(2_pInt*faceID(1)-1_pInt,iGr3P) ! identify the interface ID of the grain + intFaceP = getInterface(2*faceID(1)-1,iGr3P) ! identify the interface ID of the grain normP = interfaceNormal(intFaceP,instance,of) !-------------------------------------------------------------------------------------------------- ! compute the residual stress (contribution of mismatch and volume penalties) from perturbed state ! at all interfaces - do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt + 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) & + (pD(i,j,iGrP) - D(i,j,iGrP))*normP(j) & @@ -665,10 +665,10 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) enddo #ifdef DEBUG - if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt) then + if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0) then write(6,'(1x,a30)')'Jacobian matrix of penalty' - do i = 1_pInt,3_pInt*nIntFaceTot - write(6,'(1x,100(e11.4,1x))')(pmatrix(i,j), j = 1_pInt,3_pInt*nIntFaceTot) + do i = 1,3*nIntFaceTot + write(6,'(1x,100(e11.4,1x))')(pmatrix(i,j), j = 1,3*nIntFaceTot) enddo write(6,*)' ' flush(6) @@ -678,15 +678,15 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) !-------------------------------------------------------------------------------------------------- ! ... of the numerical viscosity traction "rmatrix" allocate(rmatrix(3*nIntFaceTot,3*nIntFaceTot),source=0.0_pReal) - forall (i=1_pInt:3_pInt*nIntFaceTot) & + forall (i=1:3*nIntFaceTot) & rmatrix(i,i) = viscModus_RGC*viscPower_RGC/(refRelaxRate_RGC*dt)* & ! tangent due to numerical viscosity traction appears (abs(drelax(i))/(refRelaxRate_RGC*dt))**(viscPower_RGC - 1.0_pReal) ! only in the main diagonal term #ifdef DEBUG - if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt) then + if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0) then write(6,'(1x,a30)')'Jacobian matrix of penalty' - do i = 1_pInt,3_pInt*nIntFaceTot - write(6,'(1x,100(e11.4,1x))')(rmatrix(i,j), j = 1_pInt,3_pInt*nIntFaceTot) + do i = 1,3*nIntFaceTot + write(6,'(1x,100(e11.4,1x))')(rmatrix(i,j), j = 1,3*nIntFaceTot) enddo write(6,*)' ' flush(6) @@ -698,10 +698,10 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) allocate(jmatrix(3*nIntFaceTot,3*nIntFaceTot)); jmatrix = smatrix + pmatrix + rmatrix #ifdef DEBUG - if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt) then + if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0) then write(6,'(1x,a30)')'Jacobian matrix (total)' - do i = 1_pInt,3_pInt*nIntFaceTot - write(6,'(1x,100(e11.4,1x))')(jmatrix(i,j), j = 1_pInt,3_pInt*nIntFaceTot) + do i = 1,3*nIntFaceTot + write(6,'(1x,100(e11.4,1x))')(jmatrix(i,j), j = 1,3*nIntFaceTot) enddo write(6,*)' ' flush(6) @@ -710,14 +710,14 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) !-------------------------------------------------------------------------------------------------- ! computing the update of the state variable (relaxation vectors) using the Jacobian matrix - allocate(jnverse(3_pInt*nIntFaceTot,3_pInt*nIntFaceTot),source=0.0_pReal) + allocate(jnverse(3*nIntFaceTot,3*nIntFaceTot),source=0.0_pReal) call math_invert2(jnverse,error,jmatrix) #ifdef DEBUG - if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt) then + if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0) then write(6,'(1x,a30)')'Jacobian inverse' - do i = 1_pInt,3_pInt*nIntFaceTot - write(6,'(1x,100(e11.4,1x))')(jnverse(i,j), j = 1_pInt,3_pInt*nIntFaceTot) + do i = 1,3*nIntFaceTot + write(6,'(1x,100(e11.4,1x))')(jnverse(i,j), j = 1,3*nIntFaceTot) enddo write(6,*)' ' flush(6) @@ -727,7 +727,7 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) !-------------------------------------------------------------------------------------------------- ! calculate the state update (global relaxation vectors) for the next Newton-Raphson iteration drelax = 0.0_pReal - do i = 1_pInt,3_pInt*nIntFaceTot;do j = 1_pInt,3_pInt*nIntFaceTot + do i = 1,3*nIntFaceTot;do j = 1,3*nIntFaceTot drelax(i) = drelax(i) - jnverse(i,j)*resid(j) ! Calculate the correction for the state variable enddo; enddo stt%relaxationVector(:,of) = relax + drelax ! Updateing the state variable for the next iteration @@ -741,9 +741,9 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) endif #ifdef DEBUG - if (iand(debug_homogenization, debug_levelExtensive) > 0_pInt) then + if (iand(debug_homogenization, debug_levelExtensive) > 0) then write(6,'(1x,a30)')'Returned state: ' - do i = 1_pInt,size(stt%relaxationVector(:,of)) + do i = 1,size(stt%relaxationVector(:,of)) write(6,'(1x,2(e15.8,1x))') stt%relaxationVector(i,of) enddo write(6,*)' ' @@ -769,14 +769,14 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) real(pReal), dimension (:,:,:), intent(in) :: fDef !< deformation gradients real(pReal), dimension (3,3), intent(in) :: avgF !< initial effective stretch tensor - integer(pInt), intent(in) :: ip,el,instance,of + integer, intent(in) :: ip,el,instance,of - integer(pInt), dimension (4) :: intFace - integer(pInt), dimension (3) :: iGrain3,iGNghb3,nGDim + integer, dimension (4) :: intFace + integer, dimension (3) :: iGrain3,iGNghb3,nGDim real(pReal), dimension (3,3) :: gDef,nDef real(pReal), dimension (3) :: nVect,surfCorr real(pReal), dimension (2) :: Gmoduli - integer(pInt) :: iGrain,iGNghb,iFace,i,j,k,l + integer :: iGrain,iGNghb,iFace,i,j,k,l real(pReal) :: muGrain,muGNghb,nDefNorm,bgGrain,bgGNghb real(pReal), parameter :: nDefToler = 1.0e-10_pReal #ifdef DEBUG @@ -796,7 +796,7 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) associate(prm => param(instance)) #ifdef DEBUG - debugActive = iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & + debugActive = iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0 & .and. prm%of_debug == of if (debugActive) then @@ -807,20 +807,20 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) !-------------------------------------------------------------------------------------------------- ! computing the mismatch and penalty stress tensor of all grains - grainLoop: do iGrain = 1_pInt,product(prm%Nconstituents) + grainLoop: do iGrain = 1,product(prm%Nconstituents) Gmoduli = equivalentModuli(iGrain,ip,el) muGrain = Gmoduli(1) ! collecting the equivalent shear modulus of grain bgGrain = Gmoduli(2) ! and the lengthh of Burgers vector iGrain3 = grain1to3(iGrain,prm%Nconstituents) ! get the grain ID in local 3-dimensional index (x,y,z)-position - interfaceLoop: do iFace = 1_pInt,6_pInt + interfaceLoop: do iFace = 1,6 intFace = getInterface(iFace,iGrain3) ! get the 4-dimensional index of the interface in local numbering system of the grain nVect = interfaceNormal(intFace,instance,of) iGNghb3 = iGrain3 ! identify the neighboring grain across the interface iGNghb3(abs(intFace(1))) = iGNghb3(abs(intFace(1))) & + int(real(intFace(1),pReal)/real(abs(intFace(1)),pReal),pInt) where(iGNghb3 < 1) iGNghb3 = nGDim - where(iGNghb3 >nGDim) iGNghb3 = 1_pInt + where(iGNghb3 >nGDim) iGNghb3 = 1 iGNghb = grain3to1(iGNghb3,prm%Nconstituents) ! get the ID of the neighboring grain Gmoduli = equivalentModuli(iGNghb,ip,el) ! collect the shear modulus and Burgers vector of the neighbor muGNghb = Gmoduli(1) @@ -831,8 +831,8 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) ! compute the mismatch tensor of all interfaces nDefNorm = 0.0_pReal nDef = 0.0_pReal - do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt - do k = 1_pInt,3_pInt; do l = 1_pInt,3_pInt + do 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 from the jump of deformation gradient enddo; enddo nDefNorm = nDefNorm + nDef(i,j)**2.0_pReal ! compute the norm of the mismatch tensor @@ -849,7 +849,7 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) !-------------------------------------------------------------------------------------------------- ! compute the stress penalty of all interfaces - do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt; do k = 1_pInt,3_pInt; do l = 1_pInt,3_pInt + do 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*bgGrain + muGNghb*bgGNghb)*prm%xiAlpha & *surfCorr(abs(intFace(1)))/prm%dAlpha(abs(intFace(1))) & *cosh(prm%ciAlpha*nDefNorm) & @@ -889,18 +889,18 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) real(pReal), dimension (:,:,:), intent(in) :: fDef ! deformation gradients real(pReal), dimension (3,3), intent(in) :: fAvg ! overall deformation gradient - integer(pInt), intent(in) :: & + integer, intent(in) :: & Ngrain, & instance, & of real(pReal), dimension(size(vPen,3)) :: gVol - integer(pInt) :: i + integer :: i !-------------------------------------------------------------------------------------------------- ! compute the volumes of grains and of cluster vDiscrep = math_det33(fAvg) ! compute the volume of the cluster - do i = 1_pInt,nGrain + do i = 1,nGrain gVol(i) = math_det33(fDef(1:3,1:3,i)) ! compute the volume of individual grains vDiscrep = vDiscrep - gVol(i)/real(nGrain,pReal) ! calculate the difference/dicrepancy between ! the volume of the cluster and the the total volume of grains @@ -909,13 +909,13 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) !-------------------------------------------------------------------------------------------------- ! calculate the stress and penalty due to volume discrepancy vPen = 0.0_pReal - do i = 1_pInt,nGrain + do i = 1,nGrain vPen(:,:,i) = -1.0_pReal/real(nGrain,pReal)*volDiscrMod_RGC*volDiscrPow_RGC/maxVolDiscr_RGC* & sign((abs(vDiscrep)/maxVolDiscr_RGC)**(volDiscrPow_RGC - 1.0),vDiscrep)* & gVol(i)*transpose(math_inv33(fDef(:,:,i))) #ifdef DEBUG - if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & + if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0 & .and. param(instance)%of_debug == of) then write(6,'(1x,a30,i2)')'Volume penalty of grain: ',i write(6,*) transpose(vPen(:,:,i)) @@ -936,22 +936,23 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) implicit none real(pReal), dimension(3) :: surfaceCorrection + real(pReal), dimension(3,3), intent(in) :: avgF !< average F - integer(pInt), intent(in) :: & + integer, intent(in) :: & instance, & of real(pReal), dimension(3,3) :: invC real(pReal), dimension(3) :: nVect real(pReal) :: detF - integer(pInt) :: i,j,iBase + integer :: i,j,iBase logical :: error call math_invert33(matmul(transpose(avgF),avgF),invC,detF,error) surfaceCorrection = 0.0_pReal - do iBase = 1_pInt,3_pInt - nVect = interfaceNormal([iBase,1_pInt,1_pInt,1_pInt],instance,of) - do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt + do iBase = 1,3 + nVect = interfaceNormal([iBase,1,1,1],instance,of) + do i = 1,3; do j = 1,3 surfaceCorrection(iBase) = surfaceCorrection(iBase) + invC(i,j)*nVect(i)*nVect(j) ! compute the component of (the inverse of) the stretch in the direction of the normal enddo; enddo surfaceCorrection(iBase) = sqrt(surfaceCorrection(iBase))*detF ! get the surface correction factor (area contraction/enlargement) @@ -969,7 +970,8 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) implicit none real(pReal), dimension(2) :: equivalentModuli - integer(pInt), intent(in) :: & + + integer, intent(in) :: & grainID,& ip, & !< integration point number el !< element number @@ -1003,17 +1005,17 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) subroutine grainDeformation(F, avgF, instance, of) implicit none - real(pReal), dimension (:,:,:), intent(out) :: F !< partioned F per grain + real(pReal), dimension(:,:,:), intent(out) :: F !< partioned F per grain - real(pReal), dimension (:,:), intent(in) :: avgF !< averaged F - integer(pInt), intent(in) :: & + real(pReal), dimension(:,:), intent(in) :: avgF !< averaged F + integer, intent(in) :: & instance, & of - real(pReal), dimension (3) :: aVect,nVect - integer(pInt), dimension (4) :: intFace - integer(pInt), dimension (3) :: iGrain3 - integer(pInt) :: iGrain,iFace,i,j + real(pReal), dimension(3) :: aVect,nVect + integer, dimension(4) :: intFace + integer, dimension(3) :: iGrain3 + integer :: iGrain,iFace,i,j !------------------------------------------------------------------------------------------------- ! compute the deformation gradient of individual grains due to relaxations @@ -1021,13 +1023,13 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) associate(prm => param(instance)) F = 0.0_pReal - do iGrain = 1_pInt,product(prm%Nconstituents) + do iGrain = 1,product(prm%Nconstituents) iGrain3 = grain1to3(iGrain,prm%Nconstituents) - do iFace = 1_pInt,6_pInt + do iFace = 1,6 intFace = getInterface(iFace,iGrain3) aVect = relaxationVector(intFace,instance,of) nVect = interfaceNormal(intFace,instance,of) - forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) & + forall (i=1:3,j=1:3) & F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations enddo F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF ! relaxed deformation gradient @@ -1046,12 +1048,12 @@ end function homogenization_RGC_updateState subroutine homogenization_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) implicit none - real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point - real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point + real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point + real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point - real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses - real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses - integer(pInt), intent(in) :: instance + real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses + real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses + integer, intent(in) :: instance avgP = sum(P,3) /real(product(param(instance)%Nconstituents),pReal) dAvgPdAvgF = sum(dPdF,5)/real(product(param(instance)%Nconstituents),pReal) @@ -1065,40 +1067,40 @@ end subroutine homogenization_RGC_averageStressAndItsTangent pure function homogenization_RGC_postResults(instance,of) result(postResults) implicit none - integer(pInt), intent(in) :: & + integer, intent(in) :: & instance, & of - integer(pInt) :: & + integer :: & o,c real(pReal), dimension(sum(homogenization_RGC_sizePostResult(:,instance))) :: & postResults associate(stt => state(instance), dst => dependentState(instance), prm => param(instance)) - c = 0_pInt + c = 0 - outputsLoop: do o = 1_pInt,size(prm%outputID) + outputsLoop: do o = 1,size(prm%outputID) select case(prm%outputID(o)) case (constitutivework_ID) postResults(c+1) = stt%work(of) - c = c + 1_pInt + c = c + 1 case (magnitudemismatch_ID) postResults(c+1:c+3) = dst%mismatch(1:3,of) - c = c + 3_pInt + c = c + 3 case (penaltyenergy_ID) postResults(c+1) = stt%penaltyEnergy(of) - c = c + 1_pInt + c = c + 1 case (volumediscrepancy_ID) postResults(c+1) = dst%volumeDiscrepancy(of) - c = c + 1_pInt + c = c + 1 case (averagerelaxrate_ID) postResults(c+1) = dst%relaxationrate_avg(of) - c = c + 1_pInt + c = c + 1 case (maximumrelaxrate_ID) postResults(c+1) = dst%relaxationrate_max(of) - c = c + 1_pInt + c = c + 1 end select enddo outputsLoop @@ -1114,18 +1116,20 @@ end function homogenization_RGC_postResults pure function relaxationVector(intFace,instance,of) implicit none - integer(pInt), intent(in) :: instance,of + real(pReal), dimension (3) :: relaxationVector + + integer, intent(in) :: instance,of + integer, dimension(4), intent(in) :: intFace !< set of interface ID in 4D array (normal and position) - real(pReal), dimension (3) :: relaxationVector - integer(pInt), dimension (4), intent(in) :: intFace !< set of interface ID in 4D array (normal and position) - integer(pInt) :: & - iNum + integer :: iNum + + !-------------------------------------------------------------------------------------------------- ! collect the interface relaxation vector from the global state array iNum = interface4to1(intFace,param(instance)%Nconstituents) ! identify the position of the interface in global state array - if (iNum > 0_pInt) then + if (iNum > 0) then relaxationVector = state(instance)%relaxationVector((3*iNum-2):(3*iNum),of) else relaxationVector = 0.0_pReal @@ -1140,12 +1144,14 @@ end function relaxationVector pure function interfaceNormal(intFace,instance,of) implicit none - real(pReal), dimension (3) :: interfaceNormal - integer(pInt), dimension (4), intent(in) :: intFace !< interface ID in 4D array (normal and position) - integer(pInt), intent(in) :: & + real(pReal), dimension(3) :: interfaceNormal + + integer, dimension(4), intent(in) :: intFace !< interface ID in 4D array (normal and position) + integer, intent(in) :: & instance, & of - integer(pInt) :: nPos + + integer :: nPos !-------------------------------------------------------------------------------------------------- ! get the normal of the interface, identified from the value of intFace(1) @@ -1153,7 +1159,7 @@ pure function interfaceNormal(intFace,instance,of) nPos = abs(intFace(1)) ! identify the position of the interface in global state array interfaceNormal(nPos) = real(intFace(1)/abs(intFace(1)),pReal) ! get the normal vector w.r.t. cluster axis - interfaceNormal = matmul(dependentState(instance)%orientation(1:3,1:3,of),interfaceNormal) ! map the normal vector into sample coordinate system (basis) + interfaceNormal = matmul(dependentState(instance)%orientation(1:3,1:3,of),interfaceNormal) ! map the normal vector into sample coordinate system (basis) end function interfaceNormal @@ -1164,19 +1170,21 @@ end function interfaceNormal pure function getInterface(iFace,iGrain3) implicit none - integer(pInt), dimension (4) :: getInterface - integer(pInt), dimension (3), intent(in) :: iGrain3 !< grain ID in 3D array - integer(pInt), intent(in) :: iFace !< face index (1..6) mapped like (-e1,-e2,-e3,+e1,+e2,+e3) or iDir = (-1,-2,-3,1,2,3) - integer(pInt) :: iDir + integer, dimension(4) :: getInterface + + integer, dimension(3), intent(in) :: iGrain3 !< grain ID in 3D array + integer, intent(in) :: iFace !< face index (1..6) mapped like (-e1,-e2,-e3,+e1,+e2,+e3) or iDir = (-1,-2,-3,1,2,3) + + integer :: iDir !* Direction of interface normal - iDir = (int(real(iFace-1_pInt,pReal)/2.0_pReal,pInt)+1_pInt)*(-1_pInt)**iFace + iDir = (int(real(iFace-1,pReal)/2.0_pReal)+1)*(-1)**iFace getInterface(1) = iDir !-------------------------------------------------------------------------------------------------- ! identify the interface position by the direction of its normal getInterface(2:4) = iGrain3 - if (iDir < 0_pInt) getInterface(1_pInt-iDir) = getInterface(1_pInt-iDir)-1_pInt ! to have a correlation with coordinate/position in real space + if (iDir < 0) getInterface(1-iDir) = getInterface(1-iDir)-1 ! to have a correlation with coordinate/position in real space end function getInterface @@ -1187,13 +1195,14 @@ end function getInterface pure function grain1to3(grain1,nGDim) implicit none - integer(pInt), dimension(3) :: grain1to3 - integer(pInt), intent(in) :: grain1 !< grain ID in 1D array - integer(pInt), dimension(3), intent(in) :: nGDim + integer, dimension(3) :: grain1to3 - grain1to3 = 1_pInt + [mod((grain1-1_pInt),nGDim(1)), & - mod((grain1-1_pInt)/nGDim(1),nGDim(2)), & - (grain1-1_pInt)/(nGDim(1)*nGDim(2))] + integer, intent(in) :: grain1 !< grain ID in 1D array + integer, dimension(3), intent(in) :: nGDim + + grain1to3 = 1 + [mod((grain1-1),nGDim(1)), & + mod((grain1-1)/nGDim(1),nGDim(2)), & + (grain1-1)/(nGDim(1)*nGDim(2))] end function grain1to3 @@ -1201,15 +1210,15 @@ end function grain1to3 !-------------------------------------------------------------------------------------------------- !> @brief map grain ID from in 3D (local position) to in 1D (global array) !-------------------------------------------------------------------------------------------------- -integer(pInt) pure function grain3to1(grain3,nGDim) +integer pure function grain3to1(grain3,nGDim) implicit none - integer(pInt), dimension(3), intent(in) :: grain3 !< grain ID in 3D array (pos.x,pos.y,pos.z) - integer(pInt), dimension(3), intent(in) :: nGDim + integer, dimension(3), intent(in) :: grain3 !< grain ID in 3D array (pos.x,pos.y,pos.z) + integer, dimension(3), intent(in) :: nGDim grain3to1 = grain3(1) & - + nGDim(1)*(grain3(2)-1_pInt) & - + nGDim(1)*nGDim(2)*(grain3(3)-1_pInt) + + nGDim(1)*(grain3(2)-1) & + + nGDim(1)*nGDim(2)*(grain3(3)-1) end function grain3to1 @@ -1217,44 +1226,44 @@ end function grain3to1 !-------------------------------------------------------------------------------------------------- !> @brief maps interface ID from 4D (normal and local position) into 1D (global array) !-------------------------------------------------------------------------------------------------- -integer(pInt) pure function interface4to1(iFace4D, nGDim) +integer pure function interface4to1(iFace4D, nGDim) implicit none - integer(pInt), dimension(4), intent(in) :: iFace4D !< interface ID in 4D array (n.dir,pos.x,pos.y,pos.z) - integer(pInt), dimension(3), intent(in) :: nGDim + integer, dimension(4), intent(in) :: iFace4D !< interface ID in 4D array (n.dir,pos.x,pos.y,pos.z) + integer, dimension(3), intent(in) :: nGDim select case(abs(iFace4D(1))) - case(1_pInt) - if ((iFace4D(2) == 0_pInt) .or. (iFace4D(2) == nGDim(1))) then - interface4to1 = 0_pInt + case(1) + if ((iFace4D(2) == 0) .or. (iFace4D(2) == nGDim(1))) then + interface4to1 = 0 else - interface4to1 = iFace4D(3) + nGDim(2)*(iFace4D(4)-1_pInt) & - + nGDim(2)*nGDim(3)*(iFace4D(2)-1_pInt) + interface4to1 = iFace4D(3) + nGDim(2)*(iFace4D(4)-1) & + + nGDim(2)*nGDim(3)*(iFace4D(2)-1) endif - case(2_pInt) - if ((iFace4D(3) == 0_pInt) .or. (iFace4D(3) == nGDim(2))) then - interface4to1 = 0_pInt + case(2) + if ((iFace4D(3) == 0) .or. (iFace4D(3) == nGDim(2))) then + interface4to1 = 0 else - interface4to1 = iFace4D(4) + nGDim(3)*(iFace4D(2)-1_pInt) & - + nGDim(3)*nGDim(1)*(iFace4D(3)-1_pInt) & - + (nGDim(1)-1_pInt)*nGDim(2)*nGDim(3) ! total number of interfaces normal //e1 + interface4to1 = iFace4D(4) + nGDim(3)*(iFace4D(2)-1) & + + nGDim(3)*nGDim(1)*(iFace4D(3)-1) & + + (nGDim(1)-1)*nGDim(2)*nGDim(3) ! total number of interfaces normal //e1 endif - case(3_pInt) - if ((iFace4D(4) == 0_pInt) .or. (iFace4D(4) == nGDim(3))) then - interface4to1 = 0_pInt + case(3) + if ((iFace4D(4) == 0) .or. (iFace4D(4) == nGDim(3))) then + interface4to1 = 0 else - interface4to1 = iFace4D(2) + nGDim(1)*(iFace4D(3)-1_pInt) & - + nGDim(1)*nGDim(2)*(iFace4D(4)-1_pInt) & - + (nGDim(1)-1_pInt)*nGDim(2)*nGDim(3) & ! total number of interfaces normal //e1 - + nGDim(1)*(nGDim(2)-1_pInt)*nGDim(3) ! total number of interfaces normal //e2 + interface4to1 = iFace4D(2) + nGDim(1)*(iFace4D(3)-1) & + + nGDim(1)*nGDim(2)*(iFace4D(4)-1) & + + (nGDim(1)-1)*nGDim(2)*nGDim(3) & ! total number of interfaces normal //e1 + + nGDim(1)*(nGDim(2)-1)*nGDim(3) ! total number of interfaces normal //e2 endif case default - interface4to1 = -1_pInt + interface4to1 = -1 end select @@ -1267,61 +1276,35 @@ end function interface4to1 pure function interface1to4(iFace1D, nGDim) implicit none - integer(pInt), dimension(4) :: interface1to4 - integer(pInt), intent(in) :: iFace1D !< interface ID in 1D array - integer(pInt), dimension(3), intent(in) :: nGDim - integer(pInt), dimension(3) :: nIntFace + integer, dimension(4) :: interface1to4 + + integer, intent(in) :: iFace1D !< interface ID in 1D array + integer, dimension(3), intent(in) :: nGDim + integer, dimension(3) :: nIntFace !-------------------------------------------------------------------------------------------------- ! compute the total number of interfaces, which ... - nIntFace = [(nGDim(1)-1_pInt)*nGDim(2)*nGDim(3), & ! ... normal //e1 - nGDim(1)*(nGDim(2)-1_pInt)*nGDim(3), & ! ... normal //e2 - nGDim(1)*nGDim(2)*(nGDim(3)-1_pInt)] ! ... normal //e3 + nIntFace = [(nGDim(1)-1)*nGDim(2)*nGDim(3), & ! ... normal //e1 + nGDim(1)*(nGDim(2)-1)*nGDim(3), & ! ... normal //e2 + nGDim(1)*nGDim(2)*(nGDim(3)-1)] ! ... normal //e3 !-------------------------------------------------------------------------------------------------- ! get the corresponding interface ID in 4D (normal and local position) if (iFace1D > 0 .and. iFace1D <= nIntFace(1)) then ! interface with normal //e1 - interface1to4(1) = 1_pInt - interface1to4(3) = mod((iFace1D-1_pInt),nGDim(2))+1_pInt - interface1to4(4) = mod(& - int(& - real(iFace1D-1_pInt,pReal)/& - real(nGDim(2),pReal)& - ,pInt)& - ,nGDim(3))+1_pInt - interface1to4(2) = int(& - real(iFace1D-1_pInt,pReal)/& - real(nGDim(2),pReal)/& - real(nGDim(3),pReal)& - ,pInt)+1_pInt + interface1to4(1) = 1 + interface1to4(3) = mod((iFace1D-1),nGDim(2))+1 + interface1to4(4) = mod(int(real(iFace1D-1,pReal)/real(nGDim(2),pReal)),nGDim(3))+1 + interface1to4(2) = int(real(iFace1D-1,pReal)/real(nGDim(2),pReal)/real(nGDim(3),pReal))+1 elseif (iFace1D > nIntFace(1) .and. iFace1D <= (nIntFace(2) + nIntFace(1))) then ! interface with normal //e2 - interface1to4(1) = 2_pInt - interface1to4(4) = mod((iFace1D-nIntFace(1)-1_pInt),nGDim(3))+1_pInt - interface1to4(2) = mod(& - int(& - real(iFace1D-nIntFace(1)-1_pInt,pReal)/& - real(nGDim(3),pReal)& - ,pInt)& - ,nGDim(1))+1_pInt - interface1to4(3) = int(& - real(iFace1D-nIntFace(1)-1_pInt,pReal)/& - real(nGDim(3),pReal)/& - real(nGDim(1),pReal)& - ,pInt)+1_pInt + interface1to4(1) = 2 + interface1to4(4) = mod((iFace1D-nIntFace(1)-1),nGDim(3))+1 + interface1to4(2) = mod(int(real(iFace1D-nIntFace(1)-1,pReal)/real(nGDim(3),pReal)),nGDim(1))+1 + interface1to4(3) = int(real(iFace1D-nIntFace(1)-1,pReal)/real(nGDim(3),pReal)/real(nGDim(1),pReal))+1 elseif (iFace1D > nIntFace(2) + nIntFace(1) .and. iFace1D <= (nIntFace(3) + nIntFace(2) + nIntFace(1))) then ! interface with normal //e3 - interface1to4(1) = 3_pInt - interface1to4(2) = mod((iFace1D-nIntFace(2)-nIntFace(1)-1_pInt),nGDim(1))+1_pInt - interface1to4(3) = mod(& - int(& - real(iFace1D-nIntFace(2)-nIntFace(1)-1_pInt,pReal)/& - real(nGDim(1),pReal)& - ,pInt)& - ,nGDim(2))+1_pInt - interface1to4(4) = int(& - real(iFace1D-nIntFace(2)-nIntFace(1)-1_pInt,pReal)/& - real(nGDim(1),pReal)/& - real(nGDim(2),pReal)& - ,pInt)+1_pInt + interface1to4(1) = 3 + interface1to4(2) = mod((iFace1D-nIntFace(2)-nIntFace(1)-1),nGDim(1))+1 + interface1to4(3) = mod(int(real(iFace1D-nIntFace(2)-nIntFace(1)-1,pReal)/real(nGDim(1),pReal)),nGDim(2))+1 + interface1to4(4) = int(real(iFace1D-nIntFace(2)-nIntFace(1)-1,pReal)/real(nGDim(1),pReal)/real(nGDim(2),pReal))+1 endif end function interface1to4 diff --git a/src/homogenization_isostrain.f90 b/src/homogenization_isostrain.f90 index 366d76b59..f668f2396 100644 --- a/src/homogenization_isostrain.f90 +++ b/src/homogenization_isostrain.f90 @@ -17,7 +17,7 @@ module homogenization_isostrain end enum type, private :: tParameters !< container type for internal constitutive parameters - integer(pInt) :: & + integer :: & Nconstituents integer(kind(average_ID)) :: & mapping @@ -53,7 +53,7 @@ subroutine homogenization_isostrain_init() config_homogenization implicit none - integer(pInt) :: & + integer :: & Ninstance, & h, & NofMyHomog @@ -63,12 +63,12 @@ subroutine homogenization_isostrain_init() write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_ISOSTRAIN_label//' init -+>>>' Ninstance = int(count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID),pInt) - if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0_pInt) & + if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance allocate(param(Ninstance)) ! one container of parameters per instance - do h = 1_pInt, size(homogenization_type) + do h = 1, size(homogenization_type) if (homogenization_type(h) /= HOMOGENIZATION_ISOSTRAIN_ID) cycle associate(prm => param(homogenization_typeInstance(h)),& @@ -82,15 +82,15 @@ subroutine homogenization_isostrain_init() case ('avg') prm%mapping = average_ID case default - call IO_error(211_pInt,ext_msg=trim(tag)//' ('//HOMOGENIZATION_isostrain_label//')') + call IO_error(211,ext_msg=trim(tag)//' ('//HOMOGENIZATION_isostrain_label//')') end select NofMyHomog = count(material_homogenizationAt == h) - homogState(h)%sizeState = 0_pInt - homogState(h)%sizePostResults = 0_pInt - allocate(homogState(h)%state0 (0_pInt,NofMyHomog)) - allocate(homogState(h)%subState0(0_pInt,NofMyHomog)) - allocate(homogState(h)%state (0_pInt,NofMyHomog)) + homogState(h)%sizeState = 0 + homogState(h)%sizePostResults = 0 + allocate(homogState(h)%state0 (0,NofMyHomog)) + allocate(homogState(h)%subState0(0,NofMyHomog)) + allocate(homogState(h)%state (0,NofMyHomog)) end associate @@ -129,7 +129,7 @@ subroutine homogenization_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses - integer(pInt), intent(in) :: instance + integer, intent(in) :: instance associate(prm => param(instance))