DAMASK_EICMD/src/homogenization_RGC.f90

!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Denny Tjahjanto, Max-Planck-Institut für Eisenforschung GmbH
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Relaxed grain cluster (RGC) homogenization scheme
!> Nconstituents is defined as p x q x r (cluster)
!--------------------------------------------------------------------------------------------------
module homogenization_RGC
 use prec, only: &
   pReal, &
   pInt

 implicit none
 private
 integer(pInt),              dimension(:,:),     allocatable,target, public :: &
   homogenization_RGC_sizePostResult
 character(len=64),          dimension(:,:),     allocatable,target, public :: &
   homogenization_RGC_output                                                                        ! name of each post result output

 enum, bind(c) 
   enumerator :: &
     undefined_ID, &
     constitutivework_ID, &
     penaltyenergy_ID, &
     volumediscrepancy_ID, &
     averagerelaxrate_ID,&
     maximumrelaxrate_ID,&
     magnitudemismatch_ID
 end enum
 
 type, private :: tParameters
   integer(pInt), dimension(:), allocatable :: &
     Nconstituents
   real(pReal) :: &
     xiAlpha, &
     ciAlpha
   real(pReal), dimension(:), allocatable :: &
     dAlpha, &
     angles
   integer(pInt) :: &
     of_debug = 0_pInt
   integer(kind(undefined_ID)),         dimension(:),   allocatable   :: &
     outputID
 end type

 type, private :: tRGCstate
   real(pReal), pointer,     dimension(:) :: &
     work, &
     penaltyEnergy
   real(pReal), pointer,     dimension(:,:) :: &
     relaxationVector
 end type tRGCstate

 type, private :: tRGCdependentState
   real(pReal), allocatable,     dimension(:) :: &
     volumeDiscrepancy, &
     relaxationRate_avg, &
     relaxationRate_max
   real(pReal), allocatable,     dimension(:,:) :: &
     mismatch
   real(pReal), allocatable,     dimension(:,:,:) :: &
     orientation
 end type tRGCdependentState

 type(tparameters),          dimension(:), allocatable, private :: &
   param
 type(tRGCstate),            dimension(:), allocatable, private :: &
   state, &
   state0
 type(tRGCdependentState),   dimension(:), allocatable, private :: &
   dependentState

 public :: &
   homogenization_RGC_init, &
   homogenization_RGC_partitionDeformation, &
   homogenization_RGC_averageStressAndItsTangent, &
   homogenization_RGC_updateState, &
   homogenization_RGC_postResults
 private :: &
   relaxationVector, &
   interfaceNormal, &
   getInterface, &
   grain1to3, &
   grain3to1, &
   interface4to1, &
   interface1to4

contains

!--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields, reads information from material configuration file
!--------------------------------------------------------------------------------------------------
subroutine homogenization_RGC_init()
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
 use, intrinsic :: iso_fortran_env, only: &
   compiler_version, &
   compiler_options
#endif
 use prec, only: &
   pReal, &
   pInt 
 use debug, only: &
#ifdef DEBUG
   debug_i, &
   debug_e, &
#endif
   debug_level, &
   debug_homogenization, &
   debug_levelBasic
 use math, only: &
   math_EulerToR,&
   INRAD
 use IO, only: &
   IO_error, &
   IO_timeStamp
 use material, only: &
   homogenization_type, &
   material_homog, &
   homogState, &
   HOMOGENIZATION_RGC_ID, &
   HOMOGENIZATION_RGC_LABEL, &
   homogenization_typeInstance, &
   homogenization_Noutput, &
   homogenization_Ngrains
 use config, only: &
   config_homogenization

 implicit none
 integer(pInt) :: &
   Ninstance, &
   h, i, &
   NofMyHomog, outputSize, &
   sizeState, nIntFaceTot

 character(len=65536),   dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
 
 integer(kind(undefined_ID)) :: &
   outputID
   
 character(len=65536), dimension(:), allocatable :: &
   outputs
 
 write(6,'(/,a)')   ' <<<+-  homogenization_'//HOMOGENIZATION_RGC_label//' init  -+>>>'
 write(6,'(/,a)')   ' Tjahjanto et al., International Journal of Material Forming, 2(1):939–942, 2009'
 write(6,'(a)')     ' https://doi.org/10.1007/s12289-009-0619-1'
 write(6,'(/,a)')   ' Tjahjanto et al., Modelling and Simulation in Materials Science and Engineering, 18:015006, 2010'
 write(6,'(a)')     ' https://doi.org/10.1088/0965-0393/18/1/015006'
 write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"

 Ninstance = int(count(homogenization_type == HOMOGENIZATION_RGC_ID),pInt)
 if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0_pInt) &
   write(6,'(a16,1x,i5,/)') '# instances:',Ninstance

 allocate(param(Ninstance))
 allocate(state(Ninstance))
 allocate(state0(Ninstance))
 allocate(dependentState(Ninstance))
 
 allocate(homogenization_RGC_sizePostResult(maxval(homogenization_Noutput),Ninstance),source=0_pInt)
 allocate(homogenization_RGC_output(maxval(homogenization_Noutput),Ninstance))
          homogenization_RGC_output=''

 do h = 1_pInt, size(homogenization_type)
   if (homogenization_type(h) /= HOMOGENIZATION_RGC_ID) cycle
   associate(prm => param(homogenization_typeInstance(h)), &
             stt => state(homogenization_typeInstance(h)), &
             st0 => state0(homogenization_typeInstance(h)), &
             dst => dependentState(homogenization_typeInstance(h)), &
             config => config_homogenization(h))
             
#ifdef DEBUG
   !if  (h==material_homogenizationAt(debug_e)) then
   !   prm%of_debug = mappingHomogenization(1,debug_i,debug_e)
   !endif
#endif

   prm%Nconstituents = config%getInts('clustersize',requiredShape=[3])
   if (homogenization_Ngrains(h) /= product(prm%Nconstituents)) &
     call IO_error(211_pInt,ext_msg='clustersize ('//HOMOGENIZATION_RGC_label//')')
   prm%xiAlpha = config%getFloat('scalingparameter')
   prm%ciAlpha = config%getFloat('overproportionality')
   prm%dAlpha  = config%getFloats('grainsize',requiredShape=[3])
   prm%angles  = config%getFloats('clusterorientation',requiredShape=[3])

   outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
   allocate(prm%outputID(0))

   do i=1_pInt, size(outputs)
     outputID = undefined_ID
     select case(outputs(i))
     
       case('constitutivework')
         outputID = constitutivework_ID
         outputSize = 1_pInt
       case('penaltyenergy')
         outputID = penaltyenergy_ID
         outputSize = 1_pInt
       case('volumediscrepancy')
         outputID = volumediscrepancy_ID
         outputSize = 1_pInt
       case('averagerelaxrate')
         outputID = averagerelaxrate_ID
         outputSize = 1_pInt
       case('maximumrelaxrate')
         outputID = maximumrelaxrate_ID
         outputSize = 1_pInt
       case('magnitudemismatch')
         outputID = magnitudemismatch_ID
         outputSize = 3_pInt

     end select
     
     if (outputID /= undefined_ID) then
       homogenization_RGC_output(i,homogenization_typeInstance(h)) = outputs(i)
       homogenization_RGC_sizePostResult(i,homogenization_typeInstance(h)) = outputSize
       prm%outputID = [prm%outputID , outputID]
     endif
     
   enddo

   NofMyHomog = count(material_homog == 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))
   sizeState = nIntFaceTot &
             + size(['avg constitutive work ','average penalty energy'])

   homogState(h)%sizeState = sizeState
   homogState(h)%sizePostResults = sum(homogenization_RGC_sizePostResult(:,homogenization_typeInstance(h)))
   allocate(homogState(h)%state0   (sizeState,NofMyHomog), source=0.0_pReal)
   allocate(homogState(h)%subState0(sizeState,NofMyHomog), source=0.0_pReal)
   allocate(homogState(h)%state    (sizeState,NofMyHomog), source=0.0_pReal)

   stt%relaxationVector   => homogState(h)%state(1:nIntFaceTot,:)
   st0%relaxationVector   => homogState(h)%state0(1:nIntFaceTot,:)
   stt%work               => homogState(h)%state(nIntFaceTot+1,:)
   stt%penaltyEnergy      => homogState(h)%state(nIntFaceTot+2,:)

   allocate(dst%volumeDiscrepancy(     NofMyHomog))
   allocate(dst%relaxationRate_avg(    NofMyHomog))
   allocate(dst%relaxationRate_max(    NofMyHomog))
   allocate(dst%mismatch(          3,  NofMyHomog))
   allocate(dst%orientation(       3,3,NofMyHomog))

!--------------------------------------------------------------------------------------------------
! assigning cluster orientations
   dst%orientation = spread(math_EulerToR(prm%angles*inRad),3,NofMyHomog)

   end associate
   
 enddo  
 
end subroutine homogenization_RGC_init


!--------------------------------------------------------------------------------------------------
!> @brief partitions the deformation gradient onto the constituents
!--------------------------------------------------------------------------------------------------
subroutine homogenization_RGC_partitionDeformation(F,avgF,instance,of)
#ifdef DEBUG
 use debug, only: &
   debug_level, &
   debug_homogenization, &
   debug_levelExtensive
#endif

 implicit none
 real(pReal),   dimension (:,:,:), intent(out) :: F                                                 !< partioned F  per grain
 
 real(pReal),   dimension (:,:),   intent(in)  :: avgF                                              !< averaged F
 integer(pInt),                    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

 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)
   iGrain3 = grain1to3(iGrain,prm%Nconstituents)
   do iFace = 1_pInt,6_pInt
     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) &
       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
     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)
     enddo
     write(6,*)' '
     flush(6)
   endif
#endif
 enddo
 
 end associate

end subroutine homogenization_RGC_partitionDeformation


!--------------------------------------------------------------------------------------------------
!> @brief update the internal state of the homogenization scheme and tell whether "done" and 
! "happy" with result
!--------------------------------------------------------------------------------------------------
function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
 use prec, only: &
   dEq0
#ifdef DEBUG
 use debug, only: &
   debug_level, &
   debug_homogenization,&
   debug_levelExtensive, &
   debug_e, &
   debug_i
#endif
 use math, only: &
   math_invert2
 use material, only: &
   material_homogenizationAt, &
   homogenization_typeInstance, &
   mappingHomogenization, &
   homogenization_maxNgrains
 use numerics, only: &
   absTol_RGC, &
   relTol_RGC, &
   absMax_RGC, &
   relMax_RGC, &
   pPert_RGC, &
   maxdRelax_RGC, &
   viscPower_RGC, &
   viscModus_RGC, &
   refRelaxRate_RGC

 implicit none

 real(pReal), dimension (3,3,homogenization_maxNgrains),     intent(in)    :: & 
   P,&                                                                                              !< array of P
   F,&                                                                                              !< array of F
   F0                                                                                               !< array of initial F
 real(pReal), dimension (3,3,3,3,homogenization_maxNgrains), 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) :: &
   ip, &                                                                                            !< integration point number
   el                                                                                               !< element number

 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) instance,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ipert,iGrain,nGrain, of
 real(pReal), dimension (3,3,homogenization_maxNgrains) :: R,pF,pR,D,pD
 real(pReal), dimension (3,homogenization_maxNgrains)   :: NN,devNull33
 real(pReal), dimension (3)                             :: normP,normN,mornP,mornN
 real(pReal) :: residMax,stresMax,devNull
 logical error

 real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix,rmatrix
 real(pReal), dimension(:), allocatable   :: resid,relax,p_relax,p_resid,drelax
 
 zeroTimeStep: if(dEq0(dt)) then
   homogenization_RGC_updateState = .true.                                                          ! pretend everything is fine and return
   return                                                                    
 endif zeroTimeStep

 instance  = homogenization_typeInstance(material_homogenizationAt(el))
 of = mappingHomogenization(1,ip,el)
 
 associate(stt => state(instance), st0 => state0(instance), dst => dependentState(instance), prm => param(instance))
 
!--------------------------------------------------------------------------------------------------
! 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)

!--------------------------------------------------------------------------------------------------
! 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(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
   write(6,'(1x,a30)')'Obtained state: '
   do i = 1_pInt,size(stt%relaxationVector(:,of))
     write(6,'(1x,2(e15.8,1x))') stt%relaxationVector(i,of)
   enddo
   write(6,*)' '
 endif
#endif

!--------------------------------------------------------------------------------------------------
! computing interface mismatch and stress penalty tensor for all interfaces of all grains
 call stressPenalty(R,NN,avgF,F,ip,el,instance,of)

!--------------------------------------------------------------------------------------------------
! calculating volume discrepancy and stress penalty related to overall volume discrepancy 
 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
     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)
     enddo
     write(6,*)' '
   enddo
 endif
#endif

!------------------------------------------------------------------------------------------------
! 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)

!--------------------------------------------------------------------------------------------------
! 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)
   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)
   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
       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
     enddo
   enddo
   
#ifdef DEBUG
   if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) 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,*)' '
   endif
#endif
 enddo

!--------------------------------------------------------------------------------------------------
! convergence check for stress residual
 stresMax = maxval(abs(P))                                                                          ! get the maximum of first Piola-Kirchhoff (material) stress
 stresLoc = int(maxloc(abs(P)),pInt)                                                                ! get the location of the maximum stress
 residMax = maxval(abs(tract))                                                                      ! get the maximum of the residual
 residLoc = int(maxloc(abs(tract)),pInt)                                                            ! get the position of the maximum residual

#ifdef DEBUG
 if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
     .and. debug_e == el .and. debug_i == ip) then
   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, &
              '@ grain',stresLoc(3),'in component',stresLoc(1),stresLoc(2)
   write(6,'(1x,a15,1x,e15.8,1x,a7,i3,1x,a12,i2)')'Max residual: ',residMax, &
              '@ iface',residLoc(1),'in direction',residLoc(2)
   flush(6)
 endif
#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.
#ifdef DEBUG
    if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
        .and. debug_e == el .and. debug_i == ip) then 
     write(6,'(1x,a55,/)')'... done and happy'
     flush(6)
   endif
#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
       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) &
                             + R(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal)
     enddo; enddo
   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_max(of) = maxval(abs(drelax))/dt
   
#ifdef DEBUG
   if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
        .and. debug_e == el .and. debug_i == ip) 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), &
                                                           dst%mismatch(2,of), &
                                                           dst%mismatch(3,of)
     write(6,'(1x,a30,1x,e15.8)')   'Penalty energy: ',          stt%penaltyEnergy(of)
     write(6,'(1x,a30,1x,e15.8,/)') 'Volume discrepancy: ',      dst%volumeDiscrepancy(of)
     write(6,'(1x,a30,1x,e15.8)')   'Maximum relaxation rate: ', dst%relaxationRate_max(of)
     write(6,'(1x,a30,1x,e15.8,/)') 'Average relaxation rate: ', dst%relaxationRate_avg(of)
     flush(6)
   endif
#endif
   
   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
   
#ifdef DEBUG
   if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
       .and. debug_e == el .and. debug_i == ip) then
     write(6,'(1x,a55,/)')'... broken'
     flush(6)
   endif
#endif

   return
   
 else                                                                                               ! proceed with computing the Jacobian and state update
#ifdef DEBUG
   if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
     .and. debug_e == el .and. debug_i == ip) then
     write(6,'(1x,a55,/)')'... not yet done'
     flush(6)
   endif
#endif
   
 endif

!---------------------------------------------------------------------------------------------------
! construct the global Jacobian matrix for updating the global relaxation vector array when 
! convergence is not yet reached ...

!--------------------------------------------------------------------------------------------------
! ... 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
   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
   normN = interfaceNormal(intFaceN,instance,of)
   do iFace = 1_pInt,6_pInt
     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
         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
! projecting the material tangent dPdF into the interface
! to obtain the Jacobian matrix contribution of dPdF
     endif
   enddo
   
!--------------------------------------------------------------------------------------------------
! 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
   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
   normP = interfaceNormal(intFaceP,instance,of)
   do iFace = 1_pInt,6_pInt
     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
         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
     endif
   enddo
 enddo
 
#ifdef DEBUG
 if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) 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)
   enddo
   write(6,*)' '
   flush(6)
 endif
#endif
 
!--------------------------------------------------------------------------------------------------
! ... of the stress penalty tangent (mismatch penalty and volume penalty, computed using numerical 
! perturbation method) "pmatrix"
 allocate(pmatrix(3*nIntFaceTot,3*nIntFaceTot), source=0.0_pReal)
 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
   p_relax = relax
   p_relax(ipert) = relax(ipert) + pPert_RGC                                                        ! perturb the relaxation vector
   stt%relaxationVector(:,of) = p_relax
   call grainDeformation(pF,avgF,instance,of)                                                       ! rain deformation from perturbed state
   call stressPenalty(pR,DevNull33, avgF,pF,ip,el,instance,of)                                      ! stress penalty due to interface mismatch from perturbed state
   call volumePenalty(pD,devNull,   avgF,pF,nGrain,instance,of)                                     ! stress penalty due to volume discrepancy from perturbed state

!--------------------------------------------------------------------------------------------------
! computing the global stress residual array from the perturbed state
   p_resid = 0.0_pReal
   do iNum = 1_pInt,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
     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)
     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
     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
       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) &
                                                     + (pD(i,j,iGrN) - D(i,j,iGrN))*normN(j)
     enddo; enddo
   enddo
   pmatrix(:,ipert) = p_resid/pPert_RGC
 enddo
 
