unifying code style

src/source_damage_anisoBrittle.f90

@@ -9,8 +9,8 @@ module source_damage_anisoBrittle
   use debug
   use IO
   use math
-  use material
   use discretization
+  use material
   use config
   use lattice
   use results
@@ -58,7 +58,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine source_damage_anisoBrittle_init
 
-  integer :: Ninstance,instance,source,sourceOffset,NofMyPhase,p
+  integer :: Ninstance,source,sourceOffset,NofMyPhase,p
   character(len=pStringLen) :: &
     extmsg = ''
 
@@ -72,7 +72,6 @@ subroutine source_damage_anisoBrittle_init
   allocate(source_damage_anisoBrittle_instance(size(config_phase)), source=0)
   allocate(param(Ninstance))
 
-
   do p = 1, size(config_phase)
     source_damage_anisoBrittle_instance(p) = count(phase_source(:,1:p) == SOURCE_DAMAGE_ANISOBRITTLE_ID)
     do source = 1, phase_Nsources(p)
@@ -84,21 +83,24 @@ subroutine source_damage_anisoBrittle_init
     associate(prm => param(source_damage_anisoBrittle_instance(p)), &
               config => config_phase(p))
 
+    prm%Ncleavage = config%getInts('ncleavage',defaultVal=emptyIntArray)
+    prm%totalNcleavage = sum(prm%Ncleavage)
+
     prm%aTol      = config%getFloat('anisobrittle_atol',defaultVal = 1.0e-3_pReal)
     prm%N         = config%getFloat('anisobrittle_ratesensitivity')
     prm%sdot_0    = config%getFloat('anisobrittle_sdot0')
-    prm%Ncleavage = config%getInts('ncleavage',defaultVal=emptyIntArray)
+
-    prm%totalNcleavage = sum(prm%Ncleavage)
     prm%critDisp  = config%getFloats('anisobrittle_criticaldisplacement',requiredSize=size(prm%Ncleavage))
     prm%critLoad  = config%getFloats('anisobrittle_criticalload',        requiredSize=size(prm%Ncleavage))
 
-    prm%cleavage_systems  = lattice_SchmidMatrix_cleavage (prm%Ncleavage,config%getString('lattice_structure'),&
+    prm%cleavage_systems = lattice_SchmidMatrix_cleavage(prm%Ncleavage,config%getString('lattice_structure'),&
                                                          config%getFloat('c/a',defaultVal=0.0_pReal))
 
     ! expand: family => system
     prm%critDisp = math_expand(prm%critDisp, prm%Ncleavage)
     prm%critLoad = math_expand(prm%critLoad, prm%Ncleavage)
 
+    ! sanity checks
     if (prm%aTol         <  0.0_pReal)  extmsg = trim(extmsg)//' anisobrittle_atol'
     if (prm%N            <= 0.0_pReal)  extmsg = trim(extmsg)//' anisobrittle_ratesensitivity'
     if (prm%sdot_0       <= 0.0_pReal)  extmsg = trim(extmsg)//' anisobrittle_sdot0'
@@ -107,19 +109,17 @@ subroutine source_damage_anisoBrittle_init
 
 !--------------------------------------------------------------------------------------------------
 !  exit if any parameter is out of range
-    if (extmsg /= '') &
+    if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//')')
-      call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//')')
 
 !--------------------------------------------------------------------------------------------------
 !  output pararameters
     prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
 
-    NofMyPhase=count(material_phaseAt==p) * discretization_nIP
+    NofMyPhase = count(material_phaseAt==p) * discretization_nIP
-    instance = source_damage_anisoBrittle_instance(p)
     sourceOffset = source_damage_anisoBrittle_offset(p)
 
     call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0)
-    sourceState(p)%p(sourceOffset)%aTolState=param(instance)%aTol
+    sourceState(p)%p(sourceOffset)%aTolState=prm%aTol
 
     end associate
   enddo
@@ -138,44 +138,43 @@ subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el)
     el                                                                                              !< element
   real(pReal),  intent(in), dimension(3,3) :: &
     S
+
   integer :: &
     phase, &
     constituent, &
-    instance, &
     sourceOffset, &
     damageOffset, &
     homog, &
-    index
+    i
   real(pReal) :: &
     traction_d, traction_t, traction_n, traction_crit
 
   phase = material_phaseAt(ipc,el)
   constituent = material_phasememberAt(ipc,ip,el)
-  instance = source_damage_anisoBrittle_instance(phase)
   sourceOffset = source_damage_anisoBrittle_offset(phase)
   homog = material_homogenizationAt(el)
   damageOffset = damageMapping(homog)%p(ip,el)
 
