diff --git a/src/CPFEM.f90 b/src/CPFEM.f90
index dfa8dd09d..be11c0107 100644
--- a/src/CPFEM.f90
+++ b/src/CPFEM.f90
@@ -262,7 +262,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
      if (debug_e <= discretization_nElem .and. debug_i <=discretization_nIP) then
        write(6,'(a,1x,i8,1x,i2,1x,i4,/,(12x,6(e20.8,1x)),/)') &
              '<< CPFEM >> aged state of elFE ip grain',debug_e, debug_i, 1, &
-              plasticState(phaseAt(1,debug_i,debug_e))%state(:,phasememberAt(1,debug_i,debug_e))
+              plasticState(material_phaseAt(1,debug_e))%state(:,material_phasememberAt(1,debug_i,debug_e))
        endif
    endif
 
diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90
index 51d9af152..c9b497f9d 100644
--- a/src/CPFEM2.f90
+++ b/src/CPFEM2.f90
@@ -89,7 +89,6 @@ subroutine CPFEM_init
 
     fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5')
    
-    call HDF5_read(fileHandle,material_phase,      'recordedPhase')
     call HDF5_read(fileHandle,crystallite_F0,      'convergedF')
     call HDF5_read(fileHandle,crystallite_Fp0,     'convergedFp')
     call HDF5_read(fileHandle,crystallite_Fi0,     'convergedFi')
@@ -158,7 +157,6 @@ subroutine CPFEM_age
     write(rankStr,'(a1,i0)')'_',worldrank
     fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','a')
     
-    call HDF5_write(fileHandle,material_phase,      'recordedPhase')
     call HDF5_write(fileHandle,crystallite_F0,      'convergedF')
     call HDF5_write(fileHandle,crystallite_Fp0,     'convergedFp')
     call HDF5_write(fileHandle,crystallite_Fi0,     'convergedFi')
diff --git a/src/constitutive.f90 b/src/constitutive.f90
index cc70e8f22..ec66ab126 100644
--- a/src/constitutive.f90
+++ b/src/constitutive.f90
@@ -251,7 +251,7 @@ function constitutive_homogenizedC(ipc,ip,el)
    ip, &                                                                                             !< integration point
    el                                                                                                !< element
 
- plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
+ plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el)))
    case (PLASTICITY_DISLOTWIN_ID) plasticityType
      constitutive_homogenizedC = plastic_dislotwin_homogenizedC(ipc,ip,el)
    case default plasticityType
@@ -280,14 +280,14 @@ subroutine constitutive_microstructure(Fe, Fp, ipc, ip, el)
  ho = material_homogenizationAt(el)
  tme = thermalMapping(ho)%p(ip,el)
 
- plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
+ plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el)))
    case (PLASTICITY_DISLOTWIN_ID) plasticityType
-     of = phasememberAt(ipc,ip,el)
-     instance = phase_plasticityInstance(material_phase(ipc,ip,el))
+     of = material_phasememberAt(ipc,ip,el)
+     instance = phase_plasticityInstance(material_phaseAt(ipc,el))
      call plastic_dislotwin_dependentState(temperature(ho)%p(tme),instance,of)
    case (PLASTICITY_DISLOUCLA_ID) plasticityType
-     of = phasememberAt(ipc,ip,el)
-     instance = phase_plasticityInstance(material_phase(ipc,ip,el))
+     of = material_phasememberAt(ipc,ip,el)
+     instance = phase_plasticityInstance(material_phaseAt(ipc,el))
      call plastic_disloUCLA_dependentState(instance,of)
    case (PLASTICITY_NONLOCAL_ID) plasticityType
      call plastic_nonlocal_dependentState (Fe,Fp,ip,el)
@@ -331,25 +331,25 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
 
  Mp  = matmul(matmul(transpose(Fi),Fi),S)
 
- plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
+ plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el)))
 
    case (PLASTICITY_NONE_ID) plasticityType
      Lp = 0.0_pReal
      dLp_dMp = 0.0_pReal
 
    case (PLASTICITY_ISOTROPIC_ID) plasticityType
-     of = phasememberAt(ipc,ip,el)
-     instance = phase_plasticityInstance(material_phase(ipc,ip,el))
+     of = material_phasememberAt(ipc,ip,el)
+     instance = phase_plasticityInstance(material_phaseAt(ipc,el))
      call plastic_isotropic_LpAndItsTangent       (Lp,dLp_dMp,Mp,instance,of)
 
    case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
-     of = phasememberAt(ipc,ip,el)
-     instance = phase_plasticityInstance(material_phase(ipc,ip,el))
+     of = material_phasememberAt(ipc,ip,el)
+     instance = phase_plasticityInstance(material_phaseAt(ipc,el))
      call plastic_phenopowerlaw_LpAndItsTangent   (Lp,dLp_dMp,Mp,instance,of)
 
    case (PLASTICITY_KINEHARDENING_ID) plasticityType
-     of = phasememberAt(ipc,ip,el)
-     instance = phase_plasticityInstance(material_phase(ipc,ip,el))
+     of = material_phasememberAt(ipc,ip,el)
+     instance = phase_plasticityInstance(material_phaseAt(ipc,el))
      call plastic_kinehardening_LpAndItsTangent   (Lp,dLp_dMp, Mp,instance,of)
 
    case (PLASTICITY_NONLOCAL_ID) plasticityType
@@ -357,13 +357,13 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
                                                    temperature(ho)%p(tme),geometry_plastic_nonlocal_IPvolume0(ip,el),ip,el)
 
    case (PLASTICITY_DISLOTWIN_ID) plasticityType
-     of = phasememberAt(ipc,ip,el)
-     instance = phase_plasticityInstance(material_phase(ipc,ip,el))
+     of = material_phasememberAt(ipc,ip,el)
+     instance = phase_plasticityInstance(material_phaseAt(ipc,el))
      call plastic_dislotwin_LpAndItsTangent       (Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of)
 
    case (PLASTICITY_DISLOUCLA_ID) plasticityType
-     of = phasememberAt(ipc,ip,el)
-     instance = phase_plasticityInstance(material_phase(ipc,ip,el))
+     of = material_phasememberAt(ipc,ip,el)
+     instance = phase_plasticityInstance(material_phaseAt(ipc,el))
      call plastic_disloucla_LpAndItsTangent       (Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of)
 
  end select plasticityType
@@ -414,10 +414,10 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
  dLi_dS  = 0.0_pReal
  dLi_dFi = 0.0_pReal
 
- plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
+ plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el)))
    case (PLASTICITY_isotropic_ID) plasticityType
-     of = phasememberAt(ipc,ip,el)
-     instance = phase_plasticityInstance(material_phase(ipc,ip,el))
+     of = material_phasememberAt(ipc,ip,el)
+     instance = phase_plasticityInstance(material_phaseAt(ipc,el))
      call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,instance,of)
    case default plasticityType
      my_Li = 0.0_pReal
@@ -427,8 +427,8 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
  Li = Li + my_Li
  dLi_dS = dLi_dS + my_dLi_dS
 
- KinematicsLoop: do k = 1, phase_Nkinematics(material_phase(ipc,ip,el))
-   kinematicsType: select case (phase_kinematics(k,material_phase(ipc,ip,el)))
+ KinematicsLoop: do k = 1, phase_Nkinematics(material_phaseAt(ipc,el))
+   kinematicsType: select case (phase_kinematics(k,material_phaseAt(ipc,el)))
      case (KINEMATICS_cleavage_opening_ID) kinematicsType
        call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ipc, ip, el)
      case (KINEMATICS_slipplane_opening_ID) kinematicsType
@@ -475,7 +475,7 @@ pure function constitutive_initialFi(ipc, ip, el)
    homog, offset
 
  constitutive_initialFi = math_I3
- phase = material_phase(ipc,ip,el)
+ phase = material_phaseAt(ipc,el)
 
  KinematicsLoop: do k = 1, phase_Nkinematics(phase)                                            !< Warning: small initial strain assumption
    kinematicsType: select case (phase_kinematics(k,phase))
@@ -546,8 +546,8 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
  ho = material_homogenizationAt(el)
  C = math_66toSym3333(constitutive_homogenizedC(ipc,ip,el))
 
- DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phase(ipc,ip,el))
-   degradationType: select case(phase_stiffnessDegradation(d,material_phase(ipc,ip,el)))
+ DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phaseAt(ipc,el))
+   degradationType: select case(phase_stiffnessDegradation(d,material_phaseAt(ipc,el)))
      case (STIFFNESS_DEGRADATION_damage_ID) degradationType
        C = C * damage(ho)%p(damageMapping(ho)%p(ip,el))**2
    end select degradationType
@@ -597,31 +597,31 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip,
 
  Mp  = matmul(matmul(transpose(Fi),Fi),S)
 
- plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
+ plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el)))
 
    case (PLASTICITY_ISOTROPIC_ID) plasticityType
