Merge branch 'clean-damage' into 'development'

Clean damage

See merge request damask/DAMASK!336

PRIVATE

@@ -1 +1 @@
-Subproject commit e74cf00628285a587ced1e551cc9837c1011ca1c
+Subproject commit 2ed5cd4ba97b44ad9c8b61ced94060aee57a2dd8

src/commercialFEM_fileList.f90

@@ -41,8 +41,6 @@
 #include "phase_damage_isoductile.f90"
 #include "phase_damage_anisobrittle.f90"
 #include "phase_damage_anisoductile.f90"
-#include "damage_none.f90"
-#include "damage_nonlocal.f90"
 #include "homogenization.f90"
 #include "homogenization_mechanical.f90"
 #include "homogenization_mechanical_pass.f90"

src/config.f90

@@ -5,16 +5,10 @@
 !! precedence over material.yaml.
 !--------------------------------------------------------------------------------------------------
 module config
-  use prec
-  use DAMASK_interface
   use IO
   use YAML_parse
   use YAML_types
 
-#ifdef PETSc
-#include <petsc/finclude/petscsys.h>
-   use petscsys
-#endif
 
   implicit none
   private
@@ -50,17 +44,12 @@ end subroutine config_init
 subroutine parse_material
 
   logical :: fileExists
-  character(len=:), allocatable :: fname
 
-  fname = getSolverJobName()//'.yaml'
+
-  inquire(file=fname,exist=fileExists)
+  inquire(file='material.yaml',exist=fileExists)
-  if(.not. fileExists) then
+  if(.not. fileExists) call IO_error(100,ext_msg='material.yaml')
-    fname = 'material.yaml'
+  print*, 'reading material.yaml'; flush(IO_STDOUT)
-    inquire(file=fname,exist=fileExists)
+  config_material => YAML_parse_file('material.yaml')
-    if(.not. fileExists) call IO_error(100,ext_msg=fname)
-  endif
-  print*, 'reading '//fname; flush(IO_STDOUT)
-  config_material => YAML_parse_file(fname)
 
 end subroutine parse_material
 
@@ -72,6 +61,7 @@ subroutine parse_numerics
 
   logical :: fexist
 
+
   config_numerics => emptyDict
   inquire(file='numerics.yaml', exist=fexist)
   if (fexist) then
@@ -89,6 +79,7 @@ subroutine parse_debug
 
   logical :: fexist
 
+
   config_debug => emptyDict
   inquire(file='debug.yaml', exist=fexist)
   fileExists: if (fexist) then

src/damage_none.f90

@@ -1,38 +0,0 @@
-!--------------------------------------------------------------------------------------------------
-!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief material subroutine for constant damage field
-!--------------------------------------------------------------------------------------------------
-module damage_none
-  use prec
-  use config
-  use material
-
-  implicit none
-  public
-
-contains
-
-!--------------------------------------------------------------------------------------------------
-!> @brief allocates all neccessary fields, reads information from material configuration file
-!--------------------------------------------------------------------------------------------------
-subroutine damage_none_init
-
-  integer :: h,Nmaterialpoints
-
-  print'(/,a)', ' <<<+-  damage_none init  -+>>>'; flush(6)
-
-  do h = 1, size(material_name_homogenization)
-    if (damage_type(h) /= DAMAGE_NONE_ID) cycle
-
-    Nmaterialpoints = count(material_homogenizationAt == h)
-    damageState_h(h)%sizeState = 0
-    allocate(damageState_h(h)%state0   (0,Nmaterialpoints))
-    allocate(damageState_h(h)%state    (0,Nmaterialpoints))
-
-    allocate  (damage(h)%p(Nmaterialpoints), source=1.0_pReal)
-
-  enddo
-
-end subroutine damage_none_init
-
-end module damage_none

src/damage_nonlocal.f90

@@ -1,94 +0,0 @@
-!--------------------------------------------------------------------------------------------------
-!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief material subroutine for non-locally evolving damage field
-!--------------------------------------------------------------------------------------------------
-module damage_nonlocal
-  use prec
-  use material
-  use config
-  use YAML_types
-  use lattice
-  use phase
-  use results
-
-  implicit none
-  private
-
-  type, private :: tNumerics
-    real(pReal) :: &
-    charLength                                                                                      !< characteristic length scale for gradient problems
-  end type tNumerics
-
-  type(tNumerics), private :: &
-    num
-
-  public :: &
-    damage_nonlocal_init, &
-    damage_nonlocal_getDiffusion
-
-contains
-
-!--------------------------------------------------------------------------------------------------
-!> @brief module initialization
-!> @details reads in material parameters, allocates arrays, and does sanity checks
-!--------------------------------------------------------------------------------------------------
-subroutine damage_nonlocal_init
-
-  integer :: Ninstances,Nmaterialpoints,h
-  class(tNode), pointer :: &
-    num_generic, &
-    material_homogenization
-
-  print'(/,a)', ' <<<+-  damage_nonlocal init  -+>>>'; flush(6)
-
-!------------------------------------------------------------------------------------
-! read numerics parameter
-  num_generic => config_numerics%get('generic',defaultVal= emptyDict)
-  num%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal)
-
-  Ninstances = count(damage_type == DAMAGE_nonlocal_ID)
-
-  material_homogenization => config_material%get('homogenization')
-  do h = 1, material_homogenization%length
-    if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle
-
-    Nmaterialpoints = count(material_homogenizationAt == h)
-    damageState_h(h)%sizeState = 1
-    allocate(damageState_h(h)%state0   (1,Nmaterialpoints), source=1.0_pReal)
-    allocate(damageState_h(h)%state    (1,Nmaterialpoints), source=1.0_pReal)
-
-    damage(h)%p => damageState_h(h)%state(1,:)
-
-  enddo
-
-end subroutine damage_nonlocal_init
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief returns homogenized non local damage diffusion tensor in reference configuration
-!--------------------------------------------------------------------------------------------------
-function damage_nonlocal_getDiffusion(ip,el)
-
-  integer, intent(in) :: &
-    ip, &                                                                                           !< integration point number
-    el                                                                                              !< element number
-  real(pReal), dimension(3,3) :: &
-    damage_nonlocal_getDiffusion
-  integer :: &
-    homog, &
-    grain
-
-  homog  = material_homogenizationAt(el)
-  damage_nonlocal_getDiffusion = 0.0_pReal
-  do grain = 1, homogenization_Nconstituents(homog)
-    damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + &
-      crystallite_push33ToRef(grain,ip,el,lattice_D(1:3,1:3,material_phaseAt(grain,el)))
-  enddo
-
-  damage_nonlocal_getDiffusion = &
-    num%charLength**2*damage_nonlocal_getDiffusion/real(homogenization_Nconstituents(homog),pReal)
-
-end function damage_nonlocal_getDiffusion
-
-
-end module damage_nonlocal

src/grid/grid_damage_spectral.f90

@@ -15,7 +15,6 @@ module grid_damage_spectral
   use IO
   use spectral_utilities
   use discretization_grid
-  use damage_nonlocal
   use YAML_types
   use homogenization
   use config

src/homogenization.f90

@@ -12,14 +12,53 @@ module homogenization
   use material
   use phase
   use discretization
-  use damage_none
-  use damage_nonlocal
   use HDF5_utilities
   use results
+  use lattice
 
   implicit none
   private
 
+
+  enum, bind(c); enumerator :: &
+    THERMAL_ISOTHERMAL_ID, &
+    THERMAL_CONDUCTION_ID, &
+    DAMAGE_NONE_ID, &
+    DAMAGE_NONLOCAL_ID, &
+    HOMOGENIZATION_UNDEFINED_ID, &
+    HOMOGENIZATION_NONE_ID, &
+    HOMOGENIZATION_ISOSTRAIN_ID, &
+    HOMOGENIZATION_RGC_ID
+  end enum
+
+    type(tState),        allocatable, dimension(:), public :: &
+    homogState, &
+    damageState_h
+
+  integer, dimension(:), allocatable, public, protected :: &
+    homogenization_typeInstance, &                                                                  !< instance of particular type of each homogenization
+    thermal_typeInstance, &                                                                         !< instance of particular type of each thermal transport
+    damage_typeInstance                                                                             !< instance of particular type of each nonlocal damage
+
+  real(pReal), dimension(:), allocatable, public, protected :: &
+    thermal_initialT
+
+  integer(kind(THERMAL_isothermal_ID)),       dimension(:),   allocatable, public, protected :: &
+    thermal_type                                                                                    !< thermal transport model
+  integer(kind(DAMAGE_none_ID)),              dimension(:),   allocatable, public, protected :: &
+    damage_type                                                                                     !< nonlocal damage model
+  integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:),   allocatable, public, protected :: &
+    homogenization_type                                                                             !< type of each homogenization
+
+  type, private :: tNumerics_damage
+    real(pReal) :: &
+    charLength                                                                                      !< characteristic length scale for gradient problems
+  end type tNumerics_damage
+
+  type(tNumerics_damage), private :: &
+    num_damage
+
+
   logical, public :: &
     terminallyIll = .false.                                                                         !< at least one material point is terminally ill
 
@@ -193,7 +232,13 @@ module homogenization
     homogenization_forward, &
     homogenization_results, &
     homogenization_restartRead, &
-    homogenization_restartWrite
+    homogenization_restartWrite, &
+    THERMAL_CONDUCTION_ID, &
+    DAMAGE_NONLOCAL_ID
+
+  public :: &
+    damage_nonlocal_init, &
+    damage_nonlocal_getDiffusion
 
 contains
 
@@ -209,6 +254,12 @@ subroutine homogenization_init()
 
   print'(/,a)', ' <<<+-  homogenization init  -+>>>'; flush(IO_STDOUT)
 
+
+  allocate(homogState      (size(material_name_homogenization)))
+  allocate(damageState_h   (size(material_name_homogenization)))
+  call material_parseHomogenization
+
+
   num_homog        => config_numerics%get('homogenization',defaultVal=emptyDict)
   num_homogGeneric => num_homog%get('generic',defaultVal=emptyDict)
 
@@ -220,8 +271,7 @@ subroutine homogenization_init()
   call thermal_init()
   call damage_init()
 
-  if (any(damage_type == DAMAGE_none_ID))      call damage_none_init
+  call damage_nonlocal_init
-  if (any(damage_type == DAMAGE_nonlocal_ID))  call damage_nonlocal_init
 
 
 end subroutine homogenization_init
@@ -257,6 +307,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
 
       if(homogState(ho)%sizeState > 0)  homogState(ho)%State(:,me) = homogState(ho)%State0(:,me)
       if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%State(:,me) = damageState_h(ho)%State0(:,me)
+      call damage_partition(ce)
 
       doneAndHappy = [.false.,.true.]
 
@@ -311,26 +362,6 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
     enddo
     !$OMP END DO
 
-!    !$OMP DO PRIVATE(ho,ph,ce)
-!    do el = FEsolving_execElem(1),FEsolving_execElem(2)
-!      if (terminallyIll) continue
-!      ho = material_homogenizationAt(el)
-!      do ip = FEsolving_execIP(1),FEsolving_execIP(2)
-!        ce = (el-1)*discretization_nIPs + ip
-!        call damage_partition(ce)
-!        do co = 1, homogenization_Nconstituents(ho)
-!          ph = material_phaseAt(co,el)
-!          if (.not. thermal_stress(dt,ph,material_phaseMemberAt(co,ip,el))) then
-!            if (.not. terminallyIll) &                                                           ! so first signals terminally ill...
-!              print*, ' Integration point ', ip,' at element ', el, ' terminally ill'
-!            terminallyIll = .true.                                                                  ! ...and kills all others
-!         endif
-!         call thermal_homogenize(ip,el)
-!        enddo
-!      enddo
-!    enddo
-!    !$OMP END DO
-
     !$OMP DO PRIVATE(ho)
     elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2)
       ho = material_homogenizationAt(el)
@@ -397,7 +428,8 @@ subroutine homogenization_forward
 
   do ho = 1, size(material_name_homogenization)
     homogState (ho)%state0 = homogState (ho)%state
-    damageState_h(ho)%state0 = damageState_h(ho)%state
+    if(damageState_h(ho)%sizeState > 0) &
+      damageState_h(ho)%state0 = damageState_h(ho)%state
   enddo
 
 end subroutine homogenization_forward
@@ -457,4 +489,143 @@ subroutine homogenization_restartRead(fileHandle)
 end subroutine homogenization_restartRead
 
 
+
+!--------------------------------------------------------------------------------------------------
+!> @brief module initialization
+!> @details reads in material parameters, allocates arrays, and does sanity checks
+!--------------------------------------------------------------------------------------------------
+subroutine damage_nonlocal_init
+
+  integer :: Ninstances,Nmaterialpoints,h
+  class(tNode), pointer :: &
+    num_generic, &
+    material_homogenization
+
+  print'(/,a)', ' <<<+-  damage_nonlocal init  -+>>>'; flush(6)
+
+!------------------------------------------------------------------------------------
+! read numerics parameter
+  num_generic => config_numerics%get('generic',defaultVal= emptyDict)
+  num_damage%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal)
+
+  Ninstances = count(damage_type == DAMAGE_nonlocal_ID)
+
+  material_homogenization => config_material%get('homogenization')
+  do h = 1, material_homogenization%length
+    if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle
+
+    Nmaterialpoints = count(material_homogenizationAt == h)
+    damageState_h(h)%sizeState = 1
+    allocate(damageState_h(h)%state0   (1,Nmaterialpoints), source=1.0_pReal)
+    allocate(damageState_h(h)%state    (1,Nmaterialpoints), source=1.0_pReal)
+
+  enddo
+
+end subroutine damage_nonlocal_init
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief returns homogenized non local damage diffusion tensor in reference configuration
+!--------------------------------------------------------------------------------------------------
+function damage_nonlocal_getDiffusion(ip,el)
+
+  integer, intent(in) :: &
+    ip, &                                                                                           !< integration point number
+    el                                                                                              !< element number
+  real(pReal), dimension(3,3) :: &
+    damage_nonlocal_getDiffusion
+  integer :: &
+    homog, &
+    grain
+
+  homog  = material_homogenizationAt(el)
+  damage_nonlocal_getDiffusion = 0.0_pReal
+  do grain = 1, homogenization_Nconstituents(homog)
+    damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + &
+      crystallite_push33ToRef(grain,ip,el,lattice_D(1:3,1:3,material_phaseAt(grain,el)))
+  enddo
+
+  damage_nonlocal_getDiffusion = &
+    num_damage%charLength**2*damage_nonlocal_getDiffusion/real(homogenization_Nconstituents(homog),pReal)
+
+end function damage_nonlocal_getDiffusion
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief parses the homogenization part from the material configuration
+! ToDo: This should be done in homogenization
+!--------------------------------------------------------------------------------------------------
+subroutine material_parseHomogenization
+
+  class(tNode), pointer :: &
+    material_homogenization, &
+    homog, &
+    homogMech, &
+    homogThermal, &
+    homogDamage
+
+  integer :: h
+
+  material_homogenization => config_material%get('homogenization')
+
+  allocate(homogenization_type(size(material_name_homogenization)),           source=HOMOGENIZATION_undefined_ID)
+  allocate(thermal_type(size(material_name_homogenization)),                  source=THERMAL_isothermal_ID)
+  allocate(damage_type (size(material_name_homogenization)),                  source=DAMAGE_none_ID)
+  allocate(homogenization_typeInstance(size(material_name_homogenization)),   source=0)
+  allocate(thermal_typeInstance(size(material_name_homogenization)),          source=0)
+  allocate(damage_typeInstance(size(material_name_homogenization)),           source=0)
+  allocate(thermal_initialT(size(material_name_homogenization)),              source=300.0_pReal)
+
+  do h=1, size(material_name_homogenization)
+    homog => material_homogenization%get(h)
+    homogMech => homog%get('mechanics')
+    select case (homogMech%get_asString('type'))
+      case('pass')
+        homogenization_type(h) = HOMOGENIZATION_NONE_ID
+      case('isostrain')
+        homogenization_type(h) = HOMOGENIZATION_ISOSTRAIN_ID
+      case('RGC')
+        homogenization_type(h) = HOMOGENIZATION_RGC_ID
+      case default
+        call IO_error(500,ext_msg=homogMech%get_asString('type'))
+    end select
+
+    homogenization_typeInstance(h) = count(homogenization_type==homogenization_type(h))
+
+    if(homog%contains('thermal')) then
+      homogThermal => homog%get('thermal')
+        thermal_initialT(h) =  homogThermal%get_asFloat('T_0',defaultVal=300.0_pReal)
+
+        select case (homogThermal%get_asString('type'))
+          case('isothermal')
+            thermal_type(h) = THERMAL_isothermal_ID
+          case('conduction')
+            thermal_type(h) = THERMAL_conduction_ID
+          case default
+            call IO_error(500,ext_msg=homogThermal%get_asString('type'))
+        end select
+    endif
+
+    if(homog%contains('damage')) then
+      homogDamage => homog%get('damage')
+        select case (homogDamage%get_asString('type'))
+          case('none')
+            damage_type(h) = DAMAGE_none_ID
+          case('nonlocal')
+            damage_type(h) = DAMAGE_nonlocal_ID
+          case default
+            call IO_error(500,ext_msg=homogDamage%get_asString('type'))
+        end select
+    endif
+  enddo
+
+  do h=1, size(material_name_homogenization)
+    homogenization_typeInstance(h)  = count(homogenization_type(1:h) == homogenization_type(h))
+    thermal_typeInstance(h)         = count(thermal_type       (1:h) == thermal_type       (h))
+    damage_typeInstance(h)          = count(damage_type        (1:h) == damage_type        (h))
+  enddo
+
+end subroutine material_parseHomogenization
+
+
 end module homogenization

