DAMASK_EICMD/src/homogenization_mechanical_R...

!--------------------------------------------------------------------------------------------------
!> @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
!> N_constituents is defined as p x q x r (cluster)
!--------------------------------------------------------------------------------------------------
submodule(homogenization:mechanical) RGC
  use rotations
  use lattice

  type :: tParameters
    integer, dimension(:), allocatable :: &
      N_constituents
    real(pReal) :: &
      xi_alpha, &
      c_Alpha
    real(pReal), dimension(:), allocatable :: &
      D_alpha, &
      a_g
    character(len=pStringLen), allocatable, dimension(:) :: &
      output
  end type tParameters

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

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

  type :: tNumerics_RGC
    real(pReal) :: &
      atol, &                                                                                       !< absolute tolerance of RGC residuum
      rtol, &                                                                                       !< relative tolerance of RGC residuum
      absMax, &                                                                                     !< absolute maximum of RGC residuum
      relMax, &                                                                                     !< relative maximum of RGC residuum
      pPert, &                                                                                      !< perturbation for computing RGC penalty tangent
      xSmoo, &                                                                                      !< RGC penalty smoothing parameter (hyperbolic tangent)
      viscPower, &                                                                                  !< power (sensitivity rate) of numerical viscosity in RGC scheme, Default 1.0e0: Newton viscosity (linear model)
      viscModus, &                                                                                  !< stress modulus of RGC numerical viscosity, Default 0.0e0: No viscosity is applied
      refRelaxRate, &                                                                               !< reference relaxation rate in RGC viscosity
      maxdRelax, &                                                                                  !< threshold of maximum relaxation vector increment (if exceed this then cutback)
      maxVolDiscr, &                                                                                !< threshold of maximum volume discrepancy allowed
      volDiscrMod, &                                                                                !< stiffness of RGC volume discrepancy (zero = without volume discrepancy constraint)
      volDiscrPow                                                                                   !< powerlaw penalty for volume discrepancy
  end type tNumerics_RGC

  type(tparameters),          dimension(:), allocatable :: &
    param
  type(tRGCstate),            dimension(:), allocatable :: &
    state, &
    state0
  type(tRGCdependentState),   dimension(:), allocatable :: &
    dependentState
  type(tNumerics_RGC) :: &
    num                                                                                             ! numerics parameters. Better name?

contains

!--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields, reads information from material configuration file
!--------------------------------------------------------------------------------------------------
module subroutine RGC_init(num_homogMech)

  class(tNode), pointer, intent(in) :: &
    num_homogMech                                                                                   !< pointer to mechanical homogenization numerics data

  integer :: &
    ho, &
    Nmembers, &
    sizeState, nIntFaceTot

  class (tNode), pointer :: &
    num_RGC, &                                                                                      ! pointer to RGC numerics data
    material_homogenization, &
    homog, &
    homogMech

  print'(/,1x,a)', '<<<+-  homogenization:mechanical:RGC init  -+>>>'

  print'(/,a,i0)', ' # homogenizations: ',count(homogenization_type == HOMOGENIZATION_RGC_ID)
  flush(IO_STDOUT)

  print'(/,1x,a)', 'D.D. Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009'
  print'(  1x,a)', 'https://doi.org/10.1007/s12289-009-0619-1'//IO_EOL

  print'(/,1x,a)', 'D.D. Tjahjanto et al., Modelling and Simulation in Materials Science and Engineering 18:015006, 2010'
  print'(  1x,a)', 'https://doi.org/10.1088/0965-0393/18/1/015006'//IO_EOL


  material_homogenization => config_material%get('homogenization')
  allocate(param(material_homogenization%length))
  allocate(state(material_homogenization%length))
  allocate(state0(material_homogenization%length))
  allocate(dependentState(material_homogenization%length))

  num_RGC => num_homogMech%get('RGC',defaultVal=emptyDict)

  num%atol         =  num_RGC%get_asFloat('atol',              defaultVal=1.0e+4_pReal)
  num%rtol         =  num_RGC%get_asFloat('rtol',              defaultVal=1.0e-3_pReal)
  num%absMax       =  num_RGC%get_asFloat('amax',              defaultVal=1.0e+10_pReal)
  num%relMax       =  num_RGC%get_asFloat('rmax',              defaultVal=1.0e+2_pReal)
  num%pPert        =  num_RGC%get_asFloat('perturbpenalty',    defaultVal=1.0e-7_pReal)
  num%xSmoo        =  num_RGC%get_asFloat('relvantmismatch',   defaultVal=1.0e-5_pReal)
  num%viscPower    =  num_RGC%get_asFloat('viscositypower',    defaultVal=1.0e+0_pReal)
  num%viscModus    =  num_RGC%get_asFloat('viscositymodulus',  defaultVal=0.0e+0_pReal)
  num%refRelaxRate =  num_RGC%get_asFloat('refrelaxationrate', defaultVal=1.0e-3_pReal)
  num%maxdRelax    =  num_RGC%get_asFloat('maxrelaxationrate', defaultVal=1.0e+0_pReal)
  num%maxVolDiscr  =  num_RGC%get_asFloat('maxvoldiscrepancy', defaultVal=1.0e-5_pReal)
  num%volDiscrMod  =  num_RGC%get_asFloat('voldiscrepancymod', defaultVal=1.0e+12_pReal)
  num%volDiscrPow  =  num_RGC%get_asFloat('dicrepancypower',   defaultVal=5.0_pReal)

  if (num%atol <= 0.0_pReal)         call IO_error(301,ext_msg='absTol_RGC')
  if (num%rtol <= 0.0_pReal)         call IO_error(301,ext_msg='relTol_RGC')
  if (num%absMax <= 0.0_pReal)       call IO_error(301,ext_msg='absMax_RGC')
  if (num%relMax <= 0.0_pReal)       call IO_error(301,ext_msg='relMax_RGC')
  if (num%pPert <= 0.0_pReal)        call IO_error(301,ext_msg='pPert_RGC')
  if (num%xSmoo <= 0.0_pReal)        call IO_error(301,ext_msg='xSmoo_RGC')
  if (num%viscPower < 0.0_pReal)     call IO_error(301,ext_msg='viscPower_RGC')
  if (num%viscModus < 0.0_pReal)     call IO_error(301,ext_msg='viscModus_RGC')
  if (num%refRelaxRate <= 0.0_pReal) call IO_error(301,ext_msg='refRelaxRate_RGC')
  if (num%maxdRelax <= 0.0_pReal)    call IO_error(301,ext_msg='maxdRelax_RGC')
  if (num%maxVolDiscr <= 0.0_pReal)  call IO_error(301,ext_msg='maxVolDiscr_RGC')
  if (num%volDiscrMod < 0.0_pReal)   call IO_error(301,ext_msg='volDiscrMod_RGC')
  if (num%volDiscrPow <= 0.0_pReal)  call IO_error(301,ext_msg='volDiscrPw_RGC')


  do ho = 1, size(homogenization_type)
    if (homogenization_type(ho) /= HOMOGENIZATION_RGC_ID) cycle
    homog => material_homogenization%get(ho)
    homogMech => homog%get('mechanical')
    associate(prm => param(ho), &
              stt => state(ho), &
              st0 => state0(ho), &
              dst => dependentState(ho))

#if defined (__GFORTRAN__)
    prm%output = output_as1dString(homogMech)
#else
    prm%output = homogMech%get_as1dString('output',defaultVal=emptyStringArray)
#endif

    prm%N_constituents = homogMech%get_as1dInt('cluster_size',requiredSize=3)
    if (homogenization_Nconstituents(ho) /= product(prm%N_constituents)) &
      call IO_error(211,ext_msg='N_constituents (RGC)')

    prm%xi_alpha = homogMech%get_asFloat('xi_alpha')
    prm%c_alpha  = homogMech%get_asFloat('c_alpha')

    prm%D_alpha  = homogMech%get_as1dFloat('D_alpha', requiredSize=3)
    prm%a_g      = homogMech%get_as1dFloat('a_g',     requiredSize=3)

    Nmembers = count(material_homogenizationID == ho)
    nIntFaceTot = 3*(  (prm%N_constituents(1)-1)*prm%N_constituents(2)*prm%N_constituents(3) &
                      + prm%N_constituents(1)*(prm%N_constituents(2)-1)*prm%N_constituents(3) &
                      + prm%N_constituents(1)*prm%N_constituents(2)*(prm%N_constituents(3)-1))
    sizeState = nIntFaceTot

    homogState(ho)%sizeState = sizeState
    allocate(homogState(ho)%state0   (sizeState,Nmembers), source=0.0_pReal)
    allocate(homogState(ho)%state    (sizeState,Nmembers), source=0.0_pReal)

    stt%relaxationVector   => homogState(ho)%state(1:nIntFaceTot,:)
    st0%relaxationVector   => homogState(ho)%state0(1:nIntFaceTot,:)

    allocate(dst%volumeDiscrepancy(   Nmembers), source=0.0_pReal)
    allocate(dst%relaxationRate_avg(  Nmembers), source=0.0_pReal)
    allocate(dst%relaxationRate_max(  Nmembers), source=0.0_pReal)
    allocate(dst%mismatch(          3,Nmembers), source=0.0_pReal)

!--------------------------------------------------------------------------------------------------
! assigning cluster orientations
    dependentState(ho)%orientation = spread(eu2om(prm%a_g*inRad),3,Nmembers)
    !dst%orientation = spread(eu2om(prm%a_g*inRad),3,Nmembers) ifort version 18.0.1 crashes (for whatever reason)

    end associate

  end do

end subroutine RGC_init


!--------------------------------------------------------------------------------------------------
!> @brief partitions the deformation gradient onto the constituents
!--------------------------------------------------------------------------------------------------
module subroutine RGC_partitionDeformation(F,avgF,ce)

  real(pReal),   dimension (:,:,:), intent(out) :: F                                                !< partitioned F  per grain

  real(pReal),   dimension (3,3),   intent(in)  :: avgF                                             !< averaged F
  integer,                          intent(in)  :: &
    ce

  real(pReal), dimension(3) :: aVect,nVect
  integer,     dimension(4) :: intFace
  integer,     dimension(3) :: iGrain3
  integer ::  iGrain,iFace,i,j,ho,en

  associate(prm => param(material_homogenizationID(ce)))

  ho = material_homogenizationID(ce)
  en = material_homogenizationEntry(ce)
!--------------------------------------------------------------------------------------------------
! compute the deformation gradient of individual grains due to relaxations
  F = 0.0_pReal
  do iGrain = 1,product(prm%N_constituents)
    iGrain3 = grain1to3(iGrain,prm%N_constituents)
    do iFace = 1,6
      intFace = getInterface(iFace,iGrain3)                                                         ! identifying 6 interfaces of each grain
      aVect = relaxationVector(intFace,ho,en)                                                       ! get the relaxation vectors for each interface from global relaxation vector array
      nVect = interfaceNormal(intFace,ho,en)
      forall (i=1:3,j=1:3) &
        F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j)                                           ! calculating deformation relaxations due to interface relaxation
    end do
    F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF                                                    ! resulting relaxed deformation gradient
  end do

  end associate

end subroutine RGC_partitionDeformation


!--------------------------------------------------------------------------------------------------
!> @brief update the internal state of the homogenization scheme and tell whether "done" and
! "happy" with result
!--------------------------------------------------------------------------------------------------
module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
      logical, dimension(2) :: doneAndHappy
      real(pReal), dimension(:,:,:),     intent(in)    :: &
        P,&                                                                                         !< partitioned stresses
        F                                                                                           !< partitioned deformation gradients
      real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF                                         !< partitioned stiffnesses
      real(pReal), dimension(3,3),       intent(in) :: avgF                                         !< average F
      real(pReal),                       intent(in) :: dt                                           !< time increment
      integer,                           intent(in) :: &
        ce                                                                                          !< cell

  integer, dimension(4) :: intFaceN,intFaceP,faceID
  integer, dimension(3) :: nGDim,iGr3N,iGr3P
  integer :: ho,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ipert,nGrain, en
  real(pReal), dimension(3,3,size(P,3)) :: R,pF,pR,D,pD
  real(pReal), dimension(3,size(P,3))   :: NN,devNull
  real(pReal), dimension(3)             :: normP,normN,mornP,mornN
  real(pReal) :: residMax,stresMax
  logical :: error
  real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix,rmatrix
  real(pReal), dimension(:), allocatable   :: resid,relax,p_relax,p_resid,drelax

  zeroTimeStep: if (dEq0(dt)) then
    doneAndHappy = .true.                                                                   ! pretend everything is fine and return
    return
  end if zeroTimeStep

  ho  = material_homogenizationID(ce)
  en = material_homogenizationEntry(ce)

  associate(stt => state(ho), st0 => state0(ho), dst => dependentState(ho), prm => param(ho))

