!-------------------------------------------------------------------------------------------------- !> @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) homogenization_mech_RGC use rotations type :: tParameters integer, dimension(:), allocatable :: & N_constituents real(pReal) :: & xi_alpha, & c_Alpha real(pReal), dimension(:), allocatable :: & D_alpha, & a_g integer :: & of_debug = 0 character(len=pStringLen), allocatable, dimension(:) :: & output end type tParameters type :: tRGCstate real(pReal), pointer, dimension(:) :: & work, & penaltyEnergy 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 mech_RGC_init(num_homogMech) class(tNode), pointer, intent(in) :: & num_homogMech !< pointer to mechanical homogenization numerics data integer :: & Ninstance, & h, & NofMyHomog, & sizeState, nIntFaceTot class (tNode), pointer :: & num_RGC, & ! pointer to RGC numerics data material_homogenization, & homog, & homogMech print'(/,a)', ' <<<+- homogenization_mech_rgc init -+>>>' Ninstance = count(homogenization_type == HOMOGENIZATION_RGC_ID) print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print*, 'Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009' print*, 'https://doi.org/10.1007/s12289-009-0619-1'//IO_EOL print*, 'Tjahjanto et al., Modelling and Simulation in Materials Science and Engineering 18:015006, 2010' print*, 'https://doi.org/10.1088/0965-0393/18/1/015006'//IO_EOL allocate(param(Ninstance)) allocate(state(Ninstance)) allocate(state0(Ninstance)) allocate(dependentState(Ninstance)) num_RGC => num_homogMech%get('RGC',defaultVal=emptyDict) num%atol = num_RGC%get_asFloat('atol', defaultVal=1.0e+4_pReal) num%rtol = num_RGC%get_asFloat('rtol', defaultVal=1.0e-3_pReal) num%absMax = num_RGC%get_asFloat('amax', defaultVal=1.0e+10_pReal) num%relMax = num_RGC%get_asFloat('rmax', defaultVal=1.0e+2_pReal) num%pPert = num_RGC%get_asFloat('perturbpenalty', defaultVal=1.0e-7_pReal) num%xSmoo = num_RGC%get_asFloat('relvantmismatch', defaultVal=1.0e-5_pReal) num%viscPower = num_RGC%get_asFloat('viscositypower', defaultVal=1.0e+0_pReal) num%viscModus = num_RGC%get_asFloat('viscositymodulus', defaultVal=0.0e+0_pReal) num%refRelaxRate = num_RGC%get_asFloat('refrelaxationrate', defaultVal=1.0e-3_pReal) num%maxdRelax = num_RGC%get_asFloat('maxrelaxationrate', defaultVal=1.0e+0_pReal) num%maxVolDiscr = num_RGC%get_asFloat('maxvoldiscrepancy', defaultVal=1.0e-5_pReal) num%volDiscrMod = num_RGC%get_asFloat('voldiscrepancymod', defaultVal=1.0e+12_pReal) num%volDiscrPow = num_RGC%get_asFloat('dicrepancypower', defaultVal=5.0_pReal) if (num%atol <= 0.0_pReal) call IO_error(301,ext_msg='absTol_RGC') if (num%rtol <= 0.0_pReal) call IO_error(301,ext_msg='relTol_RGC') if (num%absMax <= 0.0_pReal) call IO_error(301,ext_msg='absMax_RGC') if (num%relMax <= 0.0_pReal) call IO_error(301,ext_msg='relMax_RGC') if (num%pPert <= 0.0_pReal) call IO_error(301,ext_msg='pPert_RGC') if (num%xSmoo <= 0.0_pReal) call IO_error(301,ext_msg='xSmoo_RGC') if (num%viscPower < 0.0_pReal) call IO_error(301,ext_msg='viscPower_RGC') if (num%viscModus < 0.0_pReal) call IO_error(301,ext_msg='viscModus_RGC') if (num%refRelaxRate <= 0.0_pReal) call IO_error(301,ext_msg='refRelaxRate_RGC') if (num%maxdRelax <= 0.0_pReal) call IO_error(301,ext_msg='maxdRelax_RGC') if (num%maxVolDiscr <= 0.0_pReal) call IO_error(301,ext_msg='maxVolDiscr_RGC') if (num%volDiscrMod < 0.0_pReal) call IO_error(301,ext_msg='volDiscrMod_RGC') if (num%volDiscrPow <= 0.0_pReal) call IO_error(301,ext_msg='volDiscrPw_RGC') material_homogenization => config_material%get('homogenization') do h = 1, size(homogenization_type) if (homogenization_type(h) /= HOMOGENIZATION_RGC_ID) cycle homog => material_homogenization%get(h) homogMech => homog%get('mech') associate(prm => param(homogenization_typeInstance(h)), & stt => state(homogenization_typeInstance(h)), & st0 => state0(homogenization_typeInstance(h)), & dst => dependentState(homogenization_typeInstance(h))) #ifdef DEBUG if (h==material_homogenizationAt(debugHomog%element)) then prm%of_debug = material_homogenizationMemberAt(debugHomog%ip,debugHomog%element) endif #endif #if defined (__GFORTRAN__) prm%output = output_asStrings(homogMech) #else prm%output = homogMech%get_asStrings('output',defaultVal=emptyStringArray) #endif prm%N_constituents = homogMech%get_asInts('cluster_size',requiredSize=3) if (homogenization_Nconstituent(h) /= product(prm%N_constituents)) & call IO_error(211,ext_msg='N_constituents (mech_RGC)') prm%xi_alpha = homogMech%get_asFloat('xi_alpha') prm%c_alpha = homogMech%get_asFloat('c_alpha') prm%D_alpha = homogMech%get_asFloats('D_alpha', requiredSize=3) prm%a_g = homogMech%get_asFloats('a_g', requiredSize=3) NofMyHomog = count(material_homogenizationAt == h) 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 & + size(['avg constitutive work ','average penalty energy']) homogState(h)%sizeState = sizeState allocate(homogState(h)%state0 (sizeState,NofMyHomog), source=0.0_pReal) allocate(homogState(h)%subState0(sizeState,NofMyHomog), source=0.0_pReal) allocate(homogState(h)%state (sizeState,NofMyHomog), source=0.0_pReal) stt%relaxationVector => homogState(h)%state(1:nIntFaceTot,:) st0%relaxationVector => homogState(h)%state0(1:nIntFaceTot,:) stt%work => homogState(h)%state(nIntFaceTot+1,:) stt%penaltyEnergy => homogState(h)%state(nIntFaceTot+2,:) allocate(dst%volumeDiscrepancy( NofMyHomog), source=0.0_pReal) allocate(dst%relaxationRate_avg( NofMyHomog), source=0.0_pReal) allocate(dst%relaxationRate_max( NofMyHomog), source=0.0_pReal) allocate(dst%mismatch( 3,NofMyHomog), source=0.0_pReal) !-------------------------------------------------------------------------------------------------- ! assigning cluster orientations dependentState(homogenization_typeInstance(h))%orientation = spread(eu2om(prm%a_g*inRad),3,NofMyHomog) !dst%orientation = spread(eu2om(prm%a_g*inRad),3,NofMyHomog) ifort version 18.0.1 crashes (for whatever reason) end associate enddo end subroutine mech_RGC_init !-------------------------------------------------------------------------------------------------- !> @brief partitions the deformation gradient onto the constituents !-------------------------------------------------------------------------------------------------- module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of) real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned F per grain real(pReal), dimension (3,3), intent(in) :: avgF !< averaged F integer, intent(in) :: & instance, & of real(pReal), dimension(3) :: aVect,nVect integer, dimension(4) :: intFace integer, dimension(3) :: iGrain3 integer :: iGrain,iFace,i,j associate(prm => param(instance)) !-------------------------------------------------------------------------------------------------- ! compute the deformation gradient of individual grains due to relaxations F = 0.0_pReal do iGrain = 1,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,instance,of) ! get the relaxation vectors for each interface from global relaxation vector array nVect = interfaceNormal(intFace,instance,of) forall (i=1:3,j=1:3) & F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! calculating deformation relaxations due to interface relaxation enddo F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF ! resulting relaxed deformation gradient #ifdef DEBUG if (debugHomog%extensive) then print'(a,i3)',' Deformation gradient of grain: ',iGrain do i = 1,3 print'(1x,3(e15.8,1x))',(F(i,j,iGrain), j = 1,3) enddo print*,' ' flush(IO_STDOUT) endif #endif enddo end associate end subroutine mech_RGC_partitionDeformation !-------------------------------------------------------------------------------------------------- !> @brief update the internal state of the homogenization scheme and tell whether "done" and ! "happy" with result !-------------------------------------------------------------------------------------------------- module procedure mech_RGC_updateState integer, dimension(4) :: intFaceN,intFaceP,faceID integer, dimension(3) :: nGDim,iGr3N,iGr3P integer :: instance,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ipert,iGrain,nGrain, of real(pReal), dimension(3,3,size(P,3)) :: R,pF,pR,D,pD real(pReal), dimension(3,size(P,3)) :: NN,devNull real(pReal), dimension(3) :: normP,normN,mornP,mornN real(pReal) :: residMax,stresMax logical :: error real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix,rmatrix real(pReal), dimension(:), allocatable :: resid,relax,p_relax,p_resid,drelax #ifdef DEBUG integer, dimension(3) :: stresLoc integer, dimension(2) :: residLoc #endif zeroTimeStep: if(dEq0(dt)) then mech_RGC_updateState = .true. ! pretend everything is fine and return return endif zeroTimeStep instance = homogenization_typeInstance(material_homogenizationAt(el)) of = material_homogenizationMemberAt(ip,el) associate(stt => state(instance), st0 => state0(instance), dst => dependentState(instance), prm => param(instance)) !-------------------------------------------------------------------------------------------------- ! get the dimension of the cluster (grains and interfaces) nGDim = prm%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(:,of) drelax = stt%relaxationVector(:,of) - st0%relaxationVector(:,of) #ifdef DEBUG if (debugHomog%extensive) then print*, 'Obtained state: ' do i = 1,size(stt%relaxationVector(:,of)) print'(1x,2(e15.8,1x))', stt%relaxationVector(i,of) enddo print*,' ' endif #endif !-------------------------------------------------------------------------------------------------- ! computing interface mismatch and stress penalty tensor for all interfaces of all grains call stressPenalty(R,NN,avgF,F,ip,el,instance,of) !-------------------------------------------------------------------------------------------------- ! calculating volume discrepancy and stress penalty related to overall volume discrepancy call volumePenalty(D,dst%volumeDiscrepancy(of),avgF,F,nGrain,instance,of) !------------------------------------------------------------------------------------------------ ! computing the residual stress from the balance of traction at all (interior) interfaces do iNum = 1,nIntFaceTot faceID = interface1to4(iNum,param(instance)%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(instance)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) intFaceN = getInterface(2*faceID(1),iGr3N) normN = interfaceNormal(intFaceN,instance,of) !-------------------------------------------------------------------------------------------------- ! identify the right/up/front grain (+|P) iGr3P = iGr3N iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate system (3-dimensional index) iGrP = grain3to1(iGr3P,param(instance)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) intFaceP = getInterface(2*faceID(1)-1,iGr3P) normP = interfaceNormal(intFaceP,instance,of) !-------------------------------------------------------------------------------------------------- ! compute the residual of traction at the interface (in local system, 4-dimensional index) do i = 1,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 enddo enddo #ifdef DEBUG if (debugHomog%extensive) then print'(a,i3)',' Traction at interface: ',iNum print'(1x,3(e15.8,1x))',(tract(iNum,j), j = 1,3) print*,' ' endif #endif enddo !-------------------------------------------------------------------------------------------------- ! 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 #ifdef DEBUG if (debugHomog%extensive .and. prm%of_debug == of) then stresLoc = maxloc(abs(P)) residLoc = maxloc(abs(tract)) print'(a,i2,1x,i4)',' RGC residual check ... ',ip,el print'(a,e15.8,a,i3,a,i2,i2)', ' Max stress: ',stresMax, & '@ grain ',stresLoc(3),' in component ',stresLoc(1),stresLoc(2) print'(a,e15.8,a,i3,a,i2)',' Max residual: ',residMax, & ' @ iface ',residLoc(1),' in direction ',residLoc(2) flush(IO_STDOUT) endif #endif mech_RGC_updateState = .false. !-------------------------------------------------------------------------------------------------- ! If convergence reached => done and happy if (residMax < num%rtol*stresMax .or. residMax < num%atol) then mech_RGC_updateState = .true. #ifdef DEBUG if (debugHomog%extensive .and. prm%of_debug == of) & print*, '... done and happy'; flush(IO_STDOUT) #endif !-------------------------------------------------------------------------------------------------- ! compute/update the state for postResult, i.e., all energy densities computed by time-integration do iGrain = 1,product(prm%N_constituents) do i = 1,3;do j = 1,3 stt%work(of) = stt%work(of) & + P(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal) stt%penaltyEnergy(of) = stt%penaltyEnergy(of) & + R(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal) enddo; enddo enddo dst%mismatch(1:3,of) = sum(NN,2)/real(nGrain,pReal) dst%relaxationRate_avg(of) = sum(abs(drelax))/dt/real(3*nIntFaceTot,pReal) dst%relaxationRate_max(of) = maxval(abs(drelax))/dt #ifdef DEBUG if (debugHomog%extensive .and. prm%of_debug == of) then print'(a,e15.8)', ' Constitutive work: ',stt%work(of) print'(a,3(1x,e15.8))', ' Magnitude mismatch: ',dst%mismatch(1,of), & dst%mismatch(2,of), & dst%mismatch(3,of) print'(a,e15.8)', ' Penalty energy: ', stt%penaltyEnergy(of) print'(a,e15.8,/)', ' Volume discrepancy: ', dst%volumeDiscrepancy(of) print'(a,e15.8)', ' Maximum relaxation rate: ', dst%relaxationRate_max(of) print'(a,e15.8,/)', ' Average relaxation rate: ', dst%relaxationRate_avg(of) flush(IO_STDOUT) endif #endif 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 mech_RGC_updateState = [.true.,.false.] ! with direct cut-back #ifdef DEBUG if (debugHomog%extensive .and. prm%of_debug == of) & print'(a,/)', ' ... broken'; flush(IO_STDOUT) #endif return else ! proceed with computing the Jacobian and state update #ifdef DEBUG if (debugHomog%extensive .and. prm%of_debug == of) & print'(a,/)', ' ... not yet done'; flush(IO_STDOUT) #endif endif !--------------------------------------------------------------------------------------------------- ! construct the global Jacobian matrix for updating the global relaxation vector array when ! convergence is not yet reached ... !-------------------------------------------------------------------------------------------------- ! ... of the constitutive stress tangent, assembled from dPdF or material constitutive model "smatrix" allocate(smatrix(3*nIntFaceTot,3*nIntFaceTot), source=0.0_pReal) do iNum = 1,nIntFaceTot faceID = interface1to4(iNum,param(instance)%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(instance)%N_constituents) ! translate into global grain ID intFaceN = getInterface(2*faceID(1),iGr3N) ! identifying the connecting interface in local coordinate system normN = interfaceNormal(intFaceN,instance,of) do iFace = 1,6 intFaceN = getInterface(iFace,iGr3N) ! identifying all interfaces that influence relaxation of the above interface mornN = interfaceNormal(intFaceN,instance,of) iMun = interface4to1(intFaceN,param(instance)%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) enddo;enddo;enddo;enddo ! projecting the material tangent dPdF into the interface ! to obtain the Jacobian matrix contribution of dPdF endif enddo !-------------------------------------------------------------------------------------------------- ! identify the right/up/front grain (+|P) iGr3P = iGr3N iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate sytem iGrP = grain3to1(iGr3P,param(instance)%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,instance,of) do iFace = 1,6 intFaceP = getInterface(iFace,iGr3P) ! identifying all interfaces that influence relaxation of the above interface mornP = interfaceNormal(intFaceP,instance,of) iMun = interface4to1(intFaceP,param(instance)%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) enddo;enddo;enddo;enddo endif enddo enddo #ifdef DEBUG if (debugHomog%extensive) then print*, 'Jacobian matrix of stress' do i = 1,3*nIntFaceTot print'(1x,100(e11.4,1x))',(smatrix(i,j), j = 1,3*nIntFaceTot) enddo print*,' ' flush(IO_STDOUT) endif #endif !-------------------------------------------------------------------------------------------------- ! ... of the stress penalty tangent (mismatch penalty and volume penalty, computed using numerical ! perturbation method) "pmatrix" allocate(pmatrix(3*nIntFaceTot,3*nIntFaceTot), source=0.0_pReal) allocate(p_relax(3*nIntFaceTot), source=0.0_pReal) allocate(p_resid(3*nIntFaceTot), source=0.0_pReal) do ipert = 1,3*nIntFaceTot p_relax = relax p_relax(ipert) = relax(ipert) + num%pPert ! perturb the relaxation vector stt%relaxationVector(:,of) = p_relax call grainDeformation(pF,avgF,instance,of) ! rain deformation from perturbed state call stressPenalty(pR,DevNull, avgF,pF,ip,el,instance,of) ! stress penalty due to interface mismatch from perturbed state call volumePenalty(pD,devNull(1,1), avgF,pF,nGrain,instance,of) ! stress penalty due to volume discrepancy from perturbed state !-------------------------------------------------------------------------------------------------- ! computing the global stress residual array from the perturbed state p_resid = 0.0_pReal do iNum = 1,nIntFaceTot faceID = interface1to4(iNum,param(instance)%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(instance)%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,instance,of) !-------------------------------------------------------------------------------------------------- ! 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(instance)%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,instance,of) !-------------------------------------------------------------------------------------------------- ! 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) enddo; enddo enddo pmatrix(:,ipert) = p_resid/num%pPert enddo #ifdef DEBUG if (debugHomog%extensive) then print*, 'Jacobian matrix of penalty' do i = 1,3*nIntFaceTot print'(1x,100(e11.4,1x))',(pmatrix(i,j), j = 1,3*nIntFaceTot) enddo print*,' ' flush(IO_STDOUT) endif #endif !-------------------------------------------------------------------------------------------------- ! ... 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 enddo #ifdef DEBUG if (debugHomog%extensive) then print*, 'Jacobian matrix of penalty' do i = 1,3*nIntFaceTot print'(1x,100(e11.4,1x))',(rmatrix(i,j), j = 1,3*nIntFaceTot) enddo print*,' ' flush(IO_STDOUT) endif #endif !-------------------------------------------------------------------------------------------------- ! The overall Jacobian matrix summarizing contributions of smatrix, pmatrix, rmatrix allocate(jmatrix(3*nIntFaceTot,3*nIntFaceTot)); jmatrix = smatrix + pmatrix + rmatrix #ifdef DEBUG if (debugHomog%extensive) then print*, 'Jacobian matrix (total)' do i = 1,3*nIntFaceTot print'(1x,100(e11.4,1x))',(jmatrix(i,j), j = 1,3*nIntFaceTot) enddo print*,' ' flush(IO_STDOUT) endif #endif !-------------------------------------------------------------------------------------------------- ! 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) #ifdef DEBUG if (debugHomog%extensive) then print*, 'Jacobian inverse' do i = 1,3*nIntFaceTot print'(1x,100(e11.4,1x))',(jnverse(i,j), j = 1,3*nIntFaceTot) enddo print*,' ' flush(IO_STDOUT) endif #endif !-------------------------------------------------------------------------------------------------- ! 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 enddo; enddo stt%relaxationVector(:,of) = 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 mech_RGC_updateState = [.true.,.false.] !$OMP CRITICAL (write2out) print'(a,i3,a,i3,a)',' RGC_updateState: ip ',ip,' | el ',el,' enforces cutback' print'(a,e15.8)',' due to large relaxation change = ',maxval(abs(drelax)) flush(IO_STDOUT) !$OMP END CRITICAL (write2out) endif #ifdef DEBUG if (debugHomog%extensive) then print*, 'Returned state: ' do i = 1,size(stt%relaxationVector(:,of)) print'(1x,2(e15.8,1x))', stt%relaxationVector(i,of) enddo print*,' ' flush(IO_STDOUT) endif #endif end associate contains !------------------------------------------------------------------------------------------------ !> @brief calculate stress-like penalty due to deformation mismatch !------------------------------------------------------------------------------------------------ subroutine stressPenalty(rPen,nMis,avgF,fDef,ip,el,instance,of) 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) :: ip,el,instance,of integer, dimension (4) :: intFace integer, dimension (3) :: iGrain3,iGNghb3,nGDim real(pReal), dimension (3,3) :: gDef,nDef real(pReal), dimension (3) :: nVect,surfCorr real(pReal), dimension (2) :: Gmoduli integer :: iGrain,iGNghb,iFace,i,j,k,l real(pReal) :: muGrain,muGNghb,nDefNorm,bgGrain,bgGNghb real(pReal), parameter :: nDefToler = 1.0e-10_pReal nGDim = param(instance)%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,instance,of) associate(prm => param(instance)) #ifdef DEBUG if (debugHomog%extensive .and. prm%of_debug == of) then print'(a,2(1x,i3))', ' Correction factor: ',ip,el print*, surfCorr endif #endif !----------------------------------------------------------------------------------------------- ! computing the mismatch and penalty stress tensor of all grains grainLoop: do iGrain = 1,product(prm%N_constituents) Gmoduli = equivalentModuli(iGrain,ip,el) muGrain = Gmoduli(1) ! collecting the equivalent shear modulus of grain bgGrain = Gmoduli(2) ! and the lengthh of Burgers vector iGrain3 = grain1to3(iGrain,prm%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,instance,of) iGNghb3 = iGrain3 ! identify the neighboring grain across the interface iGNghb3(abs(intFace(1))) = iGNghb3(abs(intFace(1))) & + int(real(intFace(1),pReal)/real(abs(intFace(1)),pReal)) where(iGNghb3 < 1) iGNghb3 = nGDim where(iGNghb3 >nGDim) iGNghb3 = 1 iGNghb = grain3to1(iGNghb3,prm%N_constituents) ! get the ID of the neighboring grain Gmoduli = equivalentModuli(iGNghb,ip,el) ! collect the shear modulus and Burgers vector of the neighbor muGNghb = Gmoduli(1) bgGNghb = Gmoduli(2) gDef = 0.5_pReal*(fDef(1:3,1:3,iGNghb) - fDef(1:3,1:3,iGrain)) ! difference/jump in deformation gradeint across the neighbor !------------------------------------------------------------------------------------------- ! compute the mismatch tensor of all interfaces nDefNorm = 0.0_pReal nDef = 0.0_pReal do i = 1,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 enddo; enddo nDefNorm = nDefNorm + nDef(i,j)**2.0_pReal ! compute the norm of the mismatch tensor enddo; enddo nDefNorm = max(nDefToler,sqrt(nDefNorm)) ! approximation to zero mismatch if mismatch is zero (singularity) nMis(abs(intFace(1)),iGrain) = nMis(abs(intFace(1)),iGrain) + nDefNorm ! total amount of mismatch experienced by the grain (at all six interfaces) #ifdef DEBUG if (debugHomog%extensive .and. prm%of_debug == of) then print'(a,i2,a,i3)',' Mismatch to face: ',intFace(1),' neighbor grain: ',iGNghb print*, transpose(nDef) print'(a,e11.4)', ' with magnitude: ',nDefNorm endif #endif !------------------------------------------------------------------------------------------- ! 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*bgGrain + muGNghb*bgGNghb)*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) enddo; enddo;enddo; enddo enddo interfaceLoop #ifdef DEBUG if (debugHomog%extensive .and. prm%of_debug == of) then print'(a,i2)', ' Penalty of grain: ',iGrain print*, transpose(rPen(1:3,1:3,iGrain)) endif #endif enddo grainLoop end associate end subroutine stressPenalty !------------------------------------------------------------------------------------------------ !> @brief calculate stress-like penalty due to volume discrepancy !------------------------------------------------------------------------------------------------ subroutine volumePenalty(vPen,vDiscrep,fAvg,fDef,nGrain,instance,of) 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, & instance, & of 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 enddo !---------------------------------------------------------------------------------------------- ! calculate the stress and penalty due to volume discrepancy vPen = 0.0_pReal do i = 1,nGrain vPen(:,:,i) = -1.0_pReal/real(nGrain,pReal)*num%volDiscrMod*num%volDiscrPow/num%maxVolDiscr* & sign((abs(vDiscrep)/num%maxVolDiscr)**(num%volDiscrPow - 1.0),vDiscrep)* & gVol(i)*transpose(math_inv33(fDef(:,:,i))) #ifdef DEBUG if (debugHomog%extensive .and. param(instance)%of_debug == of) then print'(a,i2)',' Volume penalty of grain: ',i print*, transpose(vPen(:,:,i)) endif #endif enddo end subroutine volumePenalty !-------------------------------------------------------------------------------------------------- !> @brief compute the correction factor accouted for surface evolution (area change) due to ! deformation !-------------------------------------------------------------------------------------------------- function surfaceCorrection(avgF,instance,of) real(pReal), dimension(3) :: surfaceCorrection real(pReal), dimension(3,3), intent(in) :: avgF !< average F integer, intent(in) :: & instance, & of 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],instance,of) do i = 1,3; do j = 1,3 surfaceCorrection(iBase) = surfaceCorrection(iBase) + invC(i,j)*nVect(i)*nVect(j) ! compute the component of (the inverse of) the stretch in the direction of the normal enddo; enddo surfaceCorrection(iBase) = sqrt(surfaceCorrection(iBase))*detF ! get the surface correction factor (area contraction/enlargement) enddo end function surfaceCorrection !-------------------------------------------------------------------------------------------------- !> @brief compute the equivalent shear and bulk moduli from the elasticity tensor !-------------------------------------------------------------------------------------------------- function equivalentModuli(grainID,ip,el) real(pReal), dimension(2) :: equivalentModuli integer, intent(in) :: & grainID,& ip, & !< integration point number el !< element number real(pReal), dimension(6,6) :: elasTens real(pReal) :: & cEquiv_11, & cEquiv_12, & cEquiv_44 elasTens = constitutive_homogenizedC(grainID,ip,el) !---------------------------------------------------------------------------------------------- ! compute the equivalent shear modulus after Turterltaub and Suiker, JMPS (2005) cEquiv_11 = (elasTens(1,1) + elasTens(2,2) + elasTens(3,3))/3.0_pReal cEquiv_12 = (elasTens(1,2) + elasTens(2,3) + elasTens(3,1) + & elasTens(1,3) + elasTens(2,1) + elasTens(3,2))/6.0_pReal cEquiv_44 = (elasTens(4,4) + elasTens(5,5) + elasTens(6,6))/3.0_pReal equivalentModuli(1) = 0.2_pReal*(cEquiv_11 - cEquiv_12) + 0.6_pReal*cEquiv_44 !---------------------------------------------------------------------------------------------- ! obtain the length of Burgers vector (could be model dependend) equivalentModuli(2) = 2.5e-10_pReal end function equivalentModuli !-------------------------------------------------------------------------------------------------- !> @brief calculating the grain deformation gradient (the same with ! homogenization_RGC_partitionDeformation, but used only for perturbation scheme) !-------------------------------------------------------------------------------------------------- subroutine grainDeformation(F, avgF, instance, of) real(pReal), dimension(:,:,:), intent(out) :: F !< partitioned F per grain real(pReal), dimension(:,:), intent(in) :: avgF !< averaged F integer, intent(in) :: & instance, & of 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(instance)) 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,instance,of) nVect = interfaceNormal(intFace,instance,of) forall (i=1:3,j=1:3) & F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations enddo F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF ! relaxed deformation gradient enddo end associate end subroutine grainDeformation end procedure mech_RGC_updateState !-------------------------------------------------------------------------------------------------- !> @brief derive average stress and stiffness from constituent quantities !-------------------------------------------------------------------------------------------------- module subroutine mech_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses integer, intent(in) :: instance avgP = sum(P,3) /real(product(param(instance)%N_constituents),pReal) dAvgPdAvgF = sum(dPdF,5)/real(product(param(instance)%N_constituents),pReal) end subroutine mech_RGC_averageStressAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- module subroutine mech_RGC_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group integer :: o associate(stt => state(instance), dst => dependentState(instance), prm => param(instance)) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) case('W') call results_writeDataset(group,stt%work,trim(prm%output(o)), & 'work density','J/m³') case('M') call results_writeDataset(group,dst%mismatch,trim(prm%output(o)), & 'average mismatch tensor','1') case('R') call results_writeDataset(group,stt%penaltyEnergy,trim(prm%output(o)), & 'mismatch penalty density','J/m³') case('Delta_V') call results_writeDataset(group,dst%volumeDiscrepancy,trim(prm%output(o)), & 'volume discrepancy','m³') case('max_a_dot') call results_writeDataset(group,dst%relaxationrate_max,trim(prm%output(o)), & 'maximum relaxation rate','m/s') case('avg_a_dot') call results_writeDataset(group,dst%relaxationrate_avg,trim(prm%output(o)), & 'average relaxation rate','m/s') end select enddo outputsLoop end associate end subroutine mech_RGC_results !-------------------------------------------------------------------------------------------------- !> @brief collect relaxation vectors of an interface !-------------------------------------------------------------------------------------------------- pure function relaxationVector(intFace,instance,of) real(pReal), dimension (3) :: relaxationVector integer, intent(in) :: instance,of integer, dimension(4), intent(in) :: intFace !< set of interface ID in 4D array (normal and position) integer :: iNum !-------------------------------------------------------------------------------------------------- ! collect the interface relaxation vector from the global state array iNum = interface4to1(intFace,param(instance)%N_constituents) ! identify the position of the interface in global state array if (iNum > 0) then relaxationVector = state(instance)%relaxationVector((3*iNum-2):(3*iNum),of) else relaxationVector = 0.0_pReal endif end function relaxationVector !-------------------------------------------------------------------------------------------------- !> @brief identify the normal of an interface !-------------------------------------------------------------------------------------------------- pure function interfaceNormal(intFace,instance,of) real(pReal), dimension(3) :: interfaceNormal integer, dimension(4), intent(in) :: intFace !< interface ID in 4D array (normal and position) integer, intent(in) :: & instance, & of integer :: nPos !-------------------------------------------------------------------------------------------------- ! get the normal of the interface, identified from the value of intFace(1) interfaceNormal = 0.0_pReal nPos = abs(intFace(1)) ! identify the position of the interface in global state array interfaceNormal(nPos) = real(intFace(1)/abs(intFace(1)),pReal) ! get the normal vector w.r.t. cluster axis interfaceNormal = matmul(dependentState(instance)%orientation(1:3,1:3,of),interfaceNormal) ! map the normal vector into sample coordinate system (basis) end function interfaceNormal !-------------------------------------------------------------------------------------------------- !> @brief collect six faces of a grain in 4D (normal and position) !-------------------------------------------------------------------------------------------------- pure function getInterface(iFace,iGrain3) integer, dimension(4) :: getInterface integer, dimension(3), intent(in) :: iGrain3 !< grain ID in 3D array integer, intent(in) :: iFace !< face index (1..6) mapped like (-e1,-e2,-e3,+e1,+e2,+e3) or iDir = (-1,-2,-3,1,2,3) integer :: iDir !< direction of interface normal iDir = (int(real(iFace-1,pReal)/2.0_pReal)+1)*(-1)**iFace getInterface(1) = iDir !-------------------------------------------------------------------------------------------------- ! identify the interface position by the direction of its normal getInterface(2:4) = iGrain3 if (iDir < 0) getInterface(1-iDir) = getInterface(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) endif 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 endif 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 endif 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 endif end function interface1to4 end submodule homogenization_mech_RGC