src/homogenization_damage.f90

@@ -72,7 +72,8 @@ module subroutine damage_partition(ce)
   integer :: co
 
 
-  phi     = current(material_homogenizationAt2(ce))%phi(material_homogenizationMemberAt2(ce))
+  if(damageState_h(material_homogenizationAt2(ce))%sizeState < 1) return
+  phi     = damagestate_h(material_homogenizationAt2(ce))%state(1,material_homogenizationMemberAt2(ce))
   do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce))
     call phase_damage_set_phi(phi,co,ce)
   enddo
@@ -143,7 +144,7 @@ module subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el)
 
   homog  = material_homogenizationAt(el)
   offset = material_homogenizationMemberAt(ip,el)
-  damage(homog)%p(offset) = phi
+  damagestate_h(homog)%state(1,offset) = phi
 
 end subroutine damage_nonlocal_putNonLocalDamage
 
@@ -162,7 +163,7 @@ module subroutine damage_nonlocal_results(homog,group)
   outputsLoop: do o = 1,size(prm%output)
     select case(prm%output(o))
       case ('phi')
-        call results_writeDataset(group,damage(homog)%p,prm%output(o),&
+        call results_writeDataset(group,damagestate_h(homog)%state(1,:),prm%output(o),&
                                   'damage indicator','-')
     end select
   enddo outputsLoop

src/homogenization_mechanical.f90

@@ -2,7 +2,7 @@
 !> @author Martin Diehl, KU Leuven
 !> @brief Partition F and homogenize P/dPdF
 !--------------------------------------------------------------------------------------------------
-submodule(homogenization) mechanics
+submodule(homogenization) mechanical
 
 
   interface
@@ -86,7 +86,7 @@ module subroutine mechanical_init(num_homog)
   class(tNode), pointer :: &
     num_homogMech
 
-  print'(/,a)', ' <<<+-  homogenization:mechanics init  -+>>>'
+  print'(/,a)', ' <<<+-  homogenization:mechanical init  -+>>>'
 
   allocate(homogenization_dPdF(3,3,3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal)
   homogenization_F0 = spread(math_I3,3,discretization_nIPs*discretization_Nelems)                   ! initialize to identity
@@ -225,8 +225,6 @@ end function mechanical_updateState
 !> @brief Write results to file.
 !--------------------------------------------------------------------------------------------------
 module subroutine mechanical_results(group_base,h)
-  use material, only: &
-    material_homogenization_type => homogenization_type
 
   character(len=*), intent(in) :: group_base
   integer, intent(in)          :: h
@@ -236,7 +234,7 @@ module subroutine mechanical_results(group_base,h)
   group = trim(group_base)//'/mech'
   call results_closeGroup(results_addGroup(group))
 
-  select case(material_homogenization_type(h))
+  select case(homogenization_type(h))
 
     case(HOMOGENIZATION_rgc_ID)
       call mechanical_RGC_results(homogenization_typeInstance(h),group)
@@ -253,4 +251,4 @@ module subroutine mechanical_results(group_base,h)
 end subroutine mechanical_results
 
 
-end submodule mechanics
+end submodule mechanical

src/homogenization_mechanical_RGC.f90

@@ -6,7 +6,7 @@
 !> @brief Relaxed grain cluster (RGC) homogenization scheme
 !> N_constituents is defined as p x q x r (cluster)
 !--------------------------------------------------------------------------------------------------
-submodule(homogenization:mechanics) RGC
+submodule(homogenization:mechanical) RGC
   use rotations
   use lattice
 
@@ -88,7 +88,7 @@ module subroutine mechanical_RGC_init(num_homogMech)
     homog, &
     homogMech
 
-  print'(/,a)', ' <<<+-  homogenization:mechanics:RGC init  -+>>>'
+  print'(/,a)', ' <<<+-  homogenization:mechanical:RGC init  -+>>>'
 
   Ninstances = count(homogenization_type == HOMOGENIZATION_RGC_ID)
   print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)

src/homogenization_mechanical_isostrain.f90

@@ -4,7 +4,7 @@
 !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
 !> @brief Isostrain (full constraint Taylor assuption) homogenization scheme
 !--------------------------------------------------------------------------------------------------
-submodule(homogenization:mechanics) isostrain
+submodule(homogenization:mechanical) isostrain
 
   enum, bind(c); enumerator :: &
     parallel_ID, &
@@ -37,7 +37,7 @@ module subroutine mechanical_isostrain_init
     homog, &
     homogMech
 
-  print'(/,a)', ' <<<+-  homogenization:mechanics:isostrain init  -+>>>'
+  print'(/,a)', ' <<<+-  homogenization:mechanical:isostrain init  -+>>>'
 
   Ninstances = count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)
   print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)

src/homogenization_mechanical_pass.f90

@@ -4,7 +4,7 @@
 !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
 !> @brief dummy homogenization homogenization scheme for 1 constituent per material point
 !--------------------------------------------------------------------------------------------------
-submodule(homogenization:mechanics) none
+submodule(homogenization:mechanical) none
 
 contains
 
@@ -18,7 +18,7 @@ module subroutine mechanical_pass_init
     h, &
     Nmaterialpoints
 
-  print'(/,a)', ' <<<+-  homogenization:mechanics:pass init  -+>>>'
+  print'(/,a)', ' <<<+-  homogenization:mechanical:pass init  -+>>>'
 
   Ninstances = count(homogenization_type == HOMOGENIZATION_NONE_ID)
   print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)

src/lattice.f90

@@ -459,7 +459,8 @@ subroutine lattice_init
     phase, &
     mech, &
     elasticity, &
-    thermal
+    thermal, &
+    damage
 
   print'(/,a)', ' <<<+-  lattice init  -+>>>'; flush(IO_STDOUT)
 
@@ -535,13 +536,17 @@ subroutine lattice_init
     endif
 
 
-    lattice_D(1,1,ph) = phase%get_asFloat('D_11',defaultVal=0.0_pReal)
+    if (phase%contains('damage')) then
-    lattice_D(2,2,ph) = phase%get_asFloat('D_22',defaultVal=0.0_pReal)
+      damage => phase%get('damage')
-    lattice_D(3,3,ph) = phase%get_asFloat('D_33',defaultVal=0.0_pReal)
+      damage => damage%get(1)
-    lattice_D(1:3,1:3,ph) = lattice_applyLatticeSymmetry33(lattice_D(1:3,1:3,ph), &
+      lattice_D(1,1,ph) = damage%get_asFloat('D_11',defaultVal=0.0_pReal)
+      lattice_D(2,2,ph) = damage%get_asFloat('D_22',defaultVal=0.0_pReal)
+      lattice_D(3,3,ph) = damage%get_asFloat('D_33',defaultVal=0.0_pReal)
+      lattice_D(1:3,1:3,ph) = lattice_applyLatticeSymmetry33(lattice_D(1:3,1:3,ph), &
                                                           phase%get_asString('lattice'))
 
-    lattice_M(ph) = phase%get_asFloat('M',defaultVal=0.0_pReal)
+      lattice_M(ph) = damage%get_asFloat('M',defaultVal=0.0_pReal)
+    endif
     ! SHOULD NOT BE PART OF LATTICE END
 
     call selfTest

src/material.f90

@@ -6,50 +6,26 @@
 !--------------------------------------------------------------------------------------------------
 module material
   use prec
-  use math
   use config
   use results
   use IO
   use rotations
   use discretization
+  use YAML_types
 
   implicit none
   private
 
-  enum, bind(c); enumerator :: &
+  integer, dimension(:), allocatable, public, protected :: &
-    THERMAL_ISOTHERMAL_ID, &
+    homogenization_Nconstituents                                                                 !< number of grains in each homogenization
-    THERMAL_CONDUCTION_ID, &
-    DAMAGE_NONE_ID, &
-    DAMAGE_NONLOCAL_ID, &
-    HOMOGENIZATION_UNDEFINED_ID, &
-    HOMOGENIZATION_NONE_ID, &
-    HOMOGENIZATION_ISOSTRAIN_ID, &
-    HOMOGENIZATION_RGC_ID
-  end enum
 
   character(len=:), public, protected, allocatable, dimension(:) :: &
     material_name_phase, &                                                                          !< name of each phase
     material_name_homogenization                                                                    !< name of each homogenization
 
-  integer(kind(THERMAL_isothermal_ID)),       dimension(:),   allocatable, public, protected :: &
-    thermal_type                                                                                    !< thermal transport model
-  integer(kind(DAMAGE_none_ID)),              dimension(:),   allocatable, public, protected :: &
-    damage_type                                                                                     !< nonlocal damage model
-  integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:),   allocatable, public, protected :: &
-    homogenization_type                                                                             !< type of each homogenization
-
   integer, public, protected :: &
     homogenization_maxNconstituents                                                                 !< max number of grains in any USED homogenization
 
-  integer, dimension(:), allocatable, public, protected :: &
-    homogenization_Nconstituents, &                                                                 !< number of grains in each homogenization
-    homogenization_typeInstance, &                                                                  !< instance of particular type of each homogenization
-    thermal_typeInstance, &                                                                         !< instance of particular type of each thermal transport
-    damage_typeInstance                                                                             !< instance of particular type of each nonlocal damage
-
-  real(pReal), dimension(:), allocatable, public, protected :: &
-    thermal_initialT                                                                                !< initial temperature per each homogenization
-
   integer, dimension(:),     allocatable, public, protected :: &                                    ! (elem)
     material_homogenizationAt, &                                                                    !< homogenization ID of each element
     material_homogenizationAt2, &                                                                   !< per cell
@@ -63,50 +39,28 @@ module material
   integer, dimension(:,:,:), allocatable, public, protected :: &                                    ! (constituent,IP,elem)
     material_phaseMemberAt                                                                          !< position of the element within its phase instance
 
-  type(tState),        allocatable, dimension(:), public :: &
-    homogState, &
-    damageState_h
-
   type(Rotation), dimension(:,:,:), allocatable, public, protected :: &
     material_orientation0                                                                           !< initial orientation of each grain,IP,element
 
-  type(group_float),  allocatable, dimension(:), public :: &
-    damage                                                                                       !< damage field
-
   public :: &
-    material_init, &
+    material_init
-    THERMAL_ISOTHERMAL_ID, &
-    THERMAL_CONDUCTION_ID, &
-    DAMAGE_NONE_ID, &
-    DAMAGE_NONLOCAL_ID, &
-    HOMOGENIZATION_NONE_ID, &
-    HOMOGENIZATION_ISOSTRAIN_ID, &
-    HOMOGENIZATION_RGC_ID
 
 contains
 
 !--------------------------------------------------------------------------------------------------
-!> @brief parses material configuration file
+!> @brief Parse material configuration file (material.yaml).
 !--------------------------------------------------------------------------------------------------
 subroutine material_init(restart)
 
   logical, intent(in) :: restart
 
+
   print'(/,a)', ' <<<+-  material init  -+>>>'; flush(IO_STDOUT)
 
 
   call material_parseMaterial
   print*, 'Material parsed'
 
-  call material_parseHomogenization
-  print*, 'Homogenization parsed'
-
-
-  allocate(homogState      (size(material_name_homogenization)))
-  allocate(damageState_h     (size(material_name_homogenization)))
-
-  allocate(damage          (size(material_name_homogenization)))
-
 
   if (.not. restart) then
     call results_openJobFile
@@ -118,82 +72,6 @@ subroutine material_init(restart)
 end subroutine material_init
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief parses the homogenization part from the material configuration
-! ToDo: This should be done in homogenization
-!--------------------------------------------------------------------------------------------------
-subroutine material_parseHomogenization
-
-  class(tNode), pointer :: &
-    material_homogenization, &
-    homog, &
-    homogMech, &
-    homogThermal, &
-    homogDamage
-
-  integer :: h
-
-  material_homogenization => config_material%get('homogenization')
-
-  allocate(homogenization_type(size(material_name_homogenization)),           source=HOMOGENIZATION_undefined_ID)
-  allocate(thermal_type(size(material_name_homogenization)),                  source=THERMAL_isothermal_ID)
-  allocate(damage_type (size(material_name_homogenization)),                  source=DAMAGE_none_ID)
-  allocate(homogenization_typeInstance(size(material_name_homogenization)),   source=0)
-  allocate(thermal_typeInstance(size(material_name_homogenization)),          source=0)
-  allocate(damage_typeInstance(size(material_name_homogenization)),           source=0)
-  allocate(thermal_initialT(size(material_name_homogenization)),              source=300.0_pReal)
-
-  do h=1, size(material_name_homogenization)
-    homog => material_homogenization%get(h)
-    homogMech => homog%get('mechanics')
-    select case (homogMech%get_asString('type'))
-      case('pass')
-        homogenization_type(h) = HOMOGENIZATION_NONE_ID
-      case('isostrain')
-        homogenization_type(h) = HOMOGENIZATION_ISOSTRAIN_ID
-      case('RGC')
-        homogenization_type(h) = HOMOGENIZATION_RGC_ID
-      case default
-        call IO_error(500,ext_msg=homogMech%get_asString('type'))
-    end select
-
-    homogenization_typeInstance(h) = count(homogenization_type==homogenization_type(h))
-
-    if(homog%contains('thermal')) then
-      homogThermal => homog%get('thermal')
-        thermal_initialT(h) =  homogThermal%get_asFloat('T_0',defaultVal=300.0_pReal)
-
-        select case (homogThermal%get_asString('type'))
-          case('isothermal')
-            thermal_type(h) = THERMAL_isothermal_ID
-          case('conduction')
-            thermal_type(h) = THERMAL_conduction_ID
-          case default
-            call IO_error(500,ext_msg=homogThermal%get_asString('type'))
-        end select
-    endif
-
-    if(homog%contains('damage')) then
-      homogDamage => homog%get('damage')
-        select case (homogDamage%get_asString('type'))
-          case('none')
-            damage_type(h) = DAMAGE_none_ID
-          case('nonlocal')
-            damage_type(h) = DAMAGE_nonlocal_ID
-          case default
-            call IO_error(500,ext_msg=homogDamage%get_asString('type'))
-        end select
-    endif
-  enddo
-
-  do h=1, size(material_name_homogenization)
-    homogenization_typeInstance(h)  = count(homogenization_type(1:h) == homogenization_type(h))
-    thermal_typeInstance(h)         = count(thermal_type       (1:h) == thermal_type       (h))
-    damage_typeInstance(h)          = count(damage_type        (1:h) == damage_type        (h))
-  enddo
-
-end subroutine material_parseHomogenization
-
 
 !--------------------------------------------------------------------------------------------------
 !> @brief parses the material part in the material configuration file

src/phase.f90

@@ -58,20 +58,17 @@ module phase
   type(tDebugOptions) :: debugCrystallite
 
   integer, dimension(:), allocatable, public :: &                                                   !< ToDo: should be protected (bug in Intel compiler)
-    thermal_Nsources, &
+    phase_elasticityInstance, &
-    phase_Nsources, &                                                                               !< number of source mechanisms active in each phase
-    phase_Nkinematics, &                                                                            !< number of kinematic mechanisms active in each phase
     phase_NstiffnessDegradations, &                                                                 !< number of stiffness degradation mechanisms active in each phase