#ifdef DEBUG
 if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt) 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)
   enddo
   write(6,*)' '
   flush(6)
 endif
#endif
 
!--------------------------------------------------------------------------------------------------
! ... of the numerical viscosity traction "rmatrix"
 allocate(rmatrix(3*nIntFaceTot,3*nIntFaceTot),source=0.0_pReal)
 forall (i=1_pInt:3_pInt*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
   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)
   enddo
   write(6,*)' '
   flush(6)
 endif
#endif

!--------------------------------------------------------------------------------------------------
! The overall Jacobian matrix summarizing contributions of smatrix, pmatrix, rmatrix
 allocate(jmatrix(3*nIntFaceTot,3*nIntFaceTot)); jmatrix = smatrix + pmatrix + rmatrix

#ifdef DEBUG
 if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt) 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)
   enddo
   write(6,*)' '
   flush(6)
 endif
#endif

!--------------------------------------------------------------------------------------------------
! 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)
 call math_invert2(jnverse,error,jmatrix)
 
#ifdef DEBUG
 if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt) 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)
   enddo
   write(6,*)' '
   flush(6)
 endif
#endif

!--------------------------------------------------------------------------------------------------
! 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
   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
 if (any(abs(drelax) > maxdRelax_RGC)) then                                                         ! Forcing cutback when the incremental change of relaxation vector becomes too large
   homogenization_RGC_updateState = [.true.,.false.]
   !$OMP CRITICAL (write2out)
   write(6,'(1x,a,1x,i3,1x,a,1x,i3,1x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback'
   write(6,'(1x,a,1x,e15.8)')'due to large relaxation change =',maxval(abs(drelax))
   flush(6)
   !$OMP END CRITICAL (write2out)
 endif

#ifdef DEBUG
 if (iand(debug_homogenization, debug_levelExtensive) > 0_pInt) then
   write(6,'(1x,a30)')'Returned state: '
   do i = 1_pInt,size(stt%relaxationVector(:,of))
     write(6,'(1x,2(e15.8,1x))') stt%relaxationVector(i,of)
   enddo
   write(6,*)' '
   flush(6)
 endif
#endif

 end associate

 contains
 !--------------------------------------------------------------------------------------------------
 !> @brief calculate stress-like penalty due to deformation mismatch
 !--------------------------------------------------------------------------------------------------
 subroutine stressPenalty(rPen,nMis,avgF,fDef,ip,el,instance,of)
  use math, only: &
    math_civita
  use numerics, only: &
    xSmoo_RGC
  
  implicit none
  real(pReal),   dimension (3,3,homogenization_maxNgrains), intent(out) :: rPen                      !< stress-like penalty
  real(pReal),   dimension (3,homogenization_maxNgrains),   intent(out) :: nMis                      !< total amount of mismatch
  
  real(pReal),   dimension (3,3,homogenization_maxNgrains), 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(pInt), dimension (4)   :: intFace
  integer(pInt), 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
  real(pReal)   :: muGrain,muGNghb,nDefNorm,bgGrain,bgGNghb
  real(pReal),                                               parameter  :: nDefToler = 1.0e-10_pReal
#ifdef DEBUG
  logical :: debugActive
#endif

  nGDim = param(instance)%Nconstituents
  rPen = 0.0_pReal
  nMis = 0.0_pReal
 
 !--------------------------------------------------------------------------------------------------
 ! get the correction factor the modulus of penalty stress representing the evolution of area of 
 ! the interfaces due to deformations

  surfCorr = surfaceCorrection(avgF,instance,of)
  
  associate(prm => param(instance))

#ifdef DEBUG
  debugActive = iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
                .and. prm%of_debug == of

  if (debugActive) then
    write(6,'(1x,a20,2(1x,i3))')'Correction factor: ',ip,el
    write(6,*) surfCorr
  endif
#endif
 
 !--------------------------------------------------------------------------------------------------
 ! computing the mismatch and penalty stress tensor of all grains 
  grainLoop: do iGrain = 1_pInt,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
      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
      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)
      bgGNghb = Gmoduli(2)
      gDef = 0.5_pReal*(fDef(1:3,1:3,iGNghb) - fDef(1:3,1:3,iGrain))                                ! difference/jump in deformation gradeint across the neighbor
 
 !--------------------------------------------------------------------------------------------------
 ! 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
          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
      enddo; enddo
      nDefNorm = max(nDefToler,sqrt(nDefNorm))                                                      ! approximation to zero mismatch if mismatch is zero (singularity)
      nMis(abs(intFace(1)),iGrain) = nMis(abs(intFace(1)),iGrain) + nDefNorm                        ! total amount of mismatch experienced by the grain (at all six interfaces)
#ifdef DEBUG
      if (debugActive) then
        write(6,'(1x,a20,i2,1x,a20,1x,i3)')'Mismatch to face: ',intFace(1),'neighbor grain: ',iGNghb
        write(6,*) transpose(nDef)
        write(6,'(1x,a20,e11.4)')'with magnitude: ',nDefNorm
      endif
#endif

 !--------------------------------------------------------------------------------------------------
 ! 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
        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) &
                                               *0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_civita(k,l,j) &
                                               *tanh(nDefNorm/xSmoo_RGC)
      enddo; enddo;enddo; enddo
    enddo interfaceLoop
#ifdef DEBUG
    if (debugActive) then
      write(6,'(1x,a20,i2)')'Penalty of grain: ',iGrain
      write(6,*) transpose(rPen(1:3,1:3,iGrain))
    endif
#endif
 
  enddo grainLoop
  
  end associate
 
 end subroutine stressPenalty
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief calculate stress-like penalty due to volume discrepancy 
 !--------------------------------------------------------------------------------------------------
 subroutine volumePenalty(vPen,vDiscrep,fAvg,fDef,nGrain,instance,of)
  use math, only: &
    math_det33, &
    math_inv33
  use numerics, only: &
    maxVolDiscr_RGC,&
    volDiscrMod_RGC,&
    volDiscrPow_RGC
 
  implicit none
  real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: vPen                        ! stress-like penalty due to volume
  real(pReal), intent(out)                  :: vDiscrep                                              ! total volume discrepancy
  
  real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in)  :: fDef                        ! deformation gradients
  real(pReal), dimension (3,3), intent(in)  :: fAvg                                                  ! overall deformation gradient
  integer(pInt), intent(in) :: &
    Ngrain, &
    instance, &
    of
    
  real(pReal), dimension (homogenization_maxNgrains) :: gVol
  integer(pInt) :: i
   
 !--------------------------------------------------------------------------------------------------
 ! compute the volumes of grains and of cluster
  vDiscrep = math_det33(fAvg)                                                                       ! compute the volume of the cluster
  do i = 1_pInt,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
  enddo
 
 !--------------------------------------------------------------------------------------------------
 ! calculate the stress and penalty due to volume discrepancy
  vPen      = 0.0_pReal
  do i = 1_pInt,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 &
                .and. param(instance)%of_debug == of) then
      write(6,'(1x,a30,i2)')'Volume penalty of grain: ',i
      write(6,*) transpose(vPen(:,:,i))
    endif
#endif
  enddo

 end subroutine volumePenalty


 !--------------------------------------------------------------------------------------------------
 !> @brief compute the correction factor accouted for surface evolution (area change) due to 
 ! deformation
 !--------------------------------------------------------------------------------------------------
 function surfaceCorrection(avgF,instance,of)
  use math, only: &
    math_invert33, &
    math_mul33x33
  
  implicit none
  real(pReal), dimension(3)               :: surfaceCorrection
  real(pReal), dimension(3,3), intent(in) :: avgF                                                    !< average F
  integer(pInt),               intent(in) :: &
    instance, &
    of
  real(pReal), dimension(3,3)             :: invC
  real(pReal), dimension(3)               :: nVect
  real(pReal)  :: detF
  integer(pInt) :: i,j,iBase
  logical     ::  error
 
  call math_invert33(math_mul33x33(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
      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)
  enddo
 
 end function surfaceCorrection


 !--------------------------------------------------------------------------------------------------
 !> @brief compute the equivalent shear and bulk moduli from the elasticity tensor
 !--------------------------------------------------------------------------------------------------
 function equivalentModuli(grainID,ip,el)
  use constitutive, only: &
    constitutive_homogenizedC
 
  implicit none
  real(pReal), dimension(2)    :: equivalentModuli
  integer(pInt), intent(in)    :: &
    grainID,&
    ip, &                                                                                           !< integration point number
    el                                                                                              !< element number
  real(pReal), dimension(6,6) :: elasTens
  real(pReal) :: &
    cEquiv_11, &
    cEquiv_12, &
    cEquiv_44
 
  elasTens = constitutive_homogenizedC(grainID,ip,el)
 
 !--------------------------------------------------------------------------------------------------
 ! compute the equivalent shear modulus after 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
  equivalentModuli(1) = 0.2_pReal*(cEquiv_11 - cEquiv_12) + 0.6_pReal*cEquiv_44
 
 !--------------------------------------------------------------------------------------------------
 ! obtain the length of Burgers vector (could be model dependend)
  equivalentModuli(2) = 2.5e-10_pReal
 
 end function equivalentModuli
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief calculating the grain deformation gradient (the same with 
 ! homogenization_RGC_partitionDeformation, but used only for perturbation scheme)
 !--------------------------------------------------------------------------------------------------
 subroutine grainDeformation(F, avgF, instance, of)
  
  implicit none
  real(pReal),   dimension (:,:,:), intent(out) :: F                                                 !< partioned F  per grain
 
  real(pReal),   dimension (:,:),   intent(in)  :: avgF                                              !< averaged F
  integer(pInt),                    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
 
  !-------------------------------------------------------------------------------------------------
  ! compute the deformation gradient of individual grains due to relaxations
  
  associate(prm => param(instance))
  
  F = 0.0_pReal
  do iGrain = 1_pInt,product(prm%Nconstituents)
    iGrain3 = grain1to3(iGrain,prm%Nconstituents)
    do iFace = 1_pInt,6_pInt
      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) &
        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
  enddo
  
  end associate
  
 end subroutine grainDeformation
 
end function homogenization_RGC_updateState


!--------------------------------------------------------------------------------------------------
!> @brief derive average stress and stiffness from constituent quantities 
!--------------------------------------------------------------------------------------------------
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 (:,:,:),                             intent(in)  :: P                       !< partitioned stresses
 real(pReal), dimension (:,:,:,:,:),                         intent(in)  :: dPdF                    !< partitioned stiffnesses
 integer(pInt),                                              intent(in)  :: instance

 avgP       = sum(P,3)   /real(product(param(instance)%Nconstituents),pReal)
 dAvgPdAvgF = sum(dPdF,5)/real(product(param(instance)%Nconstituents),pReal)

end subroutine homogenization_RGC_averageStressAndItsTangent


!--------------------------------------------------------------------------------------------------
!> @brief return array of homogenization results for post file inclusion 
!--------------------------------------------------------------------------------------------------
pure function homogenization_RGC_postResults(instance,of) result(postResults)

 implicit none
 integer(pInt), intent(in) :: &
   instance, &
   of
   
 integer(pInt) :: &
   o,c
 real(pReal), dimension(sum(homogenization_RGC_sizePostResult(:,instance))) :: &
   postResults

 associate(stt => state(instance), dst => dependentState(instance), prm => param(instance))

 c = 0_pInt

 outputsLoop: do o = 1_pInt,size(prm%outputID)
   select case(prm%outputID(o))
   
     case (constitutivework_ID)
       postResults(c+1) = stt%work(of)
       c = c + 1_pInt
     case (magnitudemismatch_ID)
       postResults(c+1:c+3) = dst%mismatch(1:3,of)
       c = c + 3_pInt
     case (penaltyenergy_ID)
       postResults(c+1) = stt%penaltyEnergy(of)
       c = c + 1_pInt
     case (volumediscrepancy_ID)
       postResults(c+1) = dst%volumeDiscrepancy(of)
       c = c + 1_pInt
     case (averagerelaxrate_ID)
       postResults(c+1) = dst%relaxationrate_avg(of)
       c = c + 1_pInt
     case (maximumrelaxrate_ID)
       postResults(c+1) = dst%relaxationrate_max(of)
       c = c + 1_pInt
   end select
   
 enddo outputsLoop
 
 end associate
 
end function homogenization_RGC_postResults


!--------------------------------------------------------------------------------------------------
!> @brief collect relaxation vectors of an interface
!--------------------------------------------------------------------------------------------------
pure function relaxationVector(intFace,instance,of)

 implicit none
 integer(pInt),                intent(in) :: instance,of
 
 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

!--------------------------------------------------------------------------------------------------
! 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
   relaxationVector = state(instance)%relaxationVector((3*iNum-2):(3*iNum),of)
 else
   relaxationVector = 0.0_pReal
 endif

end function relaxationVector


!--------------------------------------------------------------------------------------------------
!> @brief identify the normal of an interface 
!--------------------------------------------------------------------------------------------------
pure function interfaceNormal(intFace,instance,of)
 use math, only: &
   math_mul33x3
 
 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) :: &
   instance, &
   of
 integer(pInt) ::                            nPos

!--------------------------------------------------------------------------------------------------
! get the normal of the interface, identified from the value of intFace(1)
 interfaceNormal = 0.0_pReal
 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 = math_mul33x3(dependentState(instance)%orientation(1:3,1:3,of),interfaceNormal)   ! map the normal vector into sample coordinate system (basis)

end function interfaceNormal


!--------------------------------------------------------------------------------------------------
!> @brief collect six faces of a grain in 4D (normal and position)
!--------------------------------------------------------------------------------------------------
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
 
!* Direction of interface normal
 iDir = (int(real(iFace-1_pInt,pReal)/2.0_pReal,pInt)+1_pInt)*(-1_pInt)**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

end function getInterface


!--------------------------------------------------------------------------------------------------
!> @brief map grain ID from in 1D (global array) to in 3D (local position)
!--------------------------------------------------------------------------------------------------
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

 grain1to3 = 1_pInt + [mod((grain1-1_pInt),nGDim(1)), &
                       mod((grain1-1_pInt)/nGDim(1),nGDim(2)), &
                       (grain1-1_pInt)/(nGDim(1)*nGDim(2))]

end function grain1to3


!--------------------------------------------------------------------------------------------------
!> @brief map grain ID from in 3D (local position) to in 1D (global array)
!--------------------------------------------------------------------------------------------------
integer(pInt) 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

 grain3to1 = grain3(1) &
           + nGDim(1)*(grain3(2)-1_pInt) &
           + nGDim(1)*nGDim(2)*(grain3(3)-1_pInt)

end function grain3to1


!--------------------------------------------------------------------------------------------------
!> @brief maps interface ID from 4D (normal and local position) into 1D (global array)
!--------------------------------------------------------------------------------------------------
integer(pInt) 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

 
 select case(abs(iFace4D(1)))

   case(1_pInt)
     if ((iFace4D(2) == 0_pInt) .or. (iFace4D(2) == nGDim(1))) then
       interface4to1 = 0_pInt
     else
       interface4to1 = iFace4D(3) + nGDim(2)*(iFace4D(4)-1_pInt) &
                     + nGDim(2)*nGDim(3)*(iFace4D(2)-1_pInt)
     endif

   case(2_pInt)
     if ((iFace4D(3) == 0_pInt) .or. (iFace4D(3) == nGDim(2))) then
       interface4to1 = 0_pInt
     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
     endif

   case(3_pInt)
     if ((iFace4D(4) == 0_pInt) .or. (iFace4D(4) == nGDim(3))) then
       interface4to1 = 0_pInt
     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
     endif

   case default
     interface4to1 = -1_pInt
     
 end select

end function interface4to1


!--------------------------------------------------------------------------------------------------
!> @brief maps interface ID from 1D (global array) into 4D (normal and local position)
!--------------------------------------------------------------------------------------------------
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

!--------------------------------------------------------------------------------------------------
! 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

!--------------------------------------------------------------------------------------------------
! 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
 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
 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
 endif

end function interface1to4


end module homogenization_RGC
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
-												reading in parameters using new style in parallel

											
										
										
											2018-08-28 16:45:38 +05:30
+								!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!> @author Denny Tjahjanto, Max-Planck-Institut für Eisenforschung GmbH
 								!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
 								!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
 								!> @brief Relaxed grain cluster (RGC) homogenization scheme
-												reading in parameters using new style in parallel

											
										
										
											2018-08-28 16:45:38 +05:30
+								!> Nconstituents is defined as p x q x r (cluster)
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
 								module homogenization_RGC
 								 use prec, only: &
 								   pReal, &
 								   pInt
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												added new, flexible debugging scheme.
now all modules have their own debug specification.
compiles and runs, I hope nothing is broken
did a lot of polishing

											
										
										
											2012-03-09 01:55:28 +05:30
+								 implicit none
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								 private
-												homogenization: added enums and sourced allocation for RGC, some higher verbosity for isostrain

											
										
										
											2013-12-18 12:58:01 +05:30
+								 integer(pInt),              dimension(:,:),     allocatable,target, public :: &
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								   homogenization_RGC_sizePostResult
-												homogenization: added enums and sourced allocation for RGC, some higher verbosity for isostrain

											
										
										
											2013-12-18 12:58:01 +05:30