-     of = phasememberAt(ipc,ip,el)
-     instance = phase_plasticityInstance(material_phase(ipc,ip,el))
+     of = material_phasememberAt(ipc,ip,el)
+     instance = phase_plasticityInstance(material_phaseAt(ipc,el))
      call plastic_isotropic_dotState    (Mp,instance,of)
 
    case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
-     of = phasememberAt(ipc,ip,el)
-     instance = phase_plasticityInstance(material_phase(ipc,ip,el))
+     of = material_phasememberAt(ipc,ip,el)
+     instance = phase_plasticityInstance(material_phaseAt(ipc,el))
      call plastic_phenopowerlaw_dotState(Mp,instance,of)
 
    case (PLASTICITY_KINEHARDENING_ID) plasticityType
-     of = phasememberAt(ipc,ip,el)
-     instance = phase_plasticityInstance(material_phase(ipc,ip,el))
+     of = material_phasememberAt(ipc,ip,el)
+     instance = phase_plasticityInstance(material_phaseAt(ipc,el))
      call plastic_kinehardening_dotState(Mp,instance,of)
 
    case (PLASTICITY_DISLOTWIN_ID) plasticityType
-     of = phasememberAt(ipc,ip,el)
-     instance = phase_plasticityInstance(material_phase(ipc,ip,el))
+     of = material_phasememberAt(ipc,ip,el)
+     instance = phase_plasticityInstance(material_phaseAt(ipc,el))
      call plastic_dislotwin_dotState    (Mp,temperature(ho)%p(tme),instance,of)
 
    case (PLASTICITY_DISLOUCLA_ID) plasticityType
-     of = phasememberAt(ipc,ip,el)
-     instance = phase_plasticityInstance(material_phase(ipc,ip,el))
+     of = material_phasememberAt(ipc,ip,el)
+     instance = phase_plasticityInstance(material_phaseAt(ipc,el))
      call plastic_disloucla_dotState    (Mp,temperature(ho)%p(tme),instance,of)
 
    case (PLASTICITY_NONLOCAL_ID) plasticityType
@@ -629,9 +629,9 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip,
                                          subdt,ip,el)
  end select plasticityType
 
- SourceLoop: do i = 1, phase_Nsources(material_phase(ipc,ip,el))
+ SourceLoop: do i = 1, phase_Nsources(material_phaseAt(ipc,el))
 
-   sourceType: select case (phase_source(i,material_phase(ipc,ip,el)))
+   sourceType: select case (phase_source(i,material_phaseAt(ipc,el)))
 
      case (SOURCE_damage_anisoBrittle_ID) sourceType
        call source_damage_anisoBrittle_dotState (S, ipc, ip, el) !< correct stress?
@@ -643,8 +643,8 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip,
        call source_damage_anisoDuctile_dotState (         ipc, ip, el)
 
      case (SOURCE_thermal_externalheat_ID) sourceType
-       of = phasememberAt(ipc,ip,el)
-       call source_thermal_externalheat_dotState(material_phase(ipc,ip,el),of)
+       of = material_phasememberAt(ipc,ip,el)
+       call source_thermal_externalheat_dotState(material_phaseAt(ipc,el),of)
 
    end select sourceType
 
@@ -674,11 +674,11 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
 
  Mp  = matmul(matmul(transpose(Fi),Fi),S)
 
- plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
+ plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el)))
 
    case (PLASTICITY_KINEHARDENING_ID) plasticityType
-     of = phasememberAt(ipc,ip,el)
-     instance = phase_plasticityInstance(material_phase(ipc,ip,el))
+     of = material_phasememberAt(ipc,ip,el)
+     instance = phase_plasticityInstance(material_phaseAt(ipc,el))
      call plastic_kinehardening_deltaState(Mp,instance,of)
 
    case (PLASTICITY_NONLOCAL_ID) plasticityType
@@ -686,9 +686,9 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
 
  end select plasticityType
 
- sourceLoop: do i = 1, phase_Nsources(material_phase(ipc,ip,el))
+ sourceLoop: do i = 1, phase_Nsources(material_phaseAt(ipc,el))
 
-    sourceType: select case (phase_source(i,material_phase(ipc,ip,el)))
+    sourceType: select case (phase_source(i,material_phaseAt(ipc,el)))
 
      case (SOURCE_damage_isoBrittle_ID) sourceType
        call source_damage_isoBrittle_deltaState  (constitutive_homogenizedC(ipc,ip,el), Fe, &
@@ -710,8 +710,8 @@ function constitutive_postResults(S, Fi, ipc, ip, el)
    ipc, &                                                                                           !< component-ID of integration point
    ip, &                                                                                            !< integration point
    el                                                                                               !< element
- real(pReal), dimension(plasticState(material_phase(ipc,ip,el))%sizePostResults + &
-                        sum(sourceState(material_phase(ipc,ip,el))%p(:)%sizePostResults)) :: &
+ real(pReal), dimension(plasticState(material_phaseAt(ipc,el))%sizePostResults + &
+                        sum(sourceState(material_phaseAt(ipc,el))%p(:)%sizePostResults)) :: &
    constitutive_postResults
  real(pReal),  intent(in), dimension(3,3) :: &
    Fi                                                                                               !< intermediate deformation gradient
@@ -734,12 +734,12 @@ function constitutive_postResults(S, Fi, ipc, ip, el)
  tme = thermalMapping(ho)%p(ip,el)
 
  startPos = 1
- endPos = plasticState(material_phase(ipc,ip,el))%sizePostResults
+ endPos = plasticState(material_phaseAt(ipc,el))%sizePostResults
 
- of = phasememberAt(ipc,ip,el)
- instance = phase_plasticityInstance(material_phase(ipc,ip,el))
+ of = material_phasememberAt(ipc,ip,el)
+ instance = phase_plasticityInstance(material_phaseAt(ipc,el))
 
- plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
+ plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el)))
    case (PLASTICITY_ISOTROPIC_ID) plasticityType
      constitutive_postResults(startPos:endPos) = &
        plastic_isotropic_postResults(Mp,instance,of)
@@ -762,23 +762,23 @@ function constitutive_postResults(S, Fi, ipc, ip, el)
 
    case (PLASTICITY_NONLOCAL_ID) plasticityType
      constitutive_postResults(startPos:endPos) = &
-       plastic_nonlocal_postResults (material_phase(ipc,ip,el),instance,of)
+       plastic_nonlocal_postResults (material_phaseAt(ipc,el),instance,of)
 
  end select plasticityType
 
- SourceLoop: do i = 1, phase_Nsources(material_phase(ipc,ip,el))
+ SourceLoop: do i = 1, phase_Nsources(material_phaseAt(ipc,el))
    startPos = endPos + 1
-   endPos = endPos + sourceState(material_phase(ipc,ip,el))%p(i)%sizePostResults
-   of = phasememberAt(ipc,ip,el)
-   sourceType: select case (phase_source(i,material_phase(ipc,ip,el)))
+   endPos = endPos + sourceState(material_phaseAt(ipc,el))%p(i)%sizePostResults
+   of = material_phasememberAt(ipc,ip,el)
+   sourceType: select case (phase_source(i,material_phaseAt(ipc,el)))
      case (SOURCE_damage_isoBrittle_ID) sourceType
-       constitutive_postResults(startPos:endPos) = source_damage_isoBrittle_postResults(material_phase(ipc,ip,el),of)
+       constitutive_postResults(startPos:endPos) = source_damage_isoBrittle_postResults(material_phaseAt(ipc,el),of)
      case (SOURCE_damage_isoDuctile_ID) sourceType
-       constitutive_postResults(startPos:endPos) = source_damage_isoDuctile_postResults(material_phase(ipc,ip,el),of)
+       constitutive_postResults(startPos:endPos) = source_damage_isoDuctile_postResults(material_phaseAt(ipc,el),of)
      case (SOURCE_damage_anisoBrittle_ID) sourceType
-       constitutive_postResults(startPos:endPos) = source_damage_anisoBrittle_postResults(material_phase(ipc,ip,el),of)
+       constitutive_postResults(startPos:endPos) = source_damage_anisoBrittle_postResults(material_phaseAt(ipc,el),of)
      case (SOURCE_damage_anisoDuctile_ID) sourceType
-       constitutive_postResults(startPos:endPos) = source_damage_anisoDuctile_postResults(material_phase(ipc,ip,el),of)
+       constitutive_postResults(startPos:endPos) = source_damage_anisoDuctile_postResults(material_phaseAt(ipc,el),of)
    end select sourceType
 
  enddo SourceLoop
diff --git a/src/crystallite.f90 b/src/crystallite.f90
index 98070f0fa..e47eac4df 100644
--- a/src/crystallite.f90
+++ b/src/crystallite.f90
@@ -386,7 +386,7 @@ subroutine crystallite_init
       crystallite_Fp0(1:3,1:3,c,i,e) = math_EulerToR(material_EulerAngles(1:3,c,i,e))               ! plastic def gradient reflects init orientation
       crystallite_Fi0(1:3,1:3,c,i,e) = constitutive_initialFi(c,i,e)
       crystallite_F0(1:3,1:3,c,i,e)  = math_I3