-    phase_plasticInstance, &                                                                     !< instance of particular plasticity of each phase
+    phase_plasticInstance
-    phase_elasticityInstance                                                                        !< instance of particular elasticity of each phase
 
   logical, dimension(:), allocatable, public :: &                                                   ! ToDo: should be protected (bug in Intel Compiler)
     phase_localPlasticity                                                                           !< flags phases with local constitutive law
 
   type(tPlasticState), allocatable, dimension(:), public :: &
     plasticState
-  type(tSourceState),  allocatable, dimension(:), public :: &
+  type(tState),  allocatable, dimension(:), public :: &
-    damageState, thermalState
+    damageState
 
 
   integer, public, protected :: &
@@ -181,6 +178,11 @@ module phase
       real(pReal) :: dot_T
     end function thermal_dot_T
 
+    module function damage_phi(ph,me) result(phi)
+      integer, intent(in) :: ph,me
+      real(pReal) :: phi
+    end function damage_phi
+
 
     module subroutine phase_mechanical_setF(F,co,ip,el)
       real(pReal), dimension(3,3), intent(in) :: F
@@ -261,6 +263,27 @@ module phase
         el                                                                                          !< element
     end subroutine plastic_dependentState
 
+
+    module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me)
+      integer, intent(in) :: ph, me
+      real(pReal),   intent(in),  dimension(3,3) :: &
+        S
+      real(pReal),   intent(out), dimension(3,3) :: &
+        Ld                                                                                          !< damage velocity gradient
+      real(pReal),   intent(out), dimension(3,3,3,3) :: &
+        dLd_dTstar                                                                                  !< derivative of Ld with respect to Tstar (4th-order tensor)
+    end subroutine kinematics_cleavage_opening_LiAndItsTangent
+
+    module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me)
+      integer, intent(in) :: ph, me
+      real(pReal),   intent(in),  dimension(3,3) :: &
+        S
+      real(pReal),   intent(out), dimension(3,3) :: &
+        Ld                                                                                          !< damage velocity gradient
+      real(pReal),   intent(out), dimension(3,3,3,3) :: &
+        dLd_dTstar                                                                                  !< derivative of Ld with respect to Tstar (4th-order tensor)
+    end subroutine kinematics_slipplane_opening_LiAndItsTangent
+
   end interface
 
 
@@ -345,14 +368,12 @@ subroutine phase_init
   PhaseLoop2:do ph = 1,phases%length
 !--------------------------------------------------------------------------------------------------
 ! partition and initialize state
-    plasticState(ph)%state             = plasticState(ph)%state0
+    plasticState(ph)%state = plasticState(ph)%state0
-    forall(so = 1:phase_Nsources(ph))
+    if(damageState(ph)%sizeState > 0) &
-      damageState(ph)%p(so)%state = damageState(ph)%p(so)%state0
+      damageState(ph)%state  = damageState(ph)%state0
-    end forall
-
-    phase_source_maxSizeDotState   = max(phase_source_maxSizeDotState, &
-                                                maxval(damageState(ph)%p%sizeDotState))
   enddo PhaseLoop2
+
+  phase_source_maxSizeDotState     = maxval(damageState%sizeDotState)
   phase_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState)
 
 end subroutine phase_init
@@ -362,7 +383,7 @@ end subroutine phase_init
 !> @brief Allocate the components of the state structure for a given phase
 !--------------------------------------------------------------------------------------------------
 subroutine phase_allocateState(state, &
-                                  Nconstituents,sizeState,sizeDotState,sizeDeltaState)
+                               Nconstituents,sizeState,sizeDotState,sizeDeltaState)
 
   class(tState), intent(out) :: &
     state
@@ -399,15 +420,13 @@ subroutine phase_restore(ce,includeL)
   integer, intent(in) :: ce
 
   integer :: &
-    co, &                                                                                            !< constituent number
+    co
-    so
 
 
   do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce))
-    do so = 1, phase_Nsources(material_phaseAt2(co,ce))
+    if (damageState(material_phaseAt2(co,ce))%sizeState > 0) &
-      damageState(material_phaseAt2(co,ce))%p(so)%state( :,material_phasememberAt2(co,ce)) = &
+    damageState(material_phaseAt2(co,ce))%state( :,material_phasememberAt2(co,ce)) = &
-      damageState(material_phaseAt2(co,ce))%p(so)%state0(:,material_phasememberAt2(co,ce))
+      damageState(material_phaseAt2(co,ce))%state0(:,material_phasememberAt2(co,ce))
-    enddo
   enddo
 
   call mechanical_restore(ce,includeL)
@@ -421,16 +440,16 @@ end subroutine phase_restore
 !--------------------------------------------------------------------------------------------------
 subroutine phase_forward()
 
-  integer :: ph, so
+  integer :: ph
 
 
   call mechanical_forward()
   call thermal_forward()
 
   do ph = 1, size(damageState)
-    do so = 1,phase_Nsources(ph)
+    if (damageState(ph)%sizeState > 0) &
-      damageState(ph)%p(so)%state0 = damageState(ph)%p(so)%state
+      damageState(ph)%state0 = damageState(ph)%state
-  enddo; enddo
+  enddo
 
 end subroutine phase_forward
 
@@ -526,9 +545,8 @@ subroutine crystallite_init()
   phases => config_material%get('phase')
 
   do ph = 1, phases%length
-    do so = 1, phase_Nsources(ph)
+    if (damageState(ph)%sizeState > 0) &
-      allocate(damageState(ph)%p(so)%subState0,source=damageState(ph)%p(so)%state0)                 ! ToDo: hack
+      allocate(damageState(ph)%subState0,source=damageState(ph)%state0)                 ! ToDo: hack
-    enddo
   enddo
 
   print'(a42,1x,i10)', '    # of elements:                       ', eMax
@@ -574,9 +592,8 @@ subroutine phase_windForward(ip,el)
 
     call mechanical_windForward(ph,me)
 
-    do so = 1, phase_Nsources(material_phaseAt(co,el))
+    if(damageState(ph)%sizeState > 0) damageState(ph)%state0(:,me) = damageState(ph)%state(:,me)
-      damageState(ph)%p(so)%state0(:,me) = damageState(ph)%p(so)%state(:,me)
+
-    enddo
 
   enddo
 

src/phase_damage.f90

@@ -15,110 +15,92 @@ submodule(phase) damagee
     real(pReal), dimension(:), allocatable :: phi, d_phi_d_dot_phi
   end type tDataContainer
 
-  integer(kind(DAMAGE_UNDEFINED_ID)),     dimension(:,:), allocatable :: &
+  integer(kind(DAMAGE_UNDEFINED_ID)),     dimension(:), allocatable :: &
     phase_source                                                                                !< active sources mechanisms of each phase
 
+  integer, dimension(:), allocatable :: &
+    phase_Nsources
+
   type(tDataContainer), dimension(:), allocatable :: current
 
   interface
 
-    module function anisobrittle_init(source_length) result(mySources)
+    module function anisobrittle_init() result(mySources)
-      integer, intent(in) :: source_length
+      logical, dimension(:), allocatable :: mySources
-      logical, dimension(:,:), allocatable :: mySources
     end function anisobrittle_init
 
-    module function anisoductile_init(source_length) result(mySources)
+    module function anisoductile_init() result(mySources)
-      integer, intent(in) :: source_length
+      logical, dimension(:), allocatable :: mySources
-      logical, dimension(:,:), allocatable :: mySources
     end function anisoductile_init
 
-    module function isobrittle_init(source_length) result(mySources)
+    module function isobrittle_init() result(mySources)
-      integer, intent(in) :: source_length
+      logical, dimension(:), allocatable :: mySources
-      logical, dimension(:,:), allocatable :: mySources
     end function isobrittle_init
 
-    module function isoductile_init(source_length) result(mySources)
+    module function isoductile_init() result(mySources)
-      integer, intent(in) :: source_length
+      logical, dimension(:), allocatable :: mySources
-      logical, dimension(:,:), allocatable :: mySources
     end function isoductile_init
 
 
-    module subroutine source_damage_isoBrittle_deltaState(C, Fe, ph, me)
+    module subroutine isobrittle_deltaState(C, Fe, ph, me)
       integer, intent(in) :: ph,me
       real(pReal),  intent(in), dimension(3,3) :: &
         Fe
       real(pReal),  intent(in), dimension(6,6) :: &
         C
-    end subroutine source_damage_isoBrittle_deltaState
+    end subroutine isobrittle_deltaState
 
 
-    module subroutine anisobrittle_dotState(S, co, ip, el)
+    module subroutine anisobrittle_dotState(S, ph, me)
-      integer, intent(in) :: &
+      integer, intent(in) :: ph,me
-        co, &                                                                                      !< component-ID of integration point
-        ip, &                                                                                       !< integration point
-        el                                                                                          !< element
       real(pReal),  intent(in), dimension(3,3) :: &
         S
     end subroutine anisobrittle_dotState
 
-    module subroutine anisoductile_dotState(co, ip, el)
+    module subroutine anisoductile_dotState(ph,me)
-      integer, intent(in) :: &
+      integer, intent(in) :: ph,me
-        co, &                                                                                      !< component-ID of integration point
-        ip, &                                                                                       !< integration point
-        el                                                                                          !< element
     end subroutine anisoductile_dotState
 
-    module subroutine isoductile_dotState(co, ip, el)
+    module subroutine isoductile_dotState(ph,me)
-      integer, intent(in) :: &
+      integer, intent(in) :: ph,me
-        co, &                                                                                      !< component-ID of integration point
-        ip, &                                                                                       !< integration point
-        el                                                                                          !< element
     end subroutine isoductile_dotState
 
 
-    module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
+    module subroutine anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me)
-      integer, intent(in) :: &
+      integer, intent(in) :: ph,me
-        phase, &                                                                                      !< phase ID of element
-        constituent                                                                                   !< position of element within its phase instance
       real(pReal),  intent(in) :: &
         phi                                                                                           !< damage parameter
       real(pReal),  intent(out) :: &
         localphiDot, &
         dLocalphiDot_dPhi
-    end subroutine source_damage_anisoBrittle_getRateAndItsTangent
+    end subroutine anisobrittle_getRateAndItsTangent
 
-    module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
+    module subroutine anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me)
-      integer, intent(in) :: &
+      integer, intent(in) :: ph,me
-        phase, &                                                                                      !< phase ID of element
-        constituent                                                                                   !< position of element within its phase instance
       real(pReal),  intent(in) :: &
         phi                                                                                           !< damage parameter
       real(pReal),  intent(out) :: &
         localphiDot, &
         dLocalphiDot_dPhi
-    end subroutine source_damage_anisoDuctile_getRateAndItsTangent
+    end subroutine anisoductile_getRateAndItsTangent
 
-    module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
+    module subroutine isobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me)
-      integer, intent(in) :: &
+      integer, intent(in) :: ph,me
-        phase, &                                                                                      !< phase ID of element
-        constituent                                                                                   !< position of element within its phase instance
       real(pReal),  intent(in) :: &
         phi                                                                                           !< damage parameter
       real(pReal),  intent(out) :: &
         localphiDot, &
         dLocalphiDot_dPhi
-    end subroutine source_damage_isoBrittle_getRateAndItsTangent
+    end subroutine isobrittle_getRateAndItsTangent
 
-    module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
+    module subroutine isoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me)
-      integer, intent(in) :: &
+      integer, intent(in) :: ph,me
-        phase, &                                                                                      !< phase ID of element
-        constituent                                                                                   !< position of element within its phase instance
       real(pReal),  intent(in) :: &
         phi                                                                                           !< damage parameter
       real(pReal),  intent(out) :: &
         localphiDot, &
         dLocalphiDot_dPhi
-    end subroutine source_damage_isoDuctile_getRateAndItsTangent
+    end subroutine isoductile_getRateAndItsTangent
 
     module subroutine anisobrittle_results(phase,group)
       integer,          intent(in) :: phase
@@ -175,19 +157,20 @@ module subroutine damage_init
     allocate(current(ph)%d_phi_d_dot_phi(Nconstituents),source=0.0_pReal)
 
     phase => phases%get(ph)
-    sources => phase%get('source',defaultVal=emptyList)
+    sources => phase%get('damage',defaultVal=emptyList)
+    if (sources%length > 1) error stop
     phase_Nsources(ph) = sources%length
-    allocate(damageState(ph)%p(phase_Nsources(ph)))
+
   enddo
 
-  allocate(phase_source(maxval(phase_Nsources),phases%length), source = DAMAGE_UNDEFINED_ID)
+  allocate(phase_source(phases%length), source = DAMAGE_UNDEFINED_ID)
 
 ! initialize source mechanisms
   if(maxval(phase_Nsources) /= 0) then
-    where(isobrittle_init   (maxval(phase_Nsources))) phase_source = DAMAGE_ISOBRITTLE_ID
+    where(isobrittle_init()  ) phase_source = DAMAGE_ISOBRITTLE_ID
-    where(isoductile_init   (maxval(phase_Nsources))) phase_source = DAMAGE_ISODUCTILE_ID
+    where(isoductile_init()  ) phase_source = DAMAGE_ISODUCTILE_ID
-    where(anisobrittle_init (maxval(phase_Nsources))) phase_source = DAMAGE_ANISOBRITTLE_ID
+    where(anisobrittle_init()) phase_source = DAMAGE_ANISOBRITTLE_ID
-    where(anisoductile_init (maxval(phase_Nsources))) phase_source = DAMAGE_ANISODUCTILE_ID
+    where(anisoductile_init()) phase_source = DAMAGE_ANISODUCTILE_ID
   endif
 
 end subroutine damage_init
@@ -213,7 +196,6 @@ module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi,
   integer :: &
     ph, &
     co, &
-    so, &
     me
 
    phiDot = 0.0_pReal
@@ -222,19 +204,19 @@ module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi,
    do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
      ph = material_phaseAt(co,el)
      me = material_phasememberAt(co,ip,el)
-     do so = 1, phase_Nsources(ph)
+
-       select case(phase_source(so,ph))
+       select case(phase_source(ph))
          case (DAMAGE_ISOBRITTLE_ID)
-           call source_damage_isobrittle_getRateAndItsTangent  (localphiDot, dLocalphiDot_dPhi, phi, ph, me)
+           call isobrittle_getRateAndItsTangent  (localphiDot, dLocalphiDot_dPhi, phi, ph, me)
 
          case (DAMAGE_ISODUCTILE_ID)
-           call source_damage_isoductile_getRateAndItsTangent  (localphiDot, dLocalphiDot_dPhi, phi, ph, me)
+           call isoductile_getRateAndItsTangent  (localphiDot, dLocalphiDot_dPhi, phi, ph, me)
 
          case (DAMAGE_ANISOBRITTLE_ID)
-           call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me)
+           call anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me)
 
          case (DAMAGE_ANISODUCTILE_ID)
-           call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me)
+           call anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me)
 
          case default
          localphiDot = 0.0_pReal
@@ -243,7 +225,6 @@ module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi,
       end select
       phiDot = phiDot + localphiDot
       dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi
-    enddo
   enddo
 
 end subroutine phase_damage_getRateAndItsTangents
@@ -267,57 +248,52 @@ module function integrateDamageState(dt,co,ip,el) result(broken)
     NiterationState, &                                                                              !< number of iterations in state loop
     ph, &
     me, &
-    so
-  integer, dimension(maxval(phase_Nsources)) :: &
     size_so
   real(pReal) :: &
     zeta
   real(pReal), dimension(phase_source_maxSizeDotState) :: &
     r                                                                                               ! state residuum