+								 character(len=64),          dimension(:,:),     allocatable,target, public :: &
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								   homogenization_RGC_output                                                                        ! name of each post result output
-												parameters easier to identify

											
										
										
											2018-08-28 16:27:22 +05:30
-												almost done

											
										
										
											2018-08-31 04:54:28 +05:30
+								 enum, bind(c)
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								   enumerator :: &
 								     undefined_ID, &
 								     constitutivework_ID, &
 								     penaltyenergy_ID, &
 								     volumediscrepancy_ID, &
 								     averagerelaxrate_ID,&
 								     maximumrelaxrate_ID,&
 								     magnitudemismatch_ID
-												almost done

											
										
										
											2018-08-31 04:54:28 +05:30
+								 end enum
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								 type, private :: tParameters
-												parameters easier to identify

											
										
										
											2018-08-28 16:27:22 +05:30
+								   integer(pInt), dimension(:), allocatable :: &
 								     Nconstituents
 								   real(pReal) :: &
 								     xiAlpha, &
 								     ciAlpha
 								   real(pReal), dimension(:), allocatable :: &
 								     dAlpha, &
 								     angles
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								   integer(pInt) :: &
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								     of_debug = 0_pInt
-												almost done

											
										
										
											2018-08-31 04:54:28 +05:30
+								   integer(kind(undefined_ID)),         dimension(:),   allocatable   :: &
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								     outputID
-												parameters easier to identify

											
										
										
											2018-08-28 16:27:22 +05:30
+								 end type
-												cleaning

											
										
										
											2018-11-04 02:06:49 +05:30
+								 type, private :: tRGCstate
-												cleaned and unified

											
										
										
											2018-11-03 21:20:43 +05:30
+								   real(pReal), pointer,     dimension(:) :: &
-												mesh_element should not be used anymore

											
										
										
											2018-10-13 21:51:13 +05:30
+								     work, &
-												no need to be real state

											
										
										
											2019-01-13 03:52:13 +05:30
+								     penaltyEnergy
-												cleaned and unified

											
										
										
											2018-11-03 21:20:43 +05:30
+								   real(pReal), pointer,     dimension(:,:) :: &
-												no need to be real state

											
										
										
											2019-01-13 03:52:13 +05:30
+								     relaxationVector
-												cleaning

											
										
										
											2018-11-04 02:06:49 +05:30
+								 end type tRGCstate
 								 type, private :: tRGCdependentState
-												no need to be real state

											
										
										
											2019-01-13 03:52:13 +05:30
+								   real(pReal), allocatable,     dimension(:) :: &
 								     volumeDiscrepancy, &
 								     relaxationRate_avg, &
 								     relaxationRate_max
 								   real(pReal), allocatable,     dimension(:,:) :: &
 								     mismatch
-												cleaning

											
										
										
											2018-11-04 02:06:49 +05:30
+								   real(pReal), allocatable,     dimension(:,:,:) :: &
 								     orientation
 								 end type tRGCdependentState
-												no need to be real state

											
										
										
											2019-01-13 03:52:13 +05:30
+								 type(tparameters),          dimension(:), allocatable, private :: &
 								   param
 								 type(tRGCstate),            dimension(:), allocatable, private :: &
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								   state, &
 								   state0
-												no need to be real state

											
										
										
											2019-01-13 03:52:13 +05:30
+								 type(tRGCdependentState),   dimension(:), allocatable, private :: &
 								   dependentState
-												cleaning

											
										
										
											2018-11-04 02:06:49 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								 public :: &
 								   homogenization_RGC_init, &
 								   homogenization_RGC_partitionDeformation, &
 								   homogenization_RGC_averageStressAndItsTangent, &
 								   homogenization_RGC_updateState, &
 								   homogenization_RGC_postResults
 								 private :: &
-												no need to have long prefixes for local variables and functions

											
										
										
											2018-08-30 14:20:52 +05:30
+								   relaxationVector, &
 								   interfaceNormal, &
 								   getInterface, &
 								   grain1to3, &
 								   grain3to1, &
 								   interface4to1, &
 								   interface1to4
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
 								contains
 								!--------------------------------------------------------------------------------------------------
-												fixed typo

											
										
										
											2018-04-17 18:12:04 +05:30
+								!> @brief allocates all necessary fields, reads information from material configuration file
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
-												cleaned and unified

											
										
										
											2018-11-03 21:20:43 +05:30
+								subroutine homogenization_RGC_init()
-												compiler_options and compiler_version supported by new Intel compiler

											
										
										
											2018-02-02 17:06:09 +05:30
+								#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
-												suppress warning when compiling with gfortran

											
										
										
											2017-10-05 20:05:34 +05:30
+								 use, intrinsic :: iso_fortran_env, only: &
 								   compiler_version, &
 								   compiler_options
 								#endif
-												started introducing new state structure in homogenisation

											
										
										
											2014-08-21 23:18:20 +05:30
+								 use prec, only: &
 								   pReal, &
 								   pInt
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								 use debug, only: &
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								#ifdef DEBUG
 								   debug_i, &
 								   debug_e, &
 								#endif
 								   debug_level, &
 								   debug_homogenization, &
 								   debug_levelBasic
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								 use math, only: &
 								   math_EulerToR,&
 								   INRAD
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								 use IO, only: &
 								   IO_error, &
 								   IO_timeStamp
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								 use material, only: &
 								   homogenization_type, &
 								   material_homog, &
 								   homogState, &
 								   HOMOGENIZATION_RGC_ID, &
 								   HOMOGENIZATION_RGC_LABEL, &
 								   homogenization_typeInstance, &
 								   homogenization_Noutput, &
 								   homogenization_Ngrains
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								 use config, only: &
 								   config_homogenization
-												added new, flexible debugging scheme.
now all modules have their own debug specification.
compiles and runs, I hope nothing is broken
did a lot of polishing

											
										
										
											2012-03-09 01:55:28 +05:30
 								 implicit none
-												cleaned and unified

											
										
										
											2018-11-03 21:20:43 +05:30
+								 integer(pInt) :: &
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								   Ninstance, &
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								   h, i, &
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								   NofMyHomog, outputSize, &
 								   sizeState, nIntFaceTot
-												almost done

											
										
										
											2018-08-31 04:54:28 +05:30
+								 character(len=65536),   dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
 								 integer(kind(undefined_ID)) :: &
 								   outputID
 								 character(len=65536), dimension(:), allocatable :: &
 								   outputs
-												remove omp statements during initialization, "line" for reading in from material.config now has a default (empty) value to prevent conditional jump depend on this value

											
										
										
											2013-01-09 03:41:59 +05:30
-												output already disabled for worldrank !=0

											
										
										
											2017-11-21 19:38:45 +05:30
+								 write(6,'(/,a)')   ' <<<+-  homogenization_'//HOMOGENIZATION_RGC_label//' init  -+>>>'
-												"et al." instead of "et. al" abbreviates "et alii", meaning "and others"

											
										
										
											2018-04-22 13:37:49 +05:30
+								 write(6,'(/,a)')   ' Tjahjanto et al., International Journal of Material Forming, 2(1):939–942, 2009'
-												almost done

											
										
										
											2018-08-31 04:54:28 +05:30
+								 write(6,'(a)')     ' https://doi.org/10.1007/s12289-009-0619-1'
-												"et al." instead of "et. al" abbreviates "et alii", meaning "and others"

											
										
										
											2018-04-22 13:37:49 +05:30
+								 write(6,'(/,a)')   ' Tjahjanto et al., Modelling and Simulation in Materials Science and Engineering, 18:015006, 2010'
-												almost done

											
										
										
											2018-08-31 04:54:28 +05:30
+								 write(6,'(a)')     ' https://doi.org/10.1088/0965-0393/18/1/015006'
-												output already disabled for worldrank !=0

											
										
										
											2017-11-21 19:38:45 +05:30
+								 write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
-												made code standard conform to Fortran 2008 (ignoring warning concerning comments beyond character 132). Basically, changing "x" format specifier to "Nx" ("1x") plus removing $ format specifier
added compiler switches for gfortran and ifort to check for standard conformity
old gnu compilers <4.4 are not longer supported because they don't provide the c binding for fftw 

											
										
										
											2012-02-01 00:48:55 +05:30
+								#include "compilation_info.f90"
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								 Ninstance = int(count(homogenization_type == HOMOGENIZATION_RGC_ID),pInt)
-												parameters easier to identify

											
										
										
											2018-08-28 16:27:22 +05:30
+								 if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0_pInt) &
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								   write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
-												fixed two memory faults (unallocated and wrong index)

											
										
										
											2018-08-29 13:22:22 +05:30
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								 allocate(param(Ninstance))
 								 allocate(state(Ninstance))
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								 allocate(state0(Ninstance))
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								 allocate(dependentState(Ninstance))
 								 allocate(homogenization_RGC_sizePostResult(maxval(homogenization_Noutput),Ninstance),source=0_pInt)
 								 allocate(homogenization_RGC_output(maxval(homogenization_Noutput),Ninstance))
 								          homogenization_RGC_output=''
-												fixed two memory faults (unallocated and wrong index)

											
										
										
											2018-08-29 13:22:22 +05:30
-												reading in parameters using new style in parallel

											
										
										
											2018-08-28 16:45:38 +05:30
+								 do h = 1_pInt, size(homogenization_type)
 								   if (homogenization_type(h) /= HOMOGENIZATION_RGC_ID) cycle
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								   associate(prm => param(homogenization_typeInstance(h)), &
 								             stt => state(homogenization_typeInstance(h)), &
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								             st0 => state0(homogenization_typeInstance(h)), &
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								             dst => dependentState(homogenization_typeInstance(h)), &
 								             config => config_homogenization(h))
 								#ifdef DEBUG
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								   !if  (h==material_homogenizationAt(debug_e)) then
 								   !   prm%of_debug = mappingHomogenization(1,debug_i,debug_e)
 								   !endif
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								#endif
-												almost done

											
										
										
											2018-08-31 04:54:28 +05:30
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								   prm%Nconstituents = config%getInts('clustersize',requiredShape=[3])
-												fixed two memory faults (unallocated and wrong index)

											
										
										
											2018-08-29 13:22:22 +05:30
+								   if (homogenization_Ngrains(h) /= product(prm%Nconstituents)) &
-												almost done

											
										
										
											2018-08-31 04:54:28 +05:30
+								     call IO_error(211_pInt,ext_msg='clustersize ('//HOMOGENIZATION_RGC_label//')')
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								   prm%xiAlpha = config%getFloat('scalingparameter')
 								   prm%ciAlpha = config%getFloat('overproportionality')
 								   prm%dAlpha  = config%getFloats('grainsize',requiredShape=[3])
 								   prm%angles  = config%getFloats('clusterorientation',requiredShape=[3])
-												handling cluster orientation

											
										
										
											2018-08-29 17:06:46 +05:30
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								   outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
-												almost done

											
										
										
											2018-08-31 04:54:28 +05:30
+								   allocate(prm%outputID(0))
 								   do i=1_pInt, size(outputs)
 								     outputID = undefined_ID
 								     select case(outputs(i))
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
-												almost done

											
										
										
											2018-08-31 04:54:28 +05:30
+								       case('constitutivework')
 								         outputID = constitutivework_ID
 								         outputSize = 1_pInt
 								       case('penaltyenergy')
 								         outputID = penaltyenergy_ID
 								         outputSize = 1_pInt
 								       case('volumediscrepancy')
 								         outputID = volumediscrepancy_ID
 								         outputSize = 1_pInt
 								       case('averagerelaxrate')
 								         outputID = averagerelaxrate_ID
 								         outputSize = 1_pInt
 								       case('maximumrelaxrate')
 								         outputID = maximumrelaxrate_ID
 								         outputSize = 1_pInt
 								       case('magnitudemismatch')
 								         outputID = magnitudemismatch_ID
 								         outputSize = 3_pInt
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
-												almost done

											
										
										
											2018-08-31 04:54:28 +05:30
+								     end select
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
-												bugfixes

flaws introduced in 42-xx branch and detected by new test
- outputsize was wrong
- nconstituents is product along the 3 directions, not the sum

											
										
										
											2018-11-03 19:43:11 +05:30
+								     if (outputID /= undefined_ID) then
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								       homogenization_RGC_output(i,homogenization_typeInstance(h)) = outputs(i)
 								       homogenization_RGC_sizePostResult(i,homogenization_typeInstance(h)) = outputSize
-												bugfixes

flaws introduced in 42-xx branch and detected by new test
- outputsize was wrong
- nconstituents is product along the 3 directions, not the sum

											
										
										
											2018-11-03 19:43:11 +05:30
+								       prm%outputID = [prm%outputID , outputID]
 								     endif
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
-												almost done

											
										
										
											2018-08-31 04:54:28 +05:30
+								   enddo
-												handling cluster orientation

											
										
										
											2018-08-29 17:06:46 +05:30
-												requested output is stored in prm%outputID

											
										
										
											2018-08-31 10:05:30 +05:30
+								   NofMyHomog = count(material_homog == h)
-												cleaned and unified

											
										
										
											2018-11-03 21:20:43 +05:30
+								   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))
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								   sizeState = nIntFaceTot &
-												no need to be real state

											
										
										
											2019-01-13 03:52:13 +05:30
+								             + size(['avg constitutive work ','average penalty energy'])
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
 								   homogState(h)%sizeState = sizeState
 								   homogState(h)%sizePostResults = sum(homogenization_RGC_sizePostResult(:,homogenization_typeInstance(h)))
 								   allocate(homogState(h)%state0   (sizeState,NofMyHomog), source=0.0_pReal)
 								   allocate(homogState(h)%subState0(sizeState,NofMyHomog), source=0.0_pReal)
 								   allocate(homogState(h)%state    (sizeState,NofMyHomog), source=0.0_pReal)
 								   stt%relaxationVector   => homogState(h)%state(1:nIntFaceTot,:)
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								   st0%relaxationVector   => homogState(h)%state0(1:nIntFaceTot,:)
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								   stt%work               => homogState(h)%state(nIntFaceTot+1,:)
-												no need to be real state

											
										
										
											2019-01-13 03:52:13 +05:30
+								   stt%penaltyEnergy      => homogState(h)%state(nIntFaceTot+2,:)
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
-												no need to be real state

											
										
										
											2019-01-13 03:52:13 +05:30
+								   allocate(dst%volumeDiscrepancy(     NofMyHomog))
 								   allocate(dst%relaxationRate_avg(    NofMyHomog))
 								   allocate(dst%relaxationRate_max(    NofMyHomog))
 								   allocate(dst%mismatch(          3,  NofMyHomog))
 								   allocate(dst%orientation(       3,3,NofMyHomog))
-												cleaning

											
										
										
											2018-11-04 02:06:49 +05:30
 								!--------------------------------------------------------------------------------------------------
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								! assigning cluster orientations
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								   dst%orientation = spread(math_EulerToR(prm%angles*inRad),3,NofMyHomog)
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
-												reading in parameters using new style in parallel

											
										
										
											2018-08-28 16:45:38 +05:30
+								   end associate
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
-												reading in parameters using new style in parallel

											
										
										
											2018-08-28 16:45:38 +05:30
+								 enddo
-												fixed buggy initialization. Now initialising in similar way to constitutive models.

											
										
										
											2014-09-10 19:44:03 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								end subroutine homogenization_RGC_init
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
 								!> @brief partitions the deformation gradient onto the constituents
 								!--------------------------------------------------------------------------------------------------
-												cleaning

											
										
										
											2018-11-04 03:53:52 +05:30
+								subroutine homogenization_RGC_partitionDeformation(F,avgF,instance,of)
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								#ifdef DEBUG
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								 use debug, only: &
 								   debug_level, &
 								   debug_homogenization, &
 								   debug_levelExtensive
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								#endif
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
 								 implicit none
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								 real(pReal),   dimension (:,:,:), intent(out) :: F                                                 !< partioned F  per grain
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								 real(pReal),   dimension (:,:),   intent(in)  :: avgF                                              !< averaged F
 								 integer(pInt),                    intent(in)  :: &
-												cleaning

											
										
										
											2018-11-04 03:53:52 +05:30
+								   instance, &
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								   of
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								 real(pReal),   dimension (3) :: aVect,nVect
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								 integer(pInt), dimension (4) :: intFace
 								 integer(pInt), dimension (3) :: iGrain3
-												cleaning

											
										
										
											2018-11-04 03:53:52 +05:30
+								 integer(pInt) ::  iGrain,iFace,i,j
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												simplifying

											
										
										
											2018-11-04 03:43:20 +05:30
+								 associate(prm => param(instance))
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								!--------------------------------------------------------------------------------------------------
 								! compute the deformation gradient of individual grains due to relaxations
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								 F = 0.0_pReal
-												cleaning

											
										
										
											2018-11-04 03:53:52 +05:30
+								 do iGrain = 1_pInt,product(prm%Nconstituents)
-												simplified

											
										
										
											2018-11-04 12:41:35 +05:30
+								   iGrain3 = grain1to3(iGrain,prm%Nconstituents)
-												cleaning

											
										
										
											2018-11-04 02:06:49 +05:30
+								   do iFace = 1_pInt,6_pInt
-												simplifying

											
										
										
											2018-11-04 03:43:20 +05:30