!--------------------------------------------------------------------------------------------------
! get the dimension of the cluster (grains and interfaces)
  nGDim  = prm%N_constituents
  nGrain = product(nGDim)
  nIntFaceTot = (nGDim(1)-1)*nGDim(2)*nGDim(3) &
              + nGDim(1)*(nGDim(2)-1)*nGDim(3) &
              + nGDim(1)*nGDim(2)*(nGDim(3)-1)

!--------------------------------------------------------------------------------------------------
! allocate the size of the global relaxation arrays/jacobian matrices depending on the size of the cluster
  allocate(resid(3*nIntFaceTot), source=0.0_pReal)
  allocate(tract(nIntFaceTot,3), source=0.0_pReal)
  relax  = stt%relaxationVector(:,en)
  drelax = stt%relaxationVector(:,en) - st0%relaxationVector(:,en)

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

!--------------------------------------------------------------------------------------------------
! calculating volume discrepancy and stress penalty related to overall volume discrepancy
  call volumePenalty(D,dst%volumeDiscrepancy(en),avgF,F,nGrain)

!------------------------------------------------------------------------------------------------
! computing the residual stress from the balance of traction at all (interior) interfaces
  do iNum = 1,nIntFaceTot
    faceID = interface1to4(iNum,param(ho)%N_constituents)                                           ! 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(ho)%N_constituents)                                                ! translate the local grain ID into global coordinate system (1-dimensional index)
    intFaceN = getInterface(2*faceID(1),iGr3N)
    normN = interfaceNormal(intFaceN,ho,en)

!--------------------------------------------------------------------------------------------------
! identify the right/up/front grain (+|P)
    iGr3P = iGr3N
    iGr3P(faceID(1)) = iGr3N(faceID(1))+1                                                           ! identifying the grain ID in local coordinate system (3-dimensional index)
    iGrP = grain3to1(iGr3P,param(ho)%N_constituents)                                                ! translate the local grain ID into global coordinate system (1-dimensional index)
    intFaceP = getInterface(2*faceID(1)-1,iGr3P)
    normP = interfaceNormal(intFaceP,ho,en)

!--------------------------------------------------------------------------------------------------
! compute the residual of traction at the interface (in local system, 4-dimensional index)
    do i = 1,3
      tract(iNum,i) = sign(num%viscModus*(abs(drelax(i+3*(iNum-1)))/(num%refRelaxRate*dt))**num%viscPower, &
                           drelax(i+3*(iNum-1)))                                                    ! contribution from the relaxation viscosity
      do j = 1,3
        tract(iNum,i) = tract(iNum,i) + (P(i,j,iGrP) + R(i,j,iGrP) + D(i,j,iGrP))*normP(j) &        ! contribution from material stress P, mismatch penalty R, and volume penalty D projected into the interface
                                      + (P(i,j,iGrN) + R(i,j,iGrN) + D(i,j,iGrN))*normN(j)
        resid(i+3*(iNum-1)) = tract(iNum,i)                                                         ! translate the local residual into global 1-dimensional residual array
      end do
    end do

  end do

!--------------------------------------------------------------------------------------------------
! convergence check for stress residual
  stresMax = maxval(abs(P))                                                                         ! get the maximum of first Piola-Kirchhoff (material) stress
  residMax = maxval(abs(tract))                                                                     ! get the maximum of the residual

  doneAndHappy = .false.

!--------------------------------------------------------------------------------------------------
!  If convergence reached => done and happy
  if (residMax < num%rtol*stresMax .or. residMax < num%atol) then
    doneAndHappy = .true.

    dst%mismatch(1:3,en)       = sum(NN,2)/real(nGrain,pReal)
    dst%relaxationRate_avg(en) = sum(abs(drelax))/dt/real(3*nIntFaceTot,pReal)
    dst%relaxationRate_max(en) = maxval(abs(drelax))/dt

    return

!--------------------------------------------------------------------------------------------------
! if residual blows-up => done but unhappy
  elseif (residMax > num%relMax*stresMax .or. residMax > num%absMax) then                           ! try to restart when residual blows up exceeding maximum bound
    doneAndHappy = [.true.,.false.]                                                         ! with direct cut-back
   return
  end if

!---------------------------------------------------------------------------------------------------
! 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,nIntFaceTot
    faceID = interface1to4(iNum,param(ho)%N_constituents)                                           ! 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(ho)%N_constituents)                                                ! translate into global grain ID
    intFaceN = getInterface(2*faceID(1),iGr3N)                                                      ! identifying the connecting interface in local coordinate system
    normN = interfaceNormal(intFaceN,ho,en)
    do iFace = 1,6
      intFaceN = getInterface(iFace,iGr3N)                                                          ! identifying all interfaces that influence relaxation of the above interface
      mornN = interfaceNormal(intFaceN,ho,en)
      iMun = interface4to1(intFaceN,param(ho)%N_constituents)                                       ! translate the interfaces ID into local 4-dimensional index
      if (iMun > 0) then                                                                            ! get the corresponding tangent
        do i=1,3; do j=1,3; do k=1,3; do l=1,3
          smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) &
                                             + dPdF(i,k,j,l,iGrN)*normN(k)*mornN(l)
        end do;enddo;enddo;enddo
! projecting the material tangent dPdF into the interface
! to obtain the Jacobian matrix contribution of dPdF
      end if
    end do

!--------------------------------------------------------------------------------------------------
! identify the right/up/front grain (+|P)
    iGr3P = iGr3N
    iGr3P(faceID(1)) = iGr3N(faceID(1))+1                                                           ! identifying the grain ID in local coordinate sytem
    iGrP = grain3to1(iGr3P,param(ho)%N_constituents)                                                ! translate into global grain ID
    intFaceP = getInterface(2*faceID(1)-1,iGr3P)                                                    ! identifying the connecting interface in local coordinate system
    normP = interfaceNormal(intFaceP,ho,en)
    do iFace = 1,6
      intFaceP = getInterface(iFace,iGr3P)                                                          ! identifying all interfaces that influence relaxation of the above interface
      mornP = interfaceNormal(intFaceP,ho,en)
      iMun = interface4to1(intFaceP,param(ho)%N_constituents)                                       ! translate the interfaces ID into local 4-dimensional index
      if (iMun > 0) then                                                                            ! get the corresponding tangent
        do i=1,3; do j=1,3; do k=1,3; do l=1,3
          smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) &
                                             + dPdF(i,k,j,l,iGrP)*normP(k)*mornP(l)
        end do;enddo;enddo;enddo
      end if
    end do
  end do

!--------------------------------------------------------------------------------------------------
! ... 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,3*nIntFaceTot
    p_relax = relax
    p_relax(ipert) = relax(ipert) + num%pPert                                                       ! perturb the relaxation vector
    stt%relaxationVector(:,en) = p_relax
    call grainDeformation(pF,avgF,ho,en)                                                            ! rain deformation from perturbed state
    call stressPenalty(pR,DevNull,      avgF,pF,ho,en)                                              ! stress penalty due to interface mismatch from perturbed state
    call volumePenalty(pD,devNull(1,1), avgF,pF,nGrain)                                             ! 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,nIntFaceTot
      faceID = interface1to4(iNum,param(ho)%N_constituents)                                         ! 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(ho)%N_constituents)                                              ! translate the local grain ID into global coordinate system (1-dimensional index)
      intFaceN = getInterface(2*faceID(1),iGr3N)                                                    ! identify the interface ID of the grain
      normN = interfaceNormal(intFaceN,ho,en)

!--------------------------------------------------------------------------------------------------
! identify the right/up/front grain (+|P)
      iGr3P = iGr3N
      iGr3P(faceID(1)) = iGr3N(faceID(1))+1                                                         ! identify the grain ID in local coordinate system (3-dimensional index)
      iGrP = grain3to1(iGr3P,param(ho)%N_constituents)                                              ! translate the local grain ID into global coordinate system (1-dimensional index)
      intFaceP = getInterface(2*faceID(1)-1,iGr3P)                                                  ! identify the interface ID of the grain
      normP = interfaceNormal(intFaceP,ho,en)

!--------------------------------------------------------------------------------------------------
! compute the residual stress (contribution of mismatch and volume penalties) from perturbed state
! at all interfaces
      do i = 1,3; do j = 1,3
        p_resid(i+3*(iNum-1)) = p_resid(i+3*(iNum-1)) + (pR(i,j,iGrP) - R(i,j,iGrP))*normP(j) &
                                                      + (pR(i,j,iGrN) - R(i,j,iGrN))*normN(j) &
                                                      + (pD(i,j,iGrP) - D(i,j,iGrP))*normP(j) &
                                                      + (pD(i,j,iGrN) - D(i,j,iGrN))*normN(j)
      end do; end do
    end do
    pmatrix(:,ipert) = p_resid/num%pPert
  end do


!--------------------------------------------------------------------------------------------------
! ... of the numerical viscosity traction "rmatrix"
  allocate(rmatrix(3*nIntFaceTot,3*nIntFaceTot),source=0.0_pReal)
  do i=1,3*nIntFaceTot
    rmatrix(i,i) = num%viscModus*num%viscPower/(num%refRelaxRate*dt)* &                             ! tangent due to numerical viscosity traction appears
                   (abs(drelax(i))/(num%refRelaxRate*dt))**(num%viscPower - 1.0_pReal)              ! only in the main diagonal term
  end do


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

!--------------------------------------------------------------------------------------------------
! computing the update of the state variable (relaxation vectors) using the Jacobian matrix
  allocate(jnverse(3*nIntFaceTot,3*nIntFaceTot),source=0.0_pReal)
  call math_invert(jnverse,error,jmatrix)