-  real(pReal), dimension(phase_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState
+  real(pReal), dimension(phase_source_maxSizeDotState,2) :: source_dotState
   logical :: &
     converged_
 
-
   ph = material_phaseAt(co,el)
   me = material_phaseMemberAt(co,ip,el)
 
+  if (damageState(ph)%sizeState == 0) then
+    broken = .false.
+    return
+  endif
+
   converged_ = .true.
-  broken = phase_damage_collectDotState(co,ip,el,ph,me)
+  broken = phase_damage_collectDotState(ph,me)
   if(broken) return
 
-  do so = 1, phase_Nsources(ph)
+    size_so = damageState(ph)%sizeDotState
-    size_so(so) = damageState(ph)%p(so)%sizeDotState
+    damageState(ph)%state(1:size_so,me) = damageState(ph)%subState0(1:size_so,me) &
-    damageState(ph)%p(so)%state(1:size_so(so),me) = damageState(ph)%p(so)%subState0(1:size_so(so),me) &
+                                        + damageState(ph)%dotState (1:size_so,me) * dt
-                                                  + damageState(ph)%p(so)%dotState (1:size_so(so),me) * dt
+    source_dotState(1:size_so,2) = 0.0_pReal
-    source_dotState(1:size_so(so),2,so) = 0.0_pReal
-  enddo
 
   iteration: do NiterationState = 1, num%nState
 
-    do so = 1, phase_Nsources(ph)
+    if(nIterationState > 1) source_dotState(1:size_so,2) = source_dotState(1:size_so,1)
-      if(nIterationState > 1) source_dotState(1:size_so(so),2,so) = source_dotState(1:size_so(so),1,so)
+    source_dotState(1:size_so,1) = damageState(ph)%dotState(:,me)
-      source_dotState(1:size_so(so),1,so) = damageState(ph)%p(so)%dotState(:,me)
-    enddo
 
-    broken = phase_damage_collectDotState(co,ip,el,ph,me)
+    broken = phase_damage_collectDotState(ph,me)
     if(broken) exit iteration
 
-    do so = 1, phase_Nsources(ph)
+
-      zeta = damper(damageState(ph)%p(so)%dotState(:,me), &
+      zeta = damper(damageState(ph)%dotState(:,me),source_dotState(1:size_so,1),source_dotState(1:size_so,2))
-                    source_dotState(1:size_so(so),1,so),&
+      damageState(ph)%dotState(:,me) = damageState(ph)%dotState(:,me) * zeta &
-                    source_dotState(1:size_so(so),2,so))
+                                    + source_dotState(1:size_so,1)* (1.0_pReal - zeta)
-      damageState(ph)%p(so)%dotState(:,me) = damageState(ph)%p(so)%dotState(:,me) * zeta &
+      r(1:size_so) = damageState(ph)%state    (1:size_so,me)  &
-                                        + source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta)
+                   - damageState(ph)%subState0(1:size_so,me)  &
-      r(1:size_so(so)) = damageState(ph)%p(so)%state    (1:size_so(so),me)  &
+                   - damageState(ph)%dotState (1:size_so,me) * dt
-                       - damageState(ph)%p(so)%subState0(1:size_so(so),me)  &
+      damageState(ph)%state(1:size_so,me) = damageState(ph)%state(1:size_so,me) - r(1:size_so)
-                       - damageState(ph)%p(so)%dotState (1:size_so(so),me) * dt
+      converged_ = converged_  .and. converged(r(1:size_so), &
-      damageState(ph)%p(so)%state(1:size_so(so),me) = damageState(ph)%p(so)%state(1:size_so(so),me) &
+                                               damageState(ph)%state(1:size_so,me), &
-                                                - r(1:size_so(so))
+                                               damageState(ph)%atol(1:size_so))
-      converged_ = converged_  .and. converged(r(1:size_so(so)), &
+
-                                               damageState(ph)%p(so)%state(1:size_so(so),me), &
-                                               damageState(ph)%p(so)%atol(1:size_so(so)))
-    enddo
 
     if(converged_) then
       broken = phase_damage_deltaState(mechanical_F_e(ph,me),ph,me)
@@ -366,10 +342,10 @@ module subroutine damage_results(group,ph)
 
   sourceLoop: do so = 1, phase_Nsources(ph)
 
-  if (phase_source(so,ph) /= DAMAGE_UNDEFINED_ID) &
+  if (phase_source(ph) /= DAMAGE_UNDEFINED_ID) &
     call results_closeGroup(results_addGroup(group//'sources/')) ! should be 'damage'
 
-    sourceType: select case (phase_source(so,ph))
+    sourceType: select case (phase_source(ph))
 
       case (DAMAGE_ISOBRITTLE_ID) sourceType
         call isobrittle_results(ph,group//'sources/')
@@ -393,39 +369,34 @@ end subroutine damage_results
 !--------------------------------------------------------------------------------------------------
 !> @brief contains the constitutive equation for calculating the rate of change of microstructure
 !--------------------------------------------------------------------------------------------------
-function phase_damage_collectDotState(co,ip,el,ph,me) result(broken)
+function phase_damage_collectDotState(ph,me) result(broken)
 
   integer, intent(in) :: &
-    co, &                                                                                          !< component-ID me integration point
-    ip, &                                                                                           !< integration point
-    el, &                                                                                           !< element
     ph, &
-    me
+    me                                                                                         !< counter in source loop
-  integer :: &
-    so                                                                                               !< counter in source loop
   logical :: broken
 
 
   broken = .false.
 
-  SourceLoop: do so = 1, phase_Nsources(ph)
+  if (damageState(ph)%sizeState > 0) then
 
-    sourceType: select case (phase_source(so,ph))
+    sourceType: select case (phase_source(ph))
 
       case (DAMAGE_ISODUCTILE_ID) sourceType
-        call isoductile_dotState(co, ip, el)
+        call isoductile_dotState(ph,me)
 
       case (DAMAGE_ANISODUCTILE_ID) sourceType
-        call anisoductile_dotState(co, ip, el)
+        call anisoductile_dotState(ph,me)
 
       case (DAMAGE_ANISOBRITTLE_ID) sourceType
-        call anisobrittle_dotState(mechanical_S(ph,me),co, ip, el) ! correct stress?
+        call anisobrittle_dotState(mechanical_S(ph,me), ph,me) ! correct stress?
 
     end select sourceType
 
-    broken = broken .or. any(IEEE_is_NaN(damageState(ph)%p(so)%dotState(:,me)))
+    broken = broken .or. any(IEEE_is_NaN(damageState(ph)%dotState(:,me)))
 
-  enddo SourceLoop
+  endif
 
 end function phase_damage_collectDotState
 
@@ -443,7 +414,6 @@ function phase_damage_deltaState(Fe, ph, me) result(broken)
   real(pReal),   intent(in), dimension(3,3) :: &
     Fe                                                                                              !< elastic deformation gradient
   integer :: &
-    so, &
     myOffset, &
     mySize
   logical :: &
@@ -452,23 +422,22 @@ function phase_damage_deltaState(Fe, ph, me) result(broken)
 
   broken = .false.
 
-  sourceLoop: do so = 1, phase_Nsources(ph)
+  if (damageState(ph)%sizeState == 0) return
 
-     sourceType: select case (phase_source(so,ph))
+     sourceType: select case (phase_source(ph))
 
       case (DAMAGE_ISOBRITTLE_ID) sourceType
-        call source_damage_isoBrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me)
+        call isobrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me)
-        broken = any(IEEE_is_NaN(damageState(ph)%p(so)%deltaState(:,me)))
+        broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,me)))
         if(.not. broken) then
-          myOffset = damageState(ph)%p(so)%offsetDeltaState
+          myOffset = damageState(ph)%offsetDeltaState
-          mySize   = damageState(ph)%p(so)%sizeDeltaState
+          mySize   = damageState(ph)%sizeDeltaState
-          damageState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,me) = &
+          damageState(ph)%state(myOffset + 1: myOffset + mySize,me) = &
-          damageState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,me) + damageState(ph)%p(so)%deltaState(1:mySize,me)
+          damageState(ph)%state(myOffset + 1: myOffset + mySize,me) + damageState(ph)%deltaState(1:mySize,me)
         endif
 
     end select sourceType
 
-  enddo SourceLoop
 
 end function phase_damage_deltaState
 
@@ -476,28 +445,25 @@ end function phase_damage_deltaState
 !--------------------------------------------------------------------------------------------------
 !> @brief checks if a source mechanism is active or not
 !--------------------------------------------------------------------------------------------------
-function source_active(source_label,src_length)  result(active_source)
+function source_active(source_label)  result(active_source)
 
   character(len=*), intent(in)         :: source_label                                              !< name of source mechanism
-  integer,          intent(in)         :: src_length                                                !< max. number of sources in system
+  logical, dimension(:), allocatable  :: active_source
-  logical, dimension(:,:), allocatable :: active_source
 
   class(tNode), pointer :: &
     phases, &
     phase, &
     sources, &
     src
-  integer :: p,s
+  integer :: ph
 
   phases => config_material%get('phase')
-  allocate(active_source(src_length,phases%length), source = .false. )
+  allocate(active_source(phases%length))
-  do p = 1, phases%length
+  do ph = 1, phases%length
-    phase => phases%get(p)
+    phase => phases%get(ph)
-    sources => phase%get('source',defaultVal=emptyList)
+    sources => phase%get('damage',defaultVal=emptyList)
-    do s = 1, sources%length
+    src => sources%get(1)
-      src => sources%get(s)
+    active_source(ph) = src%get_asString('type',defaultVal = 'x') == source_label
-      if(src%get_asString('type') == source_label) active_source(s,p) = .true.
-    enddo
   enddo
 
 
@@ -528,4 +494,15 @@ module function phase_damage_get_phi(co,ip,el) result(phi)
 end function phase_damage_get_phi
 
 
+module function damage_phi(ph,me) result(phi)
+
+  integer, intent(in) :: ph, me
+  real(pReal) :: phi
+
+
+  phi = current(ph)%phi(me)
+
+end function damage_phi
+
+
 end submodule damagee

src/phase_damage_anisobrittle.f90

@@ -6,10 +6,6 @@
 !--------------------------------------------------------------------------------------------------
 submodule (phase:damagee) anisobrittle
 
-  integer,                       dimension(:),           allocatable :: &
-    source_damage_anisoBrittle_offset, &                                                            !< which source is my current source mechanism?
-    source_damage_anisoBrittle_instance                                                             !< instance of source mechanism
-
   type :: tParameters                                                                               !< container type for internal constitutive parameters
     real(pReal) :: &
       dot_o, &                                                                                      !< opening rate of cleavage planes
@@ -35,42 +31,38 @@ contains
 !> @brief module initialization
 !> @details reads in material parameters, allocates arrays, and does sanity checks
 !--------------------------------------------------------------------------------------------------
-module function anisobrittle_init(source_length) result(mySources)
+module function anisobrittle_init() result(mySources)
 
-  integer, intent(in)                  :: source_length
+  logical, dimension(:), allocatable :: mySources
-  logical, dimension(:,:), allocatable :: mySources
 
   class(tNode), pointer :: &
     phases, &
     phase, &
     sources, &
     src
-  integer :: Ninstances,sourceOffset,Nconstituents,p
+  integer :: Nconstituents,p
   integer, dimension(:), allocatable :: N_cl
   character(len=pStringLen) :: extmsg = ''
 
-  print'(/,a)', ' <<<+-  phase:damage:anisobrittle init  -+>>>'
 
-  mySources = source_active('damage_anisoBrittle',source_length)
+  mySources = source_active('anisobrittle')
-  Ninstances = count(mySources)
+  if(count(mySources) == 0) return
-  print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
+
-  if(Ninstances == 0) return
+  print'(/,a)', ' <<<+-  phase:damage:anisobrittle init  -+>>>'
+  print'(a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT)
+
 
   phases => config_material%get('phase')
-  allocate(param(Ninstances))
+  allocate(param(phases%length))
-  allocate(source_damage_anisoBrittle_offset  (phases%length), source=0)
+
-  allocate(source_damage_anisoBrittle_instance(phases%length), source=0)
 
   do p = 1, phases%length
+    if(mySources(p)) then
     phase => phases%get(p)
-    if(any(mySources(:,p))) source_damage_anisoBrittle_instance(p) = count(mySources(:,1:p))
+    sources => phase%get('damage')
-    if(count(mySources(:,p)) == 0) cycle
+
-    sources => phase%get('source')
+        associate(prm  => param(p))
-    do sourceOffset = 1, sources%length
+        src => sources%get(1)
-      if(mySources(sourceOffset,p)) then
-        source_damage_anisoBrittle_offset(p) = sourceOffset
-        associate(prm  => param(source_damage_anisoBrittle_instance(p)))
-        src => sources%get(sourceOffset)
 
         N_cl = src%get_asInts('N_cl',defaultVal=emptyIntArray)
         prm%sum_N_cl = sum(abs(N_cl))
@@ -101,9 +93,9 @@ module function anisobrittle_init(source_length) result(mySources)
         if (any(prm%s_crit <  0.0_pReal)) extmsg = trim(extmsg)//' s_crit'
 
         Nconstituents = count(material_phaseAt==p) * discretization_nIPs
-        call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0)
+        call phase_allocateState(damageState(p),Nconstituents,1,1,0)
-        damageState(p)%p(sourceOffset)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal)
+        damageState(p)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal)
-        if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol'
+        if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol'
 
         end associate
 
@@ -111,7 +103,7 @@ module function anisobrittle_init(source_length) result(mySources)
 !  exit if any parameter is out of range
         if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_anisoBrittle)')
       endif
-    enddo
+
   enddo
 
 end function anisobrittle_init
@@ -120,18 +112,14 @@ end function anisobrittle_init
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates derived quantities from state
 !--------------------------------------------------------------------------------------------------
-module subroutine anisobrittle_dotState(S, co, ip, el)
+module subroutine anisobrittle_dotState(S, ph,me)
 
   integer, intent(in) :: &
-    co, &                                                                                          !< component-ID of integration point
+    ph,me
-    ip, &                                                                                           !< integration point
-    el                                                                                              !< element
   real(pReal),  intent(in), dimension(3,3) :: &
     S
 
   integer :: &
-    ph, &
-    me, &
     sourceOffset, &
     damageOffset, &
     homog, &
@@ -139,28 +127,22 @@ module subroutine anisobrittle_dotState(S, co, ip, el)
   real(pReal) :: &
     traction_d, traction_t, traction_n, traction_crit
 
-  ph = material_phaseAt(co,el)
-  me = material_phasememberAt(co,ip,el)
-  sourceOffset = source_damage_anisoBrittle_offset(ph)
-  homog = material_homogenizationAt(el)
-  damageOffset = material_homogenizationMemberAt(ip,el)
 
-  associate(prm => param(source_damage_anisoBrittle_instance(ph)))
+  associate(prm => param(ph))
-  damageState(ph)%p(sourceOffset)%dotState(1,me) = 0.0_pReal
+    damageState(ph)%dotState(1,me) = 0.0_pReal
-  do i = 1, prm%sum_N_cl
+    do i = 1, prm%sum_N_cl
-    traction_d    = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i))
+      traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i))
-    traction_t    = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i))
+      traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i))
-    traction_n    = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i))
+      traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i))
 
-    traction_crit = prm%g_crit(i)*damage(homog)%p(damageOffset)**2.0_pReal
+      traction_crit = prm%g_crit(i)*damage_phi(ph,me)**2.0_pReal
 
-    damageState(ph)%p(sourceOffset)%dotState(1,me) &
+      damageState(ph)%dotState(1,me) = damageState(ph)%dotState(1,me) &
-    = damageState(ph)%p(sourceOffset)%dotState(1,me) &
+          + prm%dot_o / prm%s_crit(i) &
-    + prm%dot_o / prm%s_crit(i) &
+            * ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%q + &
-      * ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%q + &
+               (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%q + &
-         (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%q + &
+               (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**prm%q)
-         (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**prm%q)
+    enddo
-  enddo
   end associate
 
 end subroutine anisobrittle_dotState
@@ -169,28 +151,24 @@ end subroutine anisobrittle_dotState
 !--------------------------------------------------------------------------------------------------
 !> @brief returns local part of nonlocal damage driving force
 !--------------------------------------------------------------------------------------------------
-module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
+module subroutine anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me)
 
   integer, intent(in) :: &
-    phase, &
+    ph, &
-    constituent
+    me
   real(pReal),  intent(in) :: &
     phi
   real(pReal),  intent(out) :: &
     localphiDot, &
     dLocalphiDot_dPhi
 
-  integer :: &
-    sourceOffset
 
-  sourceOffset = source_damage_anisoBrittle_offset(phase)
+  dLocalphiDot_dPhi = -damageState(ph)%state(1,me)
-
-  dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent)
 
   localphiDot = 1.0_pReal &
               + dLocalphiDot_dPhi*phi
 
-end subroutine source_damage_anisoBrittle_getRateAndItsTangent
+end subroutine anisobrittle_getRateAndItsTangent
 
 
 !--------------------------------------------------------------------------------------------------
@@ -203,16 +181,78 @@ module subroutine anisobrittle_results(phase,group)
 
   integer :: o
 
-  associate(prm => param(source_damage_anisoBrittle_instance(phase)), &
+
-            stt => damageState(phase)%p(source_damage_anisoBrittle_offset(phase))%state)
+  associate(prm => param(phase), stt => damageState(phase)%state)
-  outputsLoop: do o = 1,size(prm%output)
+    outputsLoop: do o = 1,size(prm%output)
-    select case(trim(prm%output(o)))
+      select case(trim(prm%output(o)))
-      case ('f_phi')
+        case ('f_phi')
-        call results_writeDataset(group,stt,trim(prm%output(o)),'driving force','J/m³')
+          call results_writeDataset(group,stt,trim(prm%output(o)),'driving force','J/m³')
-    end select
+      end select
-  enddo outputsLoop
+    enddo outputsLoop
   end associate
 
 end subroutine anisobrittle_results
 