+								     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)
-												polishing, adding _pInt etc. where applicable
post_results now handels zero increment different (like FEM, it is always there even if it is not counted)

											
										
										
											2012-02-13 19:48:07 +05:30
+								     forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) &
-												mesh_element should not be used anymore

											
										
										
											2018-10-13 21:51:13 +05:30
+								       F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j)                                            ! calculating deformation relaxations due to interface relaxation
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								   enddo
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								   F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF                                                     ! resulting relaxed deformation gradient
-. Adding more comments
2. Changing 'subroutine' into 'function' for optimizing the computation time

											
										
										
											2010-11-26 17:20:20 +05:30
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								#ifdef DEBUG
-												changed variable name "debug_what" to "debug_level"

											
										
										
											2012-07-05 15:24:50 +05:30
+								   if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then
-												made code standard conform to Fortran 2008 (ignoring warning concerning comments beyond character 132). Basically, changing "x" format specifier to "Nx" ("1x") plus removing $ format specifier
added compiler switches for gfortran and ifort to check for standard conformity
old gnu compilers <4.4 are not longer supported because they don't provide the c binding for fftw 

											
										
										
											2012-02-01 00:48:55 +05:30
+								     write(6,'(1x,a32,1x,i3)')'Deformation gradient of grain: ',iGrain
-												polishing, adding _pInt etc. where applicable
post_results now handels zero increment different (like FEM, it is always there even if it is not counted)

											
										
										
											2012-02-13 19:48:07 +05:30
+								     do i = 1_pInt,3_pInt
-												polishing: removed variable names like 'unit' and 'data' that are keywords of fortran and ensured that integer and real precision matches independent of machine standard.
removed cut_off parameter for damask_spectral
removed outpot of derived divergence measures and added RMS output in brackets
added comments and options to the makefile

											
										
										
											2012-02-16 00:28:38 +05:30
+								       write(6,'(1x,3(e15.8,1x))')(F(i,j,iGrain), j = 1_pInt,3_pInt)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								     enddo
 								     write(6,*)' '
-												replaced call flush(6) by flush(6), slightly improved formating for output of gfortran version

											
										
										
											2012-11-06 21:20:20 +05:30
+								     flush(6)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								   endif
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								#endif
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								 enddo
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
-												simplifying

											
										
										
											2018-11-04 03:43:20 +05:30
+								 end associate
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								end subroutine homogenization_RGC_partitionDeformation
 								!--------------------------------------------------------------------------------------------------
 								!> @brief update the internal state of the homogenization scheme and tell whether "done" and
 								! "happy" with result
 								!--------------------------------------------------------------------------------------------------
-												Finished applying newstate to Homogenization.

											
										
										
											2014-09-19 23:29:06 +05:30
+								function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
-												using dEq, dNeq, cEq, dEq

											
										
										
											2016-05-29 14:15:03 +05:30
+								 use prec, only: &
-												using improved dEq0/dNeq instead of dEq/dNeq with 0.0 argument

											
										
										
											2016-10-29 13:09:08 +05:30
+								   dEq0
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#ifdef DEBUG
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								 use debug, only: &
 								   debug_level, &
 								   debug_homogenization,&
 								   debug_levelExtensive, &
 								   debug_e, &
 								   debug_i
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#endif
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								 use math, only: &
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								   math_invert2
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								 use material, only: &
-												cleaning

											
										
										
											2018-11-04 13:30:35 +05:30
+								   material_homogenizationAt, &
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								   homogenization_typeInstance, &
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								   mappingHomogenization, &
 								   homogenization_maxNgrains
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								 use numerics, only: &
 								   absTol_RGC, &
 								   relTol_RGC, &
 								   absMax_RGC, &
 								   relMax_RGC, &
 								   pPert_RGC, &
 								   maxdRelax_RGC, &
 								   viscPower_RGC, &
 								   viscModus_RGC, &
 								   refRelaxRate_RGC
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
 								 implicit none
-												Finished applying newstate to Homogenization.

											
										
										
											2014-09-19 23:29:06 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								 real(pReal), dimension (3,3,homogenization_maxNgrains),     intent(in)    :: &
 								   P,&                                                                                              !< array of P
 								   F,&                                                                                              !< array of F
 								   F0                                                                                               !< array of initial F
 								 real(pReal), dimension (3,3,3,3,homogenization_maxNgrains), 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) :: &
 								   ip, &                                                                                            !< integration point number
 								   el                                                                                               !< element number
-												cleaned and unified

											
										
										
											2018-11-03 21:20:43 +05:30
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								 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
-												starting to introduce state pointers

											
										
										
											2018-11-04 00:30:02 +05:30
+								 integer(pInt) instance,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ipert,iGrain,nGrain, of
-												Additional feature to the RGC scheme: The volumetric consistency constraint. This is implemented by penalizing any volumetric discrepancy which occurs due to severe relaxation. The penalty for volume inconsistency is described using a power-law model. The setting for this function can be done in the numerics.config via three new parameters, i.e., maxvoldiscrepancy_RGC, voldiscrepancymod_RGC and discrepancypower_RGC. For monintoring, an extra output variable in prescribed in the material.config in the RGC section, namely "(output) volumediscrepancy". 

IO.f90 and numerics.f90 has been modified accordingly to accommodate these changes

											
										
										
											2009-12-16 21:50:53 +05:30
+								 real(pReal), dimension (3,3,homogenization_maxNgrains) :: R,pF,pR,D,pD
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								 real(pReal), dimension (3,homogenization_maxNgrains)   :: NN,devNull33
-												homogenization: added enums and sourced allocation for RGC, some higher verbosity for isostrain

											
										
										
											2013-12-18 12:58:01 +05:30
+								 real(pReal), dimension (3)                             :: normP,normN,mornP,mornN
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								 real(pReal) :: residMax,stresMax,devNull
-												changed internal debug verbosity in accord with debug.config listing.


											
										
										
											2011-06-06 20:57:35 +05:30
+								 logical error
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
-												Introduction of a numerical viscosity into the RGC scheme to improve the convergent behavior and the numerical stability of the scheme (see: homogenization_RGC.f90). These changes have been incorporated in all related subroutines (homogenization.f90, IO.f90, numerics.config and numerics.f90).

											
										
										
											2009-11-17 19:12:38 +05:30
+								 real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix,rmatrix
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								 real(pReal), dimension(:), allocatable   :: resid,relax,p_relax,p_resid,drelax
-												fixed RGC homogenization in case dt==0
removed misplaced parameter statement in IO.f90

											
										
										
											2013-02-13 15:06:06 +05:30
-												using improved dEq0/dNeq instead of dEq/dNeq with 0.0 argument

											
										
										
											2016-10-29 13:09:08 +05:30
+								 zeroTimeStep: if(dEq0(dt)) then
-												using dEq, dNeq, cEq, dEq

											
										
										
											2016-05-29 14:15:03 +05:30
+								   homogenization_RGC_updateState = .true.                                                          ! pretend everything is fine and return
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								   return
-												using dEq, dNeq, cEq, dEq

											
										
										
											2016-05-29 14:15:03 +05:30
+								 endif zeroTimeStep
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												cleaning

											
										
										
											2018-11-04 13:30:35 +05:30
+								 instance  = homogenization_typeInstance(material_homogenizationAt(el))
-												starting to introduce state pointers

											
										
										
											2018-11-04 00:30:02 +05:30
+								 of = mappingHomogenization(1,ip,el)
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								 associate(stt => state(instance), st0 => state0(instance), dst => dependentState(instance), prm => param(instance))
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
 								!--------------------------------------------------------------------------------------------------
 								! get the dimension of the cluster (grains and interfaces)
 								 nGDim  = prm%Nconstituents
 								 nGrain = product(nGDim)
-												cleaning

											
										
										
											2018-11-04 13:19:40 +05:30
+								 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)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
 								! allocate the size of the global relaxation arrays/jacobian matrices depending on the size of the cluster
-												homogenization: added enums and sourced allocation for RGC, some higher verbosity for isostrain

											
										
										
											2013-12-18 12:58:01 +05:30
+								 allocate(resid(3_pInt*nIntFaceTot),   source=0.0_pReal)
 								 allocate(tract(nIntFaceTot,3),        source=0.0_pReal)
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								 relax  = stt%relaxationVector(:,of)
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								 drelax = stt%relaxationVector(:,of) - st0%relaxationVector(:,of)
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
 								#ifdef DEBUG
-												changed variable name "debug_what" to "debug_level"

											
										
										
											2012-07-05 15:24:50 +05:30
+								 if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then
-												made code standard conform to Fortran 2008 (ignoring warning concerning comments beyond character 132). Basically, changing "x" format specifier to "Nx" ("1x") plus removing $ format specifier
added compiler switches for gfortran and ifort to check for standard conformity
old gnu compilers <4.4 are not longer supported because they don't provide the c binding for fftw 

											
										
										
											2012-02-01 00:48:55 +05:30
+								   write(6,'(1x,a30)')'Obtained state: '
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								   do i = 1_pInt,size(stt%relaxationVector(:,of))
 								     write(6,'(1x,2(e15.8,1x))') stt%relaxationVector(i,of)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								   enddo
 								   write(6,*)' '
 								 endif
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								#endif
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
 								! computing interface mismatch and stress penalty tensor for all interfaces of all grains
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								 call stressPenalty(R,NN,avgF,F,ip,el,instance,of)
-												Additional feature to the RGC scheme: The volumetric consistency constraint. This is implemented by penalizing any volumetric discrepancy which occurs due to severe relaxation. The penalty for volume inconsistency is described using a power-law model. The setting for this function can be done in the numerics.config via three new parameters, i.e., maxvoldiscrepancy_RGC, voldiscrepancymod_RGC and discrepancypower_RGC. For monintoring, an extra output variable in prescribed in the material.config in the RGC section, namely "(output) volumediscrepancy". 

IO.f90 and numerics.f90 has been modified accordingly to accommodate these changes

											
										
										
											2009-12-16 21:50:53 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
 								! calculating volume discrepancy and stress penalty related to overall volume discrepancy
-												no need to be real state

											
										
										
											2019-01-13 03:52:13 +05:30
+								 call volumePenalty(D,dst%volumeDiscrepancy(of),avgF,F,nGrain,instance,of)
-												Additional feature to the RGC scheme: The volumetric consistency constraint. This is implemented by penalizing any volumetric discrepancy which occurs due to severe relaxation. The penalty for volume inconsistency is described using a power-law model. The setting for this function can be done in the numerics.config via three new parameters, i.e., maxvoldiscrepancy_RGC, voldiscrepancymod_RGC and discrepancypower_RGC. For monintoring, an extra output variable in prescribed in the material.config in the RGC section, namely "(output) volumediscrepancy". 

IO.f90 and numerics.f90 has been modified accordingly to accommodate these changes

											
										
										
											2009-12-16 21:50:53 +05:30
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								#ifdef DEBUG
-												changed variable name "debug_what" to "debug_level"

											
										
										
											2012-07-05 15:24:50 +05:30
+								 if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then
-												again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

											
										
										
											2012-02-21 22:01:37 +05:30
+								   do iGrain = 1_pInt,nGrain
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								     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
-												again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

											
										
										
											2012-02-21 22:01:37 +05:30
+								     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), &
-												mesh_element should not be used anymore

											
										
										
											2018-10-13 21:51:13 +05:30
+								                                                                (R(i,j,iGrain), j = 1_pInt,3_pInt), &
 								                                                                (D(i,j,iGrain), j = 1_pInt,3_pInt)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								     enddo
 								     write(6,*)' '
 								   enddo
 								 endif
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								#endif
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												Finished applying newstate to Homogenization.

											
										
										
											2014-09-19 23:29:06 +05:30
+								!------------------------------------------------------------------------------------------------
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								! computing the residual stress from the balance of traction at all (interior) interfaces
-												again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

											
										
										
											2012-02-21 22:01:37 +05:30
+								 do iNum = 1_pInt,nIntFaceTot
-												simplified

											
										
										
											2018-11-04 12:41:35 +05:30
+								   faceID = interface1to4(iNum,param(instance)%Nconstituents)                                         ! identifying the interface ID in local coordinate system (4-dimensional index)
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
 								!--------------------------------------------------------------------------------------------------
 								! identify the left/bottom/back grain (-|N)
 								   iGr3N = faceID(2:4)                                                                              ! identifying the grain ID in local coordinate system (3-dimensional index)
-												simplified

											
										
										
											2018-11-04 12:41:35 +05:30
+								   iGrN = grain3to1(iGr3N,param(instance)%Nconstituents)                                              ! translate the local grain ID into global coordinate system (1-dimensional index)
-												no need to have long prefixes for local variables and functions

											
										
										
											2018-08-30 14:20:52 +05:30
+								   intFaceN = getInterface(2_pInt*faceID(1),iGr3N)
-												simplifying

											
										
										
											2018-11-04 03:43:20 +05:30
+								   normN = interfaceNormal(intFaceN,instance,of)
-. Adding more comments
2. Changing 'subroutine' into 'function' for optimizing the computation time

											
										
										
											2010-11-26 17:20:20 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
 								! identify the right/up/front grain (+|P)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								   iGr3P = iGr3N
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								   iGr3P(faceID(1)) = iGr3N(faceID(1))+1_pInt                                                       ! identifying the grain ID in local coordinate system (3-dimensional index)
-												simplified

											
										
										
											2018-11-04 12:41:35 +05:30
+								   iGrP = grain3to1(iGr3P,param(instance)%Nconstituents)                                              ! translate the local grain ID into global coordinate system (1-dimensional index)
-												no need to have long prefixes for local variables and functions

											
										
										
											2018-08-30 14:20:52 +05:30
+								   intFaceP = getInterface(2_pInt*faceID(1)-1_pInt,iGr3P)
-												simplifying

											
										
										
											2018-11-04 03:43:20 +05:30
+								   normP = interfaceNormal(intFaceP,instance,of)
-. Adding more comments
2. Changing 'subroutine' into 'function' for optimizing the computation time

											
										
										
											2010-11-26 17:20:20 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
 								! compute the residual of traction at the interface (in local system, 4-dimensional index)
-												again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

											
										
										
											2012-02-21 22:01:37 +05:30
+								   do i = 1_pInt,3_pInt
-												polishing, adding _pInt etc. where applicable
post_results now handels zero increment different (like FEM, it is always there even if it is not counted)

											
										
										
											2012-02-13 19:48:07 +05:30
+								     tract(iNum,i) = sign(viscModus_RGC*(abs(drelax(i+3*(iNum-1_pInt)))/(refRelaxRate_RGC*dt))**viscPower_RGC, &
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								                          drelax(i+3*(iNum-1_pInt)))                                                ! contribution from the relaxation viscosity
-												again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

											
										
										
											2012-02-21 22:01:37 +05:30
+								     do j = 1_pInt,3_pInt
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								       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
-												Additional feature to the RGC scheme: The volumetric consistency constraint. This is implemented by penalizing any volumetric discrepancy which occurs due to severe relaxation. The penalty for volume inconsistency is described using a power-law model. The setting for this function can be done in the numerics.config via three new parameters, i.e., maxvoldiscrepancy_RGC, voldiscrepancymod_RGC and discrepancypower_RGC. For monintoring, an extra output variable in prescribed in the material.config in the RGC section, namely "(output) volumediscrepancy". 

IO.f90 and numerics.f90 has been modified accordingly to accommodate these changes

											
										
										
											2009-12-16 21:50:53 +05:30
+								                                     + (P(i,j,iGrN) + R(i,j,iGrN) + D(i,j,iGrN))*normN(j)
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								       resid(i+3_pInt*(iNum-1_pInt)) = tract(iNum,i)                                                ! translate the local residual into global 1-dimensional residual array
-												Introduction of a numerical viscosity into the RGC scheme to improve the convergent behavior and the numerical stability of the scheme (see: homogenization_RGC.f90). These changes have been incorporated in all related subroutines (homogenization.f90, IO.f90, numerics.config and numerics.f90).

											
										
										
											2009-11-17 19:12:38 +05:30
+								     enddo
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								   enddo
-. Adding more comments
2. Changing 'subroutine' into 'function' for optimizing the computation time

											
										
										
											2010-11-26 17:20:20 +05:30
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								#ifdef DEBUG
-												changed variable name "debug_what" to "debug_level"

											
										
										
											2012-07-05 15:24:50 +05:30
+								   if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then
-												made code standard conform to Fortran 2008 (ignoring warning concerning comments beyond character 132). Basically, changing "x" format specifier to "Nx" ("1x") plus removing $ format specifier
added compiler switches for gfortran and ifort to check for standard conformity
old gnu compilers <4.4 are not longer supported because they don't provide the c binding for fftw 

											
										
										
											2012-02-01 00:48:55 +05:30
+								     write(6,'(1x,a30,1x,i3)')'Traction at interface: ',iNum
-												again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

											
										
										
											2012-02-21 22:01:37 +05:30
+								     write(6,'(1x,3(e15.8,1x))')(tract(iNum,j), j = 1_pInt,3_pInt)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								     write(6,*)' '
 								   endif
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								#endif
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								 enddo
-												homogenization: added enums and sourced allocation for RGC, some higher verbosity for isostrain

											
										
										
											2013-12-18 12:58:01 +05:30