-      crystallite_localPlasticity(c,i,e) = phase_localPlasticity(material_phase(c,i,e))
+      crystallite_localPlasticity(c,i,e) = phase_localPlasticity(material_phaseAt(c,e))
       crystallite_Fe(1:3,1:3,c,i,e)  = math_inv33(matmul(crystallite_Fi0(1:3,1:3,c,i,e), &
                                                          crystallite_Fp0(1:3,1:3,c,i,e)))           ! assuming that euler angles are given in internal strain free configuration
       crystallite_Fp(1:3,1:3,c,i,e)  = crystallite_Fp0(1:3,1:3,c,i,e)
@@ -483,12 +483,12 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
   elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2)
     do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e); do c = 1,homogenization_Ngrains(material_homogenizationAt(e))
       homogenizationRequestsCalculation: if (crystallite_requested(c,i,e)) then
-        plasticState    (phaseAt(c,i,e))%subState0(      :,phasememberAt(c,i,e)) = &
-        plasticState    (phaseAt(c,i,e))%partionedState0(:,phasememberAt(c,i,e))
+        plasticState    (material_phaseAt(c,e))%subState0(      :,material_phaseMemberAt(c,i,e)) = &
+        plasticState    (material_phaseAt(c,e))%partionedState0(:,material_phaseMemberAt(c,i,e))
 
-        do s = 1, phase_Nsources(phaseAt(c,i,e))
-          sourceState(phaseAt(c,i,e))%p(s)%subState0(      :,phasememberAt(c,i,e)) = &
-          sourceState(phaseAt(c,i,e))%p(s)%partionedState0(:,phasememberAt(c,i,e))
+        do s = 1, phase_Nsources(material_phaseAt(c,e))
+          sourceState(material_phaseAt(c,e))%p(s)%subState0(      :,material_phaseMemberAt(c,i,e)) = &
+          sourceState(material_phaseAt(c,e))%p(s)%partionedState0(:,material_phaseMemberAt(c,i,e))
         enddo
         crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_partionedFp0(1:3,1:3,c,i,e)
         crystallite_subLp0(1:3,1:3,c,i,e) = crystallite_partionedLp0(1:3,1:3,c,i,e)
@@ -543,11 +543,11 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
               crystallite_subFi0(1:3,1:3,c,i,e) = crystallite_Fi  (1:3,1:3,c,i,e)
               crystallite_subS0 (1:3,1:3,c,i,e) = crystallite_S   (1:3,1:3,c,i,e)
               !if abbrevation, make c and p private in omp
-              plasticState(    phaseAt(c,i,e))%subState0(:,phasememberAt(c,i,e)) &
-                = plasticState(phaseAt(c,i,e))%state(    :,phasememberAt(c,i,e))
-              do s = 1, phase_Nsources(phaseAt(c,i,e))
-                sourceState(    phaseAt(c,i,e))%p(s)%subState0(:,phasememberAt(c,i,e)) &
-                  = sourceState(phaseAt(c,i,e))%p(s)%state(    :,phasememberAt(c,i,e))
+              plasticState(    material_phaseAt(c,e))%subState0(:,material_phaseMemberAt(c,i,e)) &
+                = plasticState(material_phaseAt(c,e))%state(    :,material_phaseMemberAt(c,i,e))
+              do s = 1, phase_Nsources(material_phaseAt(c,e))
+                sourceState(    material_phaseAt(c,e))%p(s)%subState0(:,material_phaseMemberAt(c,i,e)) &
+                  = sourceState(material_phaseAt(c,e))%p(s)%state(    :,material_phaseMemberAt(c,i,e))
               enddo
 #ifdef DEBUG
               if (iand(debug_level(debug_crystallite),debug_levelBasic) /= 0 &
@@ -572,11 +572,11 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
               crystallite_Lp (1:3,1:3,c,i,e) =            crystallite_subLp0(1:3,1:3,c,i,e)
               crystallite_Li (1:3,1:3,c,i,e) =            crystallite_subLi0(1:3,1:3,c,i,e)
             endif
-            plasticState    (phaseAt(c,i,e))%state(    :,phasememberAt(c,i,e)) &
-              = plasticState(phaseAt(c,i,e))%subState0(:,phasememberAt(c,i,e))
-            do s = 1, phase_Nsources(phaseAt(c,i,e))
-              sourceState(    phaseAt(c,i,e))%p(s)%state(    :,phasememberAt(c,i,e)) &
-                = sourceState(phaseAt(c,i,e))%p(s)%subState0(:,phasememberAt(c,i,e))
+            plasticState    (material_phaseAt(c,e))%state(    :,material_phaseMemberAt(c,i,e)) &
+              = plasticState(material_phaseAt(c,e))%subState0(:,material_phaseMemberAt(c,i,e))
+            do s = 1, phase_Nsources(material_phaseAt(c,e))
+              sourceState(    material_phaseAt(c,e))%p(s)%state(    :,material_phaseMemberAt(c,i,e)) &
+                = sourceState(material_phaseAt(c,e))%p(s)%subState0(:,material_phaseMemberAt(c,i,e))
             enddo
 
                                                                                                     ! cant restore dotState here, since not yet calculated in first cutback after initialization
@@ -839,7 +839,7 @@ subroutine crystallite_orientations
     !$OMP PARALLEL DO
     do e = FEsolving_execElem(1),FEsolving_execElem(2)
       do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
-        if (plasticState(material_phase(1,i,e))%nonLocal) &                                         ! if nonlocal model
+        if (plasticState(material_phaseAt(1,e))%nonLocal) &                                         ! if nonlocal model
           call plastic_nonlocal_updateCompatibility(crystallite_orientation,i,e)
     enddo; enddo
     !$OMP END PARALLEL DO
@@ -879,8 +879,8 @@ function crystallite_postResults(ipc, ip, el)
    ipc                           !< grain index
 
  real(pReal), dimension(1+crystallite_sizePostResults(microstructure_crystallite(discretization_microstructureAt(el))) + &
-                        1+plasticState(material_phase(ipc,ip,el))%sizePostResults + &
-                          sum(sourceState(material_phase(ipc,ip,el))%p(:)%sizePostResults)) :: &
+                        1+plasticState(material_phaseAt(ipc,el))%sizePostResults + &
+                          sum(sourceState(material_phaseAt(ipc,el))%p(:)%sizePostResults)) :: &
    crystallite_postResults
  integer :: &
    o, &
@@ -893,7 +893,7 @@ function crystallite_postResults(ipc, ip, el)
  crystID = microstructure_crystallite(discretization_microstructureAt(el))
 
  crystallite_postResults = 0.0_pReal
- crystallite_postResults(1) = real(crystallite_sizePostResults(crystID),pReal)                  ! header-like information (length)
+ crystallite_postResults(1) = real(crystallite_sizePostResults(crystID),pReal)                      ! header-like information (length)
  c = 1
 
  do o = 1,crystallite_Noutput(crystID)
@@ -901,7 +901,7 @@ function crystallite_postResults(ipc, ip, el)
    select case(crystallite_outputID(o,crystID))
      case (phase_ID)
        mySize = 1
-       crystallite_postResults(c+1) = real(material_phase(ipc,ip,el),pReal)                         ! phaseID of grain
+       crystallite_postResults(c+1) = real(material_phaseAt(ipc,el),pReal)                          ! phaseID of grain
      case (texture_ID)
        mySize = 1
        crystallite_postResults(c+1) = real(material_texture(ipc,ip,el),pReal)                       ! textureID of grain
@@ -967,7 +967,7 @@ function crystallite_postResults(ipc, ip, el)
    c = c + mySize
  enddo
 
- crystallite_postResults(c+1) = real(plasticState(material_phase(ipc,ip,el))%sizePostResults,pReal)             ! size of constitutive results
+ crystallite_postResults(c+1) = real(plasticState(material_phaseAt(ipc,el))%sizePostResults,pReal)  ! size of constitutive results
  c = c + 1
  if (size(crystallite_postResults)-c > 0) &
    crystallite_postResults(c+1:size(crystallite_postResults)) = &
@@ -1555,7 +1555,7 @@ subroutine integrateStateFPI
        do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
          do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
            if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