+
+!--------------------------------------------------------------------------------------------------
+!> @brief  contains the constitutive equation for calculating the velocity gradient
+!--------------------------------------------------------------------------------------------------
+module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me)
+
+  integer, intent(in) :: &
+    ph,me
+  real(pReal),   intent(in),  dimension(3,3) :: &
+    S
+  real(pReal),   intent(out), dimension(3,3) :: &
+    Ld                                                                                              !< damage velocity gradient
+  real(pReal),   intent(out), dimension(3,3,3,3) :: &
+    dLd_dTstar                                                                                      !< derivative of Ld with respect to Tstar (4th-order tensor)
+
+  integer :: &
+    i, k, l, m, n
+  real(pReal) :: &
+    traction_d, traction_t, traction_n, traction_crit, &
+    udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt
+
+
+  Ld = 0.0_pReal
+  dLd_dTstar = 0.0_pReal
+  associate(prm => param(ph))
+  do i = 1,prm%sum_N_cl
+    traction_crit = prm%g_crit(i)*damage_phi(ph,me)**2.0_pReal
+
+    traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i))
+    if (abs(traction_d) > traction_crit + tol_math_check) then
+      udotd = sign(1.0_pReal,traction_d)* prm%dot_o * ((abs(traction_d) - traction_crit)/traction_crit)**prm%q
+      Ld = Ld + udotd*prm%cleavage_systems(1:3,1:3,1,i)
+      dudotd_dt = sign(1.0_pReal,traction_d)*udotd*prm%q / (abs(traction_d) - traction_crit)
+      forall (k=1:3,l=1:3,m=1:3,n=1:3) &
+        dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) &
+                            + dudotd_dt*prm%cleavage_systems(k,l,1,i) * prm%cleavage_systems(m,n,1,i)
+    endif
+
+    traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i))
+    if (abs(traction_t) > traction_crit + tol_math_check) then
+      udott = sign(1.0_pReal,traction_t)* prm%dot_o * ((abs(traction_t) - traction_crit)/traction_crit)**prm%q
+      Ld = Ld + udott*prm%cleavage_systems(1:3,1:3,2,i)
+      dudott_dt = sign(1.0_pReal,traction_t)*udott*prm%q / (abs(traction_t) - traction_crit)
+      forall (k=1:3,l=1:3,m=1:3,n=1:3) &
+        dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) &
+                            + dudott_dt*prm%cleavage_systems(k,l,2,i) * prm%cleavage_systems(m,n,2,i)
+    endif
+
+    traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i))
+    if (abs(traction_n) > traction_crit + tol_math_check) then
+      udotn = sign(1.0_pReal,traction_n)* prm%dot_o * ((abs(traction_n) - traction_crit)/traction_crit)**prm%q
+      Ld = Ld + udotn*prm%cleavage_systems(1:3,1:3,3,i)
+      dudotn_dt = sign(1.0_pReal,traction_n)*udotn*prm%q / (abs(traction_n) - traction_crit)
+      forall (k=1:3,l=1:3,m=1:3,n=1:3) &
+        dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) &
+                            + dudotn_dt*prm%cleavage_systems(k,l,3,i) * prm%cleavage_systems(m,n,3,i)
+    endif
+  enddo
+  end associate
+
+end subroutine kinematics_cleavage_opening_LiAndItsTangent
+
 end submodule anisobrittle

src/phase_damage_anisoductile.f90

@@ -6,10 +6,6 @@
 !--------------------------------------------------------------------------------------------------
 submodule(phase:damagee) anisoductile
 
-  integer,                       dimension(:),           allocatable :: &
-    source_damage_anisoDuctile_offset, &                                                            !< which source is my current damage mechanism?
-    source_damage_anisoDuctile_instance                                                             !< instance of damage source mechanism
-
   type :: tParameters                                                                               !< container type for internal constitutive parameters
     real(pReal) :: &
       q                                                                                             !< damage rate sensitivity
@@ -19,7 +15,7 @@ submodule(phase:damagee) anisoductile
       output
   end type tParameters
 
-  type(tParameters), dimension(:), allocatable :: param                                             !< containers of constitutive parameters (len Ninstances)
+  type(tParameters), dimension(:), allocatable :: param                                             !< containers of constitutive parameters
 
 contains
 
@@ -28,10 +24,9 @@ contains
 !> @brief module initialization
 !> @details reads in material parameters, allocates arrays, and does sanity checks
 !--------------------------------------------------------------------------------------------------
-module function anisoductile_init(source_length) result(mySources)
+module function anisoductile_init() result(mySources)
 
-  integer, intent(in)                  :: source_length
+  logical, dimension(:), allocatable :: mySources
-  logical, dimension(:,:), allocatable :: mySources
 
   class(tNode), pointer :: &
     phases, &
@@ -40,34 +35,31 @@ module function anisoductile_init(source_length) result(mySources)
     pl, &
     sources, &
     src
-  integer :: Ninstances,sourceOffset,Nconstituents,p
+  integer :: Ninstances,Nconstituents,p
   integer, dimension(:), allocatable :: N_sl
   character(len=pStringLen) :: extmsg = ''
 
-  print'(/,a)', ' <<<+-  phase:damage:anisoductile init  -+>>>'
 
-  mySources = source_active('damage_anisoDuctile',source_length)
+  mySources = source_active('anisoductile')
-  Ninstances = count(mySources)
+  if(count(mySources) == 0) return
-  print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
+
-  if(Ninstances == 0) return
+  print'(/,a)', ' <<<+-  phase:damage:anisoductile init  -+>>>'
+  print'(a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT)
+
 
   phases => config_material%get('phase')
-  allocate(param(Ninstances))
+  allocate(param(phases%length))
-  allocate(source_damage_anisoDuctile_offset  (phases%length), source=0)
-  allocate(source_damage_anisoDuctile_instance(phases%length), source=0)
 
   do p = 1, phases%length
-    phase => phases%get(p)
+    if(mySources(p)) then
-    if(any(mySources(:,p))) source_damage_anisoDuctile_instance(p) = count(mySources(:,1:p))
+      phase => phases%get(p)
-    if(count(mySources(:,p)) == 0) cycle
+      mech  => phase%get('mechanics')
-    mech  => phase%get('mechanics')
+      pl    => mech%get('plasticity')
-    pl    => mech%get('plasticity')
+      sources => phase%get('damage')
-    sources => phase%get('source')
+
-    do sourceOffset = 1, sources%length
+
-      if(mySources(sourceOffset,p)) then
+        associate(prm  => param(p))
-        source_damage_anisoDuctile_offset(p) = sourceOffset
+        src => sources%get(1)
-        associate(prm  => param(source_damage_anisoDuctile_instance(p)))
-        src => sources%get(sourceOffset)
 
         N_sl           = pl%get_asInts('N_sl',defaultVal=emptyIntArray)
         prm%q          = src%get_asFloat('q')
@@ -86,10 +78,10 @@ module function anisoductile_init(source_length) result(mySources)
         if (prm%q              <= 0.0_pReal)  extmsg = trim(extmsg)//' q'
         if (any(prm%gamma_crit <  0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit'
 
-        Nconstituents=count(material_phaseAt==p) * discretization_nIPs
+        Nconstituents=count(material_phaseAt2==p)
-        call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0)
+        call phase_allocateState(damageState(p),Nconstituents,1,1,0)
-        damageState(p)%p(sourceOffset)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal)
+        damageState(p)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal)
-        if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol'
+        if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol'
 
         end associate
 
@@ -97,7 +89,7 @@ module function anisoductile_init(source_length) result(mySources)
 !  exit if any parameter is out of range
         if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_anisoDuctile)')
       endif
-    enddo
+
   enddo
 
 
@@ -107,29 +99,15 @@ end function anisoductile_init
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates derived quantities from state
 !--------------------------------------------------------------------------------------------------
-module subroutine anisoductile_dotState(co, ip, el)
+module subroutine anisoductile_dotState(ph,me)
 
   integer, intent(in) :: &
-    co, &                                                                                          !< component-ID of integration point
-    ip, &                                                                                           !< integration point
-    el                                                                                              !< element
-
-  integer :: &
     ph, &
-    me, &
+    me
-    sourceOffset, &
-    damageOffset, &
-    homog
 
-  ph = material_phaseAt(co,el)
-  me = material_phasememberAt(co,ip,el)
-  sourceOffset = source_damage_anisoDuctile_offset(ph)
-  homog = material_homogenizationAt(el)
-  damageOffset = material_homogenizationMemberAt(ip,el)
 
-  associate(prm => param(source_damage_anisoDuctile_instance(ph)))
+  associate(prm => param(ph))
-  damageState(ph)%p(sourceOffset)%dotState(1,me) &
+    damageState(ph)%dotState(1,me) = sum(plasticState(ph)%slipRate(:,me)/(damage_phi(ph,me)**prm%q)/prm%gamma_crit)
-    = sum(plasticState(ph)%slipRate(:,me)/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit)
   end associate
 
 end subroutine anisoductile_dotState
@@ -138,28 +116,24 @@ end subroutine anisoductile_dotState
 !--------------------------------------------------------------------------------------------------
 !> @brief returns local part of nonlocal damage driving force
 !--------------------------------------------------------------------------------------------------
-module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
+module subroutine anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me)
 
   integer, intent(in) :: &
-    phase, &
+    ph, &
-    constituent
+    me
   real(pReal),  intent(in) :: &
     phi
   real(pReal),  intent(out) :: &
     localphiDot, &
     dLocalphiDot_dPhi
 
-  integer :: &
-    sourceOffset
 
-  sourceOffset = source_damage_anisoDuctile_offset(phase)
+  dLocalphiDot_dPhi = -damageState(ph)%state(1,me)
-
-  dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent)
 
   localphiDot = 1.0_pReal &
               + dLocalphiDot_dPhi*phi
 
-end subroutine source_damage_anisoDuctile_getRateAndItsTangent
+end subroutine anisoductile_getRateAndItsTangent
 
 
 !--------------------------------------------------------------------------------------------------
@@ -172,14 +146,14 @@ module subroutine anisoductile_results(phase,group)
 
   integer :: o
 
-  associate(prm => param(source_damage_anisoDuctile_instance(phase)), &
+
-            stt => damageState(phase)%p(source_damage_anisoDuctile_offset(phase))%state)
+  associate(prm => param(phase), stt => damageState(phase)%state)
-  outputsLoop: do o = 1,size(prm%output)
+    outputsLoop: do o = 1,size(prm%output)
-    select case(trim(prm%output(o)))
+      select case(trim(prm%output(o)))
-      case ('f_phi')
+        case ('f_phi')
-        call results_writeDataset(group,stt,trim(prm%output(o)),'driving force','J/m³')
+          call results_writeDataset(group,stt,trim(prm%output(o)),'driving force','J/m³')
-    end select
+      end select
-  enddo outputsLoop
+    enddo outputsLoop
   end associate
 
 end subroutine anisoductile_results

src/phase_damage_isobrittle.f90

@@ -6,10 +6,6 @@
 !--------------------------------------------------------------------------------------------------
 submodule(phase:damagee) isobrittle
 
-  integer,                       dimension(:),           allocatable :: &
-    source_damage_isoBrittle_offset, &
-    source_damage_isoBrittle_instance
-
   type :: tParameters                                                                               !< container type for internal constitutive parameters
     real(pReal) :: &
       W_crit                                                                                        !< critical elastic strain energy
@@ -26,41 +22,36 @@ contains
 !> @brief module initialization
 !> @details reads in material parameters, allocates arrays, and does sanity checks
 !--------------------------------------------------------------------------------------------------
-module function isobrittle_init(source_length) result(mySources)
+module function isobrittle_init() result(mySources)
 
-  integer, intent(in)                  :: source_length
+  logical, dimension(:), allocatable :: mySources
-  logical, dimension(:,:), allocatable :: mySources
 
   class(tNode), pointer :: &
     phases, &
     phase, &
     sources, &
     src
-  integer :: Ninstances,sourceOffset,Nconstituents,p
+  integer :: Nconstituents,p
   character(len=pStringLen) :: extmsg = ''
 
-  print'(/,a)', ' <<<+-  phase:damage:isobrittle init  -+>>>'
 
-  mySources = source_active('damage_isoBrittle',source_length)
+  mySources = source_active('isobrittle')
-  Ninstances = count(mySources)
+  if(count(mySources) == 0) return
-  print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
+
-  if(Ninstances == 0) return
+  print'(/,a)', ' <<<+-  phase:damage:isobrittle init  -+>>>'
+  print'(a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT)
+
 
   phases => config_material%get('phase')
-  allocate(param(Ninstances))
+  allocate(param(phases%length))
-  allocate(source_damage_isoBrittle_offset  (phases%length), source=0)
-  allocate(source_damage_isoBrittle_instance(phases%length), source=0)
 
   do p = 1, phases%length
+    if(mySources(p)) then
     phase => phases%get(p)
-    if(any(mySources(:,p))) source_damage_isoBrittle_instance(p) = count(mySources(:,1:p))
+    sources => phase%get('damage')
-    if(count(mySources(:,p)) == 0) cycle
+
-    sources => phase%get('source')
+        associate(prm  => param(p))
-    do sourceOffset = 1, sources%length
+        src => sources%get(1)
-      if(mySources(sourceOffset,p)) then
-        source_damage_isoBrittle_offset(p) = sourceOffset
-        associate(prm  => param(source_damage_isoBrittle_instance(p)))
-        src => sources%get(sourceOffset)
 
         prm%W_crit = src%get_asFloat('W_crit')
 
@@ -73,10 +64,10 @@ module function isobrittle_init(source_length) result(mySources)
         ! sanity checks
         if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit'
 
-        Nconstituents = count(material_phaseAt==p) * discretization_nIPs
+        Nconstituents = count(material_phaseAt2==p)
-        call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,1)
+        call phase_allocateState(damageState(p),Nconstituents,1,1,1)
-        damageState(p)%p(sourceOffset)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal)
+        damageState(p)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal)
-        if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol'
+        if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol'
 
         end associate
 
@@ -84,7 +75,7 @@ module function isobrittle_init(source_length) result(mySources)
 !  exit if any parameter is out of range
         if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isoBrittle)')
       endif
-    enddo
+
   enddo
 
 
@@ -94,7 +85,7 @@ end function isobrittle_init
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates derived quantities from state
 !--------------------------------------------------------------------------------------------------
-module subroutine source_damage_isoBrittle_deltaState(C, Fe, ph,me)
+module subroutine isobrittle_deltaState(C, Fe, ph,me)
 
   integer, intent(in) :: ph,me
   real(pReal),  intent(in), dimension(3,3) :: &
@@ -102,60 +93,47 @@ module subroutine source_damage_isoBrittle_deltaState(C, Fe, ph,me)
   real(pReal),  intent(in), dimension(6,6) :: &
     C
 
-  integer :: &
-    sourceOffset
   real(pReal), dimension(6) :: &
     strain
   real(pReal) :: &
     strainenergy
 
-  sourceOffset = source_damage_isoBrittle_offset(ph)
 
   strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3)
 
-  associate(prm => param(source_damage_isoBrittle_instance(ph)))
+  associate(prm => param(ph))
-  strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/prm%W_crit
+    strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/prm%W_crit
-  ! ToDo: check strainenergy = 2.0_pReal*dot_product(strain,matmul(C,strain))/prm%W_crit
+    ! ToDo: check strainenergy = 2.0_pReal*dot_product(strain,matmul(C,strain))/prm%W_crit
 
-  if (strainenergy > damageState(ph)%p(sourceOffset)%subState0(1,me)) then
+    damageState(ph)%deltaState(1,me) = merge(strainenergy - damageState(ph)%state(1,me), &
-    damageState(ph)%p(sourceOffset)%deltaState(1,me) = &
+                                             damageState(ph)%subState0(1,me) - damageState(ph)%state(1,me), &
-      strainenergy - damageState(ph)%p(sourceOffset)%state(1,me)
+                                             strainenergy > damageState(ph)%subState0(1,me))
-  else
-    damageState(ph)%p(sourceOffset)%deltaState(1,me) = &
-      damageState(ph)%p(sourceOffset)%subState0(1,me) - &
-      damageState(ph)%p(sourceOffset)%state(1,me)
-  endif
   end associate
 