+								!--------------------------------------------------------------------------------------------------
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								! convergence check for stress residual
 								 stresMax = maxval(abs(P))                                                                          ! get the maximum of first Piola-Kirchhoff (material) stress
 								 stresLoc = int(maxloc(abs(P)),pInt)                                                                ! get the location of the maximum stress
 								 residMax = maxval(abs(tract))                                                                      ! get the maximum of the residual
 								 residLoc = int(maxloc(abs(tract)),pInt)                                                            ! get the position of the maximum residual
-. Adding more comments
2. Changing 'subroutine' into 'function' for optimizing the computation time

											
										
										
											2010-11-26 17:20:20 +05:30
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								#ifdef DEBUG
-												changed variable name "debug_what" to "debug_level"

											
										
										
											2012-07-05 15:24:50 +05:30
+								 if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
-												added new, flexible debugging scheme.
now all modules have their own debug specification.
compiles and runs, I hope nothing is broken
did a lot of polishing

											
										
										
											2012-03-09 01:55:28 +05:30
+								     .and. debug_e == el .and. debug_i == ip) then
-												made code standard conform to Fortran 2008 (ignoring warning concerning comments beyond character 132). Basically, changing "x" format specifier to "Nx" ("1x") plus removing $ format specifier
added compiler switches for gfortran and ifort to check for standard conformity
old gnu compilers <4.4 are not longer supported because they don't provide the c binding for fftw 

											
										
										
											2012-02-01 00:48:55 +05:30
+								   write(6,'(1x,a)')' '
 								   write(6,'(1x,a,1x,i2,1x,i4)')'RGC residual check ...',ip,el
-												polishing: removed variable names like 'unit' and 'data' that are keywords of fortran and ensured that integer and real precision matches independent of machine standard.
removed cut_off parameter for damask_spectral
removed outpot of derived divergence measures and added RMS output in brackets
added comments and options to the makefile

											
										
										
											2012-02-16 00:28:38 +05:30
+								   write(6,'(1x,a15,1x,e15.8,1x,a7,i3,1x,a12,i2,i2)')'Max stress: ',stresMax, &
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								              '@ grain',stresLoc(3),'in component',stresLoc(1),stresLoc(2)
-												polishing: removed variable names like 'unit' and 'data' that are keywords of fortran and ensured that integer and real precision matches independent of machine standard.
removed cut_off parameter for damask_spectral
removed outpot of derived divergence measures and added RMS output in brackets
added comments and options to the makefile

											
										
										
											2012-02-16 00:28:38 +05:30
+								   write(6,'(1x,a15,1x,e15.8,1x,a7,i3,1x,a12,i2)')'Max residual: ',residMax, &
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								              '@ iface',residLoc(1),'in direction',residLoc(2)
-												replaced call flush(6) by flush(6), slightly improved formating for output of gfortran version

											
										
										
											2012-11-06 21:20:20 +05:30
+								   flush(6)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								 endif
-												cleaning

											
										
										
											2019-01-13 02:34:03 +05:30
+								#endif
-. Adding more comments
2. Changing 'subroutine' into 'function' for optimizing the computation time

											
										
										
											2010-11-26 17:20:20 +05:30
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								 homogenization_RGC_updateState = .false.
-												removed temperature integration and corresponding data structures and debugging options
temperature is stored in crystallite, but homogeneous on one IP (not an component (grain) quantity and an input value parsed in by the BVP solver.
introduced heat, a component (grain) quantity which is homogenized before returned to the heat transfer solver.
went ahead with removal of dummy functions in homogenization and constitutive, this time mainly reduced function signatures to reflect actually needed quantities.


											
										
										
											2013-10-16 18:34:59 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
 								!  If convergence reached => done and happy
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								 if (residMax < relTol_RGC*stresMax .or. residMax < absTol_RGC) then
 								   homogenization_RGC_updateState = .true.
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#ifdef DEBUG
-												changed variable name "debug_what" to "debug_level"

											
										
										
											2012-07-05 15:24:50 +05:30
+								    if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
-												added new, flexible debugging scheme.
now all modules have their own debug specification.
compiles and runs, I hope nothing is broken
did a lot of polishing

											
										
										
											2012-03-09 01:55:28 +05:30
+								        .and. debug_e == el .and. debug_i == ip) then
-												removed temperature integration and corresponding data structures and debugging options
temperature is stored in crystallite, but homogeneous on one IP (not an component (grain) quantity and an input value parsed in by the BVP solver.
introduced heat, a component (grain) quantity which is homogenized before returned to the heat transfer solver.
went ahead with removal of dummy functions in homogenization and constitutive, this time mainly reduced function signatures to reflect actually needed quantities.


											
										
										
											2013-10-16 18:34:59 +05:30
+								     write(6,'(1x,a55,/)')'... done and happy'
-												replaced call flush(6) by flush(6), slightly improved formating for output of gfortran version

											
										
										
											2012-11-06 21:20:20 +05:30
+								     flush(6)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								   endif
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#endif
-. Adding more comments
2. Changing 'subroutine' into 'function' for optimizing the computation time

											
										
										
											2010-11-26 17:20:20 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
 								! compute/update the state for postResult, i.e., all energy densities computed by time-integration
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								   do iGrain = 1_pInt,product(prm%Nconstituents)
-												cleaning

											
										
										
											2018-11-04 03:53:52 +05:30
+								     do i = 1_pInt,3_pInt;do j = 1_pInt,3_pInt
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								       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) &
 								                             + R(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal)
-												cleaning

											
										
										
											2018-11-04 03:53:52 +05:30
+								     enddo; enddo
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								   enddo
-												cleaning

											
										
										
											2018-11-04 02:06:49 +05:30
-												no need to be real state

											
										
										
											2019-01-13 03:52:13 +05:30
+								   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_max(of) = maxval(abs(drelax))/dt
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
 								#ifdef DEBUG
-												changed variable name "debug_what" to "debug_level"

											
										
										
											2012-07-05 15:24:50 +05:30
+								   if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
-												added new, flexible debugging scheme.
now all modules have their own debug specification.
compiles and runs, I hope nothing is broken
did a lot of polishing

											
										
										
											2012-03-09 01:55:28 +05:30
+								        .and. debug_e == el .and. debug_i == ip) then
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								     write(6,'(1x,a30,1x,e15.8)')   'Constitutive work: ',stt%work(of)
-												no need to be real state

											
										
										
											2019-01-13 03:52:13 +05:30
+								     write(6,'(1x,a30,3(1x,e15.8))')'Magnitude mismatch: ',dst%mismatch(1,of), &
 								                                                           dst%mismatch(2,of), &
 								                                                           dst%mismatch(3,of)
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								     write(6,'(1x,a30,1x,e15.8)')   'Penalty energy: ',          stt%penaltyEnergy(of)
-												no need to be real state

											
										
										
											2019-01-13 03:52:13 +05:30
+								     write(6,'(1x,a30,1x,e15.8,/)') 'Volume discrepancy: ',      dst%volumeDiscrepancy(of)
 								     write(6,'(1x,a30,1x,e15.8)')   'Maximum relaxation rate: ', dst%relaxationRate_max(of)
 								     write(6,'(1x,a30,1x,e15.8,/)') 'Average relaxation rate: ', dst%relaxationRate_avg(of)
-												replaced call flush(6) by flush(6), slightly improved formating for output of gfortran version

											
										
										
											2012-11-06 21:20:20 +05:30
+								     flush(6)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								   endif
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#endif
-. Adding more comments
2. Changing 'subroutine' into 'function' for optimizing the computation time

											
										
										
											2010-11-26 17:20:20 +05:30
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								   return
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
 								!--------------------------------------------------------------------------------------------------
 								! 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
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
 								#ifdef DEBUG
-												changed variable name "debug_what" to "debug_level"

											
										
										
											2012-07-05 15:24:50 +05:30
+								   if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
-												added new, flexible debugging scheme.
now all modules have their own debug specification.
compiles and runs, I hope nothing is broken
did a lot of polishing

											
										
										
											2012-03-09 01:55:28 +05:30
+								       .and. debug_e == el .and. debug_i == ip) then
-												removed temperature integration and corresponding data structures and debugging options
temperature is stored in crystallite, but homogeneous on one IP (not an component (grain) quantity and an input value parsed in by the BVP solver.
introduced heat, a component (grain) quantity which is homogenized before returned to the heat transfer solver.
went ahead with removal of dummy functions in homogenization and constitutive, this time mainly reduced function signatures to reflect actually needed quantities.


											
										
										
											2013-10-16 18:34:59 +05:30
+								     write(6,'(1x,a55,/)')'... broken'
-												replaced call flush(6) by flush(6), slightly improved formating for output of gfortran version

											
										
										
											2012-11-06 21:20:20 +05:30
+								     flush(6)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								   endif
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#endif
 								   return
-. Adding more comments
2. Changing 'subroutine' into 'function' for optimizing the computation time

											
										
										
											2010-11-26 17:20:20 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								 else                                                                                               ! proceed with computing the Jacobian and state update
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#ifdef DEBUG
-												changed variable name "debug_what" to "debug_level"

											
										
										
											2012-07-05 15:24:50 +05:30
+								   if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
-												added new, flexible debugging scheme.
now all modules have their own debug specification.
compiles and runs, I hope nothing is broken
did a lot of polishing

											
										
										
											2012-03-09 01:55:28 +05:30
+								     .and. debug_e == el .and. debug_i == ip) then
-												removed temperature integration and corresponding data structures and debugging options
temperature is stored in crystallite, but homogeneous on one IP (not an component (grain) quantity and an input value parsed in by the BVP solver.
introduced heat, a component (grain) quantity which is homogenized before returned to the heat transfer solver.
went ahead with removal of dummy functions in homogenization and constitutive, this time mainly reduced function signatures to reflect actually needed quantities.


											
										
										
											2013-10-16 18:34:59 +05:30
+								     write(6,'(1x,a55,/)')'... not yet done'
-												replaced call flush(6) by flush(6), slightly improved formating for output of gfortran version

											
										
										
											2012-11-06 21:20:20 +05:30
+								     flush(6)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								   endif
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#endif
-. Adding more comments
2. Changing 'subroutine' into 'function' for optimizing the computation time

											
										
										
											2010-11-26 17:20:20 +05:30
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								 endif
-												Finished applying newstate to Homogenization.

											
										
										
											2014-09-19 23:29:06 +05:30
+								!---------------------------------------------------------------------------------------------------
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								! construct the global Jacobian matrix for updating the global relaxation vector array when
 								! convergence is not yet reached ...
 								!--------------------------------------------------------------------------------------------------
 								! ... of the constitutive stress tangent, assembled from dPdF or material constitutive model "smatrix"
-												homogenization: added enums and sourced allocation for RGC, some higher verbosity for isostrain

											
										
										
											2013-12-18 12:58:01 +05:30
+								 allocate(smatrix(3*nIntFaceTot,3*nIntFaceTot), source=0.0_pReal)
-												again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

											
										
										
											2012-02-21 22:01:37 +05:30
+								 do iNum = 1_pInt,nIntFaceTot
-												cleaning

											
										
										
											2018-12-18 22:47:06 +05:30
+								   faceID = interface1to4(iNum,param(instance)%Nconstituents)                                       ! assembling of local dPdF into global Jacobian matrix
-. Adding more comments
2. Changing 'subroutine' into 'function' for optimizing the computation time

											
										
										
											2010-11-26 17:20:20 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
 								! identify the left/bottom/back grain (-|N)
 								   iGr3N = faceID(2:4)                                                                              ! identifying the grain ID in local coordinate sytem
-												cleaning

											
										
										
											2018-12-18 22:47:06 +05:30
+								   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
-												simplifying

											
										
										
											2018-11-04 03:43:20 +05:30
+								   normN = interfaceNormal(intFaceN,instance,of)
-												cleaning

											
										
										
											2018-11-04 02:06:49 +05:30
+								   do iFace = 1_pInt,6_pInt
-												cleaning

											
										
										
											2018-12-18 22:47:06 +05:30
+								     intFaceN = getInterface(iFace,iGr3N)                                                           ! identifying all interfaces that influence relaxation of the above interface
-												simplifying

											
										
										
											2018-11-04 03:43:20 +05:30
+								     mornN = interfaceNormal(intFaceN,instance,of)
-												simplified

											
										
										
											2018-11-04 12:26:27 +05:30
+								     iMun = interface4to1(intFaceN,param(instance)%Nconstituents)                                     ! translate the interfaces ID into local 4-dimensional index
-												homogenization: added enums and sourced allocation for RGC, some higher verbosity for isostrain

											
										
										
											2013-12-18 12:58:01 +05:30
+								     if (iMun > 0) then                                                                             ! get the corresponding tangent
-												again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

											
										
										
											2012-02-21 22:01:37 +05:30
+								       do i=1_pInt,3_pInt; do j=1_pInt,3_pInt; do k=1_pInt,3_pInt; do l=1_pInt,3_pInt
-												cleaning

											
										
										
											2018-11-04 03:53:52 +05:30
+								         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)
-												exchanged forall with nested do loops... Bye, bye smatrix warning, will miss you!!

											
										
										
											2011-08-03 23:28:16 +05:30
+								       enddo;enddo;enddo;enddo
-												started introducing new state structure in homogenisation

											
										
										
											2014-08-21 23:18:20 +05:30
+								! projecting the material tangent dPdF into the interface
 								! to obtain the Jacobian matrix contribution of dPdF
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								     endif
 								   enddo
-. Adding more comments
2. Changing 'subroutine' into 'function' for optimizing the computation time

											
										
										
											2010-11-26 17:20:20 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
 								! identify the right/up/front grain (+|P)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								   iGr3P = iGr3N
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								   iGr3P(faceID(1)) = iGr3N(faceID(1))+1_pInt                                                       ! identifying the grain ID in local coordinate sytem
-												cleaning

											
										
										
											2018-12-18 22:47:06 +05:30
+								   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
-												simplifying

											
										
										
											2018-11-04 03:43:20 +05:30
+								   normP = interfaceNormal(intFaceP,instance,of)
-												cleaning

											
										
										
											2018-11-04 02:06:49 +05:30
+								   do iFace = 1_pInt,6_pInt
-												cleaning

											
										
										
											2018-12-18 22:47:06 +05:30
+								     intFaceP = getInterface(iFace,iGr3P)                                                           ! identifying all interfaces that influence relaxation of the above interface
-												simplifying

											
										
										
											2018-11-04 03:43:20 +05:30
+								     mornP = interfaceNormal(intFaceP,instance,of)
-												cleaning

											
										
										
											2018-12-18 22:47:06 +05:30
+								     iMun = interface4to1(intFaceP,param(instance)%Nconstituents)                                   ! translate the interfaces ID into local 4-dimensional index
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								     if (iMun > 0_pInt) then                                                                        ! get the corresponding tangent
-												again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

											
										
										
											2012-02-21 22:01:37 +05:30
+								       do i=1_pInt,3_pInt; do j=1_pInt,3_pInt; do k=1_pInt,3_pInt; do l=1_pInt,3_pInt
-												cleaning

											
										
										
											2018-11-04 13:19:40 +05:30
+								         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)
-												exchanged forall with nested do loops... Bye, bye smatrix warning, will miss you!!

											
										
										
											2011-08-03 23:28:16 +05:30
+								       enddo;enddo;enddo;enddo
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								     endif
 								   enddo
 								 enddo
-. Adding more comments
2. Changing 'subroutine' into 'function' for optimizing the computation time

											
										
										
											2010-11-26 17:20:20 +05:30
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#ifdef DEBUG
-												changed variable name "debug_what" to "debug_level"

											
										
										
											2012-07-05 15:24:50 +05:30
+								 if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then
-												made code standard conform to Fortran 2008 (ignoring warning concerning comments beyond character 132). Basically, changing "x" format specifier to "Nx" ("1x") plus removing $ format specifier
added compiler switches for gfortran and ifort to check for standard conformity
old gnu compilers <4.4 are not longer supported because they don't provide the c binding for fftw 

											
										
										
											2012-02-01 00:48:55 +05:30
+								   write(6,'(1x,a30)')'Jacobian matrix of stress'
-												again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

											
										
										
											2012-02-21 22:01:37 +05:30
+								   do i = 1_pInt,3_pInt*nIntFaceTot
 								     write(6,'(1x,100(e11.4,1x))')(smatrix(i,j), j = 1_pInt,3_pInt*nIntFaceTot)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								   enddo
 								   write(6,*)' '
-												replaced call flush(6) by flush(6), slightly improved formating for output of gfortran version

											
										
										
											2012-11-06 21:20:20 +05:30
+								   flush(6)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								 endif
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#endif
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
 								!--------------------------------------------------------------------------------------------------
 								! ... of the stress penalty tangent (mismatch penalty and volume penalty, computed using numerical
 								! perturbation method) "pmatrix"
-												homogenization: added enums and sourced allocation for RGC, some higher verbosity for isostrain

											
										
										
											2013-12-18 12:58:01 +05:30
+								 allocate(pmatrix(3*nIntFaceTot,3*nIntFaceTot), source=0.0_pReal)
 								 allocate(p_relax(3*nIntFaceTot),               source=0.0_pReal)
 								 allocate(p_resid(3*nIntFaceTot),               source=0.0_pReal)
-												cleaning

											
										
										
											2018-11-04 13:19:40 +05:30