-             p = phaseAt(g,i,e); c = phasememberAt(g,i,e)
+             p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
 
              plasticState(p)%previousDotState2(:,c) = merge(plasticState(p)%previousDotState(:,c),&
                                                             0.0_pReal,&
@@ -1583,7 +1583,7 @@ subroutine integrateStateFPI
      do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
        do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
          if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
-           p = phaseAt(g,i,e); c = phasememberAt(g,i,e)
+           p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
            sizeDotState = plasticState(p)%sizeDotState
 
            zeta = damper(plasticState(p)%dotState         (:,c), &
@@ -1746,7 +1746,7 @@ subroutine integrateStateAdaptiveEuler
    do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
      do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
        if (crystallite_todo(g,i,e)) then
-         p = phaseAt(g,i,e); c = phasememberAt(g,i,e)
+         p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
          sizeDotState = plasticState(p)%sizeDotState
          
          residuum_plastic(1:sizeDotState,g,i,e) = plasticState(p)%dotstate(1:sizeDotState,c) &
@@ -1775,7 +1775,7 @@ subroutine integrateStateAdaptiveEuler
    do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
      do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
        if (crystallite_todo(g,i,e)) then
-         p = phaseAt(g,i,e); c = phasememberAt(g,i,e)
+         p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
          sizeDotState = plasticState(p)%sizeDotState
          
          residuum_plastic(1:sizeDotState,g,i,e) = residuum_plastic(1:sizeDotState,g,i,e) &
@@ -1835,7 +1835,7 @@ subroutine integrateStateRK4
      do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
        do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
        if (crystallite_todo(g,i,e)) then
-         p = phaseAt(g,i,e); c = phasememberAt(g,i,e)
+         p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
 
          plasticState(p)%RK4dotState(:,c) = WEIGHT(n)*plasticState(p)%dotState(:,c) &
                                           + merge(plasticState(p)%RK4dotState(:,c),0.0_pReal,n>1)
@@ -1926,7 +1926,7 @@ subroutine integrateStateRKCK45
      do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
        do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
        if (crystallite_todo(g,i,e)) then
-         p = phaseAt(g,i,e); cc = phasememberAt(g,i,e)
+         p = material_phaseAt(g,e); cc = material_phaseMemberAt(g,i,e)
 
          plasticState(p)%RKCK45dotState(stage,:,cc) = plasticState(p)%dotState(:,cc)
          plasticState(p)%dotState(:,cc) = A(1,stage) * plasticState(p)%RKCK45dotState(1,:,cc)
@@ -1966,7 +1966,7 @@ subroutine integrateStateRKCK45
      do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
        do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
      if (crystallite_todo(g,i,e)) then
-       p = phaseAt(g,i,e); cc = phasememberAt(g,i,e)
+       p = material_phaseAt(g,e); cc = material_phaseMemberAt(g,i,e)
        
        sizeDotState = plasticState(p)%sizeDotState
               
@@ -2005,7 +2005,7 @@ subroutine integrateStateRKCK45
    do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
      do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
        if (crystallite_todo(g,i,e)) then
-         p  = phaseAt(g,i,e); cc = phasememberAt(g,i,e)
+         p  = material_phaseAt(g,e); cc = material_phaseMemberAt(g,i,e)
        
          sizeDotState = plasticState(p)%sizeDotState
          
@@ -2163,7 +2163,7 @@ subroutine update_state(timeFraction)
      do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
      do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
          if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
-       p = phaseAt(g,i,e); c = phasememberAt(g,i,e)
+       p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
 
        mySize = plasticState(p)%sizeDotState
        plasticState(p)%state(1:mySize,c) = plasticState(p)%subState0(1:mySize,c) &
@@ -2214,7 +2214,7 @@ subroutine update_dotState(timeFraction)
                                              crystallite_Fi(1:3,1:3,g,i,e), &
                                              crystallite_Fp, &
                                              crystallite_subdt(g,i,e)*timeFraction, g,i,e)
-           p = phaseAt(g,i,e); c = phasememberAt(g,i,e)
+           p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
            NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
            do s = 1, phase_Nsources(p)
              NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c)))
@@ -2259,7 +2259,7 @@ subroutine update_deltaState
                                             crystallite_Fe(1:3,1:3,g,i,e), &
                                             crystallite_Fi(1:3,1:3,g,i,e), &
                                             g,i,e)
-         p = phaseAt(g,i,e); c = phasememberAt(g,i,e)
+         p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
          myOffset = plasticState(p)%offsetDeltaState
          mySize   = plasticState(p)%sizeDeltaState
          NaN = any(IEEE_is_NaN(plasticState(p)%deltaState(1:mySize,c)))
@@ -2311,8 +2311,8 @@ logical function stateJump(ipc,ip,el)
    myOffset, &
    mySize
 
- c = phasememberAt(ipc,ip,el)
- p = phaseAt(ipc,ip,el)
+ c = material_phaseMemberAt(ipc,ip,el)
+ p = material_phaseAt(ipc,el)
 
  call constitutive_collectDeltaState(crystallite_S(1:3,1:3,ipc,ip,el), &
                                      crystallite_Fe(1:3,1:3,ipc,ip,el), &
diff --git a/src/damage_local.f90 b/src/damage_local.f90
index 2764bfcb0..74ad47c9b 100644
--- a/src/damage_local.f90
+++ b/src/damage_local.f90
@@ -178,8 +178,8 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el
   phiDot = 0.0_pReal
   dPhiDot_dPhi = 0.0_pReal
   do grain = 1, homogenization_Ngrains(material_homogenizationAt(el))
-    phase = phaseAt(grain,ip,el)
-    constituent = phasememberAt(grain,ip,el)
+    phase = material_phaseAt(grain,el)
+    constituent = material_phasememberAt(grain,ip,el)
     do source = 1, phase_Nsources(phase)
       select case(phase_source(source,phase))                                                   
         case (SOURCE_damage_isoBrittle_ID)
diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90
index d8ab8bf1b..dce88bc61 100644
--- a/src/damage_nonlocal.f90
+++ b/src/damage_nonlocal.f90
@@ -144,8 +144,8 @@ subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip,
   phiDot = 0.0_pReal
   dPhiDot_dPhi = 0.0_pReal
   do grain = 1, homogenization_Ngrains(material_homogenizationAt(el))
-    phase = phaseAt(grain,ip,el)
-    constituent = phasememberAt(grain,ip,el)
+    phase = material_phaseAt(grain,el)
+    constituent = material_phasememberAt(grain,ip,el)
     do source = 1, phase_Nsources(phase)
       select case(phase_source(source,phase))                                                   
         case (SOURCE_damage_isoBrittle_ID)
diff --git a/src/homogenization.f90 b/src/homogenization.f90
index 35e2e9e03..9f7890131 100644
--- a/src/homogenization.f90
+++ b/src/homogenization.f90
@@ -326,11 +326,11 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
    do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e);
      do g = 1,myNgrains
 
-       plasticState    (phaseAt(g,i,e))%partionedState0(:,phasememberAt(g,i,e)) = &
-       plasticState    (phaseAt(g,i,e))%state0(         :,phasememberAt(g,i,e))
-       do mySource = 1, phase_Nsources(phaseAt(g,i,e))
-         sourceState(phaseAt(g,i,e))%p(mySource)%partionedState0(:,phasememberAt(g,i,e)) = &
-         sourceState(phaseAt(g,i,e))%p(mySource)%state0(         :,phasememberAt(g,i,e))
+       plasticState    (material_phaseAt(g,e))%partionedState0(:,material_phasememberAt(g,i,e)) = &
+       plasticState    (material_phaseAt(g,e))%state0(         :,material_phasememberAt(g,i,e))
+       do mySource = 1, phase_Nsources(material_phaseAt(g,e))
+         sourceState(material_phaseAt(g,e))%p(mySource)%partionedState0(:,material_phasememberAt(g,i,e)) = &
+         sourceState(material_phaseAt(g,e))%p(mySource)%state0(         :,material_phasememberAt(g,i,e))
        enddo
 
        crystallite_partionedFp0(1:3,1:3,g,i,e) = crystallite_Fp0(1:3,1:3,g,i,e)
@@ -412,11 +412,11 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
              crystallite_S          (1:3,1:3,1:myNgrains,i,e)
 
            do g = 1,myNgrains
-             plasticState    (phaseAt(g,i,e))%partionedState0(:,phasememberAt(g,i,e)) = &
-             plasticState    (phaseAt(g,i,e))%state          (:,phasememberAt(g,i,e))
-             do mySource = 1, phase_Nsources(phaseAt(g,i,e))
-               sourceState(phaseAt(g,i,e))%p(mySource)%partionedState0(:,phasememberAt(g,i,e)) = &
-               sourceState(phaseAt(g,i,e))%p(mySource)%state          (:,phasememberAt(g,i,e))
+             plasticState    (material_phaseAt(g,e))%partionedState0(:,material_phasememberAt(g,i,e)) = &
+             plasticState    (material_phaseAt(g,e))%state          (:,material_phasememberAt(g,i,e))
+             do mySource = 1, phase_Nsources(material_phaseAt(g,e))
+               sourceState(material_phaseAt(g,e))%p(mySource)%partionedState0(:,material_phasememberAt(g,i,e)) = &
+               sourceState(material_phaseAt(g,e))%p(mySource)%state          (:,material_phasememberAt(g,i,e))
              enddo
            enddo
 