-end subroutine source_damage_isoBrittle_deltaState
+end subroutine isobrittle_deltaState
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief returns local part of nonlocal damage driving force
 !--------------------------------------------------------------------------------------------------
-module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
+module subroutine isobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi,  ph, me)
 
   integer, intent(in) :: &
-    phase, &
+    ph, me
-    constituent
   real(pReal),  intent(in) :: &
     phi
   real(pReal),  intent(out) :: &
     localphiDot, &
     dLocalphiDot_dPhi
 
-  integer :: &
-    sourceOffset
 
-  sourceOffset = source_damage_isoBrittle_offset(phase)
+  associate(prm => param(ph))
-
+    localphiDot = 1.0_pReal &
-  associate(prm => param(source_damage_isoBrittle_instance(phase)))
+                - phi*damageState(ph)%state(1,me)
-  localphiDot = 1.0_pReal &
+    dLocalphiDot_dPhi = - damageState(ph)%state(1,me)
-              - phi*damageState(phase)%p(sourceOffset)%state(1,constituent)
-  dLocalphiDot_dPhi = - damageState(phase)%p(sourceOffset)%state(1,constituent)
   end associate
 
-end subroutine source_damage_isoBrittle_getRateAndItsTangent
+end subroutine isobrittle_getRateAndItsTangent
 
 
 !--------------------------------------------------------------------------------------------------
@@ -168,14 +146,15 @@ module subroutine isobrittle_results(phase,group)
 
   integer :: o
 
-  associate(prm => param(source_damage_isoBrittle_instance(phase)), &
+
-            stt => damageState(phase)%p(source_damage_isoBrittle_offset(phase))%state)
+  associate(prm => param(phase), &
-  outputsLoop: do o = 1,size(prm%output)
+            stt => damageState(phase)%state)
-    select case(trim(prm%output(o)))
+    outputsLoop: do o = 1,size(prm%output)
-      case ('f_phi')
+      select case(trim(prm%output(o)))
-        call results_writeDataset(group,stt,trim(prm%output(o)),'driving force','J/m³')
+        case ('f_phi')
-    end select
+          call results_writeDataset(group,stt,trim(prm%output(o)),'driving force','J/m³')
-  enddo outputsLoop
+      end select
+    enddo outputsLoop
   end associate
 
 end subroutine isobrittle_results

src/phase_damage_isoductile.f90

@@ -6,10 +6,6 @@
 !--------------------------------------------------------------------------------------------------
 submodule(phase:damagee) isoductile
 
-  integer,                       dimension(:),           allocatable :: &
-    source_damage_isoDuctile_offset, &                                                              !< which source is my current damage mechanism?
-    source_damage_isoDuctile_instance                                                               !< instance of damage source mechanism
-
   type:: tParameters                                                                                !< container type for internal constitutive parameters
     real(pReal) :: &
       gamma_crit, &                                                                                 !< critical plastic strain
@@ -28,41 +24,36 @@ contains
 !> @brief module initialization
 !> @details reads in material parameters, allocates arrays, and does sanity checks
 !--------------------------------------------------------------------------------------------------
-module function isoductile_init(source_length) result(mySources)
+module function isoductile_init() result(mySources)
 
-  integer, intent(in)                  :: source_length
+  logical, dimension(:), allocatable :: mySources
-  logical, dimension(:,:), allocatable :: mySources
 
   class(tNode), pointer :: &
     phases, &
     phase, &
     sources, &
     src
-  integer :: Ninstances,sourceOffset,Nconstituents,p
+  integer :: Ninstances,Nconstituents,p
   character(len=pStringLen) :: extmsg = ''
 
-  print'(/,a)', ' <<<+-  phase:damage:isoductile init  -+>>>'
 
-  mySources = source_active('damage_isoDuctile',source_length)
+  mySources = source_active('isoductile')
-  Ninstances = count(mySources)
+  if(count(mySources) == 0) return
-  print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
+
-  if(Ninstances == 0) return
+  print'(/,a)', ' <<<+-  phase:damage:isoductile init  -+>>>'
+  print'(a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT)
+
 
   phases => config_material%get('phase')
-  allocate(param(Ninstances))
+  allocate(param(phases%length))
-  allocate(source_damage_isoDuctile_offset  (phases%length), source=0)
-  allocate(source_damage_isoDuctile_instance(phases%length), source=0)
 
   do p = 1, phases%length
-    phase => phases%get(p)
+    if(mySources(p)) then
-    if(count(mySources(:,p)) == 0) cycle
+      phase => phases%get(p)
-    if(any(mySources(:,p))) source_damage_isoDuctile_instance(p) = count(mySources(:,1:p))
+      sources => phase%get('damage')
-    sources => phase%get('source')
+
-    do sourceOffset = 1, sources%length
+        associate(prm  => param(p))
-      if(mySources(sourceOffset,p)) then
+        src => sources%get(1)
-        source_damage_isoDuctile_offset(p) = sourceOffset
-        associate(prm  => param(source_damage_isoDuctile_instance(p)))
-        src => sources%get(sourceOffset)
 
         prm%q          = src%get_asFloat('q')
         prm%gamma_crit = src%get_asFloat('gamma_crit')
@@ -77,10 +68,10 @@ module function isoductile_init(source_length) result(mySources)
         if (prm%q          <= 0.0_pReal) extmsg = trim(extmsg)//' q'
         if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit'
 
-        Nconstituents=count(material_phaseAt==p) * discretization_nIPs
+        Nconstituents=count(material_phaseAt2==p)
-        call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0)
+        call phase_allocateState(damageState(p),Nconstituents,1,1,0)
-        damageState(p)%p(sourceOffset)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal)
+        damageState(p)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal)
-        if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol'
+        if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol'
 
         end associate
 
@@ -88,7 +79,6 @@ module function isoductile_init(source_length) result(mySources)
 !  exit if any parameter is out of range
         if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isoDuctile)')
       endif
-    enddo
   enddo
 
 
@@ -98,29 +88,16 @@ end function isoductile_init
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates derived quantities from state
 !--------------------------------------------------------------------------------------------------
-module subroutine isoductile_dotState(co, ip, el)
+module subroutine isoductile_dotState(ph, me)
 
   integer, intent(in) :: &
-    co, &                                                                                          !< component-ID of integration point
-    ip, &                                                                                           !< integration point
-    el                                                                                              !< element
-
-  integer :: &
     ph, &
-    me, &
+    me
-    sourceOffset, &
-    damageOffset, &
-    homog
 
-  ph = material_phaseAt(co,el)
-  me = material_phasememberAt(co,ip,el)
-  sourceOffset = source_damage_isoDuctile_offset(ph)
-  homog = material_homogenizationAt(el)
-  damageOffset = material_homogenizationMemberAt(ip,el)
 
-  associate(prm => param(source_damage_isoDuctile_instance(ph)))
+  associate(prm => param(ph))
-  damageState(ph)%p(sourceOffset)%dotState(1,me) = &
+    damageState(ph)%dotState(1,me) = sum(plasticState(ph)%slipRate(:,me)) &
-    sum(plasticState(ph)%slipRate(:,me))/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit
+                                   / (prm%gamma_crit*damage_phi(ph,me)**prm%q)
   end associate
 
 end subroutine isoductile_dotState
@@ -129,28 +106,24 @@ end subroutine isoductile_dotState
 !--------------------------------------------------------------------------------------------------
 !> @brief returns local part of nonlocal damage driving force
 !--------------------------------------------------------------------------------------------------
-module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
+module subroutine isoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me)
 
   integer, intent(in) :: &
-    phase, &
+    ph, &
-    constituent
+    me
   real(pReal),  intent(in) :: &
     phi
   real(pReal),  intent(out) :: &
     localphiDot, &
     dLocalphiDot_dPhi
 
-  integer :: &
-    sourceOffset
 
-  sourceOffset = source_damage_isoDuctile_offset(phase)
+  dLocalphiDot_dPhi = -damageState(ph)%state(1,me)
-
-  dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent)
 
   localphiDot = 1.0_pReal &
               + dLocalphiDot_dPhi*phi
 
-end subroutine source_damage_isoDuctile_getRateAndItsTangent
+end subroutine isoductile_getRateAndItsTangent
 
 
 !--------------------------------------------------------------------------------------------------
@@ -163,14 +136,13 @@ module subroutine isoductile_results(phase,group)
 
   integer :: o
 
-  associate(prm => param(source_damage_isoDuctile_instance(phase)), &
+  associate(prm => param(phase), stt => damageState(phase)%state)
-            stt => damageState(phase)%p(source_damage_isoDuctile_offset(phase))%state)
+    outputsLoop: do o = 1,size(prm%output)
-  outputsLoop: do o = 1,size(prm%output)
+      select case(trim(prm%output(o)))
-    select case(trim(prm%output(o)))
+        case ('f_phi')
-      case ('f_phi')
+          call results_writeDataset(group,stt,trim(prm%output(o)),'driving force','J/m³')
-        call results_writeDataset(group,stt,trim(prm%output(o)),'driving force','J/m³')
+      end select
-    end select
+    enddo outputsLoop
-  enddo outputsLoop
   end associate
 
 end subroutine isoductile_results

src/phase_mechanical.f90

@@ -1,7 +1,8 @@
 !----------------------------------------------------------------------------------------------------
 !> @brief internal microstructure state for all plasticity constitutive models
 !----------------------------------------------------------------------------------------------------
-submodule(phase) mechanics
+submodule(phase) mechanical
+
 
   enum, bind(c); enumerator :: &
     ELASTICITY_UNDEFINED_ID, &
@@ -22,8 +23,6 @@ submodule(phase) mechanics
     KINEMATICS_THERMAL_EXPANSION_ID
   end enum
 
-  integer(kind(KINEMATICS_UNDEFINED_ID)),     dimension(:,:), allocatable :: &
-    phase_kinematics
   integer(kind(ELASTICITY_UNDEFINED_ID)), dimension(:),   allocatable :: &
     phase_elasticity                                                                                !< elasticity of each phase
   integer(kind(STIFFNESS_DEGRADATION_UNDEFINED_ID)),     dimension(:,:), allocatable :: &
@@ -98,12 +97,9 @@ submodule(phase) mechanics
     end function plastic_deltaState
 
     module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
-                                             S, Fi, co, ip, el)
+                                             S, Fi, ph,me)
-
       integer, intent(in) :: &
-        co, &                                                                                          !< component-ID of integration point
+        ph,me
-        ip, &                                                                                           !< integration point
-        el                                                                                              !< element
       real(pReal),   intent(in),  dimension(3,3) :: &
         S                                                                                               !< 2nd Piola-Kirchhoff stress
       real(pReal),   intent(in),  dimension(3,3) :: &
@@ -206,7 +202,7 @@ module subroutine mechanical_init(phases)
     elastic, &
     stiffDegradation
 
-  print'(/,a)', ' <<<+-  phase:mechanics init  -+>>>'
+  print'(/,a)', ' <<<+-  phase:mechanical init  -+>>>'
 
 !-------------------------------------------------------------------------------------------------
 ! initialize elasticity (hooke)                         !ToDO: Maybe move to elastic submodule along with function homogenizedC?
@@ -381,7 +377,7 @@ subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
   DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phaseAt(co,el))
     degradationType: select case(phase_stiffnessDegradation(d,material_phaseAt(co,el)))
       case (STIFFNESS_DEGRADATION_damage_ID) degradationType
-        C = C * damage(ho)%p(material_homogenizationMemberAt(ip,el))**2
+        C = C * phase_damage_get_phi(co,ip,el)**2
     end select degradationType
   enddo DegradationLoop
 
@@ -583,7 +579,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
     enddo LpLoop
 
     call phase_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, &
-                                       S, Fi_new, co, ip, el)
+                                       S, Fi_new, ph,me)
 
     !* update current residuum and check for convergence of loop
     atol_Li = max(num%rtol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), &        ! absolute tolerance from largest acceptable relative error
@@ -1141,7 +1137,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
   real(pReal) :: &
     formerSubStep
   integer :: &
-    so, ph, me, sizeDotState
+    ph, me, sizeDotState
   logical :: todo
   real(pReal) :: subFrac,subStep
   real(pReal), dimension(3,3) :: &
@@ -1162,10 +1158,9 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
   subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,me)
   subState0 = plasticState(ph)%State0(:,me)
 
+  if (damageState(ph)%sizeState > 0) &
+    damageState(ph)%subState0(:,me) = damageState(ph)%state0(:,me)
 
-  do so = 1, phase_Nsources(ph)
-    damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state0(:,me)
-  enddo
   subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,me)
   subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,me)
   subF0  = phase_mechanical_F0(ph)%data(1:3,1:3,me)
@@ -1191,9 +1186,9 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
         subFp0 = phase_mechanical_Fp(ph)%data(1:3,1:3,me)
         subFi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,me)
         subState0 = plasticState(ph)%state(:,me)
-        do so = 1, phase_Nsources(ph)
+        if (damageState(ph)%sizeState > 0) &
-          damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state(:,me)
+          damageState(ph)%subState0(:,me) = damageState(ph)%state(:,me)
-        enddo
+
       endif
 !--------------------------------------------------------------------------------------------------
 !  cut back (reduced time and restore)
@@ -1207,9 +1202,8 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
         phase_mechanical_Li(ph)%data(1:3,1:3,me) = subLi0
       endif
       plasticState(ph)%state(:,me) = subState0
-      do so = 1, phase_Nsources(ph)
+      if (damageState(ph)%sizeState > 0) &
-        damageState(ph)%p(so)%state(:,me) = damageState(ph)%p(so)%subState0(:,me)
+        damageState(ph)%state(:,me) = damageState(ph)%subState0(:,me)
-      enddo
 
       todo = subStep > num%subStepMinCryst                          ! still on track or already done (beyond repair)
     endif
@@ -1303,7 +1297,7 @@ module function phase_mechanical_dPdF(dt,co,ip,el) result(dPdF)
   call phase_LiAndItsTangents(devNull,dLidS,dLidFi, &
                                      phase_mechanical_S(ph)%data(1:3,1:3,me), &
                                      phase_mechanical_Fi(ph)%data(1:3,1:3,me), &
-                                     co,ip,el)
+                                     ph,me)
 
   invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,me))
   invFi = math_inv33(phase_mechanical_Fi(ph)%data(1:3,1:3,me))
@@ -1505,4 +1499,4 @@ module subroutine phase_mechanical_setF(F,co,ip,el)
 end subroutine phase_mechanical_setF
 
 
-end submodule mechanics
+end submodule mechanical

src/phase_mechanical_eigen.f90

@@ -1,50 +1,29 @@
-submodule(phase:mechanics) eigendeformation
+submodule(phase:mechanical) eigen
+
+  integer, dimension(:), allocatable :: &
+    Nmodels
+
+  integer(kind(KINEMATICS_UNDEFINED_ID)),  dimension(:,:), allocatable :: &
+    model
+  integer(kind(KINEMATICS_UNDEFINED_ID)),  dimension(:), allocatable :: &
+    model_damage
 
   interface
-    module function kinematics_cleavage_opening_init(kinematics_length) result(myKinematics)
+    module function kinematics_cleavage_opening_init() result(myKinematics)
-      integer, intent(in) :: kinematics_length
+      logical, dimension(:), allocatable :: myKinematics
-      logical, dimension(:,:), allocatable :: myKinematics
     end function kinematics_cleavage_opening_init
 
-    module function kinematics_slipplane_opening_init(kinematics_length) result(myKinematics)
+    module function kinematics_slipplane_opening_init() result(myKinematics)
-      integer, intent(in) :: kinematics_length
+      logical, dimension(:), allocatable :: myKinematics
-      logical, dimension(:,:), allocatable :: myKinematics
     end function kinematics_slipplane_opening_init
 
-    module function kinematics_thermal_expansion_init(kinematics_length) result(myKinematics)
+    module function thermalexpansion_init(kinematics_length) result(myKinematics)
       integer, intent(in) :: kinematics_length
       logical, dimension(:,:), allocatable :: myKinematics