-												again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

											
										
										
											2012-02-21 22:01:37 +05:30
+								 do ipert = 1_pInt,3_pInt*nIntFaceTot
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								   p_relax = relax
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								   p_relax(ipert) = relax(ipert) + pPert_RGC                                                        ! perturb the relaxation vector
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								   stt%relaxationVector(:,of) = p_relax
-												cleaning

											
										
										
											2018-12-18 22:47:06 +05:30
+								   call grainDeformation(pF,avgF,instance,of)                                                       ! rain deformation from perturbed state
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								   call stressPenalty(pR,DevNull33, avgF,pF,ip,el,instance,of)                                      ! stress penalty due to interface mismatch from perturbed state
 								   call volumePenalty(pD,devNull,   avgF,pF,nGrain,instance,of)                                     ! stress penalty due to volume discrepancy from perturbed state
-. Adding more comments
2. Changing 'subroutine' into 'function' for optimizing the computation time

											
										
										
											2010-11-26 17:20:20 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
 								! computing the global stress residual array from the perturbed state
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								   p_resid = 0.0_pReal
-												again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

											
										
										
											2012-02-21 22:01:37 +05:30
+								   do iNum = 1_pInt,nIntFaceTot
-												simplified

											
										
										
											2018-11-04 12:41:35 +05:30
+								     faceID = interface1to4(iNum,param(instance)%Nconstituents)                                     ! identifying the interface ID in local coordinate system (4-dimensional index)
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
 								!--------------------------------------------------------------------------------------------------
 								! identify the left/bottom/back grain (-|N)
-												cleaning

											
										
										
											2018-12-18 22:47:06 +05:30
+								     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
-												simplifying

											
										
										
											2018-11-04 03:43:20 +05:30
+								     normN = interfaceNormal(intFaceN,instance,of)
-. Adding more comments
2. Changing 'subroutine' into 'function' for optimizing the computation time

											
										
										
											2010-11-26 17:20:20 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
 								! identify the right/up/front grain (+|P)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								     iGr3P = iGr3N
-												cleaning

											
										
										
											2018-12-18 22:47:06 +05:30
+								     iGr3P(faceID(1)) = iGr3N(faceID(1))+1_pInt                                                     ! 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
-												simplifying

											
										
										
											2018-11-04 03:43:20 +05:30
+								     normP = interfaceNormal(intFaceP,instance,of)
-. Adding more comments
2. Changing 'subroutine' into 'function' for optimizing the computation time

											
										
										
											2010-11-26 17:20:20 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
 								! 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
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								       p_resid(i+3*(iNum-1)) = p_resid(i+3*(iNum-1)) + (pR(i,j,iGrP) - R(i,j,iGrP))*normP(j) &
-												Additional feature to the RGC scheme: The volumetric consistency constraint. This is implemented by penalizing any volumetric discrepancy which occurs due to severe relaxation. The penalty for volume inconsistency is described using a power-law model. The setting for this function can be done in the numerics.config via three new parameters, i.e., maxvoldiscrepancy_RGC, voldiscrepancymod_RGC and discrepancypower_RGC. For monintoring, an extra output variable in prescribed in the material.config in the RGC section, namely "(output) volumediscrepancy". 

IO.f90 and numerics.f90 has been modified accordingly to accommodate these changes

											
										
										
											2009-12-16 21:50:53 +05:30
+								                                                     + (pR(i,j,iGrN) - R(i,j,iGrN))*normN(j) &
 								                                                     + (pD(i,j,iGrP) - D(i,j,iGrP))*normP(j) &
 								                                                     + (pD(i,j,iGrN) - D(i,j,iGrN))*normN(j)
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								     enddo; enddo
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								   enddo
 								   pmatrix(:,ipert) = p_resid/pPert_RGC
 								 enddo
-. Adding more comments
2. Changing 'subroutine' into 'function' for optimizing the computation time

											
										
										
											2010-11-26 17:20:20 +05:30
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#ifdef DEBUG
-												changed variable name "debug_what" to "debug_level"

											
										
										
											2012-07-05 15:24:50 +05:30
+								 if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt) then
-												made code standard conform to Fortran 2008 (ignoring warning concerning comments beyond character 132). Basically, changing "x" format specifier to "Nx" ("1x") plus removing $ format specifier
added compiler switches for gfortran and ifort to check for standard conformity
old gnu compilers <4.4 are not longer supported because they don't provide the c binding for fftw 

											
										
										
											2012-02-01 00:48:55 +05:30
+								   write(6,'(1x,a30)')'Jacobian matrix of penalty'
-												again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

											
										
										
											2012-02-21 22:01:37 +05:30
+								   do i = 1_pInt,3_pInt*nIntFaceTot
 								     write(6,'(1x,100(e11.4,1x))')(pmatrix(i,j), j = 1_pInt,3_pInt*nIntFaceTot)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								   enddo
 								   write(6,*)' '
-												replaced call flush(6) by flush(6), slightly improved formating for output of gfortran version

											
										
										
											2012-11-06 21:20:20 +05:30
+								   flush(6)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								 endif
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#endif
-												Introduction of a numerical viscosity into the RGC scheme to improve the convergent behavior and the numerical stability of the scheme (see: homogenization_RGC.f90). These changes have been incorporated in all related subroutines (homogenization.f90, IO.f90, numerics.config and numerics.f90).

											
										
										
											2009-11-17 19:12:38 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
 								! ... of the numerical viscosity traction "rmatrix"
-												homogenization: added enums and sourced allocation for RGC, some higher verbosity for isostrain

											
										
										
											2013-12-18 12:58:01 +05:30
+								 allocate(rmatrix(3*nIntFaceTot,3*nIntFaceTot),source=0.0_pReal)
-												again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

											
										
										
											2012-02-21 22:01:37 +05:30
+								 forall (i=1_pInt:3_pInt*nIntFaceTot) &
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								   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
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
 								#ifdef DEBUG
-												changed variable name "debug_what" to "debug_level"

											
										
										
											2012-07-05 15:24:50 +05:30
+								 if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt) then
-												made code standard conform to Fortran 2008 (ignoring warning concerning comments beyond character 132). Basically, changing "x" format specifier to "Nx" ("1x") plus removing $ format specifier
added compiler switches for gfortran and ifort to check for standard conformity
old gnu compilers <4.4 are not longer supported because they don't provide the c binding for fftw 

											
										
										
											2012-02-01 00:48:55 +05:30
+								   write(6,'(1x,a30)')'Jacobian matrix of penalty'
-												again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

											
										
										
											2012-02-21 22:01:37 +05:30
+								   do i = 1_pInt,3_pInt*nIntFaceTot
 								     write(6,'(1x,100(e11.4,1x))')(rmatrix(i,j), j = 1_pInt,3_pInt*nIntFaceTot)
-												Introduction of a numerical viscosity into the RGC scheme to improve the convergent behavior and the numerical stability of the scheme (see: homogenization_RGC.f90). These changes have been incorporated in all related subroutines (homogenization.f90, IO.f90, numerics.config and numerics.f90).

											
										
										
											2009-11-17 19:12:38 +05:30
+								   enddo
 								   write(6,*)' '
-												replaced call flush(6) by flush(6), slightly improved formating for output of gfortran version

											
										
										
											2012-11-06 21:20:20 +05:30
+								   flush(6)
-												Introduction of a numerical viscosity into the RGC scheme to improve the convergent behavior and the numerical stability of the scheme (see: homogenization_RGC.f90). These changes have been incorporated in all related subroutines (homogenization.f90, IO.f90, numerics.config and numerics.f90).

											
										
										
											2009-11-17 19:12:38 +05:30
+								 endif
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#endif
-												Introduction of a numerical viscosity into the RGC scheme to improve the convergent behavior and the numerical stability of the scheme (see: homogenization_RGC.f90). These changes have been incorporated in all related subroutines (homogenization.f90, IO.f90, numerics.config and numerics.f90).

											
										
										
											2009-11-17 19:12:38 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
 								! The overall Jacobian matrix summarizing contributions of smatrix, pmatrix, rmatrix
-												Introduction of a numerical viscosity into the RGC scheme to improve the convergent behavior and the numerical stability of the scheme (see: homogenization_RGC.f90). These changes have been incorporated in all related subroutines (homogenization.f90, IO.f90, numerics.config and numerics.f90).

											
										
										
											2009-11-17 19:12:38 +05:30
+								 allocate(jmatrix(3*nIntFaceTot,3*nIntFaceTot)); jmatrix = smatrix + pmatrix + rmatrix
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
 								#ifdef DEBUG
-												changed variable name "debug_what" to "debug_level"

											
										
										
											2012-07-05 15:24:50 +05:30
+								 if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt) then
-												made code standard conform to Fortran 2008 (ignoring warning concerning comments beyond character 132). Basically, changing "x" format specifier to "Nx" ("1x") plus removing $ format specifier
added compiler switches for gfortran and ifort to check for standard conformity
old gnu compilers <4.4 are not longer supported because they don't provide the c binding for fftw 

											
										
										
											2012-02-01 00:48:55 +05:30
+								   write(6,'(1x,a30)')'Jacobian matrix (total)'
-												again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

											
										
										
											2012-02-21 22:01:37 +05:30
+								   do i = 1_pInt,3_pInt*nIntFaceTot
 								     write(6,'(1x,100(e11.4,1x))')(jmatrix(i,j), j = 1_pInt,3_pInt*nIntFaceTot)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								   enddo
 								   write(6,*)' '
-												replaced call flush(6) by flush(6), slightly improved formating for output of gfortran version

											
										
										
											2012-11-06 21:20:20 +05:30
+								   flush(6)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								 endif
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#endif
-. Adding more comments
2. Changing 'subroutine' into 'function' for optimizing the computation time

											
										
										
											2010-11-26 17:20:20 +05:30
-												homogenization: added enums and sourced allocation for RGC, some higher verbosity for isostrain

											
										
										
											2013-12-18 12:58:01 +05:30
+								!--------------------------------------------------------------------------------------------------
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								! computing the update of the state variable (relaxation vectors) using the Jacobian matrix
-												homogenization: added enums and sourced allocation for RGC, some higher verbosity for isostrain

											
										
										
											2013-12-18 12:58:01 +05:30
+								 allocate(jnverse(3_pInt*nIntFaceTot,3_pInt*nIntFaceTot),source=0.0_pReal)
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								 call math_invert2(jnverse,error,jmatrix)
-. Adding more comments
2. Changing 'subroutine' into 'function' for optimizing the computation time

											
										
										
											2010-11-26 17:20:20 +05:30
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#ifdef DEBUG
-												changed variable name "debug_what" to "debug_level"

											
										
										
											2012-07-05 15:24:50 +05:30
+								 if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt) then
-												made code standard conform to Fortran 2008 (ignoring warning concerning comments beyond character 132). Basically, changing "x" format specifier to "Nx" ("1x") plus removing $ format specifier
added compiler switches for gfortran and ifort to check for standard conformity
old gnu compilers <4.4 are not longer supported because they don't provide the c binding for fftw 

											
										
										
											2012-02-01 00:48:55 +05:30
+								   write(6,'(1x,a30)')'Jacobian inverse'
-												again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

											
										
										
											2012-02-21 22:01:37 +05:30
+								   do i = 1_pInt,3_pInt*nIntFaceTot
-												polishing: removed variable names like 'unit' and 'data' that are keywords of fortran and ensured that integer and real precision matches independent of machine standard.
removed cut_off parameter for damask_spectral
removed outpot of derived divergence measures and added RMS output in brackets
added comments and options to the makefile

											
										
										
											2012-02-16 00:28:38 +05:30
+								     write(6,'(1x,100(e11.4,1x))')(jnverse(i,j), j = 1_pInt,3_pInt*nIntFaceTot)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								   enddo
 								   write(6,*)' '
-												replaced call flush(6) by flush(6), slightly improved formating for output of gfortran version

											
										
										
											2012-11-06 21:20:20 +05:30
+								   flush(6)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								 endif
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#endif
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
 								! calculate the state update (global relaxation vectors) for the next Newton-Raphson iteration
-												Introduction of a numerical viscosity into the RGC scheme to improve the convergent behavior and the numerical stability of the scheme (see: homogenization_RGC.f90). These changes have been incorporated in all related subroutines (homogenization.f90, IO.f90, numerics.config and numerics.f90).

											
										
										
											2009-11-17 19:12:38 +05:30
+								 drelax = 0.0_pReal
-												cleaning

											
										
										
											2018-11-04 13:19:40 +05:30
+								 do i = 1_pInt,3_pInt*nIntFaceTot;do j = 1_pInt,3_pInt*nIntFaceTot
 								   drelax(i) = drelax(i) - jnverse(i,j)*resid(j)                                                    ! Calculate the correction for the state variable
 								 enddo; enddo
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								 stt%relaxationVector(:,of) = relax + drelax                                                        ! Updateing the state variable for the next iteration
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								 if (any(abs(drelax) > maxdRelax_RGC)) then                                                         ! Forcing cutback when the incremental change of relaxation vector becomes too large
 								   homogenization_RGC_updateState = [.true.,.false.]
-												openmp parallelization working again (at least for j2 and nonlocal constitutive model). 
In order to keep it like that, please follow these simple rules:

DON'T use implicit array subscripts:
example:    real, dimension(3,3) :: A,B
                  A(:,2) = B(:,1)               <--- DON'T USE
                  A(1:3,2) = B(1:3,1)       <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to  prevent memory leaks.

Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)

Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.


											
										
										
											2011-03-17 16:16:17 +05:30
+								   !$OMP CRITICAL (write2out)
-												made code standard conform to Fortran 2008 (ignoring warning concerning comments beyond character 132). Basically, changing "x" format specifier to "Nx" ("1x") plus removing $ format specifier
added compiler switches for gfortran and ifort to check for standard conformity
old gnu compilers <4.4 are not longer supported because they don't provide the c binding for fftw 

											
										
										
											2012-02-01 00:48:55 +05:30
+								   write(6,'(1x,a,1x,i3,1x,a,1x,i3,1x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback'
-												polishing: removed variable names like 'unit' and 'data' that are keywords of fortran and ensured that integer and real precision matches independent of machine standard.
removed cut_off parameter for damask_spectral
removed outpot of derived divergence measures and added RMS output in brackets
added comments and options to the makefile

											
										
										
											2012-02-16 00:28:38 +05:30
+								   write(6,'(1x,a,1x,e15.8)')'due to large relaxation change =',maxval(abs(drelax))
-												replaced call flush(6) by flush(6), slightly improved formating for output of gfortran version

											
										
										
											2012-11-06 21:20:20 +05:30
+								   flush(6)
-												openmp parallelization working again (at least for j2 and nonlocal constitutive model). 
In order to keep it like that, please follow these simple rules:

DON'T use implicit array subscripts:
example:    real, dimension(3,3) :: A,B
                  A(:,2) = B(:,1)               <--- DON'T USE
                  A(1:3,2) = B(1:3,1)       <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to  prevent memory leaks.

Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)

Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.


											
										
										
											2011-03-17 16:16:17 +05:30
+								   !$OMP END CRITICAL (write2out)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								 endif
-. Adding more comments
2. Changing 'subroutine' into 'function' for optimizing the computation time

											
										
										
											2010-11-26 17:20:20 +05:30
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#ifdef DEBUG
-												added new, flexible debugging scheme.
now all modules have their own debug specification.
compiles and runs, I hope nothing is broken
did a lot of polishing

											
										
										
											2012-03-09 01:55:28 +05:30
+								 if (iand(debug_homogenization, debug_levelExtensive) > 0_pInt) then
-												made code standard conform to Fortran 2008 (ignoring warning concerning comments beyond character 132). Basically, changing "x" format specifier to "Nx" ("1x") plus removing $ format specifier
added compiler switches for gfortran and ifort to check for standard conformity
old gnu compilers <4.4 are not longer supported because they don't provide the c binding for fftw 

											
										
										
											2012-02-01 00:48:55 +05:30
+								   write(6,'(1x,a30)')'Returned state: '
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								   do i = 1_pInt,size(stt%relaxationVector(:,of))
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								     write(6,'(1x,2(e15.8,1x))') stt%relaxationVector(i,of)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								   enddo
 								   write(6,*)' '
-												replaced call flush(6) by flush(6), slightly improved formating for output of gfortran version

											
										
										
											2012-11-06 21:20:20 +05:30
+								   flush(6)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								 endif
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#endif
 								 end associate
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
 								 contains
 								 !--------------------------------------------------------------------------------------------------
 								 !> @brief calculate stress-like penalty due to deformation mismatch
 								 !--------------------------------------------------------------------------------------------------
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								 subroutine stressPenalty(rPen,nMis,avgF,fDef,ip,el,instance,of)
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								  use math, only: &
 								    math_civita
 								  use numerics, only: &
 								    xSmoo_RGC
 								  implicit none
 								  real(pReal),   dimension (3,3,homogenization_maxNgrains), intent(out) :: rPen                      !< stress-like penalty
 								  real(pReal),   dimension (3,homogenization_maxNgrains),   intent(out) :: nMis                      !< total amount of mismatch
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								  real(pReal),   dimension (3,3,homogenization_maxNgrains), intent(in)  :: fDef                      !< deformation gradients
 								  real(pReal),   dimension (3,3),                           intent(in)  :: avgF                      !< initial effective stretch tensor
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								  integer(pInt),                                            intent(in)  :: ip,el,instance,of
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								  integer(pInt), dimension (4)   :: intFace
 								  integer(pInt), dimension (3)   :: iGrain3,iGNghb3,nGDim
 								  real(pReal),   dimension (3,3) :: gDef,nDef
 								  real(pReal),   dimension (3)   :: nVect,surfCorr
 								  real(pReal),   dimension (2)   :: Gmoduli
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								  integer(pInt) :: iGrain,iGNghb,iFace,i,j,k,l
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								  real(pReal)   :: muGrain,muGNghb,nDefNorm,bgGrain,bgGNghb
 								  real(pReal),                                               parameter  :: nDefToler = 1.0e-10_pReal
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#ifdef DEBUG
-												cleaning

											
										
										
											2018-11-04 03:53:52 +05:30