!--------------------------------------------------------------------------------------------------
! calculate the state update (global relaxation vectors) for the next Newton-Raphson iteration
  drelax = 0.0_pReal
  do i = 1,3*nIntFaceTot;do j = 1,3*nIntFaceTot
    drelax(i) = drelax(i) - jnverse(i,j)*resid(j)                                                   ! Calculate the correction for the state variable
  end do; end do
  stt%relaxationVector(:,en) = relax + drelax                                                       ! Updateing the state variable for the next iteration
  if (any(abs(drelax) > num%maxdRelax)) then                                                        ! Forcing cutback when the incremental change of relaxation vector becomes too large
    doneAndHappy = [.true.,.false.]
    !$OMP CRITICAL (write2out)
    print'(a,i3,a,i3,a)',' RGC_updateState: enforces cutback'
    print'(a,e15.8)',' due to large relaxation change = ',maxval(abs(drelax))
    flush(IO_STDOUT)
    !$OMP END CRITICAL (write2out)
  end if

  end associate

  contains
  !------------------------------------------------------------------------------------------------
  !> @brief calculate stress-like penalty due to deformation mismatch
  !------------------------------------------------------------------------------------------------
  subroutine stressPenalty(rPen,nMis,avgF,fDef,ho,en)

    real(pReal),   dimension (:,:,:), intent(out) :: rPen                                           !< stress-like penalty
    real(pReal),   dimension (:,:),   intent(out) :: nMis                                           !< total amount of mismatch

    real(pReal),   dimension (:,:,:), intent(in)  :: fDef                                           !< deformation gradients
    real(pReal),   dimension (3,3),   intent(in)  :: avgF                                           !< initial effective stretch tensor
    integer,                          intent(in)  :: ho, en

    integer, dimension (4)   :: intFace
    integer, dimension (3)   :: iGrain3,iGNghb3,nGDim
    real(pReal),   dimension (3,3) :: gDef,nDef
    real(pReal),   dimension (3)   :: nVect,surfCorr
    integer :: iGrain,iGNghb,iFace,i,j,k,l
    real(pReal) :: muGrain,muGNghb,nDefNorm
    real(pReal), parameter  :: &
      nDefToler = 1.0e-10_pReal, &
      b = 2.5e-10_pReal                                                                             ! Length of Burgers vector

    nGDim = param(ho)%N_constituents
    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,ho,en)

    associate(prm => param(ho))

   !-----------------------------------------------------------------------------------------------
   ! computing the mismatch and penalty stress tensor of all grains
   grainLoop: do iGrain = 1,product(prm%N_constituents)
     muGrain = equivalentMu(iGrain,ce)
     iGrain3 = grain1to3(iGrain,prm%N_constituents)                                                 ! get the grain ID in local 3-dimensional index (x,y,z)-position

     interfaceLoop: do iFace = 1,6
       intFace = getInterface(iFace,iGrain3)                                                        ! get the 4-dimensional index of the interface in local numbering system of the grain
       nVect = interfaceNormal(intFace,ho,en)
       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))
       where(iGNghb3 < 1)    iGNghb3 = nGDim
       where(iGNghb3 >nGDim) iGNghb3 = 1
       iGNghb  = grain3to1(iGNghb3,prm%N_constituents)                                              ! get the ID of the neighboring grain
       muGNghb = equivalentMu(iGNghb,ce)
       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,3; do j = 1,3
         do k = 1,3; do l = 1,3
           nDef(i,j) = nDef(i,j) - nVect(k)*gDef(i,l)*math_LeviCivita(j,k,l)                        ! compute the interface mismatch tensor from the jump of deformation gradient
         end do; end do
         nDefNorm = nDefNorm + nDef(i,j)**2.0_pReal                                                 ! compute the norm of the mismatch tensor
       end do; end do
       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)


       !-------------------------------------------------------------------------------------------
       ! compute the stress penalty of all interfaces
       do i = 1,3; do j = 1,3; do k = 1,3; do l = 1,3
         rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*b + muGNghb*b)*prm%xi_alpha &
                                                *surfCorr(abs(intFace(1)))/prm%D_alpha(abs(intFace(1))) &
                                                *cosh(prm%c_alpha*nDefNorm) &
                                                *0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_LeviCivita(k,l,j) &
                                                *tanh(nDefNorm/num%xSmoo)
       end do; end do;enddo; end do
     end do interfaceLoop


   end do grainLoop

   end associate

  end subroutine stressPenalty


  !------------------------------------------------------------------------------------------------
  !> @brief calculate stress-like penalty due to volume discrepancy
  !------------------------------------------------------------------------------------------------
  subroutine volumePenalty(vPen,vDiscrep,fAvg,fDef,nGrain)

    real(pReal), dimension (:,:,:), intent(out) :: vPen                                             ! stress-like penalty due to volume
    real(pReal),                    intent(out) :: vDiscrep                                         ! total volume discrepancy

    real(pReal), dimension (:,:,:), intent(in)  :: fDef                                             ! deformation gradients
    real(pReal), dimension (3,3),   intent(in)  :: fAvg                                             ! overall deformation gradient
    integer,                        intent(in) :: &
      Ngrain

    real(pReal), dimension(size(vPen,3)) :: gVol
    integer :: i

    !----------------------------------------------------------------------------------------------
    ! compute the volumes of grains and of cluster
    vDiscrep = math_det33(fAvg)                                                                     ! compute the volume of the cluster
    do i = 1,nGrain
      gVol(i) = math_det33(fDef(1:3,1:3,i))                                                         ! compute the volume of individual grains
      vDiscrep     = vDiscrep - gVol(i)/real(nGrain,pReal)                                          ! calculate the difference/dicrepancy between
                                                                                                    ! the volume of the cluster and the the total volume of grains
    end do

    !----------------------------------------------------------------------------------------------
    ! calculate the stress and penalty due to volume discrepancy
    vPen      = 0.0_pReal
    do i = 1,nGrain
      vPen(:,:,i) = -1.0_pReal/real(nGrain,pReal)*num%volDiscrMod*num%volDiscrPow/num%maxVolDiscr* &
                         sign((abs(vDiscrep)/num%maxVolDiscr)**(num%volDiscrPow - 1.0),vDiscrep)* &
                         gVol(i)*transpose(math_inv33(fDef(:,:,i)))
    end do

  end subroutine volumePenalty


  !--------------------------------------------------------------------------------------------------
  !> @brief compute the correction factor accouted for surface evolution (area change) due to
  ! deformation
  !--------------------------------------------------------------------------------------------------
  function surfaceCorrection(avgF,ho,en)

    real(pReal), dimension(3)               :: surfaceCorrection

    real(pReal), dimension(3,3), intent(in) :: avgF                                                 !< average F
    integer,                     intent(in) :: &
      ho, &
      en
    real(pReal), dimension(3,3)             :: invC
    real(pReal), dimension(3)               :: nVect
    real(pReal)  :: detF
    integer :: i,j,iBase
    logical     ::  error

    call math_invert33(invC,detF,error,matmul(transpose(avgF),avgF))

    surfaceCorrection = 0.0_pReal
    do iBase = 1,3
      nVect = interfaceNormal([iBase,1,1,1],ho,en)
      do i = 1,3; do j = 1,3
        surfaceCorrection(iBase) = surfaceCorrection(iBase) + invC(i,j)*nVect(i)*nVect(j)           ! compute the component of (the inverse of) the stretch in the direction of the normal
      end do; end do
      surfaceCorrection(iBase) = sqrt(surfaceCorrection(iBase))*detF                                ! get the surface correction factor (area contraction/enlargement)
    end do

  end function surfaceCorrection


  !-------------------------------------------------------------------------------------------------
  !> @brief compute the equivalent shear and bulk moduli from the elasticity tensor
  !-------------------------------------------------------------------------------------------------
  real(pReal) function equivalentMu(grainID,ce)

    integer, intent(in)    :: &
      grainID,&
      ce

    real(pReal), dimension(6,6) :: C


    C = phase_homogenizedC(material_phaseID(grainID,ce),material_phaseEntry(grainID,ce))
    equivalentMu = lattice_equivalent_mu(C,'voigt')

  end function equivalentMu


  !-------------------------------------------------------------------------------------------------
  !> @brief calculating the grain deformation gradient (the same with
  ! homogenization_RGC_partitionDeformation, but used only for perturbation scheme)
  !-------------------------------------------------------------------------------------------------
  subroutine grainDeformation(F, avgF, ho, en)

    real(pReal),   dimension(:,:,:), intent(out) :: F                                               !< partitioned F  per grain

    real(pReal),   dimension(:,:),   intent(in)  :: avgF                                            !< averaged F
    integer,                          intent(in)  :: &
      ho, &
      en

    real(pReal),   dimension(3) :: aVect,nVect
    integer,       dimension(4) :: intFace
    integer,       dimension(3) :: iGrain3
    integer :: iGrain,iFace,i,j

    !-----------------------------------------------------------------------------------------------
    ! compute the deformation gradient of individual grains due to relaxations

    associate (prm => param(ho))

    F = 0.0_pReal
    do iGrain = 1,product(prm%N_constituents)
      iGrain3 = grain1to3(iGrain,prm%N_constituents)
      do iFace = 1,6
        intFace = getInterface(iFace,iGrain3)
        aVect   = relaxationVector(intFace,ho,en)
        nVect   = interfaceNormal(intFace,ho,en)
        forall (i=1:3,j=1:3) &
          F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j)                                         ! effective relaxations
      end do
      F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF                                                  ! relaxed deformation gradient
    end do

    end associate

  end subroutine grainDeformation

end function RGC_updateState


!--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file
!--------------------------------------------------------------------------------------------------
module subroutine RGC_results(ho,group)

  integer,          intent(in) :: ho
  character(len=*), intent(in) :: group

  integer :: o

  associate(stt => state(ho), dst => dependentState(ho), prm => param(ho))
  outputsLoop: do o = 1,size(prm%output)
    select case(trim(prm%output(o)))
      case('M')
        call results_writeDataset(dst%mismatch,group,trim(prm%output(o)), &
                                  'average mismatch tensor','1')
      case('Delta_V')
        call results_writeDataset(dst%volumeDiscrepancy,group,trim(prm%output(o)), &
                                  'volume discrepancy','m³')
      case('max_dot_a')
        call results_writeDataset(dst%relaxationrate_max,group,trim(prm%output(o)), &
                                  'maximum relaxation rate','m/s')
      case('avg_dot_a')
        call results_writeDataset(dst%relaxationrate_avg,group,trim(prm%output(o)), &
                                  'average relaxation rate','m/s')
    end select
  end do outputsLoop
  end associate

end subroutine RGC_results


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

  real(pReal), dimension (3)            :: relaxationVector

  integer,                   intent(in) :: ho,en
  integer,     dimension(4), intent(in) :: intFace                                                  !< set of interface ID in 4D array (normal and position)

  integer :: iNum

!--------------------------------------------------------------------------------------------------
! collect the interface relaxation vector from the global state array

  associate (prm => param(ho), &
             stt => state(ho))

  iNum = interface4to1(intFace,prm%N_constituents)                                                  ! identify the position of the interface in global state array
  if (iNum > 0) then
    relaxationVector = stt%relaxationVector((3*iNum-2):(3*iNum),en)
  else
    relaxationVector = 0.0_pReal
  end if

  end associate

end function relaxationVector


!--------------------------------------------------------------------------------------------------
!> @brief identify the normal of an interface
!--------------------------------------------------------------------------------------------------
pure function interfaceNormal(intFace,ho,en) result(n)

  real(pReal), dimension(3)             :: n
  integer,     dimension(4), intent(in) :: intFace                                                  !< interface ID in 4D array (normal and position)
  integer,                   intent(in) :: &
    ho, &
    en


  associate (dst => dependentState(ho))

    n = 0.0_pReal
    n(abs(intFace(1))) = real(intFace(1)/abs(intFace(1)),pReal)                                     ! get the normal vector w.r.t. cluster axis

    n = matmul(dst%orientation(1:3,1:3,en),n)                                                       ! map the normal vector into sample coordinate system (basis)

  end associate

end function interfaceNormal


!--------------------------------------------------------------------------------------------------
!> @brief collect six faces of a grain in 4D (normal and position)
!--------------------------------------------------------------------------------------------------
pure function getInterface(iFace,iGrain3) result(i)

  integer, dimension(4)             :: i
  integer, dimension(3), intent(in) :: iGrain3                                                      !< grain ID in 3D array
  integer,               intent(in) :: iFace                                                        !< face index (1..6) mapped like (-e1,-e2,-e3,+e1,+e2,+e3) or iDir = (-1,-2,-3,1,2,3)

  integer :: iDir                                                                                   !< direction of interface normal


  iDir = (int(real(iFace-1,pReal)/2.0_pReal)+1)*(-1)**iFace
  i = [iDir,iGrain3]
  if (iDir < 0) i(1-iDir) = i(1-iDir)-1                                       ! 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)

  integer, dimension(3)             :: grain1to3

  integer,               intent(in) :: grain1                                                       !< grain ID in 1D array
  integer, dimension(3), intent(in) :: nGDim

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

end function grain1to3


!--------------------------------------------------------------------------------------------------
!> @brief map grain ID from in 3D (local position) to in 1D (global array)
!--------------------------------------------------------------------------------------------------
integer pure function grain3to1(grain3,nGDim)

  integer, dimension(3), intent(in) :: grain3                                                       !< grain ID in 3D array (pos.x,pos.y,pos.z)
  integer, dimension(3), intent(in) :: nGDim

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

end function grain3to1


!--------------------------------------------------------------------------------------------------
!> @brief maps interface ID from 4D (normal and local position) into 1D (global array)
!--------------------------------------------------------------------------------------------------
integer pure function interface4to1(iFace4D, nGDim)

  integer, dimension(4), intent(in) :: iFace4D                                                      !< interface ID in 4D array (n.dir,pos.x,pos.y,pos.z)
  integer, dimension(3), intent(in) :: nGDim


  select case(abs(iFace4D(1)))

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

    case(2)
      if ((iFace4D(3) == 0) .or. (iFace4D(3) == nGDim(2))) then
        interface4to1 = 0
      else
        interface4to1 = iFace4D(4) + nGDim(3)*(iFace4D(2)-1) &
                      + nGDim(3)*nGDim(1)*(iFace4D(3)-1) &
                      + (nGDim(1)-1)*nGDim(2)*nGDim(3)                                              ! total # of interfaces normal || e1
      end if

    case(3)
      if ((iFace4D(4) == 0) .or. (iFace4D(4) == nGDim(3))) then
        interface4to1 = 0
      else
        interface4to1 = iFace4D(2) + nGDim(1)*(iFace4D(3)-1) &
                      + nGDim(1)*nGDim(2)*(iFace4D(4)-1) &
                      + (nGDim(1)-1)*nGDim(2)*nGDim(3) &                                            ! total # of interfaces normal || e1
                      + nGDim(1)*(nGDim(2)-1)*nGDim(3)                                              ! total # of interfaces normal || e2
      end if

    case default
      interface4to1 = -1

  end select

end function interface4to1


!--------------------------------------------------------------------------------------------------
!> @brief maps interface ID from 1D (global array) into 4D (normal and local position)
!--------------------------------------------------------------------------------------------------
pure function interface1to4(iFace1D, nGDim)

  integer, dimension(4)             :: interface1to4

  integer,               intent(in) :: iFace1D                                                      !< interface ID in 1D array
  integer, dimension(3), intent(in) :: nGDim
  integer, dimension(3)             :: nIntFace

!--------------------------------------------------------------------------------------------------
! compute the total number of interfaces, which ...
  nIntFace = [(nGDim(1)-1)*nGDim(2)*nGDim(3), &                                                     ! ... normal || e1
              nGDim(1)*(nGDim(2)-1)*nGDim(3), &                                                     ! ... normal || e2
              nGDim(1)*nGDim(2)*(nGDim(3)-1)]                                                       ! ... normal || e3

