diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index 2f6959d05..af05ab679 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -86,6 +86,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt cycleCounter, & theInc, & theTime, & + theDelta, & FEsolving_execElem, & FEsolving_execIP use math, only: math_init, & @@ -98,7 +99,9 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt use mesh, only: mesh_init, & mesh_FEasCP, & mesh_NcpElems, & - mesh_maxNips + mesh_maxNips, & + mesh_element, & + FE_Nips use lattice, only: lattice_init use material, only: material_init, & homogenization_maxNgrains, & @@ -107,8 +110,6 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt constitutive_state0,constitutive_state use crystallite, only: crystallite_init, & crystallite_F0, & - crystallite_rexParm, & - crystallite_critVal, & crystallite_partionedF, & crystallite_Fp0, & crystallite_Fp, & @@ -163,7 +164,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt k, & l, & m, & - n + n, & + e logical updateJaco ! flag indicating if JAcobian has to be updated !*** global variables ***! diff --git a/code/constitutive.f90 b/code/constitutive.f90 index f6f10f924..b838c6910 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -310,7 +310,7 @@ return endfunction -subroutine constitutive_microstructure(Temperature,Fe,Fp,ipc,ip,el) +subroutine constitutive_microstructure(Temperature,Tstar_v,Fe,Fp,ipc,ip,el) !********************************************************************* !* This function calculates from state needed variables * !* INPUT: * @@ -579,4 +579,4 @@ return endfunction -END MODULE +END MODULE \ No newline at end of file diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index 2592e12f2..709e2c2a6 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -73,10 +73,10 @@ real(pReal), dimension(:,:), allocatable :: constitutive_nonlocal_ constitutive_nonlocal_lambda0PerSlipSystem, & ! mean free path prefactor for each slip system and instance constitutive_nonlocal_burgersPerSlipFamily, & ! absolute length of burgers vector [m] for each family and instance constitutive_nonlocal_burgersPerSlipSystem, & ! absolute length of burgers vector [m] for each slip system and instance - constitutive_nonlocal_dLowerEdgePerSlipFamily, & ! minimum stable edge dipole height for each family and instance - constitutive_nonlocal_dLowerEdgePerSlipSystem, & ! minimum stable edge dipole height for each slip system and instance - constitutive_nonlocal_dLowerScrewPerSlipFamily, & ! minimum stable screw dipole height for each family and instance - constitutive_nonlocal_dLowerScrewPerSlipSystem, & ! minimum stable screw dipole height for each slip system and instance + constitutive_nonlocal_dLowerEdgePerSlipFamily, & ! minimum stable edge dipole height for each family and instance + constitutive_nonlocal_dLowerEdgePerSlipSystem, & ! minimum stable edge dipole height for each slip system and instance + constitutive_nonlocal_dLowerScrewPerSlipFamily, & ! minimum stable screw dipole height for each family and instance + constitutive_nonlocal_dLowerScrewPerSlipSystem, & ! minimum stable screw dipole height for each slip system and instance constitutive_nonlocal_interactionSlipSlip ! coefficients for slip-slip interaction for each interaction type and instance real(pReal), dimension(:,:,:,:,:), allocatable :: constitutive_nonlocal_v, & ! dislocation velocity constitutive_nonlocal_rhoDotFlux ! dislocation convection term @@ -1230,9 +1230,17 @@ integer(pInt) myInstance, & ! current neighboring_ip, & ! integration point of my neighbor c, & ! character of dislocation n, & ! index of my current neighbor + opposite_n, & ! index of my opposite neighbor + opposite_ip, & ! ip of my opposite neighbor + opposite_el, & ! element index of my opposite neighbor t, & ! type of dislocation s, & ! index of my current slip system sLattice ! index of my current slip system according to lattice order +real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),8) :: & + rhoSgl, & ! current single dislocation densities (positive/negative screw and edge without dipoles) + previousRhoSgl, & ! previous single dislocation densities (positive/negative screw and edge without dipoles) + totalRhoDotSgl, & ! total rate of change of single dislocation densities + thisRhoDotSgl ! rate of change of single dislocation densities for this mechanism real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: & fluxdensity, & ! flux density at central material point neighboring_fluxdensity, &! flux density at neighbroing material point diff --git a/code/crystallite.f90 b/code/crystallite.f90 index d3c721cf3..3e4ff6d9e 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -291,8 +291,7 @@ subroutine crystallite_init(Temperature) write(6,'(a35,x,7(i5,x))') 'crystallite_subTstar0_v: ', shape(crystallite_subTstar0_v) write(6,'(a35,x,7(i5,x))') 'crystallite_dPdF: ', shape(crystallite_dPdF) write(6,'(a35,x,7(i5,x))') 'crystallite_fallbackdPdF: ', shape(crystallite_fallbackdPdF) - write(6,'(a35,x,7(i5,x))') 'crystallite_R: ', shape(crystallite_R) - write(6,'(a35,x,7(i5,x))') 'crystallite_eulerangles: ', shape(crystallite_eulerangles) + write(6,'(a35,x,7(i5,x))') 'crystallite_orientation: ', shape(crystallite_orientation) write(6,'(a35,x,7(i5,x))') 'crystallite_misorientation: ', shape(crystallite_misorientation) write(6,'(a35,x,7(i5,x))') 'crystallite_dt: ', shape(crystallite_dt) write(6,'(a35,x,7(i5,x))') 'crystallite_subdt: ', shape(crystallite_subdt) @@ -300,11 +299,12 @@ subroutine crystallite_init(Temperature) write(6,'(a35,x,7(i5,x))') 'crystallite_subStep: ', shape(crystallite_subStep) write(6,'(a35,x,7(i5,x))') 'crystallite_localConstitution: ', shape(crystallite_localConstitution) write(6,'(a35,x,7(i5,x))') 'crystallite_requested: ', shape(crystallite_requested) - write(6,'(a35,x,7(i5,x))') 'crystallite_onTrack: ', shape(crystallite_onTrack) + write(6,'(a35,x,7(i5,x))') 'crystallite_todo: ', shape(crystallite_todo) write(6,'(a35,x,7(i5,x))') 'crystallite_converged: ', shape(crystallite_converged) write(6,'(a35,x,7(i5,x))') 'crystallite_stateConverged: ', shape(crystallite_stateConverged) write(6,'(a35,x,7(i5,x))') 'crystallite_temperatureConverged: ', shape(crystallite_temperatureConverged) - write(6,'(a35,x,7(i5,x))') 'crystallite_todo: ', shape(crystallite_todo) + write(6,'(a35,x,7(i5,x))') 'crystallite_sizePostResults: ', shape(crystallite_sizePostResults) + write(6,'(a35,x,7(i5,x))') 'crystallite_sizePostResult: ', shape(crystallite_sizePostResult) write(6,*) write(6,*) 'Number of nonlocal grains: ',count(.not. crystallite_localConstitution) call flush(6) @@ -333,8 +333,14 @@ subroutine crystallite_stressAndItsTangent(updateJaco) pert_Fg, & pert_method, & nState, & - nCryst + nCryst, & + iJacoStiffness use debug, only: debugger, & + selectiveDebugger, & + verboseDebugger, & + debug_e, & + debug_i, & + debug_g, & debug_CrystalliteLoopDistribution, & debug_CrystalliteStateLoopDistribution, & debug_StiffnessStateLoopDistribution @@ -345,8 +351,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) math_Mandel6to33, & math_Mandel33to6, & math_I3, & - math_Plain3333to99, & - math_EulerToR + math_Plain3333to99 use FEsolving, only: FEsolving_execElem, & FEsolving_execIP, & theInc, & @@ -380,8 +385,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !*** local variables ***! real(pReal) myTemperature, & ! local copy of the temperature - myPert, & ! perturbation with correct sign - rexCritStrain ! perturbation with correct sign + myPert ! perturbation with correct sign real(pReal), dimension(3,3) :: invFp, & ! inverse of the plastic deformation gradient Fe_guess, & ! guess for elastic deformation gradient Tstar ! 2nd Piola-Kirchhoff stress tensor @@ -397,7 +401,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) myNgrains, & mySizeState, & mySizeDotState - integer(pInt), dimension(2,9) :: kl + integer(pInt), dimension(2,9,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & + kl logical onTrack, & ! flag indicating whether we are still on track temperatureConverged, & ! flag indicating if temperature converged stateConverged, & ! flag indicating if state converged @@ -407,7 +412,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) real(pReal), dimension(constitutive_maxSizeDotState) :: delta_dotState1, & ! difference between current and previous dotstate delta_dotState2 ! difference between previousDotState and previousDotState2 real(pReal) dot_prod12, & - dot_prod22 + dot_prod22, & + formerSubStep real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & storedF, & storedFp, & @@ -426,6 +432,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco) logical, dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & storedConvergenceFlag logical, dimension(3,3) :: mask + logical forceLocalStiffnessCalculation ! flag indicating that stiffness calculation is always done locally + forceLocalStiffnessCalculation = .false. + ! ------ initialize to starting condition ------ @@ -442,17 +451,17 @@ subroutine crystallite_stressAndItsTangent(updateJaco) crystallite_subStep = 0.0_pReal !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed + do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains - if (crystallite_requested(g,i,e)) then ! initialize restoration point of ... - crystallite_subTemperature0(g,i,e) = crystallite_partionedTemperature0(g,i,e) ! ...temperature - constitutive_subState0(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructure - crystallite_subFp0(:,:,g,i,e) = crystallite_partionedFp0(:,:,g,i,e) ! ...plastic def grad - crystallite_subLp0(:,:,g,i,e) = crystallite_partionedLp0(:,:,g,i,e) ! ...plastic velocity grad - crystallite_subF0(:,:,g,i,e) = crystallite_partionedF0(:,:,g,i,e) ! ...def grad - crystallite_subTstar0_v(:,g,i,e) = crystallite_partionedTstar0_v(:,g,i,e) ! ...2nd PK stress + if (crystallite_requested(g,i,e)) then ! initialize restoration point of ... + crystallite_subTemperature0(g,i,e) = crystallite_partionedTemperature0(g,i,e) ! ...temperature + constitutive_subState0(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructure + crystallite_subFp0(:,:,g,i,e) = crystallite_partionedFp0(:,:,g,i,e) ! ...plastic def grad + crystallite_subLp0(:,:,g,i,e) = crystallite_partionedLp0(:,:,g,i,e) ! ...plastic velocity grad + crystallite_subF0(:,:,g,i,e) = crystallite_partionedF0(:,:,g,i,e) ! ...def grad + crystallite_subTstar0_v(:,g,i,e) = crystallite_partionedTstar0_v(:,g,i,e) !...2nd PK stress crystallite_subFrac(g,i,e) = 0.0_pReal crystallite_subStep(g,i,e) = 1.0_pReal/subStepSizeCryst @@ -463,6 +472,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo enddo !$OMPEND PARALLEL DO + ! --+>> crystallite loop <<+-- @@ -543,7 +553,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains if (crystallite_todo(g,i,e)) then ! all undone crystallites call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & @@ -573,7 +583,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) if (crystallite_todo(g,i,e)) then ! all undone crystallites @@ -602,9 +612,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco) ! to account for substepping within _integrateStress ! results in crystallite_Fp,.._Lp !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed + do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) if (crystallite_todo(g,i,e)) then ! all undone crystallites @@ -635,10 +645,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains - if (crystallite_todo(g,i,e)) then ! all undone crystallites + if (crystallite_todo(g,i,e)) then ! all undone crystallites constitutive_previousDotState2(g,i,e)%p = constitutive_previousDotState(g,i,e)%p constitutive_previousDotState(g,i,e)%p = constitutive_dotState(g,i,e)%p - constitutive_dotState(g,i,e)%p = 0.0_pReal ! zero out dotState + constitutive_dotState(g,i,e)%p = 0.0_pReal ! zero out dotState endif enddo; enddo; enddo !$OMPEND PARALLEL DO @@ -648,7 +658,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) if (crystallite_todo(g,i,e)) then ! all undone crystallites @@ -669,9 +679,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !$OMPEND PARALLEL DO !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed + do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) if (crystallite_todo(g,i,e)) then ! all undone crystallites @@ -718,18 +728,18 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !$OMPEND CRITICAL (write2out) endif - enddo ! crystallite convergence loop + enddo ! crystallite convergence loop NiterationCrystallite = NiterationCrystallite + 1 - enddo ! cutback loop + enddo ! cutback loop ! ------ check for non-converged crystallites ------ !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed + do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains - if (.not. crystallite_converged(g,i,e)) then ! respond fully elastically (might be not required due to becoming terminally ill anyway) + if (.not. crystallite_converged(g,i,e)) then ! respond fully elastically (might be not required due to becoming terminally ill anyway) ! call IO_warning(600,e,i,g) invFp = math_inv3x3(crystallite_partionedFp0(:,:,g,i,e)) Fe_guess = math_mul33x33(crystallite_partionedF(:,:,g,i,e),invFp) @@ -851,11 +861,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo ! perturbation direction select case(pert_method) case (1) - crystallite_dPdF(:,:,:,:,g,i,e) = dPdF_perturbation(:,:,:,:,1) + crystallite_dPdF(:,:,:,:,g,i,e) = dPdF_perturbation(:,:,:,:,1) case (2) - crystallite_dPdF(:,:,:,:,g,i,e) = dPdF_perturbation(:,:,:,:,2) + crystallite_dPdF(:,:,:,:,g,i,e) = dPdF_perturbation(:,:,:,:,2) case (3) - crystallite_dPdF(:,:,:,:,g,i,e) = 0.5_pReal*(dPdF_perturbation(:,:,:,:,1)+dPdF_perturbation(:,:,:,:,2)) + crystallite_dPdF(:,:,:,:,g,i,e) = 0.5_pReal*(dPdF_perturbation(:,:,:,:,1)+dPdF_perturbation(:,:,:,:,2)) end select else ! grain did not converge crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! use (elastic) fallback @@ -866,7 +876,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo ! element loop !$OMPEND PARALLEL DO - elseif (any(.not. crystallite_localConstitution)) then ! if any nonlocal grain present, we have to do a full loop over all grains after each perturbance + elseif (any(.not. crystallite_localConstitution)) then ! if any nonlocal grain present, we have to do a full loop over all grains after each perturbance do e = FEsolving_execElem(1),FEsolving_execElem(2) myNgrains = homogenization_Ngrains(mesh_element(3,e)) diff --git a/code/homogenization_RGC.f90 b/code/homogenization_RGC.f90 index 42ee91a48..8b22e5da0 100644 --- a/code/homogenization_RGC.f90 +++ b/code/homogenization_RGC.f90 @@ -20,14 +20,14 @@ MODULE homogenization_RGC integer(pInt), dimension(:), allocatable :: homogenization_RGC_sizeState, & homogenization_RGC_sizePostResults - integer(pInt), dimension(:,:), allocatable,target :: homogenization_RGC_sizePostResult + integer(pInt), dimension(:,:), allocatable,target :: homogenization_RGC_sizePostResult integer(pInt), dimension(:,:), allocatable :: homogenization_RGC_Ngrains real(pReal), dimension(:,:), allocatable :: homogenization_RGC_dAlpha, & homogenization_RGC_angles real(pReal), dimension(:,:,:,:), allocatable :: homogenization_RGC_orientation real(pReal), dimension(:), allocatable :: homogenization_RGC_xiAlpha, & homogenization_RGC_ciAlpha - character(len=64), dimension(:,:), allocatable,target :: homogenization_RGC_output ! name of each post result output + character(len=64), dimension(:,:), allocatable,target :: homogenization_RGC_output ! name of each post result output CONTAINS !**************************************** @@ -48,17 +48,17 @@ subroutine homogenization_RGC_init(& ) use prec, only: pInt, pReal - use math, only: math_Mandel3333to66, math_Voigt66to3333 - use mesh, only: mesh_maxNips,mesh_NcpElems + use math, only: math_Mandel3333to66, math_Voigt66to3333,math_I3,math_sampleRandomOri,math_EulerToR,inRad + use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips use IO use material integer(pInt), intent(in) :: file integer(pInt), parameter :: maxNchunks = 4 integer(pInt), dimension(1+2*maxNchunks) :: positions - integer(pInt) section, maxNinstance, i,j,k,l, output, mySize + integer(pInt) section, maxNinstance, i,j,k,l,e, output, mySize, myInstance character(len=64) tag character(len=1024) line - + write(6,*) write(6,'(a20,a20,a12)') '<<<+- homogenization',homogenization_RGC_label,' init -+>>>' write(6,*) '$Id$' @@ -69,14 +69,14 @@ subroutine homogenization_RGC_init(& allocate(homogenization_RGC_sizeState(maxNinstance)); homogenization_RGC_sizeState = 0_pInt allocate(homogenization_RGC_sizePostResults(maxNinstance)); homogenization_RGC_sizePostResults = 0_pInt - allocate(homogenization_RGC_sizePostResult(maxval(homogenization_Noutput), & - maxNinstance)); homogenization_RGC_sizePostResult = 0_pInt allocate(homogenization_RGC_Ngrains(3,maxNinstance)); homogenization_RGC_Ngrains = 0_pInt allocate(homogenization_RGC_ciAlpha(maxNinstance)); homogenization_RGC_ciAlpha = 0.0_pReal allocate(homogenization_RGC_xiAlpha(maxNinstance)); homogenization_RGC_xiAlpha = 0.0_pReal allocate(homogenization_RGC_dAlpha(3,maxNinstance)); homogenization_RGC_dAlpha = 0.0_pReal allocate(homogenization_RGC_angles(3,maxNinstance)); homogenization_RGC_angles = 400.0_pReal allocate(homogenization_RGC_output(maxval(homogenization_Noutput),maxNinstance)); homogenization_RGC_output = '' + allocate(homogenization_RGC_sizePostResult(maxval(homogenization_Noutput),maxNinstance)) + homogenization_RGC_sizePostResult = 0_pInt allocate(homogenization_RGC_orientation(3,3,mesh_maxNips,mesh_NcpElems)) forall (i = 1:mesh_maxNips,e = 1:mesh_NcpElems) homogenization_RGC_orientation(:,:,i,e) = math_I3 @@ -148,10 +148,6 @@ subroutine homogenization_RGC_init(& enddo 100 do i = 1,maxNinstance ! sanity checks - write(6,*) - write(6,*) '<<<+- homogenization_RGC init -+>>>' - write(6,*) '$Id$' - write(6,*) write(6,'(a15,x,i4)') 'instance: ', i write(6,*) write(6,'(a25,3(x,i8))') 'cluster size: ',(homogenization_RGC_Ngrains(j,i),j=1,3) @@ -173,7 +169,7 @@ subroutine homogenization_RGC_init(& case('volumediscrepancy') mySize = 1 case default - mySize = 0 + mySize = 0 case('averagerelaxrate') mySize = 1 case('maximumrelaxrate') @@ -187,6 +183,8 @@ subroutine homogenization_RGC_init(& endif enddo + + homogenization_RGC_sizeState(i) & = 3*(homogenization_RGC_Ngrains(1,i)-1)*homogenization_RGC_Ngrains(2,i)*homogenization_RGC_Ngrains(3,i) & + 3*homogenization_RGC_Ngrains(1,i)*(homogenization_RGC_Ngrains(2,i)-1)*homogenization_RGC_Ngrains(3,i) & diff --git a/code/numerics.f90 b/code/numerics.f90 index 437a27a83..ced456a35 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -253,7 +253,7 @@ subroutine numerics_init() write(6,'(a24,x,e8.1)') 'rMax_RGC: ',relMax_RGC write(6,'(a24,x,e8.1)') 'perturbPenalty_RGC: ',pPert_RGC write(6,'(a24,x,e8.1)') 'relevantMismatch_RGC: ',xSmoo_RGC - write(6,'(a24,x,e8.1)') 'viscosityrate_RGC: ',ratePower_RGC + write(6,'(a24,x,e8.1)') 'viscosityrate_RGC: ',viscPower_RGC write(6,'(a24,x,e8.1)') 'viscositymodulus_RGC: ',viscModus_RGC write(6,'(a24,x,e8.1)') 'maxrelaxation_RGC: ',maxdRelax_RGC write(6,'(a24,x,e8.1)') 'maxVolDiscrepancy_RGC:',maxVolDiscr_RGC