@@ -475,11 +475,11 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
            crystallite_S(1:3,1:3,1:myNgrains,i,e) = &
               crystallite_partionedS0(1:3,1:3,1:myNgrains,i,e)
            do g = 1, myNgrains
-             plasticState    (phaseAt(g,i,e))%state(          :,phasememberAt(g,i,e)) = &
-             plasticState    (phaseAt(g,i,e))%partionedState0(:,phasememberAt(g,i,e))
-             do mySource = 1, phase_Nsources(phaseAt(g,i,e))
-               sourceState(phaseAt(g,i,e))%p(mySource)%state(          :,phasememberAt(g,i,e)) = &
-               sourceState(phaseAt(g,i,e))%p(mySource)%partionedState0(:,phasememberAt(g,i,e))
+             plasticState    (material_phaseAt(g,e))%state(          :,material_phasememberAt(g,i,e)) = &
+             plasticState    (material_phaseAt(g,e))%partionedState0(:,material_phasememberAt(g,i,e))
+             do mySource = 1, phase_Nsources(material_phaseAt(g,e))
+               sourceState(material_phaseAt(g,e))%p(mySource)%state(          :,material_phasememberAt(g,i,e)) = &
+               sourceState(material_phaseAt(g,e))%p(mySource)%partionedState0(:,material_phasememberAt(g,i,e))
              enddo
            enddo
            if(homogState(material_homogenizationAt(e))%sizeState > 0) &
diff --git a/src/kinematics_cleavage_opening.f90 b/src/kinematics_cleavage_opening.f90
index 39bfbf340..42761f8df 100644
--- a/src/kinematics_cleavage_opening.f90
+++ b/src/kinematics_cleavage_opening.f90
@@ -138,7 +138,7 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, i
     traction_d, traction_t, traction_n, traction_crit, &
     udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt
 
-  phase = material_phase(ipc,ip,el)
+  phase = material_phaseAt(ipc,el)
   instance = kinematics_cleavage_opening_instance(phase)
   homog = material_homogenizationAt(el)
   damageOffset = damageMapping(homog)%p(ip,el)
diff --git a/src/kinematics_slipplane_opening.f90 b/src/kinematics_slipplane_opening.f90
index 3e37e4c0d..a8aaa0e82 100644
--- a/src/kinematics_slipplane_opening.f90
+++ b/src/kinematics_slipplane_opening.f90
@@ -124,7 +124,7 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc,
    traction_d, traction_t, traction_n, traction_crit, &
    udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt
    
- phase = material_phase(ipc,ip,el)
+ phase = material_phaseAt(ipc,el)
  instance = kinematics_slipplane_opening_instance(phase)
  homog = material_homogenizationAt(el)
  damageOffset = damageMapping(homog)%p(ip,el)
diff --git a/src/kinematics_thermal_expansion.f90 b/src/kinematics_thermal_expansion.f90
index b4f23dfa7..4294321a5 100644
--- a/src/kinematics_thermal_expansion.f90
+++ b/src/kinematics_thermal_expansion.f90
@@ -112,7 +112,7 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip,
   real(pReal) :: &
     T, TRef, TDot  
     
-  phase = material_phase(ipc,ip,el)
+  phase = material_phaseAt(ipc,el)
   homog = material_homogenizationAt(el)
   offset = thermalMapping(homog)%p(ip,el)
   T = temperature(homog)%p(offset)
diff --git a/src/material.f90 b/src/material.f90
index f4a700229..f4014ec84 100644
--- a/src/material.f90
+++ b/src/material.f90
@@ -180,9 +180,6 @@ module material
    homogenization_active
 
 ! BEGIN DEPRECATED
- integer, dimension(:,:,:), allocatable, public :: phaseAt                                          !< phase ID of every material point (ipc,ip,el)
- integer, dimension(:,:,:), allocatable, public :: phasememberAt                                    !< memberID of given phase at every material point (ipc,ip,el)
-
  integer, dimension(:,:,:), allocatable, public, target :: mappingHomogenization                    !< mapping from material points to offset in heterogenous state/field
  integer, dimension(:,:),   allocatable, private, target :: mappingHomogenizationConst               !< mapping from material points to offset in constant state/field
 ! END DEPRECATED
@@ -247,11 +244,7 @@ contains
 subroutine material_init
 
  integer, parameter :: FILEUNIT = 210
- integer            :: m,c,h, myDebug, myPhase, myHomog
- integer :: &
-  g, &                                                                                              !< grain number
-  i, &                                                                                              !< integration point number
-  e                                                                                                 !< element number
+ integer            :: i,e,m,c,h, myDebug, myPhase, myHomog
  integer, dimension(:), allocatable :: &
   CounterPhase, &
   CounterHomogenization
@@ -375,26 +368,15 @@ subroutine material_init
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 ! BEGIN DEPRECATED
- allocate(phaseAt                   (  homogenization_maxNgrains,discretization_nIP,discretization_nElem),source=0)
- allocate(phasememberAt             (  homogenization_maxNgrains,discretization_nIP,discretization_nElem),source=0)
  allocate(mappingHomogenization     (2,                          discretization_nIP,discretization_nElem),source=0)
  allocate(mappingHomogenizationConst(                            discretization_nIP,discretization_nElem),source=1)
  
  CounterHomogenization=0
- CounterPhase         =0
-
-
  do e = 1,discretization_nElem
- myHomog = discretization_homogenizationAt(e)
+   myHomog = discretization_homogenizationAt(e)
    do i = 1, discretization_nIP
      CounterHomogenization(myHomog) = CounterHomogenization(myHomog) + 1
      mappingHomogenization(1:2,i,e) = [CounterHomogenization(myHomog),huge(1)]
-     do g = 1,homogenization_Ngrains(myHomog)
-       myPhase = material_phase(g,i,e)
-       CounterPhase(myPhase) = CounterPhase(myPhase)+1                             ! not distinguishing between instances of same phase
-       phaseAt(g,i,e)              = myPhase
-       phasememberAt(g,i,e)        = CounterPhase(myPhase)
-     enddo
    enddo
  enddo
 ! END DEPRECATED
diff --git a/src/plastic_disloUCLA.f90 b/src/plastic_disloUCLA.f90
index c8ef6fa15..28e1227aa 100644
--- a/src/plastic_disloUCLA.f90
+++ b/src/plastic_disloUCLA.f90
@@ -13,6 +13,7 @@ module plastic_disloUCLA
   use material
   use config
   use lattice
+  use discretization
   use results
  
   implicit none
@@ -295,7 +296,7 @@ subroutine plastic_disloUCLA_init()
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
-    NipcMyPhase = count(material_phase == p)
+    NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
     sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl
     sizeState = sizeDotState
  
diff --git a/src/plastic_dislotwin.f90 b/src/plastic_dislotwin.f90
index ad0ba41b6..f6696cdcc 100644
--- a/src/plastic_dislotwin.f90
+++ b/src/plastic_dislotwin.f90
@@ -15,6 +15,7 @@ module plastic_dislotwin
   use material
   use config
   use lattice
+  use discretization
 #if defined(PETSc) || defined(DAMASK_HDF5)
   use results
 #endif
@@ -494,7 +495,7 @@ subroutine plastic_dislotwin_init
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
-   NipcMyPhase = count(material_phase == p)
+   NipcMyPhase  = count(material_phaseAt == p) * discretization_nIP
    sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl &
                 + size(['f_tw'])                           * prm%sum_N_tw &
                 + size(['f_tr'])                           * prm%sum_N_tr
@@ -581,9 +582,9 @@ function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC)
             of
  real(pReal) :: f_unrotated
 
- of = phasememberAt(ipc,ip,el)
- associate(prm => param(phase_plasticityInstance(material_phase(ipc,ip,el))),&
-           stt => state(phase_plasticityInstance(material_phase(ipc,ip,el))))
+ of = material_phasememberAt(ipc,ip,el)
+ associate(prm => param(phase_plasticityInstance(material_phaseAt(ipc,el))),&
+           stt => state(phase_plasticityInstance(material_phaseAT(ipc,el))))
 
  f_unrotated = 1.0_pReal &
              - sum(stt%f_tw(1:prm%sum_N_tw,of)) &
diff --git a/src/plastic_isotropic.f90 b/src/plastic_isotropic.f90
index 46d0905dc..3b5ec6e5a 100644
--- a/src/plastic_isotropic.f90
+++ b/src/plastic_isotropic.f90
@@ -14,9 +14,8 @@ module plastic_isotropic
   use IO
   use material
   use config
-#if defined(PETSc) || defined(DAMASK_HDF5)
+  use discretization
   use results
-#endif
  
   implicit none
   private
@@ -127,8 +126,8 @@ subroutine plastic_isotropic_init
               config => config_phase(p))
  
 #ifdef DEBUG