!--------------------------------------------------------------------------------------------------
! get the corresponding interface ID in 4D (normal and local position)
  if (iFace1D > 0 .and. iFace1D <= nIntFace(1)) then                                                ! interface with normal || e1
    interface1to4(1) = 1
    interface1to4(3) = mod((iFace1D-1),nGDim(2))+1
    interface1to4(4) = mod(int(real(iFace1D-1,pReal)/real(nGDim(2),pReal)),nGDim(3))+1
    interface1to4(2) = int(real(iFace1D-1,pReal)/real(nGDim(2),pReal)/real(nGDim(3),pReal))+1
  elseif (iFace1D > nIntFace(1) .and. iFace1D <= (nIntFace(2) + nIntFace(1))) then                  ! interface with normal || e2
    interface1to4(1) = 2
    interface1to4(4) = mod((iFace1D-nIntFace(1)-1),nGDim(3))+1
    interface1to4(2) = mod(int(real(iFace1D-nIntFace(1)-1,pReal)/real(nGDim(3),pReal)),nGDim(1))+1
    interface1to4(3) = int(real(iFace1D-nIntFace(1)-1,pReal)/real(nGDim(3),pReal)/real(nGDim(1),pReal))+1
  elseif (iFace1D > nIntFace(2) + nIntFace(1) .and. iFace1D <= (nIntFace(3) + nIntFace(2) + nIntFace(1))) then ! interface with normal || e3
    interface1to4(1) = 3
    interface1to4(2) = mod((iFace1D-nIntFace(2)-nIntFace(1)-1),nGDim(1))+1
    interface1to4(3) = mod(int(real(iFace1D-nIntFace(2)-nIntFace(1)-1,pReal)/real(nGDim(1),pReal)),nGDim(2))+1
    interface1to4(4) = int(real(iFace1D-nIntFace(2)-nIntFace(1)-1,pReal)/real(nGDim(1),pReal)/real(nGDim(2),pReal))+1
  end if

end function interface1to4


end submodule 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
-												homogenization modules made consistent

											
										
										
											2020-09-23 05:03:19 +05:30
+								!> N_constituents 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
+								!--------------------------------------------------------------------------------------------------
-												missing renames

											
										
										
											2021-02-13 23:27:41 +05:30
+								submodule(homogenization:mechanical) RGC
-												Merge branch 'MiscImprovements' of magit1.mpie.de:/damask/DAMASK into MiscImprovements

											
										
										
											2020-01-29 17:46:49 +05:30
+								  use rotations
-												using central functionality

											
										
										
											2020-12-28 15:33:29 +05:30
+								  use lattice
-												parameters easier to identify

											
										
										
											2018-08-28 16:27:22 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  type :: tParameters
 								    integer, dimension(:), allocatable :: &
-												homogenization modules made consistent

											
										
										
											2020-09-23 05:03:19 +05:30
+								      N_constituents
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    real(pReal) :: &
-												homogenization modules made consistent

											
										
										
											2020-09-23 05:03:19 +05:30
+								      xi_alpha, &
 								      c_Alpha
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    real(pReal), dimension(:), allocatable :: &
-												homogenization modules made consistent

											
										
										
											2020-09-23 05:03:19 +05:30
+								      D_alpha, &
 								      a_g
-												less writing, same functionality

											
										
										
											2020-02-15 03:51:58 +05:30
+								    character(len=pStringLen), allocatable, dimension(:) :: &
 								      output
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  end type tParameters
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  type :: tRGCstate
 								    real(pReal), pointer,     dimension(:,:) :: &
 								      relaxationVector
 								  end type tRGCstate
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  type :: tRGCdependentState
 								    real(pReal), allocatable,     dimension(:) :: &
 								      volumeDiscrepancy, &
 								      relaxationRate_avg, &
 								      relaxationRate_max
 								    real(pReal), allocatable,     dimension(:,:) :: &
 								      mismatch
 								    real(pReal), allocatable,     dimension(:,:,:) :: &
 								      orientation
 								  end type tRGCdependentState
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												better store data locally

											
										
										
											2020-03-17 05:09:32 +05:30
+								  type :: tNumerics_RGC
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
+								    real(pReal) :: &
-												minor polishing

											
										
										
											2020-06-26 15:52:33 +05:30
+								      atol, &                                                                                       !< absolute tolerance of RGC residuum
 								      rtol, &                                                                                       !< relative tolerance of RGC residuum
 								      absMax, &                                                                                     !< absolute maximum of RGC residuum
 								      relMax, &                                                                                     !< relative maximum of RGC residuum
 								      pPert, &                                                                                      !< perturbation for computing RGC penalty tangent
 								      xSmoo, &                                                                                      !< RGC penalty smoothing parameter (hyperbolic tangent)
 								      viscPower, &                                                                                  !< power (sensitivity rate) of numerical viscosity in RGC scheme, Default 1.0e0: Newton viscosity (linear model)
 								      viscModus, &                                                                                  !< stress modulus of RGC numerical viscosity, Default 0.0e0: No viscosity is applied
 								      refRelaxRate, &                                                                               !< reference relaxation rate in RGC viscosity
 								      maxdRelax, &                                                                                  !< threshold of maximum relaxation vector increment (if exceed this then cutback)
 								      maxVolDiscr, &                                                                                !< threshold of maximum volume discrepancy allowed
 								      volDiscrMod, &                                                                                !< stiffness of RGC volume discrepancy (zero = without volume discrepancy constraint)
 								      volDiscrPow                                                                                   !< powerlaw penalty for volume discrepancy
-												better store data locally

											
										
										
											2020-03-17 05:09:32 +05:30
+								  end type tNumerics_RGC
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  type(tparameters),          dimension(:), allocatable :: &
 								    param
 								  type(tRGCstate),            dimension(:), allocatable :: &
 								    state, &
 								    state0
 								  type(tRGCdependentState),   dimension(:), allocatable :: &
 								    dependentState
-												better store data locally

											
										
										
											2020-03-17 05:09:32 +05:30
+								  type(tNumerics_RGC) :: &
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
+								    num                                                                                             ! numerics parameters. Better name?
-												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
+								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
+								!--------------------------------------------------------------------------------------------------
-												shorter names

need to prefix 'pass' to avoid name clashes that result in errors during
compilation

											
										
										
											2021-04-06 15:35:47 +05:30
+								module subroutine RGC_init(num_homogMech)
-												Read numerics.yaml once

											
										
										
											2020-06-16 22:17:19 +05:30
 								  class(tNode), pointer, intent(in) :: &
-												use submodule property; simplifying

											
										
										
											2020-07-02 01:50:22 +05:30
+								    num_homogMech                                                                                   !< pointer to mechanical homogenization numerics data
-												cleaning

											
										
										
											2020-07-03 20:15:11 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  integer :: &
-												cleaning homogenization

											
										
										
											2021-02-26 02:12:40 +05:30
+								    ho, &
-												new names and mappings

											
										
										
											2021-05-23 13:40:25 +05:30
+								    Nmembers, &
-												use submodule property; simplifying

											
										
										
											2020-07-02 01:50:22 +05:30
+								    sizeState, nIntFaceTot
-												ensure defined initital state

											
										
										
											2020-09-07 16:50:00 +05:30
-												Read numerics.yaml once

											
										
										
											2020-06-16 22:17:19 +05:30
+								  class (tNode), pointer :: &
-												Restructuring for material.yaml

											
										
										
											2020-08-15 19:32:10 +05:30
+								    num_RGC, &                                                                                      ! pointer to RGC numerics data
 								    material_homogenization, &
 								    homog, &
 								    homogMech
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								  print'(/,1x,a)', '<<<+-  homogenization:mechanical:RGC init  -+>>>'
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								  print'(/,a,i0)', ' # homogenizations: ',count(homogenization_type == HOMOGENIZATION_RGC_ID)
-												avoid repetition

											
										
										
											2021-04-11 15:48:26 +05:30
+								  flush(IO_STDOUT)
-												simplified print and format strings

											
										
										
											2020-09-18 02:27:56 +05:30
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								  print'(/,1x,a)', 'D.D. Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009'
 								  print'(  1x,a)', 'https://doi.org/10.1007/s12289-009-0619-1'//IO_EOL
-												simplified print and format strings

											
										
										
											2020-09-18 02:27:56 +05:30
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								  print'(/,1x,a)', 'D.D. Tjahjanto et al., Modelling and Simulation in Materials Science and Engineering 18:015006, 2010'
 								  print'(  1x,a)', 'https://doi.org/10.1088/0965-0393/18/1/015006'//IO_EOL
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												storing by instance just complicates the code

											
										
										
											2021-02-23 17:47:51 +05:30
+								  material_homogenization => config_material%get('homogenization')
 								  allocate(param(material_homogenization%length))
 								  allocate(state(material_homogenization%length))
 								  allocate(state0(material_homogenization%length))
 								  allocate(dependentState(material_homogenization%length))
-												ensure defined initital state

											
										
										
											2020-09-07 16:50:00 +05:30
-												Read numerics.yaml once

											
										
										
											2020-06-16 22:17:19 +05:30
+								  num_RGC => num_homogMech%get('RGC',defaultVal=emptyDict)
-												[skip ci] polishing

											
										
										
											2020-06-29 18:39:13 +05:30
-												minor polishing

											
										
										
											2020-06-26 15:52:33 +05:30
+								  num%atol         =  num_RGC%get_asFloat('atol',              defaultVal=1.0e+4_pReal)
 								  num%rtol         =  num_RGC%get_asFloat('rtol',              defaultVal=1.0e-3_pReal)
 								  num%absMax       =  num_RGC%get_asFloat('amax',              defaultVal=1.0e+10_pReal)
 								  num%relMax       =  num_RGC%get_asFloat('rmax',              defaultVal=1.0e+2_pReal)
 								  num%pPert        =  num_RGC%get_asFloat('perturbpenalty',    defaultVal=1.0e-7_pReal)
 								  num%xSmoo        =  num_RGC%get_asFloat('relvantmismatch',   defaultVal=1.0e-5_pReal)
 								  num%viscPower    =  num_RGC%get_asFloat('viscositypower',    defaultVal=1.0e+0_pReal)
 								  num%viscModus    =  num_RGC%get_asFloat('viscositymodulus',  defaultVal=0.0e+0_pReal)
 								  num%refRelaxRate =  num_RGC%get_asFloat('refrelaxationrate', defaultVal=1.0e-3_pReal)
 								  num%maxdRelax    =  num_RGC%get_asFloat('maxrelaxationrate', defaultVal=1.0e+0_pReal)
 								  num%maxVolDiscr  =  num_RGC%get_asFloat('maxvoldiscrepancy', defaultVal=1.0e-5_pReal)
 								  num%volDiscrMod  =  num_RGC%get_asFloat('voldiscrepancymod', defaultVal=1.0e+12_pReal)
 								  num%volDiscrPow  =  num_RGC%get_asFloat('dicrepancypower',   defaultVal=5.0_pReal)
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
 								  if (num%atol <= 0.0_pReal)         call IO_error(301,ext_msg='absTol_RGC')
 								  if (num%rtol <= 0.0_pReal)         call IO_error(301,ext_msg='relTol_RGC')
 								  if (num%absMax <= 0.0_pReal)       call IO_error(301,ext_msg='absMax_RGC')
 								  if (num%relMax <= 0.0_pReal)       call IO_error(301,ext_msg='relMax_RGC')
 								  if (num%pPert <= 0.0_pReal)        call IO_error(301,ext_msg='pPert_RGC')
 								  if (num%xSmoo <= 0.0_pReal)        call IO_error(301,ext_msg='xSmoo_RGC')
 								  if (num%viscPower < 0.0_pReal)     call IO_error(301,ext_msg='viscPower_RGC')
 								  if (num%viscModus < 0.0_pReal)     call IO_error(301,ext_msg='viscModus_RGC')
 								  if (num%refRelaxRate <= 0.0_pReal) call IO_error(301,ext_msg='refRelaxRate_RGC')
 								  if (num%maxdRelax <= 0.0_pReal)    call IO_error(301,ext_msg='maxdRelax_RGC')
 								  if (num%maxVolDiscr <= 0.0_pReal)  call IO_error(301,ext_msg='maxVolDiscr_RGC')
 								  if (num%volDiscrMod < 0.0_pReal)   call IO_error(301,ext_msg='volDiscrMod_RGC')
 								  if (num%volDiscrPow <= 0.0_pReal)  call IO_error(301,ext_msg='volDiscrPw_RGC')
-												Restructuring for material.yaml

											
										
										
											2020-08-15 19:32:10 +05:30
-												cleaning homogenization

											
										
										
											2021-02-26 02:12:40 +05:30
+								  do ho = 1, size(homogenization_type)
 								    if (homogenization_type(ho) /= HOMOGENIZATION_RGC_ID) cycle
 								    homog => material_homogenization%get(ho)
-												consistent naming in HDF5 output

											
										
										
											2021-03-25 23:52:59 +05:30
+								    homogMech => homog%get('mechanical')
-												cleaning homogenization

											
										
										
											2021-02-26 02:12:40 +05:30
+								    associate(prm => param(ho), &
 								              stt => state(ho), &
 								              st0 => state0(ho), &
 								              dst => dependentState(ho))
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												Restructuring for material.yaml

											
										
										
											2020-08-15 19:32:10 +05:30
+								#if defined (__GFORTRAN__)
-												get_asXXs --> get_as1dXX

											
										
										
											2021-03-11 22:30:07 +05:30
+								    prm%output = output_as1dString(homogMech)
-												Restructuring for material.yaml

											
										
										
											2020-08-15 19:32:10 +05:30
