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

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

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

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

  print'(/,a,i0)', ' # homogenizations: ',count(mechanical_type == MECHANICAL_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_dict('homogenization')
  allocate(param(material_homogenization%length))
  allocate(state(material_homogenization%length))
  allocate(state0(material_homogenization%length))
  allocate(dependentState(material_homogenization%length))

  num_homogenization => config_numerics%get_dict('homogenization',defaultVal=emptyDict)
  num_mechanical => num_homogenization%get_dict('mechanical',defaultVal=emptyDict)
  num_RGC => num_mechanical%get_dict('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(mechanical_type)
    if (mechanical_type(ho) /= MECHANICAL_RGC_ID) cycle
    homog => material_homogenization%get_dict(ho)
    homogMech => homog%get_dict('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;end do;end do;end do
! 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;end do;end do;end do
      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                                                         ! 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;end do; 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) = -real(nGrain,pReal)**(-1)*num%volDiscrMod*num%volDiscrPow/num%maxVolDiscr &
                  * sign((abs(vDiscrep)/num%maxVolDiscr)**(num%volDiscrPow - 1.0_pReal),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(co,ce)

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

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


    C = phase_homogenizedC66(material_phaseID(co,ce),material_phaseEntry(co,ce))                    ! damage not included!
    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
+								!--------------------------------------------------------------------------------------------------
-												polishing

RGC numerics is still annoying (and was probably never used in the last
10 years)

											
										
										
											2022-02-19 18:26:24 +05:30
+								module subroutine RGC_init()
-												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
-												Re-written YAML types
Strict typing for YAML

New access pattern requires to specify the expected type, i.e. 'scalar', 'list', or 'dict'. This ensures that the node offers the expected functionality instead of polluting 'tNode' with dummy functions which throw error messages if not overwritten.

The restructuring of the code allows to hierarchically construct methods without much code duplication.

Some aspects of the error messaging system have been improved.

											
										
										
											2022-10-25 21:39:36 +05:30
+								  class(tDict), pointer :: &
-												polishing

RGC numerics is still annoying (and was probably never used in the last
10 years)

											
										
										
											2022-02-19 18:26:24 +05:30
+								    num_homogenization, &
 								    num_mechanical, &
-												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
-												check N_constituents only if active

											
										
										
											2022-02-24 16:39:47 +05:30
+								  print'(/,a,i0)', ' # homogenizations: ',count(mechanical_type == MECHANICAL_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
-												Re-written YAML types
Strict typing for YAML

New access pattern requires to specify the expected type, i.e. 'scalar', 'list', or 'dict'. This ensures that the node offers the expected functionality instead of polluting 'tNode' with dummy functions which throw error messages if not overwritten.

The restructuring of the code allows to hierarchically construct methods without much code duplication.

Some aspects of the error messaging system have been improved.

											
										
										
											2022-10-25 21:39:36 +05:30
+								  material_homogenization => config_material%get_dict('homogenization')
-												storing by instance just complicates the code

											
										
										
											2021-02-23 17:47:51 +05:30
+								  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
-												Re-written YAML types
Strict typing for YAML

New access pattern requires to specify the expected type, i.e. 'scalar', 'list', or 'dict'. This ensures that the node offers the expected functionality instead of polluting 'tNode' with dummy functions which throw error messages if not overwritten.

The restructuring of the code allows to hierarchically construct methods without much code duplication.

Some aspects of the error messaging system have been improved.

											
										
										
											2022-10-25 21:39:36 +05:30
+								  num_homogenization => config_numerics%get_dict('homogenization',defaultVal=emptyDict)
 								  num_mechanical => num_homogenization%get_dict('mechanical',defaultVal=emptyDict)
 								  num_RGC => num_mechanical%get_dict('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
-												check N_constituents only if active

											
										
										
											2022-02-24 16:39:47 +05:30
+								  do ho = 1, size(mechanical_type)
 								    if (mechanical_type(ho) /= MECHANICAL_RGC_ID) cycle
-												Re-written YAML types
Strict typing for YAML

New access pattern requires to specify the expected type, i.e. 'scalar', 'list', or 'dict'. This ensures that the node offers the expected functionality instead of polluting 'tNode' with dummy functions which throw error messages if not overwritten.

The restructuring of the code allows to hierarchically construct methods without much code duplication.

Some aspects of the error messaging system have been improved.

											
										
										
											2022-10-25 21:39:36 +05:30
+								    homog => material_homogenization%get_dict(ho)
 								    homogMech => homog%get_dict('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
-												polishing

RGC numerics is still annoying (and was probably never used in the last
10 years)

											
										
										
											2022-02-19 18:26:24 +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)
-												more consistent formatting of loops

											
										
										
											2021-11-17 13:41:30 +05:30
+								        end do;end do;end do;end do
-												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)
-												more consistent formatting of loops

											
										
										
											2021-11-17 13:41:30 +05:30
+								        end do;end do;end do;end do
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								      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
-												integer exponents: faster and shorter

											
										
										
											2021-11-26 01:22:22 +05:30
+								         nDefNorm = nDefNorm + nDef(i,j)**2                                                         ! 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 variants including space

											
										
										
											2022-06-09 02:36:01 +05:30
+								       end do; end do;end do; end do
-												consistent indentation and line-spacings in reporting

											
										
										
											2021-11-15 23:05:44 +05:30
+								     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
-												Marc precision handling

											
										
										
											2022-06-27 14:07:41 +05:30
+								      vPen(:,:,i) = -real(nGrain,pReal)**(-1)*num%volDiscrMod*num%volDiscrPow/num%maxVolDiscr &
 								                  * sign((abs(vDiscrep)/num%maxVolDiscr)**(num%volDiscrPow - 1.0_pReal),vDiscrep) &
 								                  * 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
+								  !-------------------------------------------------------------------------------------------------
-												avoid conversions

											
										
										
											2021-11-19 02:29:09 +05:30
+								  real(pReal) function equivalentMu(co,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)    :: &
-												avoid conversions

											
										
										
											2021-11-19 02:29:09 +05:30
+								      co,&
-												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
-												using Voigt notation

This is a bugfix with a change of behavior

											
										
										
											2021-11-19 14:27:40 +05:30
+								    C = phase_homogenizedC66(material_phaseID(co,ce),material_phaseEntry(co,ce))                    ! damage not included!
 								    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