+  associate(prm => param(source_damage_anisoBrittle_instance(phase)))
   sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal
+  do i = 1, prm%totalNcleavage
 
-  do index = 1, param(instance)%totalNcleavage
+    traction_d    = math_mul33xx33(S,prm%cleavage_systems(1:3,1:3,1,i))
+    traction_t    = math_mul33xx33(S,prm%cleavage_systems(1:3,1:3,2,i))
+    traction_n    = math_mul33xx33(S,prm%cleavage_systems(1:3,1:3,3,i))
 
-      traction_d    = math_mul33xx33(S,param(instance)%cleavage_systems(1:3,1:3,1,index))
+    traction_crit = prm%critLoad(i)*damage(homog)%p(damageOffset)**2.0_pReal
-      traction_t    = math_mul33xx33(S,param(instance)%cleavage_systems(1:3,1:3,2,index))
-      traction_n    = math_mul33xx33(S,param(instance)%cleavage_systems(1:3,1:3,3,index))
 
-      traction_crit = param(instance)%critLoad(index)* &
+    sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &
-                      damage(homog)%p(damageOffset)*damage(homog)%p(damageOffset)
+    = sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &
-
+    + prm%sdot_0* &
-      sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = &
+      ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%N + &
-        sourceState(phase)%p(sourceOffset)%dotState(1,constituent) + &
+       (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%N + &
-        param(instance)%sdot_0* &
+       (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**prm%N)/ &
-        ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**param(instance)%N + &
+      prm%critDisp(i)
-         (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**param(instance)%N + &
-         (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**param(instance)%N)/ &
-        param(instance)%critDisp(index)
 
   enddo
+  end associate
 
 end subroutine source_damage_anisoBrittle_dotState
 
@@ -193,16 +192,17 @@ subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalph
   real(pReal),  intent(out) :: &
     localphiDot, &
     dLocalphiDot_dPhi
+
   integer :: &
     sourceOffset
 
   sourceOffset = source_damage_anisoBrittle_offset(phase)
 
-  localphiDot = 1.0_pReal &
-              - sourceState(phase)%p(sourceOffset)%state(1,constituent)*phi
-
   dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent)
 
+  localphiDot = 1.0_pReal &
+              + dLocalphiDot_dPhi*phi
+
 end subroutine source_damage_anisoBrittle_getRateAndItsTangent
 
 

src/source_damage_anisoDuctile.f90

@@ -53,7 +53,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine source_damage_anisoDuctile_init
 
-  integer :: Ninstance,instance,source,sourceOffset,NofMyPhase,p
+  integer :: Ninstance,source,sourceOffset,NofMyPhase,p
   character(len=pStringLen) :: &
     extmsg = ''
 
@@ -78,12 +78,15 @@ subroutine source_damage_anisoDuctile_init
     associate(prm => param(source_damage_anisoDuctile_instance(p)), &
               config => config_phase(p))
 
+    prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
+    prm%totalNslip = sum(prm%Nslip)
+
     prm%aTol              = config%getFloat('anisoductile_atol',defaultVal = 1.0e-3_pReal)
     prm%N                 = config%getFloat('anisoductile_ratesensitivity')
-    prm%Nslip             = config%getInts('nslip',defaultVal=emptyIntArray)
+
     prm%critPlasticStrain = config%getFloats('anisoductile_criticalplasticstrain',requiredSize=size(prm%Nslip))
 
-    prm%totalNslip = sum(prm%Nslip)      ! expand: family => system
+    ! expand: family => system
     prm%critPlasticStrain = math_expand(prm%critPlasticStrain,  prm%Nslip)
 
     ! sanity checks
@@ -100,11 +103,10 @@ subroutine source_damage_anisoDuctile_init
     prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
 
     NofMyPhase=count(material_phaseAt==p) * discretization_nIP
-    instance = source_damage_anisoDuctile_instance(p)
     sourceOffset = source_damage_anisoDuctile_offset(p)
 
     call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0)
-    sourceState(p)%p(sourceOffset)%aTolState=param(instance)%aTol
+    sourceState(p)%p(sourceOffset)%aTolState=prm%aTol
 
     end associate
   enddo
@@ -121,28 +123,28 @@ subroutine source_damage_anisoDuctile_dotState(ipc, ip, el)
     ipc, &                                                                                          !< component-ID of integration point
     ip, &                                                                                           !< integration point
     el                                                                                              !< element
+
   integer :: &
     phase, &
     constituent, &
     sourceOffset, &