+								#else
-												get_asXXs --> get_as1dXX

											
										
										
											2021-03-11 22:30:07 +05:30
+								    prm%output = homogMech%get_as1dString('output',defaultVal=emptyStringArray)
-												Restructuring for material.yaml

											
										
										
											2020-08-15 19:32:10 +05:30
+								#endif
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												get_asXXs --> get_as1dXX

											
										
										
											2021-03-11 22:30:07 +05:30
+								    prm%N_constituents = homogMech%get_as1dInt('cluster_size',requiredSize=3)
-												cleaning homogenization

											
										
										
											2021-02-26 02:12:40 +05:30
+								    if (homogenization_Nconstituents(ho) /= product(prm%N_constituents)) &
-												shorter names

need to prefix 'pass' to avoid name clashes that result in errors during
compilation

											
										
										
											2021-04-06 15:35:47 +05:30
+								      call IO_error(211,ext_msg='N_constituents (RGC)')
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												homogenization modules made consistent

											
										
										
											2020-09-23 05:03:19 +05:30
+								    prm%xi_alpha = homogMech%get_asFloat('xi_alpha')
 								    prm%c_alpha  = homogMech%get_asFloat('c_alpha')
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												get_asXXs --> get_as1dXX

											
										
										
											2021-03-11 22:30:07 +05:30
+								    prm%D_alpha  = homogMech%get_as1dFloat('D_alpha', requiredSize=3)
 								    prm%a_g      = homogMech%get_as1dFloat('a_g',     requiredSize=3)
-												less writing, same functionality

											
										
										
											2020-02-15 03:51:58 +05:30
-												new names and mappings

											
										
										
											2021-05-23 13:40:25 +05:30
+								    Nmembers = count(material_homogenizationID == ho)
-												homogenization modules made consistent

											
										
										
											2020-09-23 05:03:19 +05:30
+								    nIntFaceTot = 3*(  (prm%N_constituents(1)-1)*prm%N_constituents(2)*prm%N_constituents(3) &
 								                      + prm%N_constituents(1)*(prm%N_constituents(2)-1)*prm%N_constituents(3) &
 								                      + prm%N_constituents(1)*prm%N_constituents(2)*(prm%N_constituents(3)-1))
-												simplified

- no extra state
- no extra argument

at the cost of less output

											
										
										
											2020-12-29 19:24:58 +05:30
+								    sizeState = nIntFaceTot
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												cleaning homogenization

											
										
										
											2021-02-26 02:12:40 +05:30
+								    homogState(ho)%sizeState = sizeState
-												new names and mappings

											
										
										
											2021-05-23 13:40:25 +05:30
+								    allocate(homogState(ho)%state0   (sizeState,Nmembers), source=0.0_pReal)
 								    allocate(homogState(ho)%state    (sizeState,Nmembers), source=0.0_pReal)
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												cleaning homogenization

											
										
										
											2021-02-26 02:12:40 +05:30
+								    stt%relaxationVector   => homogState(ho)%state(1:nIntFaceTot,:)
 								    st0%relaxationVector   => homogState(ho)%state0(1:nIntFaceTot,:)
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												new names and mappings

											
										
										
											2021-05-23 13:40:25 +05:30
+								    allocate(dst%volumeDiscrepancy(   Nmembers), source=0.0_pReal)
 								    allocate(dst%relaxationRate_avg(  Nmembers), source=0.0_pReal)
 								    allocate(dst%relaxationRate_max(  Nmembers), source=0.0_pReal)
 								    allocate(dst%mismatch(          3,Nmembers), source=0.0_pReal)
-												cleaning

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

											
										
										
											2019-01-13 02:34:03 +05:30
+								! assigning cluster orientations
-												new names and mappings

											
										
										
											2021-05-23 13:40:25 +05:30
+								    dependentState(ho)%orientation = spread(eu2om(prm%a_g*inRad),3,Nmembers)
 								    !dst%orientation = spread(eu2om(prm%a_g*inRad),3,Nmembers) ifort version 18.0.1 crashes (for whatever reason)
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    end associate
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								  end do
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												shorter names

need to prefix 'pass' to avoid name clashes that result in errors during
compilation

											
										
										
											2021-04-06 15:35:47 +05:30
+								end subroutine 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
 								!--------------------------------------------------------------------------------------------------
-												shorter names

need to prefix 'pass' to avoid name clashes that result in errors during
compilation

											
										
										
											2021-04-06 15:35:47 +05:30
+								module subroutine RGC_partitionDeformation(F,avgF,ce)
-												# 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
-												partition, not partion

											
										
										
											2020-10-08 01:45:13 +05:30
+								  real(pReal),   dimension (:,:,:), intent(out) :: F                                                !< partitioned F  per grain
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  real(pReal),   dimension (3,3),   intent(in)  :: avgF                                             !< averaged F
 								  integer,                          intent(in)  :: &
-												use cell mapping

											
										
										
											2021-02-15 23:13:51 +05:30
+								    ce
-												use submodule property; simplifying

											
										
										
											2020-07-02 01:50:22 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  real(pReal), dimension(3) :: aVect,nVect
 								  integer,     dimension(4) :: intFace
 								  integer,     dimension(3) :: iGrain3
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								  integer ::  iGrain,iFace,i,j,ho,en
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								  associate(prm => param(material_homogenizationID(ce)))
-												avoid repetition

											
										
										
											2021-04-11 15:48:26 +05:30
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								  ho = material_homogenizationID(ce)
 								  en = material_homogenizationEntry(ce)
-												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
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  F = 0.0_pReal
-												homogenization modules made consistent

											
										
										
											2020-09-23 05:03:19 +05:30
+								  do iGrain = 1,product(prm%N_constituents)
 								    iGrain3 = grain1to3(iGrain,prm%N_constituents)
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    do iFace = 1,6
 								      intFace = getInterface(iFace,iGrain3)                                                         ! identifying 6 interfaces of each grain
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								      aVect = relaxationVector(intFace,ho,en)                                                       ! get the relaxation vectors for each interface from global relaxation vector array
 								      nVect = interfaceNormal(intFace,ho,en)
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								      forall (i=1:3,j=1:3) &
 								        F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j)                                           ! calculating deformation relaxations due to interface relaxation
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								    end do
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF                                                    ! resulting relaxed deformation gradient
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								  end do
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +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
-												shorter names

need to prefix 'pass' to avoid name clashes that result in errors during
compilation

											
										
										
											2021-04-06 15:35:47 +05:30
+								end subroutine RGC_partitionDeformation
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
 								!--------------------------------------------------------------------------------------------------
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
+								!> @brief update the internal state of the homogenization scheme and tell whether "done" and
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								! "happy" with result
 								!--------------------------------------------------------------------------------------------------
-												shorter names

need to prefix 'pass' to avoid name clashes that result in errors during
compilation

											
										
										
											2021-04-06 15:35:47 +05:30
+								module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
-												intended hierarchy

											
										
										
											2020-12-28 14:25:54 +05:30
+								      logical, dimension(2) :: doneAndHappy
 								      real(pReal), dimension(:,:,:),     intent(in)    :: &
 								        P,&                                                                                         !< partitioned stresses
-												simplified

- no extra state
- no extra argument

at the cost of less output

											
										
										
											2020-12-29 19:24:58 +05:30
+								        F                                                                                           !< partitioned deformation gradients
-												intended hierarchy

											
										
										
											2020-12-28 14:25:54 +05:30
+								      real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF                                         !< partitioned stiffnesses
 								      real(pReal), dimension(3,3),       intent(in) :: avgF                                         !< average F
 								      real(pReal),                       intent(in) :: dt                                           !< time increment
 								      integer,                           intent(in) :: &
-												instance less; use cell mapping

											
										
										
											2021-02-22 20:47:32 +05:30
+								        ce                                                                                          !< cell
-												cleaned and unified

											
										
										
											2018-11-03 21:20:43 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  integer, dimension(4) :: intFaceN,intFaceP,faceID
 								  integer, dimension(3) :: nGDim,iGr3N,iGr3P
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								  integer :: ho,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ipert,nGrain, en
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  real(pReal), dimension(3,3,size(P,3)) :: R,pF,pR,D,pD
 								  real(pReal), dimension(3,size(P,3))   :: NN,devNull
 								  real(pReal), dimension(3)             :: normP,normN,mornP,mornN
 								  real(pReal) :: residMax,stresMax
 								  logical :: error
 								  real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix,rmatrix
 								  real(pReal), dimension(:), allocatable   :: resid,relax,p_relax,p_resid,drelax
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								  zeroTimeStep: if (dEq0(dt)) then
-												intended hierarchy

											
										
										
											2020-12-28 14:25:54 +05:30
+								    doneAndHappy = .true.                                                                   ! pretend everything is fine and return
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
+								    return
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								  end if 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
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								  ho  = material_homogenizationID(ce)
 								  en = material_homogenizationEntry(ce)
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												storing by instance just complicates the code

											
										
										
											2021-02-23 17:47:51 +05:30
+								  associate(stt => state(ho), st0 => state0(ho), dst => dependentState(ho), prm => param(ho))
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								!--------------------------------------------------------------------------------------------------
 								! get the dimension of the cluster (grains and interfaces)
-												homogenization modules made consistent

											
										
										
											2020-09-23 05:03:19 +05:30
+								  nGDim  = prm%N_constituents
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  nGrain = product(nGDim)
 								  nIntFaceTot = (nGDim(1)-1)*nGDim(2)*nGDim(3) &
 								              + nGDim(1)*(nGDim(2)-1)*nGDim(3) &
 								              + nGDim(1)*nGDim(2)*(nGDim(3)-1)
-												# 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
-												simplified

- no extra state
- no extra argument

at the cost of less output

											
										
										
											2020-12-29 19:24:58 +05:30
+								  allocate(resid(3*nIntFaceTot), source=0.0_pReal)
 								  allocate(tract(nIntFaceTot,3), source=0.0_pReal)
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								  relax  = stt%relaxationVector(:,en)
 								  drelax = stt%relaxationVector(:,en) - st0%relaxationVector(:,en)
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09: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
+								!--------------------------------------------------------------------------------------------------
 								! computing interface mismatch and stress penalty tensor for all interfaces of all grains
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								  call stressPenalty(R,NN,avgF,F,ho,en)
-												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
+								!--------------------------------------------------------------------------------------------------
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
+								! calculating volume discrepancy and stress penalty related to overall volume discrepancy
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								  call volumePenalty(D,dst%volumeDiscrepancy(en),avgF,F,nGrain)
-												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
-												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
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  do iNum = 1,nIntFaceTot
-												storing by instance just complicates the code

											
										
										
											2021-02-23 17:47:51 +05:30
+								    faceID = interface1to4(iNum,param(ho)%N_constituents)                                           ! 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)
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    iGr3N = faceID(2:4)                                                                             ! identifying the grain ID in local coordinate system (3-dimensional index)
-												storing by instance just complicates the code

											
										
										
											2021-02-23 17:47:51 +05:30
+								    iGrN = grain3to1(iGr3N,param(ho)%N_constituents)                                                ! translate the local grain ID into global coordinate system (1-dimensional index)
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    intFaceN = getInterface(2*faceID(1),iGr3N)
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								    normN = interfaceNormal(intFaceN,ho,en)
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09: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
+								!--------------------------------------------------------------------------------------------------
 								! identify the right/up/front grain (+|P)
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    iGr3P = iGr3N
 								    iGr3P(faceID(1)) = iGr3N(faceID(1))+1                                                           ! identifying the grain ID in local coordinate system (3-dimensional index)
-												storing by instance just complicates the code

											
										
										
											2021-02-23 17:47:51 +05:30
+								    iGrP = grain3to1(iGr3P,param(ho)%N_constituents)                                                ! translate the local grain ID into global coordinate system (1-dimensional index)
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    intFaceP = getInterface(2*faceID(1)-1,iGr3P)
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								    normP = interfaceNormal(intFaceP,ho,en)
-. 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)
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    do i = 1,3
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
+								      tract(iNum,i) = sign(num%viscModus*(abs(drelax(i+3*(iNum-1)))/(num%refRelaxRate*dt))**num%viscPower, &
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								                           drelax(i+3*(iNum-1)))                                                    ! contribution from the relaxation viscosity
 								      do j = 1,3
 								        tract(iNum,i) = tract(iNum,i) + (P(i,j,iGrP) + R(i,j,iGrP) + D(i,j,iGrP))*normP(j) &        ! contribution from material stress P, mismatch penalty R, and volume penalty D projected into the interface
 								                                      + (P(i,j,iGrN) + R(i,j,iGrN) + D(i,j,iGrN))*normN(j)
 								        resid(i+3*(iNum-1)) = tract(iNum,i)                                                         ! translate the local residual into global 1-dimensional residual array
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								      end do
 								    end do
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								  end do
-												# 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: 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
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  stresMax = maxval(abs(P))                                                                         ! get the maximum of first Piola-Kirchhoff (material) stress
 								  residMax = maxval(abs(tract))                                                                     ! get the maximum of the residual