-    if  (p==material_phase(debug_g,debug_i,debug_e)) &
-      prm%of_debug = phasememberAt(debug_g,debug_i,debug_e)
+    if  (p==material_phaseAt(debug_g,debug_e)) &
+      prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e)
 #endif
  
     prm%xi_0            = config%getFloat('tau0')
@@ -190,7 +189,7 @@ subroutine plastic_isotropic_init
  
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
-    NipcMyPhase = count(material_phase == p)
+    NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
     sizeDotState = size(['xi               ','accumulated_shear'])
     sizeState = sizeDotState
  
diff --git a/src/plastic_kinematichardening.f90 b/src/plastic_kinematichardening.f90
index ab68eb176..f2183327c 100644
--- a/src/plastic_kinematichardening.f90
+++ b/src/plastic_kinematichardening.f90
@@ -13,9 +13,8 @@ module plastic_kinehardening
  use material
  use config
  use lattice
-#if defined(PETSc) || defined(DAMASK_HDF5)
+ use discretization
  use results
-#endif
 
  implicit none
  private
@@ -146,8 +145,8 @@ subroutine plastic_kinehardening_init
              config => config_phase(p))
 
 #ifdef DEBUG
-   if  (p==material_phase(debug_g,debug_i,debug_e)) then
-      prm%of_debug = phasememberAt(debug_g,debug_i,debug_e)
+   if  (p==material_phaseAt(debug_g,debug_e)) then
+      prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e)
    endif
 #endif
 
@@ -257,7 +256,7 @@ subroutine plastic_kinehardening_init
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
-   NipcMyPhase = count(material_phase == p)
+   NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
    sizeDotState   = size(['crss     ','crss_back', 'accshear ']) * prm%totalNslip
    sizeDeltaState = size(['sense ',   'chi0  ',    'gamma0'   ]) * prm%totalNslip
    sizeState = sizeDotState + sizeDeltaState
diff --git a/src/plastic_none.f90 b/src/plastic_none.f90
index 894cc9a40..fa913df4c 100644
--- a/src/plastic_none.f90
+++ b/src/plastic_none.f90
@@ -6,6 +6,7 @@
 !--------------------------------------------------------------------------------------------------
 module plastic_none
  use material
+ use discretization
  use debug
 
  implicit none
@@ -36,7 +37,7 @@ subroutine plastic_none_init
  do p = 1, size(phase_plasticity)
    if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle
 
-   NipcMyPhase = count(material_phase == p)
+   NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
    call material_allocatePlasticState(p,NipcMyPhase,0,0,0, &
                                       0,0,0)
    plasticState(p)%sizePostResults = 0
diff --git a/src/plastic_nonlocal.f90 b/src/plastic_nonlocal.f90
index baed5a066..854fc9177 100644
--- a/src/plastic_nonlocal.f90
+++ b/src/plastic_nonlocal.f90
@@ -556,7 +556,7 @@ subroutine plastic_nonlocal_init
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
-    NofMyPhase=count(material_phase==p)
+    NofMyPhase   = count(material_phaseAt==p) * discretization_nIP
     sizeDotState =     size([   'rhoSglEdgePosMobile   ','rhoSglEdgeNegMobile   ', &
                                 'rhoSglScrewPosMobile  ','rhoSglScrewNegMobile  ', &
                                 'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', &
@@ -677,7 +677,7 @@ subroutine plastic_nonlocal_init
   allocate(iD(maxval(totalNslip),2,maxNinstances),    source=0)
 
   initializeInstances: do p = 1, size(phase_plasticity)
-    NofMyPhase=count(material_phase==p)
+    NofMyPhase = count(material_phaseAt==p) * discretization_nIP
     myPhase2: if (phase_plasticity(p) == PLASTICITY_NONLOCAL_ID) then
 
       !*** determine indices to state array
@@ -766,7 +766,7 @@ subroutine plastic_nonlocal_init
     ! get the total volume of the instance
     do e = 1,discretization_nElem
       do i = 1,discretization_nIP
-        if (material_phase(1,i,e) == phase) volume(phasememberAt(1,i,e)) = IPvolume(i,e)
+        if (material_phaseAt(1,e) == phase) volume(material_phasememberAt(1,i,e)) = IPvolume(i,e)
       enddo
     enddo
     totalVolume = sum(volume)
@@ -854,29 +854,29 @@ subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el)
     invConnections
   real(pReal), dimension(3,nIPneighbors) :: &
     connection_latticeConf
-  real(pReal), dimension(2,totalNslip(phase_plasticityInstance(material_phase(1,ip,el)))) :: &
+  real(pReal), dimension(2,totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: &
     rhoExcess
-  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el)))) :: &
+  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: &
     rho_edg_delta, &
     rho_scr_delta
-  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),10) :: &
+  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),10) :: &
     rho, &
     rho_neighbor
-  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))), &
-                         totalNslip(phase_plasticityInstance(material_phase(1,ip,el)))) :: &
+  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))), &
+                         totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: &
     myInteractionMatrix                                                                             ! corrected slip interaction matrix
                                   
-  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),nIPneighbors) :: &
+  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),nIPneighbors) :: &
     rho_edg_delta_neighbor, &
     rho_scr_delta_neighbor
   real(pReal), dimension(2,maxval(totalNslip),nIPneighbors) :: &
     neighbor_rhoExcess, &                                                                           ! excess density at neighboring material point
     neighbor_rhoTotal                                                                               ! total density at neighboring material point
-  real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),2) :: &
+  real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),2) :: &
     m                                                                                               ! direction of dislocation motion
   
-  ph = phaseAt(1,ip,el)
-  of = phasememberAt(1,ip,el)
+  ph = material_phaseAt(1,el)
+  of = material_phasememberAt(1,ip,el)
   instance = phase_plasticityInstance(ph)
 
   associate(prm => param(instance),dst => microstructure(instance), stt => state(instance))
@@ -935,9 +935,9 @@ subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el)
     do n = 1,nIPneighbors
       neighbor_el = IPneighborhood(1,n,ip,el)
       neighbor_ip = IPneighborhood(2,n,ip,el)
-      no = phasememberAt(1,neighbor_ip,neighbor_el)
+      no = material_phasememberAt(1,neighbor_ip,neighbor_el)
       if (neighbor_el > 0 .and. neighbor_ip > 0) then
-        neighbor_instance = phase_plasticityInstance(material_phase(1,neighbor_ip,neighbor_el))
+        neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el))
         if (neighbor_instance == instance) then
   
             nRealNeighbors = nRealNeighbors + 1.0_pReal
@@ -1202,22 +1202,22 @@ subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, &
     of, &                                                                                           !offset
     t, &                                                                                            !< dislocation type
     s                                                                                               !< index of my current slip system
-  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),8) :: &
+  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),8) :: &
     rhoSgl                                                                                          !< single dislocation densities (including blocked)
-  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),10) :: &
+  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),10) :: &
     rho
-  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),4) :: &
+  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),4) :: &
     v, &                                                                                            !< velocity
     tauNS, &                                                                                        !< resolved shear stress including non Schmid and backstress terms
     dv_dtau, &                                                                                      !< velocity derivative with respect to the shear stress
     dv_dtauNS                                                                                       !< velocity derivative with respect to the shear stress
-  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el)))) :: &
+  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: &
     tau, &                                                                                          !< resolved shear stress including backstress terms
     gdotTotal                                                                                       !< shear rate
   
   !*** shortcut for mapping
-  ph = phaseAt(1,ip,el)
-  of = phasememberAt(1,ip,el)
+  ph = material_phaseAt(1,el)
+  of = material_phasememberAt(1,ip,el)
   
   instance = phase_plasticityInstance(ph)
   associate(prm => param(instance),dst=>microstructure(instance))
@@ -1323,23 +1323,23 @@ subroutine plastic_nonlocal_deltaState(Mp,ip,el)
     c, &                                                                                            ! character of dislocation
     t, &                                                                                            ! type of dislocation
     s                                                                                               ! index of my current slip system
-  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),10) :: &
+  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),10) :: &
     deltaRhoRemobilization, &                                                                       ! density increment by remobilization
     deltaRhoDipole2SingleStress                                                                     ! density increment by dipole dissociation (by stress change)
-  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),10) :: &
+  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),10) :: &
     rho                                                                                             ! current  dislocation densities
-  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),4) :: &
+  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),4) :: &
     v                                                                                               ! dislocation glide velocity
-  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el)))) :: &
+  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: &
     tau                                                                                             ! current resolved shear stress
-  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),2) :: &
+  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),2) :: &
     rhoDip, &                                                                                       ! current dipole dislocation densities (screw and edge dipoles)
     dUpper, &                                                                                       ! current maximum stable dipole distance for edges and screws
     dUpperOld, &                                                                                    ! old maximum stable dipole distance for edges and screws
     deltaDUpper                                                                                     ! change in maximum stable dipole distance for edges and screws
   