-    end function kinematics_thermal_expansion_init
+    end function thermalexpansion_init
-
-    module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el)
-      integer, intent(in) :: &
-        co, &                                                                                      !< grain number
-        ip, &                                                                                       !< integration point number
-        el                                                                                          !< element number
-      real(pReal),   intent(in),  dimension(3,3) :: &
-        S
-      real(pReal),   intent(out), dimension(3,3) :: &
-        Ld                                                                                          !< damage velocity gradient
-      real(pReal),   intent(out), dimension(3,3,3,3) :: &
-        dLd_dTstar                                                                                  !< derivative of Ld with respect to Tstar (4th-order tensor)
-    end subroutine kinematics_cleavage_opening_LiAndItsTangent
-
-    module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el)
-      integer, intent(in) :: &
-        co, &                                                                                      !< grain number
-        ip, &                                                                                       !< integration point number
-        el                                                                                          !< element number
-      real(pReal),   intent(in),  dimension(3,3) :: &
-        S
-      real(pReal),   intent(out), dimension(3,3) :: &
-        Ld                                                                                          !< damage velocity gradient
-      real(pReal),   intent(out), dimension(3,3,3,3) :: &
-        dLd_dTstar                                                                                  !< derivative of Ld with respect to Tstar (4th-order tensor)
-    end subroutine kinematics_slipplane_opening_LiAndItsTangent
 
     module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me)
       integer, intent(in) :: ph, me
-                                                                                       !< element number
       real(pReal),   intent(out), dimension(3,3) :: &
         Li                                                                                          !< thermal velocity gradient
       real(pReal),   intent(out), dimension(3,3,3,3) :: &
@@ -66,29 +45,36 @@ module subroutine eigendeformation_init(phases)
     ph
   class(tNode), pointer :: &
     phase, &
-   kinematics
+    kinematics, &
+    damage, &
+    mechanics
 
-  print'(/,a)', ' <<<+-  phase:mechanics:eigendeformation init  -+>>>'
+  print'(/,a)', ' <<<+-  phase:mechanical:eigen init  -+>>>'
 
 !--------------------------------------------------------------------------------------------------
-! initialize kinematic mechanisms
+! explicit eigen mechanisms
-  allocate(phase_Nkinematics(phases%length),source = 0)
+  allocate(Nmodels(phases%length),source = 0)
+
   do ph = 1,phases%length
     phase => phases%get(ph)
-    kinematics => phase%get('kinematics',defaultVal=emptyList)
+    mechanics => phase%get('mechanics')
-    phase_Nkinematics(ph) = kinematics%length
+    kinematics => mechanics%get('eigen',defaultVal=emptyList)
+    Nmodels(ph) = kinematics%length
   enddo
 
-  allocate(phase_kinematics(maxval(phase_Nkinematics),phases%length), source = KINEMATICS_undefined_ID)
+  allocate(model(maxval(Nmodels),phases%length), source = KINEMATICS_undefined_ID)
 
-  if(maxval(phase_Nkinematics) /= 0) then
+  if(maxval(Nmodels) /= 0) then
-    where(kinematics_cleavage_opening_init(maxval(phase_Nkinematics)))  phase_kinematics = KINEMATICS_cleavage_opening_ID
+    where(thermalexpansion_init(maxval(Nmodels))) model = KINEMATICS_thermal_expansion_ID
-    where(kinematics_slipplane_opening_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_slipplane_opening_ID
-    where(kinematics_thermal_expansion_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_thermal_expansion_ID
   endif
 
-end subroutine eigendeformation_init
+  allocate(model_damage(phases%length),  source = KINEMATICS_UNDEFINED_ID)
 
+  where(kinematics_cleavage_opening_init())  model_damage = KINEMATICS_cleavage_opening_ID
+  where(kinematics_slipplane_opening_init()) model_damage = KINEMATICS_slipplane_opening_ID
+
+
+end subroutine eigendeformation_init
 
 
 !--------------------------------------------------------------------------------------------------
@@ -104,17 +90,19 @@ function kinematics_active(kinematics_label,kinematics_length)  result(active_ki
     phases, &
     phase, &
     kinematics, &
-    kinematics_type
+    kinematics_type, &
+    mechanics
   integer :: p,k
 
   phases => config_material%get('phase')
   allocate(active_kinematics(kinematics_length,phases%length), source = .false. )
   do p = 1, phases%length
     phase => phases%get(p)
-    kinematics => phase%get('kinematics',defaultVal=emptyList)
+    mechanics => phase%get('mechanics')
+    kinematics => mechanics%get('eigen',defaultVal=emptyList)
     do k = 1, kinematics%length
       kinematics_type => kinematics%get(k)
-      if(kinematics_type%get_asString('type') == kinematics_label) active_kinematics(k,p) = .true.
+      active_kinematics(k,p) = kinematics_type%get_asString('type') == kinematics_label
     enddo
   enddo
 
@@ -122,17 +110,46 @@ function kinematics_active(kinematics_label,kinematics_length)  result(active_ki
 end function kinematics_active
 
 
+
+!--------------------------------------------------------------------------------------------------
+!> @brief checks if a kinematic mechanism is active or not
+!--------------------------------------------------------------------------------------------------
+function kinematics_active2(kinematics_label)  result(active_kinematics)
+
+  character(len=*), intent(in)       :: kinematics_label                                            !< name of kinematic mechanism
+  logical, dimension(:), allocatable :: active_kinematics
+
+  class(tNode), pointer :: &
+    phases, &
+    phase, &
+    kinematics, &
+    kinematics_type
+  integer :: p
+
+  phases => config_material%get('phase')
+  allocate(active_kinematics(phases%length), source = .false. )
+  do p = 1, phases%length
+    phase => phases%get(p)
+    kinematics => phase%get('damage',defaultVal=emptyList)
+    if(kinematics%length < 1) return
+    kinematics_type => kinematics%get(1)
+    if (.not. kinematics_type%contains('type')) continue
+    active_kinematics(p) = kinematics_type%get_asString('type',defaultVal='n/a') == kinematics_label
+  enddo
+
+
+end function kinematics_active2
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief  contains the constitutive equation for calculating the velocity gradient
 ! ToDo: MD: S is Mi?
 !--------------------------------------------------------------------------------------------------
 module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
-                                         S, Fi, co, ip, el)
+                                         S, Fi, ph,me)
 
   integer, intent(in) :: &
-    co, &                                                                                          !< component-ID of integration point
+    ph,me
-    ip, &                                                                                           !< integration point
-    el                                                                                              !< element
   real(pReal),   intent(in),  dimension(3,3) :: &
     S                                                                                               !< 2nd Piola-Kirchhoff stress
   real(pReal),   intent(in),  dimension(3,3) :: &
@@ -152,44 +169,49 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
   real(pReal) :: &
     detFi
   integer :: &
-    k, i, j, &
+    k, i, j
-    instance, of, me, ph
+  logical :: active
 
+  active = .false.
   Li = 0.0_pReal
   dLi_dS  = 0.0_pReal
   dLi_dFi = 0.0_pReal
 
-  plasticType: select case (phase_plasticity(material_phaseAt(co,el)))
+
+  plasticType: select case (phase_plasticity(ph))
     case (PLASTICITY_isotropic_ID) plasticType
-      of = material_phasememberAt(co,ip,el)
+      call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,phase_plasticInstance(ph),me)
-      instance = phase_plasticInstance(material_phaseAt(co,el))
+      Li = Li + my_Li
-      call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,instance,of)
+      dLi_dS = dLi_dS + my_dLi_dS
-    case default plasticType
+      active = .true.
-      my_Li = 0.0_pReal
-      my_dLi_dS = 0.0_pReal
   end select plasticType
 
-  Li = Li + my_Li
-  dLi_dS = dLi_dS + my_dLi_dS
 
-  KinematicsLoop: do k = 1, phase_Nkinematics(material_phaseAt(co,el))
+  KinematicsLoop: do k = 1, Nmodels(ph)
-    kinematicsType: select case (phase_kinematics(k,material_phaseAt(co,el)))
+    kinematicsType: select case (model(k,ph))
-      case (KINEMATICS_cleavage_opening_ID) kinematicsType
-        call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el)
-      case (KINEMATICS_slipplane_opening_ID) kinematicsType
-        call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el)
       case (KINEMATICS_thermal_expansion_ID) kinematicsType
-        me = material_phaseMemberAt(co,ip,el)
-        ph = material_phaseAt(co,el)
         call thermalexpansion_LiAndItsTangent(my_Li, my_dLi_dS, ph,me)
-      case default kinematicsType
+        Li = Li + my_Li
-        my_Li = 0.0_pReal
+        dLi_dS = dLi_dS + my_dLi_dS
-        my_dLi_dS = 0.0_pReal
+        active = .true.
     end select kinematicsType
-    Li = Li + my_Li
-    dLi_dS = dLi_dS + my_dLi_dS
   enddo KinematicsLoop
 
+  select case (model_damage(ph))
+    case (KINEMATICS_cleavage_opening_ID)
+      call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me)
+      Li = Li + my_Li
+      dLi_dS = dLi_dS + my_dLi_dS
+      active = .true.
+    case (KINEMATICS_slipplane_opening_ID)
+      call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me)
+      Li = Li + my_Li
+      dLi_dS = dLi_dS + my_dLi_dS
+      active = .true.
+  end select
+
+  if(.not. active) return
+
   FiInv = math_inv33(Fi)
   detFi = math_det33(Fi)
   Li = matmul(matmul(Fi,Li),FiInv)*detFi                                                            !< push forward to intermediate configuration
@@ -204,4 +226,4 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
 end subroutine phase_LiAndItsTangents
 
 
-end submodule eigendeformation
+end submodule eigen

src/phase_mechanical_eigen_cleavageopening.f90

@@ -4,24 +4,7 @@
 !> @brief material subroutine incorporating kinematics resulting from opening of cleavage planes
 !> @details to be done
 !--------------------------------------------------------------------------------------------------
-submodule(phase:eigendeformation) cleavageopening
+submodule(phase:eigen) cleavageopening
-
-  integer, dimension(:), allocatable :: kinematics_cleavage_opening_instance
-
-  type :: tParameters                                                                               !< container type for internal constitutive parameters
-    integer :: &
-      sum_N_cl                                                                                      !< total number of cleavage planes
-    real(pReal) :: &
-      dot_o, &                                                                                      !< opening rate of cleavage planes
-      q                                                                                             !< damage rate sensitivity
-    real(pReal),   dimension(:),   allocatable :: &
-      g_crit
-    real(pReal), dimension(:,:,:,:), allocatable :: &
-      cleavage_systems
-  end type tParameters
-
-  type(tParameters), dimension(:), allocatable :: param                                             !< containers of constitutive parameters (len Ninstances)
-
 
 contains
 
@@ -30,136 +13,20 @@ contains
 !> @brief module initialization
 !> @details reads in material parameters, allocates arrays, and does sanity checks
 !--------------------------------------------------------------------------------------------------
-module function kinematics_cleavage_opening_init(kinematics_length) result(myKinematics)
+module function kinematics_cleavage_opening_init() result(myKinematics)
 
-  integer, intent(in)                  :: kinematics_length
+  logical, dimension(:), allocatable :: myKinematics
-  logical, dimension(:,:), allocatable :: myKinematics
 
-  integer :: Ninstances,p,k
-  integer, dimension(:), allocatable :: N_cl                                                        !< active number of cleavage systems per family
-  character(len=pStringLen) :: extmsg = ''
-  class(tNode), pointer :: &
-    phases, &
-    phase, &
-    kinematics, &
-    kinematic_type
 
-  print'(/,a)', ' <<<+-  phase:mechanics:eigendeformation:cleavageopening init  -+>>>'
+  myKinematics = kinematics_active2('anisobrittle')
-
+  if(count(myKinematics) == 0) return
-  myKinematics = kinematics_active('cleavage_opening',kinematics_length)
-  Ninstances = count(myKinematics)
-  print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
-  if(Ninstances == 0) return
-
-  phases => config_material%get('phase')
-  allocate(param(Ninstances))
-  allocate(kinematics_cleavage_opening_instance(phases%length), source=0)
-
-  do p = 1, phases%length
-    if(any(myKinematics(:,p))) kinematics_cleavage_opening_instance(p) = count(myKinematics(:,1:p))
-    phase => phases%get(p)
-    if(count(myKinematics(:,p)) == 0) cycle
-    kinematics => phase%get('kinematics')
-    do k = 1, kinematics%length
-      if(myKinematics(k,p)) then
-        associate(prm  => param(kinematics_cleavage_opening_instance(p)))
-        kinematic_type => kinematics%get(k)
-
-        N_cl = kinematic_type%get_asInts('N_cl')
-        prm%sum_N_cl = sum(abs(N_cl))
-
-        prm%q       = kinematic_type%get_asFloat('q')
-        prm%dot_o   = kinematic_type%get_asFloat('dot_o')
-
-        prm%g_crit  = kinematic_type%get_asFloats('g_crit',requiredSize=size(N_cl))
-
-        prm%cleavage_systems  = lattice_SchmidMatrix_cleavage(N_cl,phase%get_asString('lattice'),&
-                                                        phase%get_asFloat('c/a',defaultVal=0.0_pReal))
-
-  ! expand: family => system
-        prm%g_crit = math_expand(prm%g_crit,N_cl)
-
-  ! sanity checks
-        if (prm%q          <= 0.0_pReal)  extmsg = trim(extmsg)//' q'
-        if (prm%dot_o      <= 0.0_pReal)  extmsg = trim(extmsg)//' dot_o'
-        if (any(prm%g_crit <  0.0_pReal)) extmsg = trim(extmsg)//' g_crit'
-
-!--------------------------------------------------------------------------------------------------
-!  exit if any parameter is out of range
-        if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(cleavage_opening)')
-        end associate
-      endif
-    enddo
-  enddo
 
+  print'(/,a)', ' <<<+-  phase:mechanical:eigen:cleavageopening init  -+>>>'
+  print'(a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT)
 
 end function kinematics_cleavage_opening_init
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief  contains the constitutive equation for calculating the velocity gradient
-!--------------------------------------------------------------------------------------------------
-module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el)
 
-  integer, intent(in) :: &
-    co, &                                                                                          !< grain number
-    ip, &                                                                                           !< integration point number
-    el                                                                                              !< element number
-  real(pReal),   intent(in),  dimension(3,3) :: &
-    S
-  real(pReal),   intent(out), dimension(3,3) :: &
-    Ld                                                                                              !< damage velocity gradient
-  real(pReal),   intent(out), dimension(3,3,3,3) :: &
-    dLd_dTstar                                                                                      !< derivative of Ld with respect to Tstar (4th-order tensor)
-
-  integer :: &
-    homog, damageOffset, &
-    i, k, l, m, n
-  real(pReal) :: &
-    traction_d, traction_t, traction_n, traction_crit, &
-    udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt
-
-  homog = material_homogenizationAt(el)
-  damageOffset = material_homogenizationMemberAt(ip,el)
-
-  Ld = 0.0_pReal
-  dLd_dTstar = 0.0_pReal
-  associate(prm => param(kinematics_cleavage_opening_instance(material_phaseAt(co,el))))
-  do i = 1,prm%sum_N_cl
-    traction_crit = prm%g_crit(i)* damage(homog)%p(damageOffset)**2.0_pReal
-
-    traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i))
-    if (abs(traction_d) > traction_crit + tol_math_check) then
-      udotd = sign(1.0_pReal,traction_d)* prm%dot_o * ((abs(traction_d) - traction_crit)/traction_crit)**prm%q
-      Ld = Ld + udotd*prm%cleavage_systems(1:3,1:3,1,i)
-      dudotd_dt = sign(1.0_pReal,traction_d)*udotd*prm%q / (abs(traction_d) - traction_crit)
-      forall (k=1:3,l=1:3,m=1:3,n=1:3) &
-        dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) &
-                            + dudotd_dt*prm%cleavage_systems(k,l,1,i) * prm%cleavage_systems(m,n,1,i)
-    endif
-
-    traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i))
-    if (abs(traction_t) > traction_crit + tol_math_check) then
-      udott = sign(1.0_pReal,traction_t)* prm%dot_o * ((abs(traction_t) - traction_crit)/traction_crit)**prm%q
-      Ld = Ld + udott*prm%cleavage_systems(1:3,1:3,2,i)
-      dudott_dt = sign(1.0_pReal,traction_t)*udott*prm%q / (abs(traction_t) - traction_crit)
-      forall (k=1:3,l=1:3,m=1:3,n=1:3) &
-        dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) &
-                            + dudott_dt*prm%cleavage_systems(k,l,2,i) * prm%cleavage_systems(m,n,2,i)
-    endif
-
-    traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i))
-    if (abs(traction_n) > traction_crit + tol_math_check) then
-      udotn = sign(1.0_pReal,traction_n)* prm%dot_o * ((abs(traction_n) - traction_crit)/traction_crit)**prm%q
-      Ld = Ld + udotn*prm%cleavage_systems(1:3,1:3,3,i)
-      dudotn_dt = sign(1.0_pReal,traction_n)*udotn*prm%q / (abs(traction_n) - traction_crit)
-      forall (k=1:3,l=1:3,m=1:3,n=1:3) &
-        dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) &
-                            + dudotn_dt*prm%cleavage_systems(k,l,3,i) * prm%cleavage_systems(m,n,3,i)
-    endif
-  enddo
-  end associate
-
-end subroutine kinematics_cleavage_opening_LiAndItsTangent
 
 end submodule cleavageopening