-. Adding more comments
2. Changing 'subroutine' into 'function' for optimizing the computation time

											
										
										
											2010-11-26 17:20:20 +05:30
-												intended hierarchy

											
										
										
											2020-12-28 14:25:54 +05:30
+								  doneAndHappy = .false.
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09: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
+								!--------------------------------------------------------------------------------------------------
 								!  If convergence reached => done and happy
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
+								  if (residMax < num%rtol*stresMax .or. residMax < num%atol) then
-												intended hierarchy

											
										
										
											2020-12-28 14:25:54 +05:30
+								    doneAndHappy = .true.
-. Adding more comments
2. Changing 'subroutine' into 'function' for optimizing the computation time

											
										
										
											2010-11-26 17:20:20 +05:30
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								    dst%mismatch(1:3,en)       = sum(NN,2)/real(nGrain,pReal)
 								    dst%relaxationRate_avg(en) = sum(abs(drelax))/dt/real(3*nIntFaceTot,pReal)
 								    dst%relaxationRate_max(en) = maxval(abs(drelax))/dt
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +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
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
+								  elseif (residMax > num%relMax*stresMax .or. residMax > num%absMax) then                           ! try to restart when residual blows up exceeding maximum bound
-												intended hierarchy

											
										
										
											2020-12-28 14:25:54 +05:30
+								    doneAndHappy = [.true.,.false.]                                                         ! with direct cut-back
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								   return
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								  end if
-												# 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
+								!---------------------------------------------------------------------------------------------------
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
+								! construct the global Jacobian matrix for updating the global relaxation vector array when
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								! convergence is not yet reached ...
 								!--------------------------------------------------------------------------------------------------
 								! ... of the constitutive stress tangent, assembled from dPdF or material constitutive model "smatrix"
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  allocate(smatrix(3*nIntFaceTot,3*nIntFaceTot), source=0.0_pReal)
 								  do iNum = 1,nIntFaceTot
-												storing by instance just complicates the code

											
										
										
											2021-02-23 17:47:51 +05:30
+								    faceID = interface1to4(iNum,param(ho)%N_constituents)                                           ! assembling of local dPdF into global Jacobian matrix
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09: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
+								!--------------------------------------------------------------------------------------------------
 								! identify the left/bottom/back grain (-|N)
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    iGr3N = faceID(2:4)                                                                             ! identifying the grain ID in local coordinate sytem
-												storing by instance just complicates the code

											
										
										
											2021-02-23 17:47:51 +05:30
+								    iGrN = grain3to1(iGr3N,param(ho)%N_constituents)                                                ! translate into global grain ID
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    intFaceN = getInterface(2*faceID(1),iGr3N)                                                      ! identifying the connecting interface in local coordinate system
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								    normN = interfaceNormal(intFaceN,ho,en)
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    do iFace = 1,6
 								      intFaceN = getInterface(iFace,iGr3N)                                                          ! identifying all interfaces that influence relaxation of the above interface
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								      mornN = interfaceNormal(intFaceN,ho,en)
-												storing by instance just complicates the code

											
										
										
											2021-02-23 17:47:51 +05:30
+								      iMun = interface4to1(intFaceN,param(ho)%N_constituents)                                       ! translate the interfaces ID into local 4-dimensional index
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								      if (iMun > 0) then                                                                            ! get the corresponding tangent
 								        do i=1,3; do j=1,3; do k=1,3; do l=1,3
 								          smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) &
 								                                             + dPdF(i,k,j,l,iGrN)*normN(k)*mornN(l)
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								        end do;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
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								      end if
 								    end do
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09: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
+								!--------------------------------------------------------------------------------------------------
 								! identify the right/up/front grain (+|P)
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    iGr3P = iGr3N
 								    iGr3P(faceID(1)) = iGr3N(faceID(1))+1                                                           ! identifying the grain ID in local coordinate sytem
-												storing by instance just complicates the code

											
										
										
											2021-02-23 17:47:51 +05:30
+								    iGrP = grain3to1(iGr3P,param(ho)%N_constituents)                                                ! translate into global grain ID
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    intFaceP = getInterface(2*faceID(1)-1,iGr3P)                                                    ! identifying the connecting interface in local coordinate system
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								    normP = interfaceNormal(intFaceP,ho,en)
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    do iFace = 1,6
 								      intFaceP = getInterface(iFace,iGr3P)                                                          ! identifying all interfaces that influence relaxation of the above interface
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								      mornP = interfaceNormal(intFaceP,ho,en)
-												storing by instance just complicates the code

											
										
										
											2021-02-23 17:47:51 +05:30
+								      iMun = interface4to1(intFaceP,param(ho)%N_constituents)                                       ! translate the interfaces ID into local 4-dimensional index
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								      if (iMun > 0) then                                                                            ! get the corresponding tangent
 								        do i=1,3; do j=1,3; do k=1,3; do l=1,3
 								          smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) &
 								                                             + dPdF(i,k,j,l,iGrP)*normP(k)*mornP(l)
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								        end do;enddo;enddo;enddo
 								      end if
 								    end do
 								  end do
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09: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
+								!--------------------------------------------------------------------------------------------------
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
+								! ... of the stress penalty tangent (mismatch penalty and volume penalty, computed using numerical
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								! perturbation method) "pmatrix"
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +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)
 								  do ipert = 1,3*nIntFaceTot
 								    p_relax = relax
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
+								    p_relax(ipert) = relax(ipert) + num%pPert                                                       ! perturb the relaxation vector
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								    stt%relaxationVector(:,en) = p_relax
 								    call grainDeformation(pF,avgF,ho,en)                                                            ! rain deformation from perturbed state
 								    call stressPenalty(pR,DevNull,      avgF,pF,ho,en)                                              ! stress penalty due to interface mismatch from perturbed state
-												cleaning homogenization

											
										
										
											2021-02-26 02:12:40 +05:30
+								    call volumePenalty(pD,devNull(1,1), avgF,pF,nGrain)                                             ! 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
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    p_resid = 0.0_pReal
 								    do iNum = 1,nIntFaceTot
-												storing by instance just complicates the code

											
										
										
											2021-02-23 17:47:51 +05:30
+								      faceID = interface1to4(iNum,param(ho)%N_constituents)                                         ! 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)
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								      iGr3N = faceID(2:4)                                                                           ! identify the grain ID in local coordinate system (3-dimensional index)
-												storing by instance just complicates the code

											
										
										
											2021-02-23 17:47:51 +05:30
+								      iGrN = grain3to1(iGr3N,param(ho)%N_constituents)                                              ! translate the local grain ID into global coordinate system (1-dimensional index)
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								      intFaceN = getInterface(2*faceID(1),iGr3N)                                                    ! identify the interface ID of the grain
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								      normN = interfaceNormal(intFaceN,ho,en)
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09: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
+								!--------------------------------------------------------------------------------------------------
 								! identify the right/up/front grain (+|P)
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								      iGr3P = iGr3N
 								      iGr3P(faceID(1)) = iGr3N(faceID(1))+1                                                         ! identify the grain ID in local coordinate system (3-dimensional index)
-												storing by instance just complicates the code

											
										
										
											2021-02-23 17:47:51 +05:30
+								      iGrP = grain3to1(iGr3P,param(ho)%N_constituents)                                              ! translate the local grain ID into global coordinate system (1-dimensional index)
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								      intFaceP = getInterface(2*faceID(1)-1,iGr3P)                                                  ! identify the interface ID of the grain
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								      normP = interfaceNormal(intFaceP,ho,en)
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09: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
+								!--------------------------------------------------------------------------------------------------
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
+								! compute the residual stress (contribution of mismatch and volume penalties) from perturbed state
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								! at all interfaces
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								      do i = 1,3; do j = 1,3
 								        p_resid(i+3*(iNum-1)) = p_resid(i+3*(iNum-1)) + (pR(i,j,iGrP) - R(i,j,iGrP))*normP(j) &
 								                                                      + (pR(i,j,iGrN) - R(i,j,iGrN))*normN(j) &
 								                                                      + (pD(i,j,iGrP) - D(i,j,iGrP))*normP(j) &
 								                                                      + (pD(i,j,iGrN) - D(i,j,iGrN))*normN(j)
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								      end do; end do
 								    end do
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
+								    pmatrix(:,ipert) = p_resid/num%pPert
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								  end do
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09: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
+								!--------------------------------------------------------------------------------------------------
 								! ... of the numerical viscosity traction "rmatrix"
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  allocate(rmatrix(3*nIntFaceTot,3*nIntFaceTot),source=0.0_pReal)
 								  do i=1,3*nIntFaceTot
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
+								    rmatrix(i,i) = num%viscModus*num%viscPower/(num%refRelaxRate*dt)* &                             ! tangent due to numerical viscosity traction appears
 								                   (abs(drelax(i))/(num%refRelaxRate*dt))**(num%viscPower - 1.0_pReal)              ! only in the main diagonal term
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								  end do
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +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
-												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
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  allocate(jmatrix(3*nIntFaceTot,3*nIntFaceTot)); jmatrix = smatrix + pmatrix + rmatrix
-												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
+								!--------------------------------------------------------------------------------------------------
-												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
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  allocate(jnverse(3*nIntFaceTot,3*nIntFaceTot),source=0.0_pReal)
-												using newer interface

											
										
										
											2019-09-21 06:50:33 +05:30
+								  call math_invert(jnverse,error,jmatrix)
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09: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
+								!--------------------------------------------------------------------------------------------------
 								! calculate the state update (global relaxation vectors) for the next Newton-Raphson iteration
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  drelax = 0.0_pReal
 								  do i = 1,3*nIntFaceTot;do j = 1,3*nIntFaceTot
 								    drelax(i) = drelax(i) - jnverse(i,j)*resid(j)                                                   ! Calculate the correction for the state variable
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								  end do; end do
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								  stt%relaxationVector(:,en) = relax + drelax                                                       ! Updateing the state variable for the next iteration
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
+								  if (any(abs(drelax) > num%maxdRelax)) then                                                        ! Forcing cutback when the incremental change of relaxation vector becomes too large
-												intended hierarchy

											
										
										
											2020-12-28 14:25:54 +05:30
+								    doneAndHappy = [.true.,.false.]
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    !$OMP CRITICAL (write2out)
-												cleaning homogenization

											
										
										
											2021-02-26 02:12:40 +05:30
+								    print'(a,i3,a,i3,a)',' RGC_updateState: enforces cutback'
-												simplified print and format strings

											
										
										
											2020-09-18 02:27:56 +05:30
+								    print'(a,e15.8)',' due to large relaxation change = ',maxval(abs(drelax))
-												standard conforming name

											
										
										
											2020-09-22 16:39:12 +05:30
+								    flush(IO_STDOUT)
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    !$OMP END CRITICAL (write2out)
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								  end if
-. Adding more comments
2. Changing 'subroutine' into 'function' for optimizing the computation time

											
										
										
											2010-11-26 17:20:20 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  end associate
-												functions only used within updatestate

											
										
										
											2018-11-04 01:41:43 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  contains
 								  !------------------------------------------------------------------------------------------------
 								  !> @brief calculate stress-like penalty due to deformation mismatch
 								  !------------------------------------------------------------------------------------------------
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								  subroutine stressPenalty(rPen,nMis,avgF,fDef,ho,en)
-												one implicit none is enough

											
										
										
											2019-05-15 02:42:32 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    real(pReal),   dimension (:,:,:), intent(out) :: rPen                                           !< stress-like penalty
 								    real(pReal),   dimension (:,:),   intent(out) :: nMis                                           !< total amount of mismatch
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    real(pReal),   dimension (:,:,:), intent(in)  :: fDef                                           !< deformation gradients
 								    real(pReal),   dimension (3,3),   intent(in)  :: avgF                                           !< initial effective stretch tensor
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								    integer,                          intent(in)  :: ho, en
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    integer, dimension (4)   :: intFace
 								    integer, dimension (3)   :: iGrain3,iGNghb3,nGDim
 								    real(pReal),   dimension (3,3) :: gDef,nDef
 								    real(pReal),   dimension (3)   :: nVect,surfCorr
 								    integer :: iGrain,iGNghb,iFace,i,j,k,l
-												using central functionality

											
										
										
											2020-12-28 15:33:29 +05:30
+								    real(pReal) :: muGrain,muGNghb,nDefNorm
 								    real(pReal), parameter  :: &
 								      nDefToler = 1.0e-10_pReal, &
 								      b = 2.5e-10_pReal                                                                             ! Length of Burgers vector
-												cleaning

											
										
										
											2018-11-04 03:53:52 +05:30
-												storing by instance just complicates the code

											
										
										
											2021-02-23 17:47:51 +05:30
+								    nGDim = param(ho)%N_constituents
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    rPen = 0.0_pReal
 								    nMis = 0.0_pReal
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    !----------------------------------------------------------------------------------------------
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
+								    ! get the correction factor the modulus of penalty stress representing the evolution of area of
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    ! the interfaces due to deformations
-												simplifying

											
										
										
											2018-11-04 03:43:20 +05:30
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								    surfCorr = surfaceCorrection(avgF,ho,en)
-												cleaning

											
										
										
											2018-11-04 03:53:52 +05:30