+								  logical :: debugActive
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#endif
-												cleaning

											
										
										
											2018-11-04 03:53:52 +05:30
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								  nGDim = param(instance)%Nconstituents
 								  rPen = 0.0_pReal
 								  nMis = 0.0_pReal
 								 !--------------------------------------------------------------------------------------------------
 								 ! get the correction factor the modulus of penalty stress representing the evolution of area of
 								 ! the interfaces due to deformations
-												simplifying

											
										
										
											2018-11-04 03:43:20 +05:30
 								  surfCorr = surfaceCorrection(avgF,instance,of)
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								  associate(prm => param(instance))
-												cleaning

											
										
										
											2018-11-04 03:53:52 +05:30
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#ifdef DEBUG
 								  debugActive = iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								                .and. prm%of_debug == of
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
-												cleaning

											
										
										
											2018-11-04 03:53:52 +05:30
+								  if (debugActive) then
 								    write(6,'(1x,a20,2(1x,i3))')'Correction factor: ',ip,el
 								    write(6,*) surfCorr
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								  endif
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#endif
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
 								 !--------------------------------------------------------------------------------------------------
 								 ! computing the mismatch and penalty stress tensor of all grains
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								  grainLoop: do iGrain = 1_pInt,product(prm%Nconstituents)
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								    Gmoduli = equivalentModuli(iGrain,ip,el)
-												cleaning

											
										
										
											2018-12-18 22:47:06 +05:30
+								    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
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								    interfaceLoop: do iFace = 1_pInt,6_pInt
-												cleaning

											
										
										
											2018-12-18 22:47:06 +05:30
+								      intFace = getInterface(iFace,iGrain3)                                                         ! get the 4-dimensional index of the interface in local numbering system of the grain
-												simplifying

											
										
										
											2018-11-04 03:43:20 +05:30
+								      nVect = interfaceNormal(intFace,instance,of)
-												cleaning

											
										
										
											2018-12-18 22:47:06 +05:30
+								      iGNghb3 = iGrain3                                                                             ! identify the neighboring grain across the interface
-												cleaning

											
										
										
											2018-11-04 02:06:49 +05:30
+								      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
-												cleaning

											
										
										
											2018-12-18 22:47:06 +05:30
+								      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
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								      muGNghb = Gmoduli(1)
 								      bgGNghb = Gmoduli(2)
-												cleaning

											
										
										
											2018-12-18 22:47:06 +05:30
+								      gDef = 0.5_pReal*(fDef(1:3,1:3,iGNghb) - fDef(1:3,1:3,iGrain))                                ! difference/jump in deformation gradeint across the neighbor
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
 								 !--------------------------------------------------------------------------------------------------
 								 ! 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
-												cleaning

											
										
										
											2018-12-18 22:47:06 +05:30
+								          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
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								        enddo; enddo
-												cleaning

											
										
										
											2018-12-18 22:47:06 +05:30
+								        nDefNorm = nDefNorm + nDef(i,j)**2.0_pReal                                                  ! compute the norm of the mismatch tensor
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								      enddo; enddo
-												cleaning

											
										
										
											2018-12-18 22:47:06 +05:30
+								      nDefNorm = max(nDefToler,sqrt(nDefNorm))                                                      ! approximation to zero mismatch if mismatch is zero (singularity)
 								      nMis(abs(intFace(1)),iGrain) = nMis(abs(intFace(1)),iGrain) + nDefNorm                        ! total amount of mismatch experienced by the grain (at all six interfaces)
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#ifdef DEBUG
-												cleaning

											
										
										
											2018-11-04 13:19:40 +05:30
+								      if (debugActive) then
 								        write(6,'(1x,a20,i2,1x,a20,1x,i3)')'Mismatch to face: ',intFace(1),'neighbor grain: ',iGNghb
 								        write(6,*) transpose(nDef)
 								        write(6,'(1x,a20,e11.4)')'with magnitude: ',nDefNorm
 								      endif
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#endif
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								 !--------------------------------------------------------------------------------------------------
 								 ! compute the stress penalty of all interfaces
-												cleaning

											
										
										
											2018-11-04 13:19:40 +05:30
+								      do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt; do k = 1_pInt,3_pInt; do l = 1_pInt,3_pInt
 								        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) &
 								                                               *0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_civita(k,l,j) &
 								                                               *tanh(nDefNorm/xSmoo_RGC)
 								      enddo; enddo;enddo; enddo
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								    enddo interfaceLoop
 								#ifdef DEBUG
-												cleaning

											
										
										
											2018-11-04 03:53:52 +05:30
+								    if (debugActive) then
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								      write(6,'(1x,a20,i2)')'Penalty of grain: ',iGrain
 								      write(6,*) transpose(rPen(1:3,1:3,iGrain))
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								    endif
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#endif
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								  enddo grainLoop
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								  end associate
 								 end subroutine stressPenalty
 								 !--------------------------------------------------------------------------------------------------
 								 !> @brief calculate stress-like penalty due to volume discrepancy
 								 !--------------------------------------------------------------------------------------------------
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								 subroutine volumePenalty(vPen,vDiscrep,fAvg,fDef,nGrain,instance,of)
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								  use math, only: &
 								    math_det33, &
 								    math_inv33
 								  use numerics, only: &
 								    maxVolDiscr_RGC,&
 								    volDiscrMod_RGC,&
 								    volDiscrPow_RGC
 								  implicit none
 								  real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: vPen                        ! stress-like penalty due to volume
 								  real(pReal), intent(out)                  :: vDiscrep                                              ! total volume discrepancy
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								  real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in)  :: fDef                        ! deformation gradients
 								  real(pReal), dimension (3,3), intent(in)  :: fAvg                                                  ! overall deformation gradient
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								  integer(pInt), intent(in) :: &
 								    Ngrain, &
 								    instance, &
 								    of
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								  real(pReal), dimension (homogenization_maxNgrains) :: gVol
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								  integer(pInt) :: i
-												cleaning

											
										
										
											2018-11-04 03:53:52 +05:30
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								 !--------------------------------------------------------------------------------------------------
 								 ! compute the volumes of grains and of cluster
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								  vDiscrep = math_det33(fAvg)                                                                       ! compute the volume of the cluster
 								  do i = 1_pInt,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
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								  enddo
 								 !--------------------------------------------------------------------------------------------------
 								 ! calculate the stress and penalty due to volume discrepancy
 								  vPen      = 0.0_pReal
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								  do i = 1_pInt,nGrain
 								    vPen(:,:,i) = -1.0_pReal/real(nGrain,pReal)*volDiscrMod_RGC*volDiscrPow_RGC/maxVolDiscr_RGC* &
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								                       sign((abs(vDiscrep)/maxVolDiscr_RGC)**(volDiscrPow_RGC - 1.0),vDiscrep)* &
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								                       gVol(i)*transpose(math_inv33(fDef(:,:,i)))
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#ifdef DEBUG
 								    if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
 								                .and. param(instance)%of_debug == of) then
 								      write(6,'(1x,a30,i2)')'Volume penalty of grain: ',i
 								      write(6,*) transpose(vPen(:,:,i))
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								    endif
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								#endif
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								  enddo
-												cleaning

											
										
										
											2018-11-04 02:06:49 +05:30
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								 end subroutine volumePenalty
-												cleaning

											
										
										
											2018-11-04 02:06:49 +05:30
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								 !--------------------------------------------------------------------------------------------------
 								 !> @brief compute the correction factor accouted for surface evolution (area change) due to
 								 ! deformation
 								 !--------------------------------------------------------------------------------------------------
-												simplifying

											
										
										
											2018-11-04 03:43:20 +05:30
+								 function surfaceCorrection(avgF,instance,of)
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								  use math, only: &
 								    math_invert33, &
 								    math_mul33x33
 								  implicit none
 								  real(pReal), dimension(3)               :: surfaceCorrection
 								  real(pReal), dimension(3,3), intent(in) :: avgF                                                    !< average F
-												simplifying

											
										
										
											2018-11-04 03:43:20 +05:30
+								  integer(pInt),               intent(in) :: &
 								    instance, &
 								    of
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								  real(pReal), dimension(3,3)             :: invC
 								  real(pReal), dimension(3)               :: nVect
 								  real(pReal)  :: detF
 								  integer(pInt) :: i,j,iBase
 								  logical     ::  error
 								  call math_invert33(math_mul33x33(transpose(avgF),avgF),invC,detF,error)
 								  surfaceCorrection = 0.0_pReal
 								  do iBase = 1_pInt,3_pInt
-												simplifying

											
										
										
											2018-11-04 03:43:20 +05:30
+								    nVect = interfaceNormal([iBase,1_pInt,1_pInt,1_pInt],instance,of)
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								    do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								      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
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								    enddo; enddo
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								    surfaceCorrection(iBase) = sqrt(surfaceCorrection(iBase))*detF                                  ! get the surface correction factor (area contraction/enlargement)
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								  enddo
 								 end function surfaceCorrection
-												cleaning

											
										
										
											2018-11-04 02:06:49 +05:30
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								 !--------------------------------------------------------------------------------------------------
 								 !> @brief compute the equivalent shear and bulk moduli from the elasticity tensor
 								 !--------------------------------------------------------------------------------------------------
 								 function equivalentModuli(grainID,ip,el)
 								  use constitutive, only: &
 								    constitutive_homogenizedC
 								  implicit none
-												simplified

											
										
										
											2018-11-04 12:41:35 +05:30
+								  real(pReal), dimension(2)    :: equivalentModuli
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								  integer(pInt), intent(in)    :: &
 								    grainID,&
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								    ip, &                                                                                           !< integration point number
 								    el                                                                                              !< element number
-												simplified

											
										
										
											2018-11-04 12:41:35 +05:30
+								  real(pReal), dimension(6,6) :: elasTens
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								  real(pReal) :: &
 								    cEquiv_11, &
 								    cEquiv_12, &
 								    cEquiv_44
 								  elasTens = constitutive_homogenizedC(grainID,ip,el)
 								 !--------------------------------------------------------------------------------------------------
 								 ! compute the equivalent shear modulus after 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
 								  equivalentModuli(1) = 0.2_pReal*(cEquiv_11 - cEquiv_12) + 0.6_pReal*cEquiv_44
 								 !--------------------------------------------------------------------------------------------------
 								 ! obtain the length of Burgers vector (could be model dependend)
 								  equivalentModuli(2) = 2.5e-10_pReal
 								 end function equivalentModuli
 								 !--------------------------------------------------------------------------------------------------
 								 !> @brief calculating the grain deformation gradient (the same with
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								 ! homogenization_RGC_partitionDeformation, but used only for perturbation scheme)
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								 !--------------------------------------------------------------------------------------------------
-												fixed array out of bounds during initialization

happened for inactive homogenization

											
										
										
											2018-11-04 11:57:25 +05:30
+								 subroutine grainDeformation(F, avgF, instance, of)
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
 								  implicit none
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								  real(pReal),   dimension (:,:,:), intent(out) :: F                                                 !< partioned F  per grain
 								  real(pReal),   dimension (:,:),   intent(in)  :: avgF                                              !< averaged F
 								  integer(pInt),                    intent(in)  :: &
-												fixed array out of bounds during initialization

happened for inactive homogenization

											
										
										
											2018-11-04 11:57:25 +05:30
+								    instance, &
 								    of
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								  real(pReal),   dimension (3) :: aVect,nVect
 								  integer(pInt), dimension (4) :: intFace
 								  integer(pInt), dimension (3) :: iGrain3
-												fixed array out of bounds during initialization

happened for inactive homogenization

											
										
										
											2018-11-04 11:57:25 +05:30
+								  integer(pInt) :: iGrain,iFace,i,j
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
+								  !-------------------------------------------------------------------------------------------------
-												simplifying

											
										
										
											2018-11-04 03:43:20 +05:30
+								  ! compute the deformation gradient of individual grains due to relaxations
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
 								  associate(prm => param(instance))
-												fixed array out of bounds during initialization

happened for inactive homogenization

											
										
										
											2018-11-04 11:57:25 +05:30
+								  F = 0.0_pReal
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								  do iGrain = 1_pInt,product(prm%Nconstituents)
 								    iGrain3 = grain1to3(iGrain,prm%Nconstituents)
-												fixed array out of bounds during initialization

happened for inactive homogenization

											
										
										
											2018-11-04 11:57:25 +05:30
+								    do iFace = 1_pInt,6_pInt
 								      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) &
 								        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
 								  enddo
-												simplifying

											
										
										
											2018-11-04 03:43:20 +05:30
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								  end associate
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
+								 end subroutine grainDeformation
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								end function homogenization_RGC_updateState
 								!--------------------------------------------------------------------------------------------------
 								!> @brief derive average stress and stiffness from constituent quantities
 								!--------------------------------------------------------------------------------------------------
-												misplaced debug statement

debugging the individual constituents should be done by the constitutive
laws

											
										
										
											2018-11-04 00:38:18 +05:30