src/phase_mechanical_eigen_slipplaneopening.f90

@@ -4,7 +4,7 @@
 !> @brief material subroutine incorporating kinematics resulting from opening of slip planes
 !> @details to be done
 !--------------------------------------------------------------------------------------------------
-submodule(phase:eigendeformation) slipplaneopening
+submodule(phase:eigen) slipplaneopening
 
   integer, dimension(:), allocatable :: kinematics_slipplane_opening_instance
 
@@ -32,12 +32,11 @@ contains
 !> @brief module initialization
 !> @details reads in material parameters, allocates arrays, and does sanity checks
 !--------------------------------------------------------------------------------------------------
-module function kinematics_slipplane_opening_init(kinematics_length) result(myKinematics)
+module function kinematics_slipplane_opening_init() result(myKinematics)
 
-  integer, intent(in)                  :: kinematics_length
+  logical, dimension(:), allocatable :: myKinematics
-  logical, dimension(:,:), allocatable :: myKinematics
 
-  integer :: Ninstances,p,i,k
+  integer :: p,i
   character(len=pStringLen) :: extmsg = ''
   integer,     dimension(:),   allocatable :: N_sl
   real(pReal), dimension(:,:), allocatable :: d,n,t
@@ -49,28 +48,26 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi
     kinematics, &
     kinematic_type
 
-  print'(/,a)', ' <<<+-  phase:mechanics:eigendeformation:slipplaneopening init  -+>>>'
 
-  myKinematics = kinematics_active('slipplane_opening',kinematics_length)
+  myKinematics = kinematics_active2('isoductile')
-  Ninstances = count(myKinematics)
+  if(count(myKinematics) == 0) return
-  print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
+  print'(/,a)', ' <<<+-  phase:mechanical:eigen:slipplaneopening init  -+>>>'
-  if(Ninstances == 0) return
+  print'(a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT)
+
 
   phases => config_material%get('phase')
-  allocate(kinematics_slipplane_opening_instance(phases%length), source=0)
+  allocate(param(phases%length))
-  allocate(param(Ninstances))
 
   do p = 1, phases%length
-    if(any(myKinematics(:,p))) kinematics_slipplane_opening_instance(p) = count(myKinematics(:,1:p))
+    if(myKinematics(p)) then
-    phase => phases%get(p)
+      phase => phases%get(p)
-    mech  => phase%get('mechanics')
+      mech  => phase%get('mechanics')
-    pl    => mech%get('plasticity')
+      pl    => mech%get('plasticity')
-    if(count(myKinematics(:,p)) == 0) cycle
+
-    kinematics => phase%get('kinematics')
+      kinematics => phase%get('damage')
-    do k = 1, kinematics%length
+
-      if(myKinematics(k,p)) then
+      associate(prm  => param(p))
-        associate(prm  => param(kinematics_slipplane_opening_instance(p)))
+        kinematic_type => kinematics%get(1)
-        kinematic_type => kinematics%get(k)
 
         prm%dot_o    = kinematic_type%get_asFloat('dot_o')
         prm%q        = kinematic_type%get_asFloat('q')
@@ -105,9 +102,8 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi
 !  exit if any parameter is out of range
         if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(slipplane_opening)')
 
-        end associate
+      end associate
-      endif
+    endif
-    enddo
   enddo
 
 
@@ -117,12 +113,10 @@ end function kinematics_slipplane_opening_init
 !--------------------------------------------------------------------------------------------------
 !> @brief  contains the constitutive equation for calculating the velocity gradient
 !--------------------------------------------------------------------------------------------------
-module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el)
+module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me)
 
   integer, intent(in) :: &
-    co, &                                                                                          !< grain number
+    ph, me
-    ip, &                                                                                           !< integration point number
-    el                                                                                              !< element number
   real(pReal),   intent(in),  dimension(3,3) :: &
     S
   real(pReal),   intent(out), dimension(3,3) :: &
@@ -131,19 +125,13 @@ module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S
     dLd_dTstar                                                                                      !< derivative of Ld with respect to Tstar (4th-order tensor)
 
   integer :: &
-    instance, phase, &
-    homog, damageOffset, &
     i, k, l, m, n
   real(pReal) :: &
     traction_d, traction_t, traction_n, traction_crit, &
     udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt
 
-  phase = material_phaseAt(co,el)
-  instance = kinematics_slipplane_opening_instance(phase)
-  homog = material_homogenizationAt(el)
-  damageOffset = material_homogenizationMemberAt(ip,el)
 
-  associate(prm => param(instance))
+  associate(prm => param(ph))
   Ld = 0.0_pReal
   dLd_dTstar = 0.0_pReal
   do i = 1, prm%sum_N_sl
@@ -152,7 +140,7 @@ module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S
     traction_t = math_tensordot(S,prm%P_t(1:3,1:3,i))
     traction_n = math_tensordot(S,prm%P_n(1:3,1:3,i))
 
-    traction_crit = prm%g_crit(i)* damage(homog)%p(damageOffset)                                  ! degrading critical load carrying capacity by damage
+    traction_crit = prm%g_crit(i)* damage_phi(ph,me)
 
     udotd = sign(1.0_pReal,traction_d)* prm%dot_o* (  abs(traction_d)/traction_crit &
                                                     - abs(traction_d)/prm%g_crit(i))**prm%q

src/phase_mechanical_eigen_thermalexpansion.f90

@@ -3,7 +3,7 @@
 !> @brief material subroutine incorporating kinematics resulting from thermal expansion
 !> @details to be done
 !--------------------------------------------------------------------------------------------------
-submodule(phase:eigendeformation) thermalexpansion
+submodule(phase:eigen) thermalexpansion
 
   integer, dimension(:), allocatable :: kinematics_thermal_expansion_instance
 
@@ -23,7 +23,7 @@ contains
 !> @brief module initialization
 !> @details reads in material parameters, allocates arrays, and does sanity checks
 !--------------------------------------------------------------------------------------------------
-module function kinematics_thermal_expansion_init(kinematics_length) result(myKinematics)
+module function thermalexpansion_init(kinematics_length) result(myKinematics)
 
   integer, intent(in)                  :: kinematics_length
   logical, dimension(:,:), allocatable :: myKinematics
@@ -36,7 +36,7 @@ module function kinematics_thermal_expansion_init(kinematics_length) result(myKi
     kinematics, &
     kinematic_type
 
-  print'(/,a)', ' <<<+-  phase:mechanics:eigendeformation:thermalexpansion init  -+>>>'
+  print'(/,a)', ' <<<+-  phase:mechanical:eigen:thermalexpansion init  -+>>>'
 
   myKinematics = kinematics_active('thermal_expansion',kinematics_length)
   Ninstances = count(myKinematics)
@@ -77,7 +77,7 @@ module function kinematics_thermal_expansion_init(kinematics_length) result(myKi
   enddo
 
 
-end function kinematics_thermal_expansion_init
+end function thermalexpansion_init
 
 
 !--------------------------------------------------------------------------------------------------

src/phase_mechanical_plastic.f90

@@ -1,4 +1,4 @@
-submodule(phase:mechanics) plastic
+submodule(phase:mechanical) plastic
 
   interface
 
@@ -226,7 +226,7 @@ contains
 module subroutine plastic_init
 
 
-  print'(/,a)', ' <<<+-  phase:mechanics:plastic init  -+>>>'
+  print'(/,a)', ' <<<+-  phase:mechanical:plastic init  -+>>>'
 
   where(plastic_none_init())              phase_plasticity = PLASTICITY_NONE_ID
   where(plastic_isotropic_init())         phase_plasticity = PLASTICITY_ISOTROPIC_ID

src/phase_mechanical_plastic_dislotungsten.f90

@@ -97,13 +97,14 @@ module function plastic_dislotungsten_init() result(myPlasticity)
     mech, &
     pl
 
-  print'(/,a)', ' <<<+-  phase:mechanics:plastic:dislotungsten init  -+>>>'
-
   myPlasticity = plastic_active('dislotungsten')
   Ninstances = count(myPlasticity)
-  print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
   if(Ninstances == 0) return
 
+  print'(/,a)', ' <<<+-  phase:mechanical:plastic:dislotungsten init  -+>>>'
+  print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT)
+
+
   print*, 'Cereceda et al., International Journal of Plasticity 78:242–256, 2016'
   print*, 'https://dx.doi.org/10.1016/j.ijplas.2015.09.002'
 

src/phase_mechanical_plastic_dislotwin.f90

@@ -144,13 +144,13 @@ module function plastic_dislotwin_init() result(myPlasticity)
     mech, &
     pl
 
-  print'(/,a)', ' <<<+-  phase:mechanics:plastic:dislotwin init  -+>>>'
-
   myPlasticity = plastic_active('dislotwin')
   Ninstances = count(myPlasticity)
-  print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
   if(Ninstances == 0) return
 
+  print'(/,a)', ' <<<+-  phase:mechanical:plastic:dislotwin init  -+>>>'
+  print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT)
+
   print*, 'Ma and Roters, Acta Materialia 52(12):3603–3612, 2004'
   print*, 'https://doi.org/10.1016/j.actamat.2004.04.012'//IO_EOL
 

src/phase_mechanical_plastic_isotropic.f90

@@ -68,13 +68,14 @@ module function plastic_isotropic_init() result(myPlasticity)
     mech, &
     pl
 
-  print'(/,a)', ' <<<+-  phase:mechanics:plastic:isotropic init  -+>>>'
-
   myPlasticity = plastic_active('isotropic')
   Ninstances = count(myPlasticity)
-  print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
   if(Ninstances == 0) return
 
+  print'(/,a)', ' <<<+-  phase:mechanical:plastic:isotropic init  -+>>>'
+  print'(a,i0)', ' # phses: ',Ninstances; flush(IO_STDOUT)
+
+
   print*, 'Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018'
   print*, 'https://doi.org/10.1016/j.scriptamat.2017.09.047'
 

src/phase_mechanical_plastic_kinehardening.f90

@@ -80,13 +80,14 @@ module function plastic_kinehardening_init() result(myPlasticity)
     mech, &
     pl
 
-  print'(/,a)', ' <<<+-  phase:mechanics:plastic:kinehardening init  -+>>>'
-
   myPlasticity = plastic_active('kinehardening')
   Ninstances = count(myPlasticity)
-  print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
   if(Ninstances == 0) return
 
+  print'(/,a)', ' <<<+-  phase:mechanical:plastic:kinehardening init  -+>>>'
+  print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT)
+
+
   allocate(param(Ninstances))
   allocate(state(Ninstances))
   allocate(dotState(Ninstances))

src/phase_mechanical_plastic_none.f90

@@ -16,32 +16,20 @@ module function plastic_none_init() result(myPlasticity)
 
   logical, dimension(:), allocatable :: myPlasticity
   integer :: &
-    Ninstances, &
     p, &
     Nconstituents
   class(tNode), pointer :: &
-    phases, &
+    phases
-    phase, &
-    mech, &
-    pl
 
-  print'(/,a)', ' <<<+-  phase:mechanics:plastic:none init  -+>>>'
+
+  myPlasticity = plastic_active('none')
+  if(count(myPlasticity) == 0) return
+
+  print'(/,a)', ' <<<+-  phase:mechanical:plastic:none init  -+>>>'
+  print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
 
   phases => config_material%get('phase')
-  allocate(myPlasticity(phases%length), source = .false.)
   do p = 1, phases%length
-    phase => phases%get(p)
-    mech  => phase%get('mechanics')
-    pl    => mech%get ('plasticity')
-    if(pl%get_asString('type') == 'none') myPlasticity(p) = .true.
-  enddo
-
-  Ninstances = count(myPlasticity)
-  print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
-  if(Ninstances == 0) return
-
-  do p = 1, phases%length
-    phase => phases%get(p)
     if(.not. myPlasticity(p)) cycle
     Nconstituents = count(material_phaseAt2 == p)
     call phase_allocateState(plasticState(p),Nconstituents,0,0,0)

src/phase_mechanical_plastic_nonlocal.f90

@@ -187,16 +187,16 @@ module function plastic_nonlocal_init() result(myPlasticity)
     mech, &
     pl
 
-  print'(/,a)', ' <<<+-  phase:mechanics:plastic:nonlocal init  -+>>>'
-
   myPlasticity = plastic_active('nonlocal')
   Ninstances = count(myPlasticity)
-  print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
   if(Ninstances == 0) then
     call geometry_plastic_nonlocal_disable
     return
   endif
 
+  print'(/,a)', ' <<<+-  phase:mechanical:plastic:nonlocal init  -+>>>'
+  print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT)
+
   print*, 'Reuber et al., Acta Materialia 71:333–348, 2014'
   print*, 'https://doi.org/10.1016/j.actamat.2014.03.012'//IO_EOL
 

src/phase_mechanical_plastic_phenopowerlaw.f90

@@ -89,13 +89,15 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
     mech, &
     pl
 
-  print'(/,a)', ' <<<+-  phase:mechanics:plastic:phenopowerlaw init  -+>>>'
 
   myPlasticity = plastic_active('phenopowerlaw')
   Ninstances = count(myPlasticity)
-  print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
   if(Ninstances == 0) return
 
+  print'(/,a)', ' <<<+-  phase:mechanical:plastic:phenopowerlaw init  -+>>>'
+  print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT)
+
+
   allocate(param(Ninstances))
   allocate(state(Ninstances))
   allocate(dotState(Ninstances))

src/phase_thermal.f90

@@ -3,6 +3,12 @@
 !----------------------------------------------------------------------------------------------------
 submodule(phase) thermal
 
+  integer, dimension(:), allocatable :: &
+    thermal_Nsources
+
+  type(tSourceState),  allocatable, dimension(:) :: &
+    thermalState
+
   enum, bind(c); enumerator :: &
     THERMAL_UNDEFINED_ID ,&
     THERMAL_DISSIPATION_ID, &
@@ -33,8 +39,6 @@ submodule(phase) thermal
     end function externalheat_init
 
 
-
-
     module subroutine externalheat_dotState(ph, me)
       integer, intent(in) :: &
         ph, &

src/phase_thermal_dissipation.f90

@@ -31,15 +31,14 @@ module function dissipation_init(source_length) result(mySources)
     phase, &
     sources, thermal, &
     src
-  integer :: Ninstances,so,Nconstituents,ph
+  integer :: so,Nconstituents,ph
 
-  print'(/,a)', ' <<<+-  phase:thermal:dissipation init  -+>>>'
 
   mySources = thermal_active('dissipation',source_length)
+  if(count(mySources) == 0) return
+  print'(/,a)', ' <<<+-  phase:thermal:dissipation init  -+>>>'
+  print'(a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT)
 
-  Ninstances = count(mySources)
-  print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
-  if(Ninstances == 0) return
 
   phases => config_material%get('phase')
   allocate(param(phases%length))

src/phase_thermal_externalheat.f90

@@ -8,7 +8,7 @@ submodule(phase:thermal) externalheat
 
 
   integer,           dimension(:),   allocatable :: &
-    source_thermal_externalheat_offset                                                           !< which source is my current thermal dissipation mechanism?
+    source_thermal_externalheat_offset                                                              !< which source is my current thermal dissipation mechanism?
 
   type :: tParameters                                                                               !< container type for internal constitutive parameters
     real(pReal), dimension(:), allocatable :: &
@@ -38,15 +38,14 @@ module function externalheat_init(source_length) result(mySources)
     phase, &
     sources, thermal, &
     src
-  integer :: Ninstances,so,Nconstituents,ph
+  integer :: so,Nconstituents,ph
 
-  print'(/,a)', ' <<<+-  phase:thermal:externalheat init  -+>>>'
 
   mySources = thermal_active('externalheat',source_length)
+  if(count(mySources) == 0) return
+  print'(/,a)', ' <<<+-  phase:thermal:externalheat init  -+>>>'
+  print'(a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT)
 
-  Ninstances = count(mySources)
-  print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
-  if(Ninstances == 0) return
 
   phases => config_material%get('phase')
   allocate(param(phases%length))