-												cleaning homogenization

											
										
										
											2021-02-26 02:12:40 +05:30
+								    associate(prm => param(ho))
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								   !-----------------------------------------------------------------------------------------------
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
+								   ! computing the mismatch and penalty stress tensor of all grains
-												homogenization modules made consistent

											
										
										
											2020-09-23 05:03:19 +05:30
+								   grainLoop: do iGrain = 1,product(prm%N_constituents)
-												instance less; use cell mapping

											
										
										
											2021-02-22 20:47:32 +05:30
+								     muGrain = equivalentMu(iGrain,ce)
-												too early (depends on IP)

											
										
										
											2020-12-28 14:57:52 +05:30
+								     iGrain3 = grain1to3(iGrain,prm%N_constituents)                                                 ! get the grain ID in local 3-dimensional index (x,y,z)-position
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								     interfaceLoop: do iFace = 1,6
 								       intFace = getInterface(iFace,iGrain3)                                                        ! get the 4-dimensional index of the interface in local numbering system of the grain
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								       nVect = interfaceNormal(intFace,ho,en)
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								       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))
 								       where(iGNghb3 < 1)    iGNghb3 = nGDim
 								       where(iGNghb3 >nGDim) iGNghb3 = 1
-												homogenization modules made consistent

											
										
										
											2020-09-23 05:03:19 +05:30
+								       iGNghb  = grain3to1(iGNghb3,prm%N_constituents)                                              ! get the ID of the neighboring grain
-												instance less; use cell mapping

											
										
										
											2021-02-22 20:47:32 +05:30
+								       muGNghb = equivalentMu(iGNghb,ce)
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +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
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								       !-------------------------------------------------------------------------------------------
 								       ! compute the mismatch tensor of all interfaces
 								       nDefNorm = 0.0_pReal
 								       nDef = 0.0_pReal
 								       do i = 1,3; do j = 1,3
 								         do k = 1,3; do l = 1,3
-												better readable, tested, and following standard notation

											
										
										
											2020-03-02 03:30:06 +05:30
+								           nDef(i,j) = nDef(i,j) - nVect(k)*gDef(i,l)*math_LeviCivita(j,k,l)                        ! compute the interface mismatch tensor from the jump of deformation gradient
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								         end do; end do
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								         nDefNorm = nDefNorm + nDef(i,j)**2.0_pReal                                                 ! compute the norm of the mismatch tensor
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								       end do; end do
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +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)
-												don't clutter with statements that are never used

											
										
										
											2020-12-23 18:40:26 +05:30
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								       !-------------------------------------------------------------------------------------------
 								       ! compute the stress penalty of all interfaces
 								       do i = 1,3; do j = 1,3; do k = 1,3; do l = 1,3
-												using central functionality

											
										
										
											2020-12-28 15:33:29 +05:30
+								         rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*b + muGNghb*b)*prm%xi_alpha &
-												homogenization modules made consistent

											
										
										
											2020-09-23 05:03:19 +05:30
+								                                                *surfCorr(abs(intFace(1)))/prm%D_alpha(abs(intFace(1))) &
 								                                                *cosh(prm%c_alpha*nDefNorm) &
-												better readable, tested, and following standard notation

											
										
										
											2020-03-02 03:30:06 +05:30
+								                                                *0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_LeviCivita(k,l,j) &
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
+								                                                *tanh(nDefNorm/num%xSmoo)
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								       end do; end do;enddo; end do
 								     end do interfaceLoop
-												don't clutter with statements that are never used

											
										
										
											2020-12-23 18:40:26 +05:30
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								   end do grainLoop
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								   end associate
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  end subroutine stressPenalty
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  !------------------------------------------------------------------------------------------------
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
+								  !> @brief calculate stress-like penalty due to volume discrepancy
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  !------------------------------------------------------------------------------------------------
-												cleaning homogenization

											
										
										
											2021-02-26 02:12:40 +05:30
+								  subroutine volumePenalty(vPen,vDiscrep,fAvg,fDef,nGrain)
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    real(pReal), dimension (:,:,:), intent(out) :: vPen                                             ! stress-like penalty due to volume
 								    real(pReal),                    intent(out) :: vDiscrep                                         ! total volume discrepancy
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    real(pReal), dimension (:,:,:), intent(in)  :: fDef                                             ! deformation gradients
 								    real(pReal), dimension (3,3),   intent(in)  :: fAvg                                             ! overall deformation gradient
 								    integer,                        intent(in) :: &
-												cleaning homogenization

											
										
										
											2021-02-26 02:12:40 +05:30
+								      Ngrain
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    real(pReal), dimension(size(vPen,3)) :: gVol
 								    integer :: i
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    !----------------------------------------------------------------------------------------------
 								    ! compute the volumes of grains and of cluster
 								    vDiscrep = math_det33(fAvg)                                                                     ! compute the volume of the cluster
 								    do i = 1,nGrain
 								      gVol(i) = math_det33(fDef(1:3,1:3,i))                                                         ! compute the volume of individual grains
 								      vDiscrep     = vDiscrep - gVol(i)/real(nGrain,pReal)                                          ! calculate the difference/dicrepancy between
-												polishing

											
										
										
											2019-01-13 03:37:35 +05:30
+								                                                                                                    ! the volume of the cluster and the the total volume of grains
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								    end do
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    !----------------------------------------------------------------------------------------------
 								    ! calculate the stress and penalty due to volume discrepancy
 								    vPen      = 0.0_pReal
 								    do i = 1,nGrain
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
+								      vPen(:,:,i) = -1.0_pReal/real(nGrain,pReal)*num%volDiscrMod*num%volDiscrPow/num%maxVolDiscr* &
 								                         sign((abs(vDiscrep)/num%maxVolDiscr)**(num%volDiscrPow - 1.0),vDiscrep)* &
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								                         gVol(i)*transpose(math_inv33(fDef(:,:,i)))
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								    end do
-												cleaning

											
										
										
											2018-11-04 02:06:49 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  end subroutine volumePenalty
-												cleaning

											
										
										
											2018-11-04 02:06:49 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  !--------------------------------------------------------------------------------------------------
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
+								  !> @brief compute the correction factor accouted for surface evolution (area change) due to
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  ! deformation
 								  !--------------------------------------------------------------------------------------------------
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								  function surfaceCorrection(avgF,ho,en)
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    real(pReal), dimension(3)               :: surfaceCorrection
 								    real(pReal), dimension(3,3), intent(in) :: avgF                                                 !< average F
 								    integer,                     intent(in) :: &
-												cleaning homogenization

											
										
										
											2021-02-26 02:12:40 +05:30
+								      ho, &
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								      en
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    real(pReal), dimension(3,3)             :: invC
 								    real(pReal), dimension(3)               :: nVect
 								    real(pReal)  :: detF
 								    integer :: i,j,iBase
 								    logical     ::  error
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												following standard: Output first

											
										
										
											2019-09-21 06:58:46 +05:30
+								    call math_invert33(invC,detF,error,matmul(transpose(avgF),avgF))
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    surfaceCorrection = 0.0_pReal
 								    do iBase = 1,3
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								      nVect = interfaceNormal([iBase,1,1,1],ho,en)
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								      do i = 1,3; do j = 1,3
 								        surfaceCorrection(iBase) = surfaceCorrection(iBase) + invC(i,j)*nVect(i)*nVect(j)           ! compute the component of (the inverse of) the stretch in the direction of the normal
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								      end do; end do
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								      surfaceCorrection(iBase) = sqrt(surfaceCorrection(iBase))*detF                                ! get the surface correction factor (area contraction/enlargement)
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								    end do
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  end function surfaceCorrection
-												cleaning

											
										
										
											2018-11-04 02:06:49 +05:30
-												using central functionality

											
										
										
											2020-12-28 15:33:29 +05:30
+								  !-------------------------------------------------------------------------------------------------
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  !> @brief compute the equivalent shear and bulk moduli from the elasticity tensor
-												using central functionality

											
										
										
											2020-12-28 15:33:29 +05:30
+								  !-------------------------------------------------------------------------------------------------
-												instance less; use cell mapping

											
										
										
											2021-02-22 20:47:32 +05:30
+								  real(pReal) function equivalentMu(grainID,ce)
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    integer, intent(in)    :: &
 								      grainID,&
-												cleaning homogenization

											
										
										
											2021-02-26 02:12:40 +05:30
+								      ce
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												consistent access via (ph)ase and phase(me)mber

											
										
										
											2021-01-19 14:44:34 +05:30
+								    real(pReal), dimension(6,6) :: C
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												consistent access via (ph)ase and phase(me)mber

											
										
										
											2021-01-19 14:44:34 +05:30
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								    C = phase_homogenizedC(material_phaseID(grainID,ce),material_phaseEntry(grainID,ce))
-												consistent access via (ph)ase and phase(me)mber

											
										
										
											2021-01-19 14:44:34 +05:30
+								    equivalentMu = lattice_equivalent_mu(C,'voigt')
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												using central functionality

											
										
										
											2020-12-28 15:33:29 +05:30
+								  end function equivalentMu
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												using central functionality

											
										
										
											2020-12-28 15:33:29 +05:30
+								  !-------------------------------------------------------------------------------------------------
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
+								  !> @brief calculating the grain deformation gradient (the same with
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  ! homogenization_RGC_partitionDeformation, but used only for perturbation scheme)
-												using central functionality

											
										
										
											2020-12-28 15:33:29 +05:30
+								  !-------------------------------------------------------------------------------------------------
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								  subroutine grainDeformation(F, avgF, ho, en)
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												partition, not partion

											
										
										
											2020-10-08 01:45:13 +05:30
+								    real(pReal),   dimension(:,:,:), intent(out) :: F                                               !< partitioned F  per grain
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    real(pReal),   dimension(:,:),   intent(in)  :: avgF                                            !< averaged F
 								    integer,                          intent(in)  :: &
-												cleaning homogenization

											
										
										
											2021-02-26 02:12:40 +05:30
+								      ho, &
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								      en
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    real(pReal),   dimension(3) :: aVect,nVect
 								    integer,       dimension(4) :: intFace
 								    integer,       dimension(3) :: iGrain3
 								    integer :: iGrain,iFace,i,j
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												using central functionality

											
										
										
											2020-12-28 15:33:29 +05:30
+								    !-----------------------------------------------------------------------------------------------
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    ! compute the deformation gradient of individual grains due to relaxations
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												cleaning homogenization

											
										
										
											2021-02-26 02:12:40 +05:30
+								    associate (prm => param(ho))
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    F = 0.0_pReal
-												homogenization modules made consistent

											
										
										
											2020-09-23 05:03:19 +05:30
+								    do iGrain = 1,product(prm%N_constituents)
 								      iGrain3 = grain1to3(iGrain,prm%N_constituents)
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								      do iFace = 1,6
 								        intFace = getInterface(iFace,iGrain3)
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								        aVect   = relaxationVector(intFace,ho,en)
 								        nVect   = interfaceNormal(intFace,ho,en)
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								        forall (i=1:3,j=1:3) &
 								          F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j)                                         ! effective relaxations
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								      end do
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								      F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF                                                  ! relaxed deformation gradient
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								    end do
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    end associate
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  end subroutine grainDeformation
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												shorter names

need to prefix 'pass' to avoid name clashes that result in errors during
compilation

											
										
										
											2021-04-06 15:35:47 +05:30
+								end function RGC_updateState
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
-												preparing output of mech_homog to HDF5

											
										
										
											2019-04-30 22:15:16 +05:30
+								!--------------------------------------------------------------------------------------------------
 								!> @brief writes results to HDF5 output file
 								!--------------------------------------------------------------------------------------------------
-												shorter names

need to prefix 'pass' to avoid name clashes that result in errors during
compilation

											
										
										
											2021-04-06 15:35:47 +05:30
+								module subroutine RGC_results(ho,group)
-												preparing output of mech_homog to HDF5

											
										
										
											2019-04-30 22:15:16 +05:30
-												storing by instance just complicates the code

											
										
										
											2021-02-23 17:47:51 +05:30
+								  integer,          intent(in) :: ho
-												RGC as submodule

submodules inherit use-associated entities and implicit none/private
statements

											
										
										
											2019-05-18 10:53:46 +05:30
+								  character(len=*), intent(in) :: group
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												preparing output of mech_homog to HDF5

											
										
										
											2019-04-30 22:15:16 +05:30
+								  integer :: o
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												storing by instance just complicates the code

											
										
										
											2021-02-23 17:47:51 +05:30
+								  associate(stt => state(ho), dst => dependentState(ho), prm => param(ho))
-												less writing, same functionality

											
										
										
											2020-02-15 03:51:58 +05:30
+								  outputsLoop: do o = 1,size(prm%output)
 								    select case(trim(prm%output(o)))
-												ensure defined initital state

											
										
										
											2020-09-07 16:50:00 +05:30