-    homog, damageOffset, &
+    damageOffset, &
-    instance, &
+    homog, &
     i
 
   phase = material_phaseAt(ipc,el)
   constituent = material_phasememberAt(ipc,ip,el)
-  instance = source_damage_anisoDuctile_instance(phase)
   sourceOffset = source_damage_anisoDuctile_offset(phase)
   homog = material_homogenizationAt(el)
   damageOffset = damageMapping(homog)%p(ip,el)
 
-
+  associate(prm => param(source_damage_anisoDuctile_instance(phase)))
-  do i = 1, param(instance)%totalNslip
+  do i = 1, prm%totalNslip
-      sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = &
+    sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &
-        sourceState(phase)%p(sourceOffset)%dotState(1,constituent) + &
+    = sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &
-        plasticState(phase)%slipRate(i,constituent)/ &
+    + plasticState(phase)%slipRate(i,constituent)/(damage(homog)%p(damageOffset)**prm%N)/prm%critPlasticStrain(i)
-        ((damage(homog)%p(damageOffset))**param(instance)%N)/param(instance)%critPlasticStrain(i)
   enddo
+  end associate
 
 end subroutine source_damage_anisoDuctile_dotState
 
@@ -160,16 +162,17 @@ subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalph
   real(pReal),  intent(out) :: &
     localphiDot, &
     dLocalphiDot_dPhi
+
   integer :: &
     sourceOffset
 
   sourceOffset = source_damage_anisoDuctile_offset(phase)
 
-  localphiDot = 1.0_pReal &
-              - sourceState(phase)%p(sourceOffset)%state(1,constituent) * phi
-
   dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent)
 
+  localphiDot = 1.0_pReal &
+              + dLocalphiDot_dPhi*phi
+
 end subroutine source_damage_anisoDuctile_getRateAndItsTangent
 
 

src/source_damage_isoBrittle.f90

@@ -48,7 +48,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine source_damage_isoBrittle_init
 
-  integer :: Ninstance,instance,source,sourceOffset,NofMyPhase,p
+  integer :: Ninstance,source,sourceOffset,NofMyPhase,p
   character(len=pStringLen) :: &
     extmsg = ''
 
@@ -92,17 +92,17 @@ subroutine source_damage_isoBrittle_init
     prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
 
     NofMyPhase = count(material_phaseAt==p) * discretization_nIP
-    instance = source_damage_isoBrittle_instance(p)
     sourceOffset = source_damage_isoBrittle_offset(p)
 
     call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,1)
-    sourceState(p)%p(sourceOffset)%aTolState=param(instance)%aTol
+    sourceState(p)%p(sourceOffset)%aTolState=prm%aTol
 
     end associate
   enddo
 
 end subroutine source_damage_isoBrittle_init
 
+
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates derived quantities from state
 !--------------------------------------------------------------------------------------------------
@@ -116,22 +116,25 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
     Fe
   real(pReal),  intent(in), dimension(6,6) :: &
     C
+
   integer :: &
-    phase, constituent, instance, sourceOffset
+    phase, &
+    constituent, &
+    sourceOffset
+  real(pReal), dimension(6) :: &
+    strain
   real(pReal) :: &
-    strain(6), &
     strainenergy
 
   phase = material_phaseAt(ipc,el)                                                                        !< phase ID at ipc,ip,el
   constituent = material_phasememberAt(ipc,ip,el)                                                            !< state array offset for phase ID at ipc,ip,el
   ! ToDo: capability for multiple instances of SAME source within given phase. Needs Ninstance loop from here on!
-  instance = source_damage_isoBrittle_instance(phase)                                               !< instance of damage_isoBrittle source
   sourceOffset = source_damage_isoBrittle_offset(phase)
 
-
   strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3)
 
-  strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/param(instance)%critStrainEnergy
+  associate(prm => param(source_damage_isoBrittle_instance(phase)))
+  strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/prm%critStrainEnergy
   ! ToDo: check strainenergy = 2.0_pReal*dot_product(strain,matmul(C,strain))/param(instance)%critStrainEnergy
 
   if (strainenergy > sourceState(phase)%p(sourceOffset)%subState0(1,constituent)) then
@@ -142,9 +145,11 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
       sourceState(phase)%p(sourceOffset)%subState0(1,constituent) - &
       sourceState(phase)%p(sourceOffset)%state(1,constituent)
   endif
+  end associate
 
 end subroutine source_damage_isoBrittle_deltaState
 
+
 !--------------------------------------------------------------------------------------------------
 !> @brief returns local part of nonlocal damage driving force
 !--------------------------------------------------------------------------------------------------