-   ph = phaseAt(1,ip,el)
-   of = phasememberAt(1,ip,el)
+   ph = material_phaseAt(1,el)
+   of = material_phasememberAt(1,ip,el)
    instance = phase_plasticityInstance(ph)
 
    associate(prm => param(instance),dst => microstructure(instance),del => deltaState(instance))
@@ -1459,7 +1459,7 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
     np,&                                                                                            !< neighbour phase shortcut
     topp, &                                                                                         !< type of dislocation with opposite sign to t
     s                                                                                               !< index of my current slip system
-  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),10) :: &
+  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),10) :: &
     rho, &
     rhoDot, &                                                                                       !< density evolution
     rhoDotMultiplication, &                                                                         !< density evolution by multiplication
@@ -1467,24 +1467,24 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
     rhoDotSingle2DipoleGlide, &                                                                     !< density evolution by dipole formation (by glide)
     rhoDotAthermalAnnihilation, &                                                                   !< density evolution by athermal annihilation
     rhoDotThermalAnnihilation                                                                       !< density evolution by thermal annihilation
-  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),8) :: &
+  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),8) :: &
     rhoSgl, &                                                                                       !< current single dislocation densities (positive/negative screw and edge without dipoles)
     neighbor_rhoSgl, &                                                                              !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles)
     my_rhoSgl                                                                                       !< single dislocation densities of central ip (positive/negative screw and edge without dipoles)
-  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),4) :: &
+  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),4) :: &
     v, &                                                                                            !< current dislocation glide velocity
     my_v, &                                                                                         !< dislocation glide velocity of central ip
     neighbor_v, &                                                                                   !< dislocation glide velocity of enighboring ip
     gdot                                                                                            !< shear rates
-  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el)))) :: &
+  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: &
     tau, &                                                                                          !< current resolved shear stress
     vClimb                                                                                          !< climb velocity of edge dipoles
                                               
-  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),2) :: &
+  real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),2) :: &
     rhoDip, &                                                                                       !< current dipole dislocation densities (screw and edge dipoles)
     dLower, &                                                                                       !< minimum stable dipole distance for edges and screws
     dUpper                                                                                          !< current maximum stable dipole distance for edges and screws
-  real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),4) :: &
+  real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),4) :: &
     m                                                                                               !< direction of dislocation motion
   real(pReal), dimension(3,3) :: &
     my_F, &                                                                                         !< my total deformation gradient
@@ -1507,15 +1507,15 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
     considerEnteringFlux, &
     considerLeavingFlux
                                               
-  p = phaseAt(1,ip,el)
-  o = phasememberAt(1,ip,el)
+  p = material_phaseAt(1,el)
+  o = material_phasememberAt(1,ip,el)
   
   if (timestep <= 0.0_pReal) then
     plasticState(p)%dotState = 0.0_pReal
     return
   endif
   
-  ph = material_phase(1,ip,el)
+  ph = material_phaseAt(1,el)
   instance = phase_plasticityInstance(ph)
   associate(prm => param(instance),dst => microstructure(instance),dot => dotState(instance),stt => state(instance))
   ns = totalNslip(instance)
@@ -1592,7 +1592,7 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
   !****************************************************************************
   !*** calculate dislocation fluxes (only for nonlocal plasticity)
   rhoDotFlux = 0.0_pReal
-  if (.not. phase_localPlasticity(material_phase(1,ip,el))) then
+  if (.not. phase_localPlasticity(material_phaseAt(1,el))) then
   
     !*** check CFL (Courant-Friedrichs-Lewy) condition for flux
     if (any( abs(gdot) > 0.0_pReal &                                                                ! any active slip system ...
@@ -1630,8 +1630,8 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
       neighbor_el = IPneighborhood(1,n,ip,el)
       neighbor_ip = IPneighborhood(2,n,ip,el)
       neighbor_n  = IPneighborhood(3,n,ip,el)
-      np = phaseAt(1,neighbor_ip,neighbor_el)
-      no = phasememberAt(1,neighbor_ip,neighbor_el)
+      np = material_phaseAt(1,neighbor_el)
+      no = material_phasememberAt(1,neighbor_ip,neighbor_el)
   
       opposite_neighbor = n + mod(n,2) - mod(n+1,2)
       opposite_el = IPneighborhood(1,opposite_neighbor,ip,el)
@@ -1639,7 +1639,7 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
       opposite_n  = IPneighborhood(3,opposite_neighbor,ip,el)
     
       if (neighbor_n > 0) then                                                                      ! if neighbor exists, average deformation gradient
-        neighbor_instance = phase_plasticityInstance(material_phase(1,neighbor_ip,neighbor_el))
+        neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el))
         neighbor_Fe = Fe(1:3,1:3,1,neighbor_ip,neighbor_el)
         neighbor_F = matmul(neighbor_Fe, Fp(1:3,1:3,1,neighbor_ip,neighbor_el))
         Favg = 0.5_pReal * (my_F + neighbor_F)
@@ -1661,7 +1661,7 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
       neighbor_v = 0.0_pReal        ! needed for check of sign change in flux density below 
       neighbor_rhoSgl = 0.0_pReal
       if (neighbor_n > 0) then
-        if (phase_plasticity(material_phase(1,neighbor_ip,neighbor_el)) == PLASTICITY_NONLOCAL_ID &
+        if (phase_plasticity(material_phaseAt(1,neighbor_el)) == PLASTICITY_NONLOCAL_ID &
             .and. any(compatibility(:,:,:,n,ip,el) > 0.0_pReal)) &
           considerEnteringFlux = .true.
       endif
@@ -1714,7 +1714,7 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
       
       considerLeavingFlux = .true.
       if (opposite_n > 0) then
-        if (phase_plasticity(material_phase(1,opposite_ip,opposite_el)) /= PLASTICITY_NONLOCAL_ID) &
+        if (phase_plasticity(material_phaseAt(1,opposite_el)) /= PLASTICITY_NONLOCAL_ID) &
           considerLeavingFlux = .false.
       endif
   
@@ -1905,20 +1905,20 @@ subroutine plastic_nonlocal_updateCompatibility(orientation,i,e)
     s2                                                                                              ! slip system index (my neighbor)
   real(pReal), dimension(4) :: &
     absoluteMisorientation                                                                          ! absolute misorientation (without symmetry) between me and my neighbor
-  real(pReal), dimension(2,totalNslip(phase_plasticityInstance(material_phase(1,i,e))),&
-                           totalNslip(phase_plasticityInstance(material_phase(1,i,e))),&
+  real(pReal), dimension(2,totalNslip(phase_plasticityInstance(material_phaseAt(1,e))),&
+                           totalNslip(phase_plasticityInstance(material_phaseAt(1,e))),&
                            nIPneighbors) :: &  
     my_compatibility                                                                                ! my_compatibility for current element and ip  
   real(pReal) :: &
     my_compatibilitySum, &
     thresholdValue, &
     nThresholdValues
-  logical, dimension(totalNslip(phase_plasticityInstance(material_phase(1,i,e)))) :: & 
+  logical, dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,e)))) :: & 
     belowThreshold
   type(rotation) :: rot
 
   Nneighbors = nIPneighbors
-  ph = material_phase(1,i,e)
+  ph = material_phaseAt(1,e)
   textureID = material_texture(1,i,e)
   instance = phase_plasticityInstance(ph)
   ns = totalNslip(instance)
@@ -1950,7 +1950,7 @@ subroutine plastic_nonlocal_updateCompatibility(orientation,i,e)
     !* we consider this to be a real "physical" phase boundary, so completely incompatible.
     !* If one of the two phases has a local plasticity law, 
     !* we do not consider this to be a phase boundary, so completely compatible.
-    neighbor_phase = material_phase(1,neighbor_i,neighbor_e)
+    neighbor_phase = material_phaseAt(1,neighbor_e)
     if (neighbor_phase /= ph) then
       if (.not. phase_localPlasticity(neighbor_phase) .and. .not. phase_localPlasticity(ph))&
         forall(s1 = 1:ns) my_compatibility(1:2,s1,s1,n) = 0.0_pReal
diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90
index a31891573..0a0052dfb 100644
--- a/src/plastic_phenopowerlaw.f90
+++ b/src/plastic_phenopowerlaw.f90
@@ -12,9 +12,8 @@ module plastic_phenopowerlaw
  use material
  use config
  use lattice
-#if defined(PETSc) || defined(DAMASK_HDF5)
+ use discretization
  use results
-#endif
 
  implicit none
  private
@@ -314,7 +313,7 @@ subroutine plastic_phenopowerlaw_init
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
-   NipcMyPhase = count(material_phase == p)
+   NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
    sizeDotState = size(['tau_slip  ','gamma_slip']) * prm%totalNslip &
                 + size(['tau_twin  ','gamma_twin']) * prm%totalNtwin
    sizeState = sizeDotState