+								subroutine homogenization_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								 implicit none
-												removed temperature integration and corresponding data structures and debugging options
temperature is stored in crystallite, but homogeneous on one IP (not an component (grain) quantity and an input value parsed in by the BVP solver.
introduced heat, a component (grain) quantity which is homogenized before returned to the heat transfer solver.
went ahead with removal of dummy functions in homogenization and constitutive, this time mainly reduced function signatures to reflect actually needed quantities.


											
										
										
											2013-10-16 18:34:59 +05:30
+								 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
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
 								 real(pReal), dimension (:,:,:),                             intent(in)  :: P                       !< partitioned stresses
 								 real(pReal), dimension (:,:,:,:,:),                         intent(in)  :: dPdF                    !< partitioned stiffnesses
-												misplaced debug statement

debugging the individual constituents should be done by the constitutive
laws

											
										
										
											2018-11-04 00:38:18 +05:30
+								 integer(pInt),                                              intent(in)  :: instance
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												misplaced debug statement

debugging the individual constituents should be done by the constitutive
laws

											
										
										
											2018-11-04 00:38:18 +05:30
+								 avgP       = sum(P,3)   /real(product(param(instance)%Nconstituents),pReal)
 								 dAvgPdAvgF = sum(dPdF,5)/real(product(param(instance)%Nconstituents),pReal)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												set useful default values to remove dummy functions in individual homogenization schemes (stateInit, averageTemperature, and updateState (isostrain only))

											
										
										
											2013-10-11 21:31:53 +05:30
+								end subroutine homogenization_RGC_averageStressAndItsTangent
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
 								!> @brief return array of homogenization results for post file inclusion
 								!--------------------------------------------------------------------------------------------------
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								pure function homogenization_RGC_postResults(instance,of) result(postResults)
-												limit size to actual Ngrains, not potential maximum

											
										
										
											2019-01-13 13:44:23 +05:30
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								 implicit none
-												introduced dummy heat calculation, available as an crystallite output
temperature, avgP and avgF are available as homogenization output. move crystallite output of ipcoordinates to homogenization

											
										
										
											2013-10-19 00:27:28 +05:30
+								 integer(pInt), intent(in) :: &
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								   instance, &
 								   of
 								 integer(pInt) :: &
 								   o,c
 								 real(pReal), dimension(sum(homogenization_RGC_sizePostResult(:,instance))) :: &
-												no need to have long prefixes for local variables and functions

											
										
										
											2018-08-30 14:20:52 +05:30
+								   postResults
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												no need to be real state

											
										
										
											2019-01-13 03:52:13 +05:30
+								 associate(stt => state(instance), dst => dependentState(instance), prm => param(instance))
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
 								 c = 0_pInt
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
-												requested output is stored in prm%outputID

											
										
										
											2018-08-31 10:05:30 +05:30
+								 outputsLoop: do o = 1_pInt,size(prm%outputID)
 								   select case(prm%outputID(o))
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
-												homogenization: added enums and sourced allocation for RGC, some higher verbosity for isostrain

											
										
										
											2013-12-18 12:58:01 +05:30
+								     case (constitutivework_ID)
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								       postResults(c+1) = stt%work(of)
-												again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

											
										
										
											2012-02-21 22:01:37 +05:30
+								       c = c + 1_pInt
-												homogenization: added enums and sourced allocation for RGC, some higher verbosity for isostrain

											
										
										
											2013-12-18 12:58:01 +05:30
+								     case (magnitudemismatch_ID)
-												no need to be real state

											
										
										
											2019-01-13 03:52:13 +05:30
+								       postResults(c+1:c+3) = dst%mismatch(1:3,of)
-												again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

											
										
										
											2012-02-21 22:01:37 +05:30
+								       c = c + 3_pInt
-												homogenization: added enums and sourced allocation for RGC, some higher verbosity for isostrain

											
										
										
											2013-12-18 12:58:01 +05:30
+								     case (penaltyenergy_ID)
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								       postResults(c+1) = stt%penaltyEnergy(of)
-												started introducing new state structure in homogenisation

											
										
										
											2014-08-21 23:18:20 +05:30
+								       c = c + 1_pInt
-												homogenization: added enums and sourced allocation for RGC, some higher verbosity for isostrain

											
										
										
											2013-12-18 12:58:01 +05:30
+								     case (volumediscrepancy_ID)
-												no need to be real state

											
										
										
											2019-01-13 03:52:13 +05:30
+								       postResults(c+1) = dst%volumeDiscrepancy(of)
-												started introducing new state structure in homogenisation

											
										
										
											2014-08-21 23:18:20 +05:30
+								       c = c + 1_pInt
-												homogenization: added enums and sourced allocation for RGC, some higher verbosity for isostrain

											
										
										
											2013-12-18 12:58:01 +05:30
+								     case (averagerelaxrate_ID)
-												no need to be real state

											
										
										
											2019-01-13 03:52:13 +05:30
+								       postResults(c+1) = dst%relaxationrate_avg(of)
-												started introducing new state structure in homogenisation

											
										
										
											2014-08-21 23:18:20 +05:30
+								       c = c + 1_pInt
-												homogenization: added enums and sourced allocation for RGC, some higher verbosity for isostrain

											
										
										
											2013-12-18 12:58:01 +05:30
+								     case (maximumrelaxrate_ID)
-												no need to be real state

											
										
										
											2019-01-13 03:52:13 +05:30
+								       postResults(c+1) = dst%relaxationrate_max(of)
-												started introducing new state structure in homogenisation

											
										
										
											2014-08-21 23:18:20 +05:30
+								       c = c + 1_pInt
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								   end select
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
-												requested output is stored in prm%outputID

											
										
										
											2018-08-31 10:05:30 +05:30
+								 enddo outputsLoop
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
-												consistent renames

											
										
										
											2018-08-30 14:16:30 +05:30
+								 end associate
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								end function homogenization_RGC_postResults
 								!--------------------------------------------------------------------------------------------------
 								!> @brief collect relaxation vectors of an interface
 								!--------------------------------------------------------------------------------------------------
-												these functions have no side effects

											
										
										
											2018-11-04 12:03:57 +05:30
+								pure function relaxationVector(intFace,instance,of)
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								 implicit none
-												simplifying

											
										
										
											2018-11-04 03:43:20 +05:30
+								 integer(pInt),                intent(in) :: instance,of
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
-												no need to have long prefixes for local variables and functions

											
										
										
											2018-08-30 14:20:52 +05:30
+								 real(pReal),   dimension (3)             :: relaxationVector
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								 integer(pInt), dimension (4), intent(in) :: intFace                                                !< set of interface ID in 4D array (normal and position)
-												polished a little bit

											
										
										
											2014-02-08 16:18:09 +05:30
+								 integer(pInt) :: &
-												simplifying

											
										
										
											2018-11-04 03:43:20 +05:30
+								   iNum
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
 								! collect the interface relaxation vector from the global state array
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
-												simplified

											
										
										
											2018-11-04 12:26:27 +05:30
+								 iNum = interface4to1(intFace,param(instance)%Nconstituents)                                        ! identify the position of the interface in global state array
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								 if (iNum > 0_pInt) then
 								   relaxationVector = state(instance)%relaxationVector((3*iNum-2):(3*iNum),of)
 								 else
 								   relaxationVector = 0.0_pReal
 								 endif
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
-												no need to have long prefixes for local variables and functions

											
										
										
											2018-08-30 14:20:52 +05:30
+								end function relaxationVector
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
 								!--------------------------------------------------------------------------------------------------
 								!> @brief identify the normal of an interface
 								!--------------------------------------------------------------------------------------------------
-												these functions have no side effects

											
										
										
											2018-11-04 12:03:57 +05:30
+								pure function interfaceNormal(intFace,instance,of)
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								 use math, only: &
 								   math_mul33x3
-												added new, flexible debugging scheme.
now all modules have their own debug specification.
compiles and runs, I hope nothing is broken
did a lot of polishing

											
										
										
											2012-03-09 01:55:28 +05:30
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								 implicit none
-												no need to have long prefixes for local variables and functions

											
										
										
											2018-08-30 14:20:52 +05:30
+								 real(pReal),   dimension (3)  ::            interfaceNormal
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								 integer(pInt), dimension (4), intent(in) :: intFace                                                !< interface ID in 4D array (normal and position)
-												polished a little bit

											
										
										
											2014-02-08 16:18:09 +05:30
+								 integer(pInt),                intent(in) :: &
-												simplifying

											
										
										
											2018-11-04 03:43:20 +05:30
+								   instance, &
 								   of
 								 integer(pInt) ::                            nPos
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
 								! get the normal of the interface, identified from the value of intFace(1)
-												no need to have long prefixes for local variables and functions

											
										
										
											2018-08-30 14:20:52 +05:30
+								 interfaceNormal = 0.0_pReal
-												polishing: removed variable names like 'unit' and 'data' that are keywords of fortran and ensured that integer and real precision matches independent of machine standard.
removed cut_off parameter for damask_spectral
removed outpot of derived divergence measures and added RMS output in brackets
added comments and options to the makefile

											
										
										
											2012-02-16 00:28:38 +05:30
+								 nPos = abs(intFace(1))                                                                             ! identify the position of the interface in global state array
-												using parameters from config.f90

											
										
										
											2018-08-30 21:21:31 +05:30
+								 interfaceNormal(nPos) = real(intFace(1)/abs(intFace(1)),pReal)                                     ! get the normal vector w.r.t. cluster axis
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												cleaning

											
										
										
											2018-11-04 03:53:52 +05:30
+								 interfaceNormal = math_mul33x3(dependentState(instance)%orientation(1:3,1:3,of),interfaceNormal)   ! map the normal vector into sample coordinate system (basis)
-												the latest RGC model + corrections for "element homogeneous" feature

											
										
										
											2010-03-24 18:50:12 +05:30
-												no need to have long prefixes for local variables and functions

											
										
										
											2018-08-30 14:20:52 +05:30
+								end function interfaceNormal
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
 								!--------------------------------------------------------------------------------------------------
 								!> @brief collect six faces of a grain in 4D (normal and position)
 								!--------------------------------------------------------------------------------------------------
-												these functions have no side effects

											
										
										
											2018-11-04 12:03:57 +05:30
+								pure function getInterface(iFace,iGrain3)
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								 implicit none
-												no need to have long prefixes for local variables and functions

											
										
										
											2018-08-30 14:20:52 +05:30
+								 integer(pInt), dimension (4) ::             getInterface
-												polished a little bit

											
										
										
											2014-02-08 16:18:09 +05:30
+								 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
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
 								!* Direction of interface normal
-												again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

											
										
										
											2012-02-21 22:01:37 +05:30
+								 iDir = (int(real(iFace-1_pInt,pReal)/2.0_pReal,pInt)+1_pInt)*(-1_pInt)**iFace
-												no need to have long prefixes for local variables and functions

											
										
										
											2018-08-30 14:20:52 +05:30
+								 getInterface(1) = iDir
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
 								! identify the interface position by the direction of its normal
-												no need to have long prefixes for local variables and functions

											
										
										
											2018-08-30 14:20:52 +05:30
+								 getInterface(2:4) = iGrain3
-												these functions have no side effects

											
										
										
											2018-11-04 12:03:57 +05:30
+								 if (iDir < 0_pInt) getInterface(1_pInt-iDir) = getInterface(1_pInt-iDir)-1_pInt                    ! to have a correlation with coordinate/position in real space
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												no need to have long prefixes for local variables and functions

											
										
										
											2018-08-30 14:20:52 +05:30
+								end function getInterface
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												simplified

											
										
										
											2018-11-04 12:41:35 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
 								!> @brief map grain ID from in 1D (global array) to in 3D (local position)
 								!--------------------------------------------------------------------------------------------------
-												simplified

											
										
										
											2018-11-04 12:41:35 +05:30
+								pure function grain1to3(grain1,nGDim)
-												added new, flexible debugging scheme.
now all modules have their own debug specification.
compiles and runs, I hope nothing is broken
did a lot of polishing

											
										
										
											2012-03-09 01:55:28 +05:30
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								 implicit none
-												simplified

											
										
										
											2018-11-04 12:41:35 +05:30
+								 integer(pInt), dimension(3)             :: grain1to3
 								 integer(pInt),               intent(in) :: grain1                                                  !< grain ID in 1D array
 								 integer(pInt), dimension(3), intent(in) :: nGDim
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												using parameters from config.f90

											
										
										
											2018-08-30 21:21:31 +05:30
+								 grain1to3 = 1_pInt + [mod((grain1-1_pInt),nGDim(1)), &
 								                       mod((grain1-1_pInt)/nGDim(1),nGDim(2)), &
 								                       (grain1-1_pInt)/(nGDim(1)*nGDim(2))]
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												no need to have long prefixes for local variables and functions

											
										
										
											2018-08-30 14:20:52 +05:30
+								end function grain1to3
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
 								!> @brief map grain ID from in 3D (local position) to in 1D (global array)
 								!--------------------------------------------------------------------------------------------------
-												better readable

											
										
										
											2018-11-04 13:49:24 +05:30
+								integer(pInt) pure function grain3to1(grain3,nGDim)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												added new, flexible debugging scheme.
now all modules have their own debug specification.
compiles and runs, I hope nothing is broken
did a lot of polishing

											
										
										
											2012-03-09 01:55:28 +05:30
+								 implicit none
-												simplified

											
										
										
											2018-11-04 12:41:35 +05:30
+								 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
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												simplified

											
										
										
											2018-11-04 12:41:35 +05:30
+								 grain3to1 = grain3(1) &
 								           + nGDim(1)*(grain3(2)-1_pInt) &
 								           + nGDim(1)*nGDim(2)*(grain3(3)-1_pInt)
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												no need to have long prefixes for local variables and functions

											
										
										
											2018-08-30 14:20:52 +05:30
+								end function grain3to1
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
 								!> @brief maps interface ID from 4D (normal and local position) into 1D (global array)
 								!--------------------------------------------------------------------------------------------------
-												simplified

											
										
										
											2018-11-04 12:26:27 +05:30
+								integer(pInt) pure function interface4to1(iFace4D, nGDim)
-												added new, flexible debugging scheme.
now all modules have their own debug specification.
compiles and runs, I hope nothing is broken
did a lot of polishing

											
										
										
											2012-03-09 01:55:28 +05:30
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								 implicit none
-												simplified

											
										
										
											2018-11-04 12:26:27 +05:30
+								 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
-. Adding more comments
2. Changing 'subroutine' into 'function' for optimizing the computation time

											
										
										
											2010-11-26 17:20:20 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
-												better readable

											
										
										
											2018-11-04 13:49:24 +05:30
+								 select case(abs(iFace4D(1)))
 								   case(1_pInt)
 								     if ((iFace4D(2) == 0_pInt) .or. (iFace4D(2) == nGDim(1))) then
 								       interface4to1 = 0_pInt
 								     else
 								       interface4to1 = iFace4D(3) + nGDim(2)*(iFace4D(4)-1_pInt) &
 								                     + nGDim(2)*nGDim(3)*(iFace4D(2)-1_pInt)
 								     endif
 								   case(2_pInt)
 								     if ((iFace4D(3) == 0_pInt) .or. (iFace4D(3) == nGDim(2))) then
 								       interface4to1 = 0_pInt
 								     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
 								     endif
 								   case(3_pInt)
 								     if ((iFace4D(4) == 0_pInt) .or. (iFace4D(4) == nGDim(3))) then
 								       interface4to1 = 0_pInt
 								     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
 								     endif
 								   case default
 								     interface4to1 = -1_pInt
 								 end select
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												no need to have long prefixes for local variables and functions

											
										
										
											2018-08-30 14:20:52 +05:30
+								end function interface4to1
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
 								!> @brief maps interface ID from 1D (global array) into 4D (normal and local position)
 								!--------------------------------------------------------------------------------------------------
-												simplified

											
										
										
											2018-11-04 12:41:35 +05:30
+								pure function interface1to4(iFace1D, nGDim)
-												added new, flexible debugging scheme.
now all modules have their own debug specification.
compiles and runs, I hope nothing is broken
did a lot of polishing

											
										
										
											2012-03-09 01:55:28 +05:30
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								 implicit none
-												simplified

											
										
										
											2018-11-04 12:41:35 +05:30
+								 integer(pInt), dimension(4)             :: interface1to4
-												using parameters from config.f90

											
										
										
											2018-08-30 21:21:31 +05:30
+								 integer(pInt),               intent(in) :: iFace1D                                                            !< interface ID in 1D array
-												simplified

											
										
										
											2018-11-04 12:41:35 +05:30
+								 integer(pInt), dimension(3), intent(in) :: nGDim
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								 integer(pInt), dimension(3)             :: nIntFace
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
 								! compute the total number of interfaces, which ...
-												better readable

											
										
										
											2018-11-04 13:49:24 +05:30
+								 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
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
 								!--------------------------------------------------------------------------------------------------
 								! get the corresponding interface ID in 4D (normal and local position)
 								 if (iFace1D > 0 .and. iFace1D <= nIntFace(1)) then                                                 ! interface with normal //e1
-												no need to have long prefixes for local variables and functions

											
										
										
											2018-08-30 14:20:52 +05:30
+								   interface1to4(1) = 1_pInt
 								   interface1to4(3) = mod((iFace1D-1_pInt),nGDim(2))+1_pInt
 								   interface1to4(4) = mod(&
-												again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

											
										
										
											2012-02-21 22:01:37 +05:30
+								                                             int(&
 								                                                 real(iFace1D-1_pInt,pReal)/&
 								                                                 real(nGDim(2),pReal)&
 								                                                 ,pInt)&
 								                                             ,nGDim(3))+1_pInt
-												no need to have long prefixes for local variables and functions

											
										
										
											2018-08-30 14:20:52 +05:30
+								   interface1to4(2) = int(&
-												again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

											
										
										
											2012-02-21 22:01:37 +05:30
+								                                             real(iFace1D-1_pInt,pReal)/&
 								                                             real(nGDim(2),pReal)/&
 								                                             real(nGDim(3),pReal)&
 								                                             ,pInt)+1_pInt
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								 elseif (iFace1D > nIntFace(1) .and. iFace1D <= (nIntFace(2) + nIntFace(1))) then                   ! interface with normal //e2
-												no need to have long prefixes for local variables and functions

											
										
										
											2018-08-30 14:20:52 +05:30
+								   interface1to4(1) = 2_pInt
 								   interface1to4(4) = mod((iFace1D-nIntFace(1)-1_pInt),nGDim(3))+1_pInt
 								   interface1to4(2) = mod(&
-												again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

											
										
										
											2012-02-21 22:01:37 +05:30
+								                                             int(&
 								                                                 real(iFace1D-nIntFace(1)-1_pInt,pReal)/&
 								                                                 real(nGDim(3),pReal)&
 								                                                 ,pInt)&
 								                                              ,nGDim(1))+1_pInt
-												no need to have long prefixes for local variables and functions

											
										
										
											2018-08-30 14:20:52 +05:30
+								   interface1to4(3) = int(&
-												again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

											
										
										
											2012-02-21 22:01:37 +05:30
+								                                             real(iFace1D-nIntFace(1)-1_pInt,pReal)/&
 								                                             real(nGDim(3),pReal)/&
 								                                             real(nGDim(1),pReal)&
 								                                             ,pInt)+1_pInt
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								 elseif (iFace1D > nIntFace(2) + nIntFace(1) .and. iFace1D <= (nIntFace(3) + nIntFace(2) + nIntFace(1))) then ! interface with normal //e3
-												no need to have long prefixes for local variables and functions

											
										
										
											2018-08-30 14:20:52 +05:30
+								   interface1to4(1) = 3_pInt
 								   interface1to4(2) = mod((iFace1D-nIntFace(2)-nIntFace(1)-1_pInt),nGDim(1))+1_pInt
 								   interface1to4(3) = mod(&
-												again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

											
										
										
											2012-02-21 22:01:37 +05:30
+								                                             int(&
 								                                                 real(iFace1D-nIntFace(2)-nIntFace(1)-1_pInt,pReal)/&
 								                                                 real(nGDim(1),pReal)&
 								                                                 ,pInt)&
 								                                             ,nGDim(2))+1_pInt
-												no need to have long prefixes for local variables and functions

											
										
										
											2018-08-30 14:20:52 +05:30
+								   interface1to4(4) = int(&
-												again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

											
										
										
											2012-02-21 22:01:37 +05:30
+								                                             real(iFace1D-nIntFace(2)-nIntFace(1)-1_pInt,pReal)/&
 								                                             real(nGDim(1),pReal)/&
 								                                             real(nGDim(2),pReal)&
 								                                             ,pInt)+1_pInt
-												# 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...


											
										
										
											2009-10-12 21:31:49 +05:30
+								 endif
-												no need to have long prefixes for local variables and functions

											
										
										
											2018-08-30 14:20:52 +05:30
+								end function interface1to4
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
 								end module homogenization_RGC