+								      case('M')
-												more meaningful order

											
										
										
											2021-06-01 20:35:13 +05:30
+								        call results_writeDataset(dst%mismatch,group,trim(prm%output(o)), &
-												[skip ci] mismatch tensor is dimensionless

											
										
										
											2019-05-13 21:03:45 +05:30
+								                                  'average mismatch tensor','1')
-												Restructuring for material.yaml

											
										
										
											2020-08-15 19:32:10 +05:30
+								      case('Delta_V')
-												more meaningful order

											
										
										
											2021-06-01 20:35:13 +05:30
+								        call results_writeDataset(dst%volumeDiscrepancy,group,trim(prm%output(o)), &
-												write all outputs to HDF5

not 100% sure about the units

											
										
										
											2019-05-12 18:44:29 +05:30
+								                                  'volume discrepancy','m³')
-												consistently use prefix as modifier

											
										
										
											2021-06-29 03:11:50 +05:30
+								      case('max_dot_a')
-												more meaningful order

											
										
										
											2021-06-01 20:35:13 +05:30
+								        call results_writeDataset(dst%relaxationrate_max,group,trim(prm%output(o)), &
-												write all outputs to HDF5

not 100% sure about the units

											
										
										
											2019-05-12 18:44:29 +05:30
+								                                  'maximum relaxation rate','m/s')
-												consistently use prefix as modifier

											
										
										
											2021-06-29 03:11:50 +05:30
+								      case('avg_dot_a')
-												more meaningful order

											
										
										
											2021-06-01 20:35:13 +05:30
+								        call results_writeDataset(dst%relaxationrate_avg,group,trim(prm%output(o)), &
-												write all outputs to HDF5

not 100% sure about the units

											
										
										
											2019-05-12 18:44:29 +05:30
+								                                  'average relaxation rate','m/s')
-												preparing output of mech_homog to HDF5

											
										
										
											2019-04-30 22:15:16 +05:30
+								    end select
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								  end do outputsLoop
-												preparing output of mech_homog to HDF5

											
										
										
											2019-04-30 22:15:16 +05:30
+								  end associate
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												shorter names

need to prefix 'pass' to avoid name clashes that result in errors during
compilation

											
										
										
											2021-04-06 15:35:47 +05:30
+								end subroutine RGC_results
-												preparing output of mech_homog to HDF5

											
										
										
											2019-04-30 22:15:16 +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 collect relaxation vectors of an interface
 								!--------------------------------------------------------------------------------------------------
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								pure function relaxationVector(intFace,ho,en)
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  real(pReal), dimension (3)            :: relaxationVector
-												no pInt

											
										
										
											2019-04-06 00:15:56 +05:30
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								  integer,                   intent(in) :: ho,en
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  integer,     dimension(4), intent(in) :: intFace                                                  !< set of interface ID in 4D array (normal and position)
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  integer :: 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
-												cleaning homogenization

											
										
										
											2021-02-26 02:12:40 +05:30
+								  associate (prm => param(ho), &
 								             stt => state(ho))
-												use cell mapping

											
										
										
											2021-02-15 23:13:51 +05:30
 								  iNum = interface4to1(intFace,prm%N_constituents)                                                  ! identify the position of the interface in global state array
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  if (iNum > 0) then
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								    relaxationVector = stt%relaxationVector((3*iNum-2):(3*iNum),en)
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  else
 								    relaxationVector = 0.0_pReal
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								  end if
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
-												use cell mapping

											
										
										
											2021-02-15 23:13:51 +05:30
+								  end associate
-												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
 								!--------------------------------------------------------------------------------------------------
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
+								!> @brief identify the normal of an interface
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
+								!--------------------------------------------------------------------------------------------------
-												lessons learned from compilation with Intel 2021.3

											
										
										
											2021-08-11 03:17:13 +05:30
+								pure function interfaceNormal(intFace,ho,en) result(n)
-												no pInt

											
										
										
											2019-04-06 00:15:56 +05:30
-												lessons learned from compilation with Intel 2021.3

											
										
										
											2021-08-11 03:17:13 +05:30
+								  real(pReal), dimension(3)             :: n
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  integer,     dimension(4), intent(in) :: intFace                                                  !< interface ID in 4D array (normal and position)
 								  integer,                   intent(in) :: &
-												cleaning homogenization

											
										
										
											2021-02-26 02:12:40 +05:30
+								    ho, &
-												new names

											
										
										
											2021-04-06 15:08:44 +05:30
+								    en
-												no pInt

											
										
										
											2019-04-06 00:15:56 +05:30
-												lessons learned from compilation with Intel 2021.3

											
										
										
											2021-08-11 03:17:13 +05:30
-												cleaning homogenization

											
										
										
											2021-02-26 02:12:40 +05:30
+								  associate (dst => dependentState(ho))
-												# 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
-												lessons learned from compilation with Intel 2021.3

											
										
										
											2021-08-11 03:17:13 +05:30
+								    n = 0.0_pReal
-												correct line continuation

											
										
										
											2021-08-11 12:15:24 +05:30
+								    n(abs(intFace(1))) = 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
-												correct line continuation

											
										
										
											2021-08-11 12:15:24 +05:30
+								    n = matmul(dst%orientation(1:3,1:3,en),n)                                                       ! map the normal vector into sample coordinate system (basis)
-												avoid repetition

											
										
										
											2021-04-11 15:48:26 +05:30
-												use cell mapping

											
										
										
											2021-02-15 23:13:51 +05:30
+								  end associate
-												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)
 								!--------------------------------------------------------------------------------------------------
-												lessons learned from compilation with Intel 2021.3

											
										
										
											2021-08-11 03:17:13 +05:30
+								pure function getInterface(iFace,iGrain3) result(i)
-												no pInt

											
										
										
											2019-04-06 00:15:56 +05:30
-												lessons learned from compilation with Intel 2021.3

											
										
										
											2021-08-11 03:17:13 +05:30
+								  integer, dimension(4)             :: i
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  integer, dimension(3), intent(in) :: iGrain3                                                      !< grain ID in 3D array
 								  integer,               intent(in) :: iFace                                                        !< face index (1..6) mapped like (-e1,-e2,-e3,+e1,+e2,+e3) or iDir = (-1,-2,-3,1,2,3)
-												no pInt

											
										
										
											2019-04-06 00:15:56 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  integer :: iDir                                                                                   !< direction of interface normal
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												lessons learned from compilation with Intel 2021.3

											
										
										
											2021-08-11 03:17:13 +05:30
+								  iDir = (int(real(iFace-1,pReal)/2.0_pReal)+1)*(-1)**iFace
 								  i = [iDir,iGrain3]
 								  if (iDir < 0) i(1-iDir) = i(1-iDir)-1                                       ! 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)
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  integer, dimension(3)             :: grain1to3
-												no pInt

											
										
										
											2019-04-06 00:15:56 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  integer,               intent(in) :: grain1                                                       !< grain ID in 1D array
 								  integer, 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
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  grain1to3 = 1 + [mod((grain1-1), nGDim(1)), &
 								                   mod((grain1-1)/ nGDim(1),nGDim(2)), &
-												no pInt

											
										
										
											2019-04-06 00:15:56 +05:30
+								                       (grain1-1)/(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)
 								!--------------------------------------------------------------------------------------------------
-												no pInt

											
										
										
											2019-04-06 00:15:56 +05:30
+								integer 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
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  integer, dimension(3), intent(in) :: grain3                                                       !< grain ID in 3D array (pos.x,pos.y,pos.z)
 								  integer, dimension(3), intent(in) :: nGDim
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  grain3to1 = grain3(1) &
 								            + nGDim(1)*(grain3(2)-1) &
 								            + nGDim(1)*nGDim(2)*(grain3(3)-1)
-												# 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)
 								!--------------------------------------------------------------------------------------------------
-												no pInt

											
										
										
											2019-04-06 00:15:56 +05:30
+								integer pure function interface4to1(iFace4D, nGDim)
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  integer, dimension(4), intent(in) :: iFace4D                                                      !< interface ID in 4D array (n.dir,pos.x,pos.y,pos.z)
 								  integer, dimension(3), intent(in) :: nGDim
-. Adding more comments
2. Changing 'subroutine' into 'function' for optimizing the computation time

											
										
										
											2010-11-26 17:20:20 +05:30
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  select case(abs(iFace4D(1)))
 								    case(1)
 								      if ((iFace4D(2) == 0) .or. (iFace4D(2) == nGDim(1))) then
 								        interface4to1 = 0
 								      else
 								        interface4to1 = iFace4D(3) + nGDim(2)*(iFace4D(4)-1) &
 								                      + nGDim(2)*nGDim(3)*(iFace4D(2)-1)
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								      end if
-												better readable

											
										
										
											2018-11-04 13:49:24 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    case(2)
 								      if ((iFace4D(3) == 0) .or. (iFace4D(3) == nGDim(2))) then
 								        interface4to1 = 0
 								      else
 								        interface4to1 = iFace4D(4) + nGDim(3)*(iFace4D(2)-1) &
 								                      + nGDim(3)*nGDim(1)*(iFace4D(3)-1) &
 								                      + (nGDim(1)-1)*nGDim(2)*nGDim(3)                                              ! total # of interfaces normal || e1
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								      end if
-												better readable

											
										
										
											2018-11-04 13:49:24 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    case(3)
 								      if ((iFace4D(4) == 0) .or. (iFace4D(4) == nGDim(3))) then
 								        interface4to1 = 0
 								      else
 								        interface4to1 = iFace4D(2) + nGDim(1)*(iFace4D(3)-1) &
 								                      + nGDim(1)*nGDim(2)*(iFace4D(4)-1) &
 								                      + (nGDim(1)-1)*nGDim(2)*nGDim(3) &                                            ! total # of interfaces normal || e1
 								                      + nGDim(1)*(nGDim(2)-1)*nGDim(3)                                              ! total # of interfaces normal || e2
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								      end if
-												better readable

											
										
										
											2018-11-04 13:49:24 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								    case default
 								      interface4to1 = -1
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  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)
-												store data where it is needed

avoid globals, even if they are read only

											
										
										
											2020-03-17 02:09:53 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  integer, dimension(4)             :: interface1to4
-												no pInt

											
										
										
											2019-04-06 00:15:56 +05:30
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  integer,               intent(in) :: iFace1D                                                      !< interface ID in 1D array
 								  integer, dimension(3), intent(in) :: nGDim
 								  integer, 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 ...
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  nIntFace = [(nGDim(1)-1)*nGDim(2)*nGDim(3), &                                                     ! ... normal || e1
 								              nGDim(1)*(nGDim(2)-1)*nGDim(3), &                                                     ! ... normal || e2
 								              nGDim(1)*nGDim(2)*(nGDim(3)-1)]                                                       ! ... normal || e3
-												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)
-												white space adjustments

											
										
										
											2019-06-15 21:57:38 +05:30
+								  if (iFace1D > 0 .and. iFace1D <= nIntFace(1)) then                                                ! interface with normal || e1
 								    interface1to4(1) = 1
 								    interface1to4(3) = mod((iFace1D-1),nGDim(2))+1
 								    interface1to4(4) = mod(int(real(iFace1D-1,pReal)/real(nGDim(2),pReal)),nGDim(3))+1
 								    interface1to4(2) = int(real(iFace1D-1,pReal)/real(nGDim(2),pReal)/real(nGDim(3),pReal))+1
 								  elseif (iFace1D > nIntFace(1) .and. iFace1D <= (nIntFace(2) + nIntFace(1))) then                  ! interface with normal || e2
 								    interface1to4(1) = 2
 								    interface1to4(4) = mod((iFace1D-nIntFace(1)-1),nGDim(3))+1
 								    interface1to4(2) = mod(int(real(iFace1D-nIntFace(1)-1,pReal)/real(nGDim(3),pReal)),nGDim(1))+1
 								    interface1to4(3) = int(real(iFace1D-nIntFace(1)-1,pReal)/real(nGDim(3),pReal)/real(nGDim(1),pReal))+1
 								  elseif (iFace1D > nIntFace(2) + nIntFace(1) .and. iFace1D <= (nIntFace(3) + nIntFace(2) + nIntFace(1))) then ! interface with normal || e3
 								    interface1to4(1) = 3
 								    interface1to4(2) = mod((iFace1D-nIntFace(2)-nIntFace(1)-1),nGDim(1))+1
 								    interface1to4(3) = mod(int(real(iFace1D-nIntFace(2)-nIntFace(1)-1,pReal)/real(nGDim(1),pReal)),nGDim(2))+1
 								    interface1to4(4) = int(real(iFace1D-nIntFace(2)-nIntFace(1)-1,pReal)/real(nGDim(1),pReal)/real(nGDim(2),pReal))+1
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								  end if
-												# 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 interface1to4
-												doxygen documentation, enabled uncommented debug statement to be turned on by switches, minor polishing

											
										
										
											2013-02-21 03:36:15 +05:30
-												hierarchical naming

											
										
										
											2021-01-26 16:11:19 +05:30
+								end submodule RGC