@@ -158,17 +163,19 @@ subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiD
   real(pReal),  intent(out) :: &
     localphiDot, &
     dLocalphiDot_dPhi
-  integer :: &
-    instance, sourceOffset
 
-  instance = source_damage_isoBrittle_instance(phase)
+  integer :: &
+    sourceOffset
+
   sourceOffset = source_damage_isoBrittle_offset(phase)
 
-  localphiDot = (1.0_pReal - phi)**(param(instance)%N - 1.0_pReal) - &
+  associate(prm => param(source_damage_isoBrittle_instance(phase)))
+  localphiDot = (1.0_pReal - phi)**(prm%n - 1.0_pReal) - &
                 phi*sourceState(phase)%p(sourceOffset)%state(1,constituent)
-  dLocalphiDot_dPhi = - (param(instance)%N - 1.0_pReal)* &
+  dLocalphiDot_dPhi = - (prm%n - 1.0_pReal)* &
-                        (1.0_pReal - phi)**max(0.0_pReal,param(instance)%N - 2.0_pReal) &
+                        (1.0_pReal - phi)**max(0.0_pReal,prm%n - 2.0_pReal) &
                       - sourceState(phase)%p(sourceOffset)%state(1,constituent)
+  end associate
 
 end subroutine source_damage_isoBrittle_getRateAndItsTangent
 

src/source_damage_isoDuctile.f90

@@ -15,6 +15,7 @@ module source_damage_isoDuctile
 
   implicit none
   private
+
   integer,                       dimension(:),           allocatable :: &
     source_damage_isoDuctile_offset, &                                                                 !< which source is my current damage mechanism?
     source_damage_isoDuctile_instance                                                                  !< instance of damage source mechanism
@@ -46,7 +47,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine source_damage_isoDuctile_init
 
-  integer :: Ninstance,instance,source,sourceOffset,NofMyPhase,p
+  integer :: Ninstance,source,sourceOffset,NofMyPhase,p
   character(len=pStringLen) :: &
     extmsg = ''
 
@@ -90,11 +91,10 @@ subroutine source_damage_isoDuctile_init
     prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
 
     NofMyPhase=count(material_phaseAt==p) * discretization_nIP
-    instance = source_damage_isoDuctile_instance(p)
     sourceOffset = source_damage_isoDuctile_offset(p)
 
     call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0)
-    sourceState(p)%p(sourceOffset)%aTolState=param(instance)%aTol
+    sourceState(p)%p(sourceOffset)%aTolState=prm%aTol
 
     end associate
   enddo
@@ -111,23 +111,28 @@ subroutine source_damage_isoDuctile_dotState(ipc, ip, el)
     ipc, &                                                                                          !< component-ID of integration point
     ip, &                                                                                           !< integration point
     el                                                                                              !< element
+
   integer :: &
-    phase, constituent, instance, homog, sourceOffset, damageOffset
+    phase, &
+    constituent, &
+    sourceOffset, &
+    damageOffset, &
+    homog
 
   phase = material_phaseAt(ipc,el)
   constituent = material_phasememberAt(ipc,ip,el)
-  instance = source_damage_isoDuctile_instance(phase)
   sourceOffset = source_damage_isoDuctile_offset(phase)
   homog = material_homogenizationAt(el)
   damageOffset = damageMapping(homog)%p(ip,el)
 
+  associate(prm => param(source_damage_isoDuctile_instance(phase)))
   sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = &
-    sum(plasticState(phase)%slipRate(:,constituent))/ &
+    sum(plasticState(phase)%slipRate(:,constituent))/(damage(homog)%p(damageOffset)**prm%N)/prm%critPlasticStrain
-    ((damage(homog)%p(damageOffset))**param(instance)%N)/ &
+  end associate
-    param(instance)%critPlasticStrain
 
 end subroutine source_damage_isoDuctile_dotState
 
+
 !--------------------------------------------------------------------------------------------------
 !> @brief returns local part of nonlocal damage driving force
 !--------------------------------------------------------------------------------------------------
@@ -141,16 +146,17 @@ subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiD
   real(pReal),  intent(out) :: &
     localphiDot, &
     dLocalphiDot_dPhi
+
   integer :: &
     sourceOffset
 
   sourceOffset = source_damage_isoDuctile_offset(phase)
 
-  localphiDot = 1.0_pReal &
-              - sourceState(phase)%p(sourceOffset)%state(1,constituent) * phi
-
   dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent)
 
+  localphiDot = 1.0_pReal &
+              + dLocalphiDot_dPhi*phi
+
 end subroutine source_damage_isoDuctile_getRateAndItsTangent