diff --git a/src/source_damage_anisoBrittle.f90 b/src/source_damage_anisoBrittle.f90
index ccad7c6b0..9997f81e5 100644
--- a/src/source_damage_anisoBrittle.f90
+++ b/src/source_damage_anisoBrittle.f90
@@ -10,6 +10,7 @@ module source_damage_anisoBrittle
   use IO
   use math
   use material
+  use discretization
   use config
   use lattice
 
@@ -164,7 +165,7 @@ subroutine source_damage_anisoBrittle_init
     end associate
     
     phase = p
-    NofMyPhase=count(material_phase==phase)
+    NofMyPhase=count(material_phaseAt==phase) * discretization_nIP
     instance = source_damage_anisoBrittle_instance(phase)
     sourceOffset = source_damage_anisoBrittle_offset(phase)
 
@@ -202,8 +203,8 @@ subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el)
   real(pReal) :: &
     traction_d, traction_t, traction_n, traction_crit
 
-  phase = phaseAt(ipc,ip,el)
-  constituent = phasememberAt(ipc,ip,el)
+  phase = material_phaseAt(ipc,el)
+  constituent = material_phasememberAt(ipc,ip,el)
   instance = source_damage_anisoBrittle_instance(phase)
   sourceOffset = source_damage_anisoBrittle_offset(phase)
   homog = material_homogenizationAt(el)
diff --git a/src/source_damage_anisoDuctile.f90 b/src/source_damage_anisoDuctile.f90
index aa8c9b8fa..409466e48 100644
--- a/src/source_damage_anisoDuctile.f90
+++ b/src/source_damage_anisoDuctile.f90
@@ -9,6 +9,7 @@ module source_damage_anisoDuctile
   use debug
   use IO
   use math
+  use discretization
   use material
   use config
  
@@ -150,7 +151,7 @@ subroutine source_damage_anisoDuctile_init
     
     phase = p
     
-    NofMyPhase=count(material_phase==phase)
+    NofMyPhase=count(material_phaseAt==phase) * discretization_nIP
     instance = source_damage_anisoDuctile_instance(phase)
     sourceOffset = source_damage_anisoDuctile_offset(phase)
  
@@ -180,8 +181,8 @@ subroutine source_damage_anisoDuctile_dotState(ipc, ip, el)
     instance, &
     i
  
-  phase = phaseAt(ipc,ip,el)
-  constituent = phasememberAt(ipc,ip,el)
+  phase = material_phaseAt(ipc,el)
+  constituent = material_phasememberAt(ipc,ip,el)
   instance = source_damage_anisoDuctile_instance(phase)
   sourceOffset = source_damage_anisoDuctile_offset(phase)
   homog = material_homogenizationAt(el)
diff --git a/src/source_damage_isoBrittle.f90 b/src/source_damage_isoBrittle.f90
index cf43fdfb8..89f5a038c 100644
--- a/src/source_damage_isoBrittle.f90
+++ b/src/source_damage_isoBrittle.f90
@@ -9,6 +9,7 @@ module source_damage_isoBrittle
   use debug
   use IO
   use math
+  use discretization
   use material
   use config
 
@@ -133,7 +134,7 @@ subroutine source_damage_isoBrittle_init
     
    phase = p
    
-   NofMyPhase=count(material_phase==phase)
+   NofMyPhase = count(material_phaseAt==phase) * discretization_nIP
    instance = source_damage_isoBrittle_instance(phase)
    sourceOffset = source_damage_isoBrittle_offset(phase)
       
@@ -164,8 +165,8 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
     strain(6), &
     strainenergy
 
-  phase = phaseAt(ipc,ip,el)                                                                        !< phase ID at ipc,ip,el
-  constituent = phasememberAt(ipc,ip,el)                                                            !< state array offset for phase ID at ipc,ip,el
+  phase = material_phaseAt(ipc,el)                                                                        !< phase ID at ipc,ip,el
+  constituent = material_phasememberAt(ipc,ip,el)                                                            !< state array offset for phase ID at ipc,ip,el
   ! ToDo: capability for multiple instances of SAME source within given phase. Needs Ninstance loop from here on!
   instance = source_damage_isoBrittle_instance(phase)                                               !< instance of damage_isoBrittle source
   sourceOffset = source_damage_isoBrittle_offset(phase)
diff --git a/src/source_damage_isoDuctile.f90 b/src/source_damage_isoDuctile.f90
index 524936077..65930cd07 100644
--- a/src/source_damage_isoDuctile.f90
+++ b/src/source_damage_isoDuctile.f90
@@ -8,6 +8,7 @@ module source_damage_isoDuctile
  use prec
  use debug
  use IO
+ use discretization
  use material
  use config
 
@@ -132,7 +133,7 @@ subroutine source_damage_isoDuctile_init
    end associate
    
    phase = p
-   NofMyPhase=count(material_phase==phase)
+   NofMyPhase=count(material_phaseAt==phase) * discretization_nIP
    instance = source_damage_isoDuctile_instance(phase)
    sourceOffset = source_damage_isoDuctile_offset(phase)
 
@@ -157,8 +158,8 @@ subroutine source_damage_isoDuctile_dotState(ipc, ip, el)
  integer :: &
    phase, constituent, instance, homog, sourceOffset, damageOffset
 
- phase = phaseAt(ipc,ip,el)
- constituent = phasememberAt(ipc,ip,el)
+ phase = material_phaseAt(ipc,el)
+ constituent = material_phasememberAt(ipc,ip,el)
  instance = source_damage_isoDuctile_instance(phase)
  sourceOffset = source_damage_isoDuctile_offset(phase)
  homog = material_homogenizationAt(el)
diff --git a/src/source_thermal_dissipation.f90 b/src/source_thermal_dissipation.f90
index e8464edd0..9b18efef5 100644
--- a/src/source_thermal_dissipation.f90
+++ b/src/source_thermal_dissipation.f90
@@ -7,6 +7,7 @@
 module source_thermal_dissipation
   use prec
   use debug
+  use discretization
   use material
   use config
  
@@ -75,7 +76,7 @@ subroutine source_thermal_dissipation_init
     if (all(phase_source(:,p) /= SOURCE_THERMAL_DISSIPATION_ID)) cycle
     instance = source_thermal_dissipation_instance(p)
     param(instance)%kappa = config_phase(p)%getFloat('dissipation_coldworkcoeff')
-    NofMyPhase=count(material_phase==p)
+    NofMyPhase = count(material_phaseAt==p) * discretization_nIP
     sourceOffset = source_thermal_dissipation_offset(p)
  
     call material_allocateSourceState(p,sourceOffset,NofMyPhase,0,0,0)
diff --git a/src/source_thermal_externalheat.f90 b/src/source_thermal_externalheat.f90
index 99d9a6f1f..1b9d03529 100644
--- a/src/source_thermal_externalheat.f90
+++ b/src/source_thermal_externalheat.f90
@@ -7,6 +7,7 @@
 module source_thermal_externalheat
   use prec
   use debug
+  use discretization
   use material
   use config
 
@@ -83,7 +84,7 @@ subroutine source_thermal_externalheat_init
     if (all(phase_source(:,p) /= SOURCE_thermal_externalheat_ID)) cycle
     instance = source_thermal_externalheat_instance(p)
     sourceOffset = source_thermal_externalheat_offset(p)
-    NofMyPhase=count(material_phase==p)
+    NofMyPhase = count(material_phaseAt==p) * discretization_nIP
  
     param(instance)%time       = config_phase(p)%getFloats('externalheat_time')
     param(instance)%nIntervals = size(param(instance)%time) - 1
diff --git a/src/thermal_adiabatic.f90 b/src/thermal_adiabatic.f90
index 526c98904..76f519913 100644
--- a/src/thermal_adiabatic.f90
+++ b/src/thermal_adiabatic.f90
@@ -167,8 +167,8 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
   Tdot = 0.0_pReal
   dTdot_dT = 0.0_pReal
   do grain = 1, homogenization_Ngrains(homog)
-    phase = phaseAt(grain,ip,el)
-    constituent = phasememberAt(grain,ip,el)
+    phase = material_phaseAt(grain,el)
+    constituent = material_phasememberAt(grain,ip,el)
     do source = 1, phase_Nsources(phase)
       select case(phase_source(source,phase))                                                   
         case (SOURCE_thermal_dissipation_ID)
diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90
index a31961dc7..c70c59b60 100644
--- a/src/thermal_conduction.f90
+++ b/src/thermal_conduction.f90
@@ -132,8 +132,8 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
   Tdot = 0.0_pReal
   dTdot_dT = 0.0_pReal
   do grain = 1, homogenization_Ngrains(homog)
-    phase = phaseAt(grain,ip,el)
-    constituent = phasememberAt(grain,ip,el)
+    phase = material_phaseAt(grain,el)
+    constituent = material_phasememberAt(grain,ip,el)
     do source = 1, phase_Nsources(phase)
       select case(phase_source(source,phase))                                                   
         case (SOURCE_thermal_dissipation_ID)