Merge branch 'development' into python-improvements

CMakeLists.txt

@@ -4,7 +4,7 @@ include (FindPkgConfig REQUIRED)
 # Dummy project to determine compiler names and version
 project (Prerequisites LANGUAGES)
 set(ENV{PKG_CONFIG_PATH} "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/pkgconfig")
-pkg_search_module (PETSC REQUIRED PETSc>3.12.0)
+pkg_check_modules (PETSC REQUIRED PETSc>=3.12.0 PETSc<3.15.0)
 pkg_get_variable (CMAKE_Fortran_COMPILER PETSc fcompiler)
 pkg_get_variable (CMAKE_C_COMPILER PETSc ccompiler)
 

VERSION

@@ -1 +1 @@
-v3.0.0-alpha2-97-g10bbeb561
+v3.0.0-alpha2-173-g584c7cc3a

examples/FEM/polyXtal/material.yaml

@@ -5,8 +5,8 @@ homogenization:
 
 phase:
   Aluminum:
+    lattice: cF
     mechanics:
-      lattice: cF
       output: [F, P, F_e, F_p, L_p]
       elasticity: {C_11: 106.75e9, C_12: 60.41e9, C_44: 28.34e9, type: hooke}
       plasticity:

python/damask/_rotation.py

@@ -367,7 +367,7 @@ class Rotation:
             """Intermediate representation supporting quaternion averaging."""
             return np.einsum('...i,...j',quat,quat)
 
-        if not weights:
+        if weights is None:
             weights = np.ones(self.shape,dtype=float)
 
         eig, vec = np.linalg.eig(np.sum(_M(self.quaternion) * weights[...,np.newaxis,np.newaxis],axis=-3) \

python/tests/test_Rotation.py

@@ -825,6 +825,14 @@ class TestRotation:
         print(f'append 3x {shape} --> {s.shape}')
         assert np.logical_and(s[0,...] == r[0,...], s[-1,...] == p[-1,...]).all()
 
+    @pytest.mark.parametrize('shape',[None,1,(1,),(4,2),(3,3,2)])
+    def test_append_list(self,shape):
+        r = Rotation.from_random(shape=shape)
+        p = Rotation.from_random(shape=shape)
+        s = r.append([r,p])
+        print(f'append 3x {shape} --> {s.shape}')
+        assert s[0,...] == r[0,...] and s[-1,...] == p[-1,...]
+
     @pytest.mark.parametrize('quat,standardized',[
                                                   ([-1,0,0,0],[1,0,0,0]),
                                                   ([-0.5,-0.5,-0.5,-0.5],[0.5,0.5,0.5,0.5]),

src/CPFEM.f90

@@ -5,7 +5,6 @@
 !--------------------------------------------------------------------------------------------------
 module CPFEM
   use prec
-  use FEsolving
   use math
   use rotations
   use YAML_types
@@ -197,11 +196,9 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
       CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_eye(6)
 
     else validCalculation
-      FEsolving_execElem = elCP
-      FEsolving_execIP   = ip
       if (debugCPFEM%extensive) &
         print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip
-      call materialpoint_stressAndItsTangent(dt)
+      call materialpoint_stressAndItsTangent(dt,[ip,ip],[elCP,elCP])
 
       terminalIllness: if (terminallyIll) then
 
@@ -260,7 +257,7 @@ end subroutine CPFEM_general
 !--------------------------------------------------------------------------------------------------
 subroutine CPFEM_forward
 
-  call crystallite_forward
+  call homogenization_forward
   call constitutive_forward
 
 end subroutine CPFEM_forward
@@ -277,7 +274,6 @@ subroutine CPFEM_results(inc,time)
   call results_openJobFile
   call results_addIncrement(inc,time)
   call constitutive_results
-  call crystallite_results
   call homogenization_results
   call discretization_results
   call results_finalizeIncrement

src/CPFEM2.f90

@@ -6,7 +6,6 @@
 module CPFEM2
   use prec
   use config
-  use FEsolving
   use math
   use rotations
   use YAML_types
@@ -97,7 +96,7 @@ end subroutine CPFEM_restartWrite
 !--------------------------------------------------------------------------------------------------
 subroutine CPFEM_forward
 
-  call crystallite_forward
+  call homogenization_forward
   call constitutive_forward
 
 end subroutine CPFEM_forward
@@ -114,7 +113,6 @@ subroutine CPFEM_results(inc,time)
   call results_openJobFile
   call results_addIncrement(inc,time)
   call constitutive_results
-  call crystallite_results
   call homogenization_results
   call discretization_results
   call results_finalizeIncrement

src/DAMASK_interface.f90

@@ -54,12 +54,6 @@ subroutine DAMASK_interface_init
 ===================================================================================================
 --  WRONG PETSc VERSION --- WRONG PETSc VERSION --- WRONG PETSc VERSION ---  WRONG PETSc VERSION --
 ===================================================================================================
-============   THIS VERSION OF DAMASK REQUIRES A DIFFERENT PETSc VERSION   ========================
-===============   THIS VERSION OF DAMASK REQUIRES A DIFFERENT PETSc VERSION   =====================
-==================   THIS VERSION OF DAMASK REQUIRES A DIFFERENT PETSc VERSION   ==================
-===================================================================================================
---  WRONG PETSc VERSION --- WRONG PETSc VERSION --- WRONG PETSc VERSION ---  WRONG PETSc VERSION --
-===================================================================================================
 #endif
 
   character(len=pPathLen*3+pStringLen) :: &

src/DAMASK_marc.f90

@@ -176,7 +176,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
   use DAMASK_interface
   use config
   use YAML_types
-  use FEsolving
   use discretization_marc
   use homogenization
   use CPFEM

src/FEsolving.f90

@@ -1,15 +0,0 @@
-!--------------------------------------------------------------------------------------------------
-!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
-!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief global variables for flow control
-!--------------------------------------------------------------------------------------------------
-module FEsolving
-
-  implicit none
-  public
-
-  integer, dimension(2) :: &
-    FEsolving_execElem, &                                                                           !< for ping-pong scheme always whole range, otherwise one specific element
-    FEsolving_execIP                                                                                !< for ping-pong scheme always range to max IP, otherwise one specific IP
-
-end module FEsolving

src/commercialFEM_fileList.f90

@@ -11,9 +11,7 @@
 #include "config.f90"
 #include "LAPACK_interface.f90"
 #include "math.f90"
-#include "quaternions.f90"
 #include "rotations.f90"
-#include "FEsolving.f90"
 #include "element.f90"
 #include "HDF5_utilities.f90"
 #include "results.f90"

src/constitutive_damage.f90

@@ -1,5 +1,5 @@
 !----------------------------------------------------------------------------------------------------
-!> @brief internal microstructure state for all damage sources and kinematics constitutive models  
+!> @brief internal microstructure state for all damage sources and kinematics constitutive models
 !----------------------------------------------------------------------------------------------------
 submodule(constitutive) constitutive_damage
 
@@ -8,7 +8,7 @@ submodule(constitutive) constitutive_damage
   module function source_damage_anisoBrittle_init(source_length) result(mySources)
     integer, intent(in) :: source_length
     logical, dimension(:,:), allocatable :: mySources
-  end function source_damage_anisoBrittle_init 
+  end function source_damage_anisoBrittle_init
 
   module function source_damage_anisoDuctile_init(source_length) result(mySources)
     integer, intent(in) :: source_length
@@ -23,7 +23,7 @@ submodule(constitutive) constitutive_damage
   module function source_damage_isoDuctile_init(source_length) result(mySources)
     integer, intent(in) :: source_length
     logical, dimension(:,:), allocatable :: mySources
-  end function source_damage_isoDuctile_init   
+  end function source_damage_isoDuctile_init
 
   module function kinematics_cleavage_opening_init(kinematics_length) result(myKinematics)
     integer, intent(in) :: kinematics_length
@@ -39,14 +39,14 @@ submodule(constitutive) constitutive_damage
   module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
     integer, intent(in) :: &
       phase, &                                                                                      !< phase ID of element
-      constituent                                                                                   !< position of element within its phase instance 
+      constituent                                                                                   !< position of element within its phase instance
     real(pReal),  intent(in) :: &
-      phi                                                                                           !< damage parameter 
+      phi                                                                                           !< damage parameter
     real(pReal),  intent(out) :: &
       localphiDot, &
       dLocalphiDot_dPhi
   end subroutine source_damage_anisoBrittle_getRateAndItsTangent
- 
+
   module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
     integer, intent(in) :: &
       phase, &                                                                                      !< phase ID of element
@@ -129,7 +129,7 @@ module subroutine damage_init
     allocate(sourceState(ph)%p(phase_Nsources(ph)))
   enddo
 
-  allocate(phase_source(maxval(phase_Nsources),phases%length), source = SOURCE_undefined_ID) 
+  allocate(phase_source(maxval(phase_Nsources),phases%length), source = SOURCE_undefined_ID)
 
 ! initialize source mechanisms
   if(maxval(phase_Nsources) /= 0) then
@@ -141,19 +141,19 @@ module subroutine damage_init
 
 !--------------------------------------------------------------------------------------------------
 ! initialize kinematic mechanisms
-  allocate(phase_Nkinematics(phases%length),source = 0)           
+  allocate(phase_Nkinematics(phases%length),source = 0)
   do ph = 1,phases%length
     phase => phases%get(ph)
     kinematics => phase%get('kinematics',defaultVal=emptyList)
     phase_Nkinematics(ph) = kinematics%length
   enddo
- 
+
-  allocate(phase_kinematics(maxval(phase_Nkinematics),phases%length), source = KINEMATICS_undefined_ID) 
+  allocate(phase_kinematics(maxval(phase_Nkinematics),phases%length), source = KINEMATICS_undefined_ID)
 
   if(maxval(phase_Nkinematics) /= 0) then
     where(kinematics_cleavage_opening_init(maxval(phase_Nkinematics)))  phase_kinematics = KINEMATICS_cleavage_opening_ID
     where(kinematics_slipplane_opening_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_slipplane_opening_ID
-  endif 
+  endif
 
 end subroutine damage_init
 
@@ -167,7 +167,7 @@ module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi
     ip, &                                                                                           !< integration point number
     el                                                                                              !< element number
   real(pReal), intent(in) :: &
-    phi                                                                                             !< damage parameter      
+    phi                                                                                             !< damage parameter
   real(pReal), intent(inout) :: &
     phiDot, &
     dPhiDot_dPhi
@@ -183,7 +183,7 @@ module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi
 
    phiDot = 0.0_pReal
    dPhiDot_dPhi = 0.0_pReal
- 
+
    do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el))
      phase = material_phaseAt(grain,el)
      constituent = material_phasememberAt(grain,ip,el)
@@ -217,32 +217,35 @@ end subroutine constitutive_damage_getRateAndItsTangents
 !----------------------------------------------------------------------------------------------
 !< @brief writes damage sources results to HDF5 output file
 !----------------------------------------------------------------------------------------------
-module subroutine damage_results
+module subroutine damage_results(group,ph)
 
-  integer :: p,i
+  character(len=*), intent(in) :: group
-  character(len=pStringLen) :: group
+  integer,          intent(in) :: ph
 
-  do p = 1, size(material_name_phase)
+  integer :: so
 
-    sourceLoop: do i = 1, phase_Nsources(p)
+  sourceLoop: do so = 1, phase_Nsources(ph)
-    group = trim('current/phase')//'/'//trim(material_name_phase(p))
-    group = trim(group)//'/sources'
-    call results_closeGroup(results_addGroup(group))
 
-      sourceType: select case (phase_source(i,p))
+  if (phase_source(so,ph) /= SOURCE_UNDEFINED_ID) &
+    call results_closeGroup(results_addGroup(group//'sources/')) ! should be 'damage'
 
-        case (SOURCE_damage_anisoBrittle_ID) sourceType
+    sourceType: select case (phase_source(so,ph))
-          call source_damage_anisoBrittle_results(p,group)
-        case (SOURCE_damage_anisoDuctile_ID) sourceType
-          call source_damage_anisoDuctile_results(p,group)
-        case (SOURCE_damage_isoBrittle_ID) sourceType
-          call source_damage_isoBrittle_results(p,group)
-        case (SOURCE_damage_isoDuctile_ID) sourceType
-          call source_damage_isoDuctile_results(p,group)
-      end select sourceType
 
-    enddo SourceLoop
+      case (SOURCE_damage_anisoBrittle_ID) sourceType
-  enddo
+        call source_damage_anisoBrittle_results(ph,group//'sources/')
+
+      case (SOURCE_damage_anisoDuctile_ID) sourceType
+        call source_damage_anisoDuctile_results(ph,group//'sources/')
+
+      case (SOURCE_damage_isoBrittle_ID) sourceType
+        call source_damage_isoBrittle_results(ph,group//'sources/')
+
+      case (SOURCE_damage_isoDuctile_ID) sourceType
+        call source_damage_isoDuctile_results(ph,group//'sources/')
+
+    end select sourceType
+
+  enddo SourceLoop
 
 end subroutine damage_results
 

src/constitutive_plastic_dislotwin.f90

@@ -485,12 +485,12 @@ end function plastic_dislotwin_init
 !--------------------------------------------------------------------------------------------------
 !> @brief Return the homogenized elasticity matrix.
 !--------------------------------------------------------------------------------------------------
-module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC)
+module function plastic_dislotwin_homogenizedC(co,ip,el) result(homogenizedC)
 
   real(pReal), dimension(6,6) :: &
     homogenizedC
   integer,     intent(in) :: &
-    ipc, &                                                                                          !< component-ID of integration point
+    co, &                                                                                          !< component-ID of integration point
     ip, &                                                                                           !< integration point
     el                                                                                              !< element
 
@@ -498,9 +498,9 @@ module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC)
              of
   real(pReal) :: f_unrotated
 
-  of = material_phasememberAt(ipc,ip,el)
+  of = material_phasememberAt(co,ip,el)
-  associate(prm => param(phase_plasticityInstance(material_phaseAt(ipc,el))),&
+  associate(prm => param(phase_plasticityInstance(material_phaseAt(co,el))),&
-            stt => state(phase_plasticityInstance(material_phaseAT(ipc,el))))
+            stt => state(phase_plasticityInstance(material_phaseAT(co,el))))
 
   f_unrotated = 1.0_pReal &
               - sum(stt%f_tw(1:prm%sum_N_tw,of)) &

src/constitutive_plastic_nonlocal.f90

@@ -552,11 +552,10 @@ end function plastic_nonlocal_init
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates quantities characterizing the microstructure
 !--------------------------------------------------------------------------------------------------
-module subroutine plastic_nonlocal_dependentState(F, Fp, instance, of, ip, el)
+module subroutine plastic_nonlocal_dependentState(F, instance, of, ip, el)
 
   real(pReal), dimension(3,3), intent(in) :: &
-    F, &
+    F
-    Fp
   integer, intent(in) :: &
     instance, &
     of, &
@@ -564,6 +563,8 @@ module subroutine plastic_nonlocal_dependentState(F, Fp, instance, of, ip, el)
     el
 
   integer :: &
+    ph, &
+    me, &
     no, &                                                                                           !< neighbor offset
     neighbor_el, &                                                                                  ! element number of neighboring material point
     neighbor_ip, &                                                                                  ! integration point of neighboring material point
@@ -643,8 +644,10 @@ module subroutine plastic_nonlocal_dependentState(F, Fp, instance, of, ip, el)
 
   rho0 = getRho0(instance,of,ip,el)
   if (.not. phase_localPlasticity(material_phaseAt(1,el)) .and. prm%shortRangeStressCorrection) then
-    invFp = math_inv33(Fp)
+    ph = material_phaseAt(1,el)
-    invFe = matmul(Fp,math_inv33(F))
+    me = material_phaseMemberAt(1,ip,el)
+    invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me))
+    invFe = matmul(constitutive_mech_Fp(ph)%data(1:3,1:3,me),math_inv33(F))
 
     rho_edg_delta = rho0(:,mob_edg_pos) - rho0(:,mob_edg_neg)
     rho_scr_delta = rho0(:,mob_scr_pos) - rho0(:,mob_scr_neg)
@@ -973,14 +976,13 @@ end subroutine plastic_nonlocal_deltaState
 !---------------------------------------------------------------------------------------------------
 !> @brief calculates the rate of change of microstructure
 !---------------------------------------------------------------------------------------------------
-module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, &
+module subroutine plastic_nonlocal_dotState(Mp, F, Temperature,timestep, &
                                             instance,of,ip,el)
 
   real(pReal), dimension(3,3), intent(in) :: &
     Mp                                                                                              !< MandelStress
   real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: &
-    F, &                                                                                            !< elastic deformation gradient
+    F                                                                                               !< Deformation gradient
-    Fp                                                                                              !< plastic deformation gradient
   real(pReal), intent(in) :: &
     Temperature, &                                                                                  !< temperature
     timestep                                                                                        !< substepped crystallite time increment
@@ -1147,7 +1149,7 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, &
                                          - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) &
                                                                   - rhoDotSingle2DipoleGlide(s,9))  ! make sure that we do not annihilate more dipoles than we have
 
-  rhoDot = rhoDotFlux(F,Fp,timestep,  instance,of,ip,el) &
+  rhoDot = rhoDotFlux(F,timestep,  instance,of,ip,el) &
          + rhoDotMultiplication &
          + rhoDotSingle2DipoleGlide &
          + rhoDotAthermalAnnihilation &
@@ -1176,11 +1178,10 @@ end subroutine plastic_nonlocal_dotState
 !---------------------------------------------------------------------------------------------------
 !> @brief calculates the rate of change of microstructure
 !---------------------------------------------------------------------------------------------------
-function rhoDotFlux(F,Fp,timestep,  instance,of,ip,el)
+function rhoDotFlux(F,timestep,  instance,of,ip,el)
 
   real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: &
-    F, &                                                                                            !< elastic deformation gradient
+    F                                                                                               !< Deformation gradient
-    Fp                                                                                              !< plastic deformation gradient
   real(pReal), intent(in) :: &
     timestep                                                                                        !< substepped crystallite time increment
   integer, intent(in) :: &
@@ -1293,7 +1294,7 @@ function rhoDotFlux(F,Fp,timestep,  instance,of,ip,el)
     m(1:3,:,4) =  prm%slip_transverse
 
     my_F = F(1:3,1:3,1,ip,el)
-    my_Fe = matmul(my_F, math_inv33(Fp(1:3,1:3,1,ip,el)))
+    my_Fe = matmul(my_F, math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,of)))
 
     neighbors: do n = 1,nIPneighbors
 
@@ -1311,7 +1312,7 @@ function rhoDotFlux(F,Fp,timestep,  instance,of,ip,el)
       if (neighbor_n > 0) then                                                                      ! if neighbor exists, average deformation gradient
         neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el))
         neighbor_F = F(1:3,1:3,1,neighbor_ip,neighbor_el)
-        neighbor_Fe = matmul(neighbor_F, math_inv33(Fp(1:3,1:3,1,neighbor_ip,neighbor_el)))
+        neighbor_Fe = matmul(neighbor_F, math_inv33(constitutive_mech_Fp(np)%data(1:3,1:3,no)))
         Favg = 0.5_pReal * (my_F + neighbor_F)
       else                                                                                          ! if no neighbor, take my value as average
         Favg = my_F

src/damage_nonlocal.f90

@@ -148,12 +148,12 @@ real(pReal) function damage_nonlocal_getMobility(ip,el)
     ip, &                                                                                           !< integration point number
     el                                                                                              !< element number
   integer :: &
-    ipc
+    co
 
   damage_nonlocal_getMobility = 0.0_pReal
 
-  do ipc = 1, homogenization_Nconstituents(material_homogenizationAt(el))
+  do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
-    damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_M(material_phaseAt(ipc,el))
+    damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_M(material_phaseAt(co,el))
   enddo
 
   damage_nonlocal_getMobility = damage_nonlocal_getMobility/&

src/grid/discretization_grid.f90

@@ -19,7 +19,6 @@ module discretization_grid
   use results
   use discretization
   use geometry_plastic_nonlocal
-  use FEsolving
 
   implicit none
   private
@@ -117,9 +116,6 @@ subroutine discretization_grid_init(restart)
                                  (grid(1)+1) * (grid(2)+1) *  grid3,&                               ! ...unless not last process
                                  worldrank+1==worldsize))
 
-  FEsolving_execElem = [1,product(myGrid)]                                                          ! parallel loop bounds set to comprise all elements
-  FEsolving_execIP   = [1,1]                                                                        ! parallel loop bounds set to comprise the only IP
-
 !--------------------------------------------------------------------------------------------------
 ! store geometry information for post processing
   if(.not. restart) then

src/grid/grid_mech_FEM.f90

@@ -18,7 +18,6 @@ module grid_mech_FEM
   use math
   use rotations
   use spectral_utilities
-  use FEsolving
   use config
   use homogenization
   use discretization

src/grid/grid_mech_spectral_basic.f90

@@ -18,7 +18,6 @@ module grid_mech_spectral_basic
   use math
   use rotations
   use spectral_utilities
-  use FEsolving
   use config
   use homogenization
   use discretization_grid

src/grid/grid_mech_spectral_polarisation.f90

@@ -18,7 +18,6 @@ module grid_mech_spectral_polarisation
   use math
   use rotations
   use spectral_utilities
-  use FEsolving
   use config
   use homogenization
   use discretization_grid

src/grid/spectral_utilities.f90

@@ -810,9 +810,9 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
   print'(/,a)', ' ... evaluating constitutive response ......................................'
   flush(IO_STDOUT)
 
-  homogenization_F  = reshape(F,[3,3,product(grid(1:2))*grid3])                                    ! set materialpoint target F to estimated field
+  homogenization_F  = reshape(F,[3,3,product(grid(1:2))*grid3])                                     ! set materialpoint target F to estimated field
 
-  call materialpoint_stressAndItsTangent(timeinc)                                                   ! calculate P field
+  call materialpoint_stressAndItsTangent(timeinc,[1,1],[1,product(grid(1:2))*grid3])                ! calculate P field
 
   P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3])
   P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt                                                   ! average of P

src/homogenization.f90

@@ -11,7 +11,6 @@ module homogenization
   use math
   use material
   use constitutive
-  use FEsolving
   use discretization
   use thermal_isothermal
   use thermal_conduction
@@ -48,20 +47,6 @@ module homogenization
 
   type(tNumerics) :: num
 
-  type :: tDebugOptions
-    logical :: &
-      basic, &
-      extensive, &
-      selective
-    integer :: &
-      element, &
-      ip, &
-      grain
-  end type tDebugOptions
-
-  type(tDebugOptions) :: debugHomog
-
-
 !--------------------------------------------------------------------------------------------------
   interface
 
@@ -85,33 +70,27 @@ module homogenization
     end subroutine mech_homogenize
 
     module subroutine mech_results(group_base,h)
-
       character(len=*), intent(in) :: group_base
       integer, intent(in)          :: h
-
     end subroutine mech_results
 
-! -------- ToDo ---------------------------------------------------------
+    module function mech_updateState(subdt,subF,ip,el) result(doneAndHappy)
-    module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
+      real(pReal), intent(in) :: &
-      logical, dimension(2) :: mech_RGC_updateState
+        subdt                                                                                           !< current time step
-      real(pReal), dimension(:,:,:),     intent(in)    :: &
+      real(pReal), intent(in), dimension(3,3) :: &
-        P,&                                                                                         !< partitioned stresses
+        subF
-        F,&                                                                                         !< partitioned deformation gradients
+      integer,     intent(in) :: &
-        F0                                                                                          !< partitioned initial deformation gradients
+        ip, &                                                                                           !< integration point
-      real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF                                         !< partitioned stiffnesses
+        el                                                                                              !< element number
-      real(pReal), dimension(3,3),       intent(in) :: avgF                                         !< average F
+      logical, dimension(2) :: doneAndHappy
-      real(pReal),                       intent(in) :: dt                                           !< time increment
+    end function mech_updateState
-      integer,                           intent(in) :: &
-        ip, &                                                                                       !< integration point number
-        el                                                                                          !< element number
-    end function mech_RGC_updateState
 
   end interface
-! -----------------------------------------------------------------------
 
   public ::  &
     homogenization_init, &
     materialpoint_stressAndItsTangent, &
+    homogenization_forward, &
     homogenization_results
 
 contains
@@ -124,24 +103,10 @@ subroutine homogenization_init
 
   class (tNode) , pointer :: &
     num_homog, &
-    num_homogGeneric, &
+    num_homogGeneric
-    debug_homogenization
 
   print'(/,a)', ' <<<+-  homogenization init  -+>>>'; flush(IO_STDOUT)
 
-  debug_homogenization => config_debug%get('homogenization', defaultVal=emptyList)
-  debugHomog%basic       =  debug_homogenization%contains('basic')
-  debugHomog%extensive   =  debug_homogenization%contains('extensive')
-  debugHomog%selective   =  debug_homogenization%contains('selective')
-  debugHomog%element     =  config_debug%get_asInt('element',defaultVal = 1)
-  debugHomog%ip          =  config_debug%get_asInt('integrationpoint',defaultVal = 1)
-  debugHomog%grain       =  config_debug%get_asInt('grain',defaultVal = 1)
-
-  if (debugHomog%grain < 1 &
-    .or. debugHomog%grain > homogenization_Nconstituents(material_homogenizationAt(debugHomog%element))) &
-    call IO_error(602,ext_msg='constituent', el=debugHomog%element, g=debugHomog%grain)
-
-
   num_homog        => config_numerics%get('homogenization',defaultVal=emptyDict)
   num_homogGeneric => num_homog%get('generic',defaultVal=emptyDict)
 
@@ -171,178 +136,130 @@ end subroutine homogenization_init
 !--------------------------------------------------------------------------------------------------
 !> @brief  parallelized calculation of stress and corresponding tangent at material points
 !--------------------------------------------------------------------------------------------------
-subroutine materialpoint_stressAndItsTangent(dt)
+subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execElem)
 
   real(pReal), intent(in) :: dt                                                                     !< time increment
+  integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP
   integer :: &
-    NiterationHomog, &
     NiterationMPstate, &
-    i, &                                                                                            !< integration point number
+    ip, &                                                                                            !< integration point number
-    e, &                                                                                            !< element number
+    el, &                                                                                            !< element number
-    myNgrains
+    myNgrains, co, ce, ho, me
-  real(pReal), dimension(discretization_nIPs,discretization_Nelems) :: &
+  real(pReal) :: &
     subFrac, &
     subStep
-  logical,     dimension(discretization_nIPs,discretization_Nelems) :: &
+  logical :: &
-    requested, &
     converged
-  logical,     dimension(2,discretization_nIPs,discretization_Nelems) :: &
+  logical, dimension(2) :: &
     doneAndHappy
-  integer :: m
 
 
+  !$OMP PARALLEL DO PRIVATE(ce,me,ho,myNgrains,NiterationMPstate,subFrac,converged,subStep,doneAndHappy)
+  do el = FEsolving_execElem(1),FEsolving_execElem(2)
+    ho = material_homogenizationAt(el)
+    myNgrains = homogenization_Nconstituents(ho)
+    do ip = FEsolving_execIP(1),FEsolving_execIP(2)
+      me = material_homogenizationMemberAt(ip,el)
 !--------------------------------------------------------------------------------------------------
 ! initialize restoration points
-  do e = FEsolving_execElem(1),FEsolving_execElem(2)
+      call constitutive_initializeRestorationPoints(ip,el)
-    do i = FEsolving_execIP(1),FEsolving_execIP(2);
 
-      call crystallite_initializeRestorationPoints(i,e)
+      subFrac = 0.0_pReal
+      converged = .false.                                                                           ! pretend failed step ...
+      subStep = 1.0_pReal/num%subStepSizeHomog                                                      ! ... larger then the requested calculation
 
-      subFrac(i,e) = 0.0_pReal
+      if (homogState(ho)%sizeState > 0) &
-      converged(i,e) = .false.                                                                      ! pretend failed step ...
+        homogState(ho)%subState0(:,me) = homogState(ho)%State0(:,me)
-      subStep(i,e) = 1.0_pReal/num%subStepSizeHomog                                                 ! ... larger then the requested calculation
+      if (damageState(ho)%sizeState > 0) &
-      requested(i,e) = .true.                                                                       ! everybody requires calculation
+        damageState(ho)%subState0(:,me) = damageState(ho)%State0(:,me)
 
-      if (homogState(material_homogenizationAt(e))%sizeState > 0) &
+      cutBackLooping: do while (.not. terminallyIll .and. subStep  > num%subStepMinHomog)
-          homogState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = &
-          homogState(material_homogenizationAt(e))%State0(   :,material_homogenizationMemberAt(i,e))
 
-      if (damageState(material_homogenizationAt(e))%sizeState > 0) &
+        if (converged) then
-          damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = &
+          subFrac = subFrac + subStep
-          damageState(material_homogenizationAt(e))%State0(   :,material_homogenizationMemberAt(i,e))
+          subStep = min(1.0_pReal-subFrac,num%stepIncreaseHomog*subStep)             ! introduce flexibility for step increase/acceleration
-    enddo
-  enddo
 
-  NiterationHomog = 0
+          steppingNeeded: if (subStep > num%subStepMinHomog) then
-
-  cutBackLooping: do while (.not. terminallyIll .and. &
-       any(subStep(FEsolving_execIP(1):FEsolving_execIP(2),&
-                   FEsolving_execElem(1):FEsolving_execElem(2)) > num%subStepMinHomog))
-
-    !$OMP PARALLEL DO PRIVATE(m)
-    elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2)
-      myNgrains = homogenization_Nconstituents(material_homogenizationAt(e))
-      IpLooping1: do i = FEsolving_execIP(1),FEsolving_execIP(2)
-
-        if (converged(i,e)) then
-          subFrac(i,e) = subFrac(i,e) + subStep(i,e)
-          subStep(i,e) = min(1.0_pReal-subFrac(i,e),num%stepIncreaseHomog*subStep(i,e))             ! introduce flexibility for step increase/acceleration
-
-          steppingNeeded: if (subStep(i,e) > num%subStepMinHomog) then
 
             ! wind forward grain starting point
-            call crystallite_windForward(i,e)
+            call constitutive_windForward(ip,el)
 
-            if(homogState(material_homogenizationAt(e))%sizeState > 0) &
+            if(homogState(ho)%sizeState > 0) &
-                homogState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = &
+              homogState(ho)%subState0(:,me) = homogState(ho)%State(:,me)
-                homogState(material_homogenizationAt(e))%State    (:,material_homogenizationMemberAt(i,e))
+            if(damageState(ho)%sizeState > 0) &
-            if(damageState(material_homogenizationAt(e))%sizeState > 0) &
+              damageState(ho)%subState0(:,me) = damageState(ho)%State(:,me)
-                damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = &
-                damageState(material_homogenizationAt(e))%State    (:,material_homogenizationMemberAt(i,e))
 
           endif steppingNeeded
-
+        elseif ( (myNgrains == 1 .and. subStep <= 1.0 ) .or. &                                   ! single grain already tried internal subStepping in crystallite
-        else
+               num%subStepSizeHomog * subStep <=  num%subStepMinHomog ) then                   ! would require too small subStep
-          if ( (myNgrains == 1 .and. subStep(i,e) <= 1.0 ) .or. &                                   ! single grain already tried internal subStepping in crystallite
-               num%subStepSizeHomog * subStep(i,e) <=  num%subStepMinHomog ) then                   ! would require too small subStep
                                                                                                     ! cutback makes no sense
-            if (.not. terminallyIll) then                                                           ! so first signals terminally ill...
+          if (.not. terminallyIll) &                                                           ! so first signals terminally ill...
-              print*, ' Integration point ', i,' at element ', e, ' terminally ill'
+            print*, ' Integration point ', ip,' at element ', el, ' terminally ill'
-            endif
+          terminallyIll = .true.                                                                  ! ...and kills all others
-            terminallyIll = .true.                                                                  ! ...and kills all others
+        else                                                                                      ! cutback makes sense
-          else                                                                                      ! cutback makes sense
+          subStep = num%subStepSizeHomog * subStep                                      ! crystallite had severe trouble, so do a significant cutback
-            subStep(i,e) = num%subStepSizeHomog * subStep(i,e)                                      ! crystallite had severe trouble, so do a significant cutback
 
-            call crystallite_restore(i,e,subStep(i,e) < 1.0_pReal)
+          call constitutive_restore(ip,el,subStep < 1.0_pReal)
-            call constitutive_restore(i,e)
 
-            if(homogState(material_homogenizationAt(e))%sizeState > 0) &
+          if(homogState(ho)%sizeState > 0) &
-                homogState(material_homogenizationAt(e))%State(    :,material_homogenizationMemberAt(i,e)) = &
+            homogState(ho)%State(:,me) = homogState(ho)%subState0(:,me)
-                homogState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e))
+          if(damageState(ho)%sizeState > 0) &
-            if(damageState(material_homogenizationAt(e))%sizeState > 0) &
+            damageState(ho)%State(:,me) = damageState(ho)%subState0(:,me)
-                damageState(material_homogenizationAt(e))%State(    :,material_homogenizationMemberAt(i,e)) = &
-                damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e))
-          endif
         endif
 
-        if (subStep(i,e) > num%subStepMinHomog) then
+        if (subStep > num%subStepMinHomog) doneAndHappy = [.false.,.true.]
-          requested(i,e) = .true.
-          doneAndHappy(1:2,i,e) = [.false.,.true.]
-        endif
-      enddo IpLooping1
-    enddo elementLooping1
-    !$OMP END PARALLEL DO
 
-    NiterationMPstate = 0
+        NiterationMPstate = 0
-
+        convergenceLooping: do while (.not. terminallyIll &
-    convergenceLooping: do while (.not. terminallyIll .and. &
+                       .and. .not. doneAndHappy(1) &
-              any(                 requested(:,FEsolving_execELem(1):FEsolving_execElem(2)) &
+                       .and. NiterationMPstate < num%nMPstate)
-                  .and. .not. doneAndHappy(1,:,FEsolving_execELem(1):FEsolving_execElem(2)) &
+          NiterationMPstate = NiterationMPstate + 1
-                 ) .and. &
-              NiterationMPstate < num%nMPstate)
-      NiterationMPstate = NiterationMPstate + 1
 
 !--------------------------------------------------------------------------------------------------
 ! deformation partitioning
-      !$OMP PARALLEL DO PRIVATE(myNgrains,m)
-      elementLooping2: do e = FEsolving_execElem(1),FEsolving_execElem(2)
-        myNgrains = homogenization_Nconstituents(material_homogenizationAt(e))
-        IpLooping2: do i = FEsolving_execIP(1),FEsolving_execIP(2)
-          if(requested(i,e) .and. .not. doneAndHappy(1,i,e)) then                                   ! requested but not yet done
-            m = (e-1)*discretization_nIPs + i
-            call mech_partition(homogenization_F0(1:3,1:3,m) &
-                                      + (homogenization_F(1:3,1:3,m)-homogenization_F0(1:3,1:3,m))&
-                                         *(subStep(i,e)+subFrac(i,e)), &
-                                      i,e)
-            crystallite_dt(1:myNgrains,i,e) = dt*subStep(i,e)                                       ! propagate materialpoint dt to grains
-            crystallite_requested(1:myNgrains,i,e) = .true.                                         ! request calculation for constituents
-          else
-            crystallite_requested(1:myNgrains,i,e) = .false.                                        ! calculation for constituents not required anymore
-          endif
-        enddo IpLooping2
-      enddo elementLooping2
-      !$OMP END PARALLEL DO
 
-!--------------------------------------------------------------------------------------------------
+          if (.not. doneAndHappy(1)) then
-! crystallite integration
+            ce = (el-1)*discretization_nIPs + ip
-      converged = crystallite_stress() !ToDo: MD not sure if that is the best logic
+            call mech_partition(homogenization_F0(1:3,1:3,ce) &
+                                      + (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce))&
+                                         *(subStep+subFrac), &
+                                      ip,el)
+            converged = .true.
+            do co = 1, myNgrains
+              converged = converged .and. crystallite_stress(dt*subStep,co,ip,el)
+            enddo
 
-!--------------------------------------------------------------------------------------------------
+            if (.not. converged) then
-! state update
+              doneAndHappy = [.true.,.false.]
-     !$OMP PARALLEL DO PRIVATE(m)
-      elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2)
-        IpLooping3: do i = FEsolving_execIP(1),FEsolving_execIP(2)
-          if (requested(i,e) .and. .not. doneAndHappy(1,i,e)) then
-            if (.not. converged(i,e)) then
-              doneAndHappy(1:2,i,e) = [.true.,.false.]
             else
-              m = (e-1)*discretization_nIPs + i
+              ce = (el-1)*discretization_nIPs + ip
-              doneAndHappy(1:2,i,e) = updateState(dt*subStep(i,e), &
+              doneAndHappy = mech_updateState(dt*subStep, &
-                                                  homogenization_F0(1:3,1:3,m) &
+                                              homogenization_F0(1:3,1:3,ce) &
-                                                  + (homogenization_F(1:3,1:3,m)-homogenization_F0(1:3,1:3,m)) &
+                                              + (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce)) &
-                                                     *(subStep(i,e)+subFrac(i,e)), &
+                                                     *(subStep+subFrac), &
-                                                 i,e)
+                                              ip,el)
-              converged(i,e) = all(doneAndHappy(1:2,i,e))                                           ! converged if done and happy
+              converged = all(doneAndHappy)
             endif
           endif
-        enddo IpLooping3
-      enddo elementLooping3
-      !$OMP END PARALLEL DO
 
-    enddo convergenceLooping
+        enddo convergenceLooping
-
+      enddo cutBackLooping
-    NiterationHomog = NiterationHomog + 1
+    enddo
-
+  enddo
-  enddo cutBackLooping
+  !$OMP END PARALLEL DO
 
   if (.not. terminallyIll ) then
-    call crystallite_orientations()                                                                 ! calculate crystal orientations
+    !$OMP PARALLEL DO PRIVATE(ho,myNgrains)
-    !$OMP PARALLEL DO
+    elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2)
-    elementLooping4: do e = FEsolving_execElem(1),FEsolving_execElem(2)
+      ho = material_homogenizationAt(el)
-      IpLooping4: do i = FEsolving_execIP(1),FEsolving_execIP(2)
+      myNgrains = homogenization_Nconstituents(ho)
-        call mech_homogenize(i,e)
+      IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2)
-      enddo IpLooping4
+        do co = 1, myNgrains
-    enddo elementLooping4
+          call crystallite_orientations(co,ip,el)
+        enddo
+        call mech_homogenize(ip,el)
+      enddo IpLooping3
+    enddo elementLooping3
     !$OMP END PARALLEL DO
   else
     print'(/,a,/)', ' << HOMOG >> Material Point terminally ill'
@@ -351,77 +268,56 @@ subroutine materialpoint_stressAndItsTangent(dt)
 end subroutine materialpoint_stressAndItsTangent
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief update the internal state of the homogenization scheme and tell whether "done" and
-!> "happy" with result
-!--------------------------------------------------------------------------------------------------
-function updateState(subdt,subF,ip,el)
-
-  real(pReal), intent(in) :: &
-    subdt                                                                                           !< current time step
-  real(pReal), intent(in), dimension(3,3) :: &
-    subF
-  integer,     intent(in) :: &
-    ip, &                                                                                           !< integration point
-    el                                                                                              !< element number
-  integer :: c
-  logical, dimension(2) :: updateState
-  real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
-
-  updateState = .true.
-  chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
-    case (HOMOGENIZATION_RGC_ID) chosenHomogenization
-      do c=1,homogenization_Nconstituents(material_homogenizationAt(el))
-        dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el)
-      enddo
-      updateState = &
-        updateState .and. &
-          mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
-                        crystallite_partitionedF(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
-                        crystallite_partitionedF0(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el),&
-                               subF,&
-                               subdt, &
-                               dPdFs, &
-                               ip, &
-                               el)
-  end select chosenHomogenization
-
-end function updateState
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief writes homogenization results to HDF5 output file
 !--------------------------------------------------------------------------------------------------
 subroutine homogenization_results
-  use material, only: &
-    material_homogenization_type => homogenization_type
 
-  integer :: p
+  integer :: ho
   character(len=:), allocatable :: group_base,group
 
 
-  do p=1,size(material_name_homogenization)
+  call results_closeGroup(results_addGroup('current/homogenization/'))
-    group_base = 'current/homogenization/'//trim(material_name_homogenization(p))
+
+  do ho=1,size(material_name_homogenization)
+    group_base = 'current/homogenization/'//trim(material_name_homogenization(ho))
     call results_closeGroup(results_addGroup(group_base))
 
-    call mech_results(group_base,p)
+    call mech_results(group_base,ho)
 
     group = trim(group_base)//'/damage'
     call results_closeGroup(results_addGroup(group))
-    select case(damage_type(p))
+    select case(damage_type(ho))
       case(DAMAGE_NONLOCAL_ID)
-        call damage_nonlocal_results(p,group)
+        call damage_nonlocal_results(ho,group)
     end select
 
     group = trim(group_base)//'/thermal'
     call results_closeGroup(results_addGroup(group))
-    select case(thermal_type(p))
+    select case(thermal_type(ho))
       case(THERMAL_CONDUCTION_ID)
-        call thermal_conduction_results(p,group)
+        call thermal_conduction_results(ho,group)
     end select
 
  enddo
 
 end subroutine homogenization_results
 
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Forward data after successful increment.
+! ToDo: Any guessing for the current states possible?
+!--------------------------------------------------------------------------------------------------
+subroutine homogenization_forward
+
+  integer :: ho
+
+
+  do ho = 1, size(material_name_homogenization)
+    homogState (ho)%state0 = homogState (ho)%state
+    damageState(ho)%state0 = damageState(ho)%state
+  enddo
+
+end subroutine homogenization_forward
+
 end module homogenization

src/homogenization_mech.f90

@@ -4,6 +4,7 @@
 !--------------------------------------------------------------------------------------------------
 submodule(homogenization) homogenization_mech
 
+
   interface
 
     module subroutine mech_none_init
@@ -51,6 +52,21 @@ submodule(homogenization) homogenization_mech
     end subroutine mech_RGC_averageStressAndItsTangent
 
 
+    module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHappy)
+      logical, dimension(2) :: doneAndHappy
+      real(pReal), dimension(:,:,:),     intent(in)    :: &
+        P,&                                                                                         !< partitioned stresses
+        F,&                                                                                         !< partitioned deformation gradients
+        F0                                                                                          !< partitioned initial deformation gradients
+      real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF                                         !< partitioned stiffnesses
+      real(pReal), dimension(3,3),       intent(in) :: avgF                                         !< average F
+      real(pReal),                       intent(in) :: dt                                           !< time increment
+      integer,                           intent(in) :: &
+        ip, &                                                                                       !< integration point number
+        el                                                                                          !< element number
+    end function mech_RGC_updateState
+
+
     module subroutine mech_RGC_results(instance,group)
       integer,          intent(in) :: instance                                                      !< homogenization instance
       character(len=*), intent(in) :: group                                                         !< group name in HDF5 file
@@ -100,16 +116,16 @@ module subroutine mech_partition(subF,ip,el)
   chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
 
     case (HOMOGENIZATION_NONE_ID) chosenHomogenization
-      crystallite_partitionedF(1:3,1:3,1,ip,el) = subF
+      crystallite_F(1:3,1:3,1,ip,el) = subF
 
     case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
       call mech_isostrain_partitionDeformation(&
-                           crystallite_partitionedF(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
+                           crystallite_F(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
                            subF)
 
     case (HOMOGENIZATION_RGC_ID) chosenHomogenization
       call mech_RGC_partitionDeformation(&
-                          crystallite_partitionedF(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
+                          crystallite_F(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
                           subF,&
                           ip, &
                           el)
@@ -127,35 +143,35 @@ module subroutine mech_homogenize(ip,el)
   integer, intent(in) :: &
        ip, &                                                                                        !< integration point
        el                                                                                           !< element number
-  integer :: c,m
+  integer :: co,ce
   real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
 
 
-  m = (el-1)* discretization_nIPs + ip
+  ce = (el-1)* discretization_nIPs + ip
   chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
 
     case (HOMOGENIZATION_NONE_ID) chosenHomogenization
-        homogenization_P(1:3,1:3,m)            = crystallite_P(1:3,1:3,1,ip,el)
+        homogenization_P(1:3,1:3,ce)            = crystallite_P(1:3,1:3,1,ip,el)
-        homogenization_dPdF(1:3,1:3,1:3,1:3,m) = crystallite_stressTangent(1,ip,el)
+        homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = crystallite_stressTangent(1,ip,el)
 
     case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
-      do c = 1, homogenization_Nconstituents(material_homogenizationAt(el))
+      do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
-        dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el)
+        dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el)
       enddo
       call mech_isostrain_averageStressAndItsTangent(&
-        homogenization_P(1:3,1:3,m), &
+        homogenization_P(1:3,1:3,ce), &
-        homogenization_dPdF(1:3,1:3,1:3,1:3,m),&
+        homogenization_dPdF(1:3,1:3,1:3,1:3,ce),&
         crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
         dPdFs, &
         homogenization_typeInstance(material_homogenizationAt(el)))
 
     case (HOMOGENIZATION_RGC_ID) chosenHomogenization
-      do c = 1, homogenization_Nconstituents(material_homogenizationAt(el))
+      do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
-        dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el)
+        dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el)
       enddo
       call mech_RGC_averageStressAndItsTangent(&
-        homogenization_P(1:3,1:3,m), &
+        homogenization_P(1:3,1:3,ce), &
-        homogenization_dPdF(1:3,1:3,1:3,1:3,m),&
+        homogenization_dPdF(1:3,1:3,1:3,1:3,ce),&
         crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
         dPdFs, &
         homogenization_typeInstance(material_homogenizationAt(el)))
@@ -165,6 +181,45 @@ module subroutine mech_homogenize(ip,el)
 end subroutine mech_homogenize
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief update the internal state of the homogenization scheme and tell whether "done" and
+!> "happy" with result
+!--------------------------------------------------------------------------------------------------
+module function mech_updateState(subdt,subF,ip,el) result(doneAndHappy)
+
+  real(pReal), intent(in) :: &
+    subdt                                                                                           !< current time step
+  real(pReal), intent(in), dimension(3,3) :: &
+    subF
+  integer,     intent(in) :: &
+    ip, &                                                                                           !< integration point
+    el                                                                                              !< element number
+  logical, dimension(2) :: doneAndHappy
+
+  integer :: co
+  real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
+
+
+  if (homogenization_type(material_homogenizationAt(el)) == HOMOGENIZATION_RGC_ID) then
+      do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
+        dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el)
+      enddo
+      doneAndHappy = &
+          mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
+                        crystallite_F(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
+                        crystallite_partitionedF0(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el),&
+                               subF,&
+                               subdt, &
+                               dPdFs, &
+                               ip, &
+                               el)
+  else
+    doneAndHappy = .true.
+  endif
+
+end function mech_updateState
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief Write results to file.
 !--------------------------------------------------------------------------------------------------

src/homogenization_mech_RGC.f90

@@ -8,6 +8,7 @@
 !--------------------------------------------------------------------------------------------------
 submodule(homogenization:homogenization_mech) homogenization_mech_RGC
   use rotations
+  use lattice
 
   type :: tParameters
     integer, dimension(:), allocatable :: &
@@ -18,8 +19,6 @@ submodule(homogenization:homogenization_mech) homogenization_mech_RGC
     real(pReal), dimension(:), allocatable :: &
       D_alpha, &
       a_g
-    integer :: &
-      of_debug = 0
     character(len=pStringLen), allocatable, dimension(:) :: &
       output
   end type tParameters
@@ -151,12 +150,6 @@ module subroutine mech_RGC_init(num_homogMech)
               st0 => state0(homogenization_typeInstance(h)), &
               dst => dependentState(homogenization_typeInstance(h)))
 
-#ifdef DEBUG
-    if  (h==material_homogenizationAt(debugHomog%element)) then
-      prm%of_debug = material_homogenizationMemberAt(debugHomog%ip,debugHomog%element)
-    endif
-#endif
-
 #if defined (__GFORTRAN__)
     prm%output = output_asStrings(homogMech)
 #else
@@ -239,17 +232,6 @@ module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of)
         F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j)                                           ! calculating deformation relaxations due to interface relaxation
     enddo
     F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF                                                    ! resulting relaxed deformation gradient
-
-#ifdef DEBUG
-    if (debugHomog%extensive) then
-      print'(a,i3)',' Deformation gradient of grain: ',iGrain
-      do i = 1,3
-        print'(1x,3(e15.8,1x))',(F(i,j,iGrain), j = 1,3)
-      enddo
-      print*,' '
-      flush(IO_STDOUT)
-    endif
-#endif
   enddo
 
   end associate
@@ -261,7 +243,18 @@ end subroutine mech_RGC_partitionDeformation
 !> @brief update the internal state of the homogenization scheme and tell whether "done" and
 ! "happy" with result
 !--------------------------------------------------------------------------------------------------
-module procedure mech_RGC_updateState
+module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHappy)
+      logical, dimension(2) :: doneAndHappy
+      real(pReal), dimension(:,:,:),     intent(in)    :: &
+        P,&                                                                                         !< partitioned stresses
+        F,&                                                                                         !< partitioned deformation gradients
+        F0                                                                                          !< partitioned initial deformation gradients
+      real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF                                         !< partitioned stiffnesses
+      real(pReal), dimension(3,3),       intent(in) :: avgF                                         !< average F
+      real(pReal),                       intent(in) :: dt                                           !< time increment
+      integer,                           intent(in) :: &
+        ip, &                                                                                       !< integration point number
+        el                                                                                          !< element number
 
   integer, dimension(4) :: intFaceN,intFaceP,faceID
   integer, dimension(3) :: nGDim,iGr3N,iGr3P
@@ -273,13 +266,9 @@ module procedure mech_RGC_updateState
   logical :: error
   real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix,rmatrix
   real(pReal), dimension(:), allocatable   :: resid,relax,p_relax,p_resid,drelax
-#ifdef DEBUG
-  integer, dimension(3) :: stresLoc
-  integer, dimension(2) :: residLoc
-#endif
 
   zeroTimeStep: if(dEq0(dt)) then
-    mech_RGC_updateState = .true.                                                                   ! pretend everything is fine and return
+    doneAndHappy = .true.                                                                   ! pretend everything is fine and return
     return
   endif zeroTimeStep
 
@@ -303,16 +292,6 @@ module procedure mech_RGC_updateState
   relax  = stt%relaxationVector(:,of)
   drelax = stt%relaxationVector(:,of) - st0%relaxationVector(:,of)
 
-#ifdef DEBUG
-  if (debugHomog%extensive) then
-    print*, 'Obtained state: '
-    do i = 1,size(stt%relaxationVector(:,of))
-      print'(1x,2(e15.8,1x))', stt%relaxationVector(i,of)
-    enddo
-    print*,' '
-  endif
-#endif
-
 !--------------------------------------------------------------------------------------------------
 ! computing interface mismatch and stress penalty tensor for all interfaces of all grains
   call stressPenalty(R,NN,avgF,F,ip,el,instance,of)
@@ -353,13 +332,6 @@ module procedure mech_RGC_updateState
       enddo
     enddo
 
-#ifdef DEBUG
-    if (debugHomog%extensive) then
-      print'(a,i3)',' Traction at interface: ',iNum
-      print'(1x,3(e15.8,1x))',(tract(iNum,j), j = 1,3)
-      print*,' '
-    endif
-#endif
   enddo
 
 !--------------------------------------------------------------------------------------------------
@@ -367,29 +339,12 @@ module procedure mech_RGC_updateState
   stresMax = maxval(abs(P))                                                                         ! get the maximum of first Piola-Kirchhoff (material) stress
   residMax = maxval(abs(tract))                                                                     ! get the maximum of the residual
 
-#ifdef DEBUG
+  doneAndHappy = .false.
-  if (debugHomog%extensive .and. prm%of_debug == of) then
-    stresLoc = maxloc(abs(P))
-    residLoc = maxloc(abs(tract))
-    print'(a,i2,1x,i4)',' RGC residual check ... ',ip,el
-    print'(a,e15.8,a,i3,a,i2,i2)', ' Max stress: ',stresMax, &
-               '@ grain ',stresLoc(3),' in component ',stresLoc(1),stresLoc(2)
-    print'(a,e15.8,a,i3,a,i2)',' Max residual: ',residMax, &
-               ' @ iface ',residLoc(1),' in direction ',residLoc(2)
-    flush(IO_STDOUT)
-  endif
-#endif
-
-  mech_RGC_updateState = .false.
 
 !--------------------------------------------------------------------------------------------------
 !  If convergence reached => done and happy
   if (residMax < num%rtol*stresMax .or. residMax < num%atol) then
-    mech_RGC_updateState = .true.
+    doneAndHappy = .true.
-#ifdef DEBUG
-    if (debugHomog%extensive .and. prm%of_debug == of) &
-      print*, '... done and happy'; flush(IO_STDOUT)
-#endif
 
 !--------------------------------------------------------------------------------------------------
 ! compute/update the state for postResult, i.e., all energy densities computed by time-integration
@@ -406,41 +361,14 @@ module procedure mech_RGC_updateState
     dst%relaxationRate_avg(of) = sum(abs(drelax))/dt/real(3*nIntFaceTot,pReal)
     dst%relaxationRate_max(of) = maxval(abs(drelax))/dt
 
-#ifdef DEBUG
-    if (debugHomog%extensive .and. prm%of_debug == of) then
-      print'(a,e15.8)',       ' Constitutive work: ',stt%work(of)
-      print'(a,3(1x,e15.8))', ' Magnitude mismatch: ',dst%mismatch(1,of), &
-                                                      dst%mismatch(2,of), &
-                                                      dst%mismatch(3,of)
-      print'(a,e15.8)',   ' Penalty energy: ',          stt%penaltyEnergy(of)
-      print'(a,e15.8,/)', ' Volume discrepancy: ',      dst%volumeDiscrepancy(of)
-      print'(a,e15.8)',   ' Maximum relaxation rate: ', dst%relaxationRate_max(of)
-      print'(a,e15.8,/)', ' Average relaxation rate: ', dst%relaxationRate_avg(of)
-      flush(IO_STDOUT)
-    endif
-#endif
-
     return
 
 !--------------------------------------------------------------------------------------------------
 ! if residual blows-up => done but unhappy
   elseif (residMax > num%relMax*stresMax .or. residMax > num%absMax) then                           ! try to restart when residual blows up exceeding maximum bound
-    mech_RGC_updateState = [.true.,.false.]                                                         ! with direct cut-back
+    doneAndHappy = [.true.,.false.]                                                         ! with direct cut-back
-
-#ifdef DEBUG
-    if (debugHomog%extensive .and. prm%of_debug == of) &
-      print'(a,/)', ' ... broken'; flush(IO_STDOUT)
-#endif
-
    return
-
+  endif
- else                                                                                               ! proceed with computing the Jacobian and state update
-#ifdef DEBUG
-   if (debugHomog%extensive .and. prm%of_debug == of) &
-     print'(a,/)', ' ... not yet done'; flush(IO_STDOUT)
-#endif
-
- endif
 
 !---------------------------------------------------------------------------------------------------
 ! construct the global Jacobian matrix for updating the global relaxation vector array when
@@ -492,17 +420,6 @@ module procedure mech_RGC_updateState
     enddo
   enddo
 
-#ifdef DEBUG
-  if (debugHomog%extensive) then
-    print*, 'Jacobian matrix of stress'
-    do i = 1,3*nIntFaceTot
-      print'(1x,100(e11.4,1x))',(smatrix(i,j), j = 1,3*nIntFaceTot)
-    enddo
-    print*,' '
-    flush(IO_STDOUT)
-  endif
-#endif
-
 !--------------------------------------------------------------------------------------------------
 ! ... of the stress penalty tangent (mismatch penalty and volume penalty, computed using numerical
 ! perturbation method) "pmatrix"
@@ -552,16 +469,6 @@ module procedure mech_RGC_updateState
     pmatrix(:,ipert) = p_resid/num%pPert
   enddo
 
-#ifdef DEBUG
-  if (debugHomog%extensive) then
-    print*, 'Jacobian matrix of penalty'
-    do i = 1,3*nIntFaceTot
-      print'(1x,100(e11.4,1x))',(pmatrix(i,j), j = 1,3*nIntFaceTot)
-    enddo
-    print*,' '
-    flush(IO_STDOUT)
-  endif
-#endif
 
 !--------------------------------------------------------------------------------------------------
 ! ... of the numerical viscosity traction "rmatrix"
@@ -571,48 +478,16 @@ module procedure mech_RGC_updateState
                    (abs(drelax(i))/(num%refRelaxRate*dt))**(num%viscPower - 1.0_pReal)              ! only in the main diagonal term
   enddo
 
-#ifdef DEBUG
-  if (debugHomog%extensive) then
-    print*, 'Jacobian matrix of penalty'
-    do i = 1,3*nIntFaceTot
-      print'(1x,100(e11.4,1x))',(rmatrix(i,j), j = 1,3*nIntFaceTot)
-    enddo
-    print*,' '
-    flush(IO_STDOUT)
-  endif
-#endif
 
 !--------------------------------------------------------------------------------------------------
 ! The overall Jacobian matrix summarizing contributions of smatrix, pmatrix, rmatrix
   allocate(jmatrix(3*nIntFaceTot,3*nIntFaceTot)); jmatrix = smatrix + pmatrix + rmatrix
 
-#ifdef DEBUG
-  if (debugHomog%extensive) then
-    print*, 'Jacobian matrix (total)'
-    do i = 1,3*nIntFaceTot
-      print'(1x,100(e11.4,1x))',(jmatrix(i,j), j = 1,3*nIntFaceTot)
-    enddo
-    print*,' '
-    flush(IO_STDOUT)
-  endif
-#endif
-
 !--------------------------------------------------------------------------------------------------
 ! computing the update of the state variable (relaxation vectors) using the Jacobian matrix
   allocate(jnverse(3*nIntFaceTot,3*nIntFaceTot),source=0.0_pReal)
   call math_invert(jnverse,error,jmatrix)
 
-#ifdef DEBUG
-  if (debugHomog%extensive) then
-    print*, 'Jacobian inverse'
-    do i = 1,3*nIntFaceTot
-      print'(1x,100(e11.4,1x))',(jnverse(i,j), j = 1,3*nIntFaceTot)
-    enddo
-    print*,' '
-    flush(IO_STDOUT)
-  endif
-#endif
-
 !--------------------------------------------------------------------------------------------------
 ! calculate the state update (global relaxation vectors) for the next Newton-Raphson iteration
   drelax = 0.0_pReal
@@ -621,7 +496,7 @@ module procedure mech_RGC_updateState
   enddo; enddo
   stt%relaxationVector(:,of) = relax + drelax                                                       ! Updateing the state variable for the next iteration
   if (any(abs(drelax) > num%maxdRelax)) then                                                        ! Forcing cutback when the incremental change of relaxation vector becomes too large
-    mech_RGC_updateState = [.true.,.false.]
+    doneAndHappy = [.true.,.false.]
     !$OMP CRITICAL (write2out)
     print'(a,i3,a,i3,a)',' RGC_updateState: ip ',ip,' | el ',el,' enforces cutback'
     print'(a,e15.8)',' due to large relaxation change = ',maxval(abs(drelax))
@@ -629,17 +504,6 @@ module procedure mech_RGC_updateState
     !$OMP END CRITICAL (write2out)
   endif
 
-#ifdef DEBUG
-  if (debugHomog%extensive) then
-    print*, 'Returned state: '
-    do i = 1,size(stt%relaxationVector(:,of))
-      print'(1x,2(e15.8,1x))', stt%relaxationVector(i,of)
-    enddo
-    print*,' '
-    flush(IO_STDOUT)
-  endif
-#endif
-
   end associate
 
   contains
@@ -661,8 +525,10 @@ module procedure mech_RGC_updateState
     real(pReal),   dimension (3)   :: nVect,surfCorr
     real(pReal),   dimension (2)   :: Gmoduli
     integer :: iGrain,iGNghb,iFace,i,j,k,l
-    real(pReal)   :: muGrain,muGNghb,nDefNorm,bgGrain,bgGNghb
+    real(pReal) :: muGrain,muGNghb,nDefNorm
-    real(pReal), parameter  :: nDefToler = 1.0e-10_pReal
+    real(pReal), parameter  :: &
+      nDefToler = 1.0e-10_pReal, &
+      b = 2.5e-10_pReal                                                                             ! Length of Burgers vector
 
     nGDim = param(instance)%N_constituents
     rPen = 0.0_pReal
@@ -676,19 +542,11 @@ module procedure mech_RGC_updateState
 
     associate(prm => param(instance))
 
-#ifdef DEBUG
-    if (debugHomog%extensive .and. prm%of_debug == of) then
-      print'(a,2(1x,i3))', ' Correction factor: ',ip,el
-      print*, surfCorr
-    endif
-#endif
 
    !-----------------------------------------------------------------------------------------------
    ! computing the mismatch and penalty stress tensor of all grains
    grainLoop: do iGrain = 1,product(prm%N_constituents)
-     Gmoduli = equivalentModuli(iGrain,ip,el)
+     muGrain = equivalentMu(iGrain,ip,el)
-     muGrain = Gmoduli(1)                                                                           ! collecting the equivalent shear modulus of grain
-     bgGrain = Gmoduli(2)                                                                           ! and the lengthh of Burgers vector
      iGrain3 = grain1to3(iGrain,prm%N_constituents)                                                 ! get the grain ID in local 3-dimensional index (x,y,z)-position
 
      interfaceLoop: do iFace = 1,6
@@ -700,9 +558,7 @@ module procedure mech_RGC_updateState
        where(iGNghb3 < 1)    iGNghb3 = nGDim
        where(iGNghb3 >nGDim) iGNghb3 = 1
        iGNghb  = grain3to1(iGNghb3,prm%N_constituents)                                              ! get the ID of the neighboring grain
-       Gmoduli = equivalentModuli(iGNghb,ip,el)                                                     ! collect the shear modulus and Burgers vector of the neighbor
+       muGNghb = equivalentMu(iGNghb,ip,el)
-       muGNghb = Gmoduli(1)
-       bgGNghb = Gmoduli(2)
        gDef = 0.5_pReal*(fDef(1:3,1:3,iGNghb) - fDef(1:3,1:3,iGrain))                               ! difference/jump in deformation gradeint across the neighbor
 
        !-------------------------------------------------------------------------------------------
@@ -717,30 +573,19 @@ module procedure mech_RGC_updateState
        enddo; enddo
        nDefNorm = max(nDefToler,sqrt(nDefNorm))                                                     ! approximation to zero mismatch if mismatch is zero (singularity)
        nMis(abs(intFace(1)),iGrain) = nMis(abs(intFace(1)),iGrain) + nDefNorm                       ! total amount of mismatch experienced by the grain (at all six interfaces)
-#ifdef DEBUG
+
-       if (debugHomog%extensive .and. prm%of_debug == of) then
-         print'(a,i2,a,i3)',' Mismatch to face: ',intFace(1),' neighbor grain: ',iGNghb
-         print*, transpose(nDef)
-         print'(a,e11.4)', ' with magnitude: ',nDefNorm
-       endif
-#endif
 
        !-------------------------------------------------------------------------------------------
        ! compute the stress penalty of all interfaces
        do i = 1,3; do j = 1,3; do k = 1,3; do l = 1,3
-         rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*bgGrain + muGNghb*bgGNghb)*prm%xi_alpha &
+         rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*b + muGNghb*b)*prm%xi_alpha &
                                                 *surfCorr(abs(intFace(1)))/prm%D_alpha(abs(intFace(1))) &
                                                 *cosh(prm%c_alpha*nDefNorm) &
                                                 *0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_LeviCivita(k,l,j) &
                                                 *tanh(nDefNorm/num%xSmoo)
        enddo; enddo;enddo; enddo
      enddo interfaceLoop
-#ifdef DEBUG
+
-     if (debugHomog%extensive .and. prm%of_debug == of) then
-       print'(a,i2)', ' Penalty of grain: ',iGrain
-       print*, transpose(rPen(1:3,1:3,iGrain))
-     endif
-#endif
 
    enddo grainLoop
 
@@ -783,13 +628,6 @@ module procedure mech_RGC_updateState
       vPen(:,:,i) = -1.0_pReal/real(nGrain,pReal)*num%volDiscrMod*num%volDiscrPow/num%maxVolDiscr* &
                          sign((abs(vDiscrep)/num%maxVolDiscr)**(num%volDiscrPow - 1.0),vDiscrep)* &
                          gVol(i)*transpose(math_inv33(fDef(:,:,i)))
-
-#ifdef DEBUG
-      if (debugHomog%extensive .and. param(instance)%of_debug == of) then
-        print'(a,i2)',' Volume penalty of grain: ',i
-        print*, transpose(vPen(:,:,i))
-      endif
-#endif
     enddo
 
   end subroutine volumePenalty
@@ -827,44 +665,26 @@ module procedure mech_RGC_updateState
   end function surfaceCorrection
 
 
-  !--------------------------------------------------------------------------------------------------
+  !-------------------------------------------------------------------------------------------------
   !> @brief compute the equivalent shear and bulk moduli from the elasticity tensor
-  !--------------------------------------------------------------------------------------------------
+  !-------------------------------------------------------------------------------------------------
-  function equivalentModuli(grainID,ip,el)
+  real(pReal) function equivalentMu(grainID,ip,el)
-
-    real(pReal), dimension(2)    :: equivalentModuli
 
     integer, intent(in)    :: &
       grainID,&
       ip, &                                                                                         !< integration point number
       el                                                                                            !< element number
-    real(pReal), dimension(6,6) :: elasTens
-    real(pReal) :: &
-      cEquiv_11, &
-      cEquiv_12, &
-      cEquiv_44
-
-    elasTens = constitutive_homogenizedC(grainID,ip,el)
-
-    !----------------------------------------------------------------------------------------------
-    ! compute the equivalent shear modulus after Turterltaub and Suiker, JMPS (2005)
-    cEquiv_11 = (elasTens(1,1) + elasTens(2,2) + elasTens(3,3))/3.0_pReal
-    cEquiv_12 = (elasTens(1,2) + elasTens(2,3) + elasTens(3,1) + &
-                 elasTens(1,3) + elasTens(2,1) + elasTens(3,2))/6.0_pReal
-    cEquiv_44 = (elasTens(4,4) + elasTens(5,5) + elasTens(6,6))/3.0_pReal
-    equivalentModuli(1) = 0.2_pReal*(cEquiv_11 - cEquiv_12) + 0.6_pReal*cEquiv_44
-
-    !----------------------------------------------------------------------------------------------
-    ! obtain the length of Burgers vector (could be model dependend)
-    equivalentModuli(2) = 2.5e-10_pReal
-
-  end function equivalentModuli
 
 
-  !--------------------------------------------------------------------------------------------------
+    equivalentMu = lattice_equivalent_mu(constitutive_homogenizedC(grainID,ip,el),'voigt')
+
+  end function equivalentMu
+
+
+  !-------------------------------------------------------------------------------------------------
   !> @brief calculating the grain deformation gradient (the same with
   ! homogenization_RGC_partitionDeformation, but used only for perturbation scheme)
-  !--------------------------------------------------------------------------------------------------
+  !-------------------------------------------------------------------------------------------------
   subroutine grainDeformation(F, avgF, instance, of)
 
     real(pReal),   dimension(:,:,:), intent(out) :: F                                               !< partitioned F  per grain
@@ -879,7 +699,7 @@ module procedure mech_RGC_updateState
     integer,       dimension(3) :: iGrain3
     integer :: iGrain,iFace,i,j
 
-    !-------------------------------------------------------------------------------------------------
+    !-----------------------------------------------------------------------------------------------
     ! compute the deformation gradient of individual grains due to relaxations
 
     associate(prm => param(instance))
@@ -901,7 +721,7 @@ module procedure mech_RGC_updateState
 
   end subroutine grainDeformation
 
-end procedure mech_RGC_updateState
+end function mech_RGC_updateState
 
 
 !--------------------------------------------------------------------------------------------------

src/kinematics_cleavage_opening.f90

@@ -99,10 +99,10 @@ 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, ipc, ip, el)
+module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el)
 
   integer, intent(in) :: &
-    ipc, &                                                                                          !< grain number
+    co, &                                                                                          !< grain number
     ip, &                                                                                           !< integration point number
     el                                                                                              !< element number
   real(pReal),   intent(in),  dimension(3,3) :: &
@@ -124,7 +124,7 @@ module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S,
 
   Ld = 0.0_pReal
   dLd_dTstar = 0.0_pReal
-  associate(prm => param(kinematics_cleavage_opening_instance(material_phaseAt(ipc,el))))
+  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
 

src/kinematics_slipplane_opening.f90

@@ -117,10 +117,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, ipc, ip, el)
+module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el)
 
   integer, intent(in) :: &
-    ipc, &                                                                                          !< grain number
+    co, &                                                                                          !< grain number
     ip, &                                                                                           !< integration point number
     el                                                                                              !< element number
   real(pReal),   intent(in),  dimension(3,3) :: &
@@ -138,7 +138,7 @@ module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S
     traction_d, traction_t, traction_n, traction_crit, &
     udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt
 
-  phase = material_phaseAt(ipc,el)
+  phase = material_phaseAt(co,el)
   instance = kinematics_slipplane_opening_instance(phase)
   homog = material_homogenizationAt(el)
   damageOffset = material_homogenizationMemberAt(ip,el)

src/kinematics_thermal_expansion.f90

@@ -84,10 +84,10 @@ end function kinematics_thermal_expansion_init
 !--------------------------------------------------------------------------------------------------
 !> @brief constitutive equation for calculating the velocity gradient
 !--------------------------------------------------------------------------------------------------
-module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el)
+module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, co, ip, el)
 
   integer, intent(in) :: &
-    ipc, &                                                                                          !< grain number
+    co, &                                                                                          !< grain number
     ip, &                                                                                           !< integration point number
     el                                                                                              !< element number
   real(pReal),   intent(out), dimension(3,3) :: &
@@ -101,7 +101,7 @@ module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, i
   real(pReal) :: &
     T, TDot
 
-  phase = material_phaseAt(ipc,el)
+  phase = material_phaseAt(co,el)
   homog = material_homogenizationAt(el)
   T = temperature(homog)%p(material_homogenizationMemberAt(ip,el))
   TDot = temperatureRate(homog)%p(material_homogenizationMemberAt(ip,el))

src/lattice.f90

@@ -421,6 +421,8 @@ module lattice
     lattice_BCT_ID, &
     lattice_HEX_ID, &
     lattice_ORT_ID, &
+    lattice_equivalent_nu, &
+    lattice_equivalent_mu, &
     lattice_applyLatticeSymmetry33, &
     lattice_SchmidMatrix_slip, &
     lattice_SchmidMatrix_twin, &
@@ -508,8 +510,8 @@ subroutine lattice_init
 
     lattice_C66(1:6,1:6,p) = applyLatticeSymmetryC66(lattice_C66(1:6,1:6,p),phase%get_asString('lattice'))
 
-    lattice_mu(p) = equivalent_mu(lattice_C66(1:6,1:6,p),'voigt')
+    lattice_nu(p) = lattice_equivalent_nu(lattice_C66(1:6,1:6,p),'voigt')
-    lattice_nu(p) = equivalent_nu(lattice_C66(1:6,1:6,p),'voigt')
+    lattice_mu(p) = lattice_equivalent_mu(lattice_C66(1:6,1:6,p),'voigt')
 
     lattice_C66(1:6,1:6,p) = math_sym3333to66(math_Voigt66to3333(lattice_C66(1:6,1:6,p)))           ! Literature data is in Voigt notation
     do i = 1, 6
@@ -2188,15 +2190,16 @@ end function getlabels
 !> @brief Equivalent Poisson's ratio (ν)
 !> @details https://doi.org/10.1143/JPSJ.20.635
 !--------------------------------------------------------------------------------------------------
-function equivalent_nu(C,assumption) result(nu)
+function lattice_equivalent_nu(C,assumption) result(nu)
 
   real(pReal), dimension(6,6), intent(in) :: C                                                      !< Stiffness tensor (Voigt notation)
   character(len=*),            intent(in) :: assumption                                             !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)
-
   real(pReal)                 :: K, mu, nu
+
   logical                     :: error
   real(pReal), dimension(6,6) :: S
 
+
   if    (IO_lc(assumption) == 'voigt') then
     K = (C(1,1)+C(2,2)+C(3,3) +2.0_pReal*(C(1,2)+C(2,3)+C(1,3))) &
       / 9.0_pReal
@@ -2210,25 +2213,26 @@ function equivalent_nu(C,assumption) result(nu)
     K = 0.0_pReal
   endif
 
-  mu = equivalent_mu(C,assumption)
+  mu = lattice_equivalent_mu(C,assumption)
   nu = (1.5_pReal*K -mu)/(3.0_pReal*K+mu)
 
-end function equivalent_nu
+end function lattice_equivalent_nu
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief Equivalent shear modulus (μ)
 !> @details https://doi.org/10.1143/JPSJ.20.635
 !--------------------------------------------------------------------------------------------------
-function equivalent_mu(C,assumption) result(mu)
+function lattice_equivalent_mu(C,assumption) result(mu)
 
   real(pReal), dimension(6,6), intent(in) :: C                                                      !< Stiffness tensor (Voigt notation)
   character(len=*),            intent(in) :: assumption                                             !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)
-
   real(pReal)                 :: mu
+
   logical                     :: error
   real(pReal), dimension(6,6) :: S
 
+
   if    (IO_lc(assumption) == 'voigt') then
     mu = (1.0_pReal*(C(1,1)+C(2,2)+C(3,3)) -1.0_pReal*(C(1,2)+C(2,3)+C(1,3)) +3.0_pReal*(C(4,4)+C(5,5)+C(6,6))) &
        / 15.0_pReal
@@ -2242,7 +2246,7 @@ function equivalent_mu(C,assumption) result(mu)
     mu = 0.0_pReal
   endif
 
-end function equivalent_mu
+end function lattice_equivalent_mu
 
 
 !--------------------------------------------------------------------------------------------------
@@ -2266,14 +2270,14 @@ subroutine selfTest
   call random_number(C)
   C(1,1) = C(1,1) + 1.0_pReal
   C = applyLatticeSymmetryC66(C,'aP')
-  if(dNeq(C(6,6),equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/voigt'
+  if(dNeq(C(6,6),lattice_equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/voigt'
-  if(dNeq(C(6,6),equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/reuss'
+  if(dNeq(C(6,6),lattice_equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/reuss'
 
   lambda = C(1,2)
-  if(dNeq(lambda*0.5_pReal/(lambda+equivalent_mu(C,'voigt')),equivalent_nu(C,'voigt'),1.0e-12_pReal)) &
+  if(dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'voigt')), &
-    error stop 'equivalent_nu/voigt'
+          lattice_equivalent_nu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_nu/voigt'
-  if(dNeq(lambda*0.5_pReal/(lambda+equivalent_mu(C,'reuss')),equivalent_nu(C,'reuss'),1.0e-12_pReal)) &
+  if(dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'reuss')), &
-    error stop 'equivalent_nu/reuss'
+          lattice_equivalent_nu(C,'reuss'),1.0e-12_pReal)) error stop 'equivalent_nu/reuss'
 
 end subroutine selfTest
 

src/marc/discretization_marc.f90

@@ -12,7 +12,6 @@ module discretization_marc
   use DAMASK_interface
   use IO
   use config
-  use FEsolving
   use element
   use discretization
   use geometry_plastic_nonlocal
@@ -89,9 +88,6 @@ subroutine discretization_marc_init
   if (debug_e < 1 .or. debug_e > nElems)    call IO_error(602,ext_msg='element')
   if (debug_i < 1 .or. debug_i > elem%nIPs) call IO_error(602,ext_msg='IP')
 
-  FEsolving_execElem = [1,nElems]
-  FEsolving_execIP   = [1,elem%nIPs]
-
   allocate(cellNodeDefinition(elem%nNodes-1))
   allocate(connectivity_cell(elem%NcellNodesPerCell,elem%nIPs,nElems))
   call buildCells(connectivity_cell,cellNodeDefinition,&

src/material.f90

@@ -17,29 +17,6 @@ module material
   private
 
   enum, bind(c); enumerator :: &
-    ELASTICITY_UNDEFINED_ID, &
-    ELASTICITY_HOOKE_ID, &
-    PLASTICITY_UNDEFINED_ID, &
-    PLASTICITY_NONE_ID, &
-    PLASTICITY_ISOTROPIC_ID, &
-    PLASTICITY_PHENOPOWERLAW_ID, &
-    PLASTICITY_KINEHARDENING_ID, &
-    PLASTICITY_DISLOTWIN_ID, &
-    PLASTICITY_DISLOTUNGSTEN_ID, &
-    PLASTICITY_NONLOCAL_ID, &
-    SOURCE_UNDEFINED_ID ,&
-    SOURCE_THERMAL_DISSIPATION_ID, &
-    SOURCE_THERMAL_EXTERNALHEAT_ID, &
-    SOURCE_DAMAGE_ISOBRITTLE_ID, &
-    SOURCE_DAMAGE_ISODUCTILE_ID, &
-    SOURCE_DAMAGE_ANISOBRITTLE_ID, &
-    SOURCE_DAMAGE_ANISODUCTILE_ID, &
-    KINEMATICS_UNDEFINED_ID ,&
-    KINEMATICS_CLEAVAGE_OPENING_ID, &
-    KINEMATICS_SLIPPLANE_OPENING_ID, &
-    KINEMATICS_THERMAL_EXPANSION_ID, &
-    STIFFNESS_DEGRADATION_UNDEFINED_ID, &
-    STIFFNESS_DEGRADATION_DAMAGE_ID, &
     THERMAL_ISOTHERMAL_ID, &
     THERMAL_CONDUCTION_ID, &
     DAMAGE_NONE_ID, &
@@ -96,29 +73,6 @@ module material
 
   public :: &
     material_init, &
-    ELASTICITY_UNDEFINED_ID, &
-    ELASTICITY_HOOKE_ID, &
-    PLASTICITY_UNDEFINED_ID, &
-    PLASTICITY_NONE_ID, &
-    PLASTICITY_ISOTROPIC_ID, &
-    PLASTICITY_PHENOPOWERLAW_ID, &
-    PLASTICITY_KINEHARDENING_ID, &
-    PLASTICITY_DISLOTWIN_ID, &
-    PLASTICITY_DISLOTUNGSTEN_ID, &
-    PLASTICITY_NONLOCAL_ID, &
-    SOURCE_UNDEFINED_ID ,&
-    SOURCE_THERMAL_DISSIPATION_ID, &
-    SOURCE_THERMAL_EXTERNALHEAT_ID, &
-    SOURCE_DAMAGE_ISOBRITTLE_ID, &
-    SOURCE_DAMAGE_ISODUCTILE_ID, &
-    SOURCE_DAMAGE_ANISOBRITTLE_ID, &
-    SOURCE_DAMAGE_ANISODUCTILE_ID, &
-    KINEMATICS_UNDEFINED_ID ,&
-    KINEMATICS_CLEAVAGE_OPENING_ID, &
-    KINEMATICS_SLIPPLANE_OPENING_ID, &
-    KINEMATICS_THERMAL_EXPANSION_ID, &
-    STIFFNESS_DEGRADATION_UNDEFINED_ID, &
-    STIFFNESS_DEGRADATION_DAMAGE_ID, &
     THERMAL_ISOTHERMAL_ID, &
     THERMAL_CONDUCTION_ID, &
     DAMAGE_NONE_ID, &

src/math.f90

@@ -279,9 +279,12 @@ real(pReal) pure function math_LeviCivita(i,j,k)
 
   integer, intent(in) :: i,j,k
 
-  if     (all([i,j,k] == [1,2,3]) .or. all([i,j,k] == [2,3,1]) .or. all([i,j,k] == [3,1,2])) then
+  integer :: o
+
+
+  if     (any([(all(cshift([i,j,k],o) == [1,2,3]),o=0,2)])) then
     math_LeviCivita = +1.0_pReal
-  elseif (all([i,j,k] == [3,2,1]) .or. all([i,j,k] == [2,1,3]) .or. all([i,j,k] == [1,3,2])) then
+  elseif (any([(all(cshift([i,j,k],o) == [3,2,1]),o=0,2)])) then
     math_LeviCivita = -1.0_pReal
   else
     math_LeviCivita =  0.0_pReal

src/mesh/DAMASK_mesh.f90

@@ -15,7 +15,6 @@ program DAMASK_mesh
   use IO
   use math
   use CPFEM2
-  use FEsolving
   use config
   use discretization_mesh
   use FEM_Utilities

src/mesh/FEM_utilities.f90

@@ -160,7 +160,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
 
   print'(/,a)', ' ... evaluating constitutive response ......................................'
 
-  call materialpoint_stressAndItsTangent(timeinc)                                                   ! calculate P field
+  call materialpoint_stressAndItsTangent(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems])                ! calculate P field
 
   cutBack = .false.                                                                                 ! reset cutBack status
   

src/mesh/discretization_mesh.f90

@@ -18,7 +18,6 @@ module discretization_mesh
   use config
   use discretization
   use results
-  use FEsolving
   use FEM_quadrature
   use YAML_types
   use prec
@@ -30,7 +29,7 @@ module discretization_mesh
     mesh_Nboundaries, &
     mesh_NcpElemsGlobal
 
-  integer :: &
+  integer, public, protected :: &
     mesh_NcpElems                                                                                   !< total number of CP elements in mesh
 
 !!!! BEGIN DEPRECATED !!!!!
@@ -174,9 +173,6 @@ subroutine discretization_mesh_init(restart)
   if (debug_element < 1 .or. debug_element > mesh_NcpElems) call IO_error(602,ext_msg='element')
   if (debug_ip < 1 .or. debug_ip > mesh_maxNips)            call IO_error(602,ext_msg='IP')
 
-  FEsolving_execElem = [1,mesh_NcpElems]                                                            ! parallel loop bounds set to comprise all DAMASK elements
-  FEsolving_execIP   = [1,mesh_maxNips]
-
   allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal)
 
   call discretization_init(materialAt,&

src/prec.f90

@@ -108,8 +108,10 @@ logical elemental pure function dEq(a,b,tol)
 
   real(pReal), intent(in)           :: a,b
   real(pReal), intent(in), optional :: tol
+
   real(pReal)                       :: eps
 
+
   if (present(tol)) then
     eps = tol
   else
@@ -132,11 +134,8 @@ logical elemental pure function dNeq(a,b,tol)
   real(pReal), intent(in)           :: a,b
   real(pReal), intent(in), optional :: tol
 
-  if (present(tol)) then
+
-    dNeq = .not. dEq(a,b,tol)
+  dNeq = .not. dEq(a,b,tol)
-  else
-    dNeq = .not. dEq(a,b)
-  endif
 
 end function dNeq
 
@@ -151,8 +150,10 @@ logical elemental pure function dEq0(a,tol)
 
   real(pReal), intent(in)           :: a
   real(pReal), intent(in), optional :: tol
+
   real(pReal)                       :: eps
 
+
   if (present(tol)) then
     eps = tol
   else
@@ -175,11 +176,8 @@ logical elemental pure function dNeq0(a,tol)
   real(pReal), intent(in)           :: a
   real(pReal), intent(in), optional :: tol
 
-  if (present(tol)) then
+
-    dNeq0 = .not. dEq0(a,tol)
+  dNeq0 = .not. dEq0(a,tol)
-  else
-    dNeq0 = .not. dEq0(a)
-  endif
 
 end function dNeq0
 
@@ -195,8 +193,10 @@ logical elemental pure function cEq(a,b,tol)
 
   complex(pReal), intent(in)           :: a,b
   real(pReal),    intent(in), optional :: tol
+
   real(pReal)                          :: eps
 
+
   if (present(tol)) then
     eps = tol
   else
@@ -220,11 +220,8 @@ logical elemental pure function cNeq(a,b,tol)
   complex(pReal), intent(in)           :: a,b
   real(pReal),    intent(in), optional :: tol
 
-  if (present(tol)) then
+
-    cNeq = .not. cEq(a,b,tol)
+  cNeq = .not. cEq(a,b,tol)
-  else
-    cNeq = .not. cEq(a,b)
-  endif
 
 end function cNeq
 
@@ -238,6 +235,7 @@ pure function prec_bytesToC_FLOAT(bytes)
   real(C_FLOAT), dimension(size(bytes,kind=pI64)/(storage_size(0._C_FLOAT,pI64)/8_pI64)) :: &
     prec_bytesToC_FLOAT
 
+
   prec_bytesToC_FLOAT = transfer(bytes,prec_bytesToC_FLOAT,size(prec_bytesToC_FLOAT))
 
 end function prec_bytesToC_FLOAT
@@ -252,6 +250,7 @@ pure function prec_bytesToC_DOUBLE(bytes)
   real(C_DOUBLE), dimension(size(bytes,kind=pI64)/(storage_size(0._C_DOUBLE,pI64)/8_pI64)) :: &
     prec_bytesToC_DOUBLE
 
+
   prec_bytesToC_DOUBLE = transfer(bytes,prec_bytesToC_DOUBLE,size(prec_bytesToC_DOUBLE))
 
 end function prec_bytesToC_DOUBLE
@@ -266,6 +265,7 @@ pure function prec_bytesToC_INT32_T(bytes)
   integer(C_INT32_T), dimension(size(bytes,kind=pI64)/(storage_size(0_C_INT32_T,pI64)/8_pI64)) :: &
     prec_bytesToC_INT32_T
 
+
   prec_bytesToC_INT32_T = transfer(bytes,prec_bytesToC_INT32_T,size(prec_bytesToC_INT32_T))
 
 end function prec_bytesToC_INT32_T
@@ -280,6 +280,7 @@ pure function prec_bytesToC_INT64_T(bytes)
   integer(C_INT64_T), dimension(size(bytes,kind=pI64)/(storage_size(0_C_INT64_T,pI64)/8_pI64)) :: &
      prec_bytesToC_INT64_T
 
+
   prec_bytesToC_INT64_T = transfer(bytes,prec_bytesToC_INT64_T,size(prec_bytesToC_INT64_T))
 
 end function prec_bytesToC_INT64_T
@@ -295,6 +296,7 @@ subroutine selfTest
   integer(pInt), dimension(1) :: i
   real(pReal),   dimension(2) :: r
 
+
   realloc_lhs_test = [1,2]
   if (any(realloc_lhs_test/=[1,2]))        error stop 'LHS allocation'
 

src/quaternions.f90

@@ -1,534 +0,0 @@
-!---------------------------------------------------------------------------------------------------
-!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
-!> @author Philip Eisenlohr, Michigan State University
-!> @brief general quaternion math, not limited to unit quaternions
-!> @details w is the real part, (x, y, z) are the imaginary parts.
-!> @details https://en.wikipedia.org/wiki/Quaternion
-!---------------------------------------------------------------------------------------------------
-module quaternions
-  use prec
-
-  implicit none
-  private
-
-  real(pReal), parameter, public :: P = -1.0_pReal                                                  !< parameter for orientation conversion.
-
-  type, public :: quaternion
-    real(pReal), private :: w = 0.0_pReal
-    real(pReal), private :: x = 0.0_pReal
-    real(pReal), private :: y = 0.0_pReal
-    real(pReal), private :: z = 0.0_pReal
-
-
-  contains
-    procedure, private :: add__
-    procedure, private :: pos__
-    generic,   public  :: operator(+) => add__,pos__
-
-    procedure, private :: sub__
-    procedure, private :: neg__
-    generic,   public  :: operator(-) => sub__,neg__
-
-    procedure, private :: mul_quat__
-    procedure, private :: mul_scal__
-    generic,   public  :: operator(*) => mul_quat__, mul_scal__
-
-    procedure, private :: div_quat__
-    procedure, private :: div_scal__
-    generic,   public  :: operator(/) => div_quat__, div_scal__
-
-    procedure, private :: eq__
-    generic,   public  :: operator(==) => eq__
-
-    procedure, private :: neq__
-    generic,   public  :: operator(/=) => neq__
-
-    procedure, private :: pow_quat__
-    procedure, private :: pow_scal__
-    generic,   public  :: operator(**) => pow_quat__, pow_scal__
-
-    procedure, public  :: abs   => abs__
-    procedure, public  :: conjg => conjg__
-    procedure, public  :: real  => real__
-    procedure, public  :: aimag => aimag__
-
-    procedure, public  :: homomorphed
-    procedure, public  :: asArray
-    procedure, public  :: inverse
-
-  end type
-
-  interface assignment (=)
-    module procedure assign_quat__
-    module procedure assign_vec__
-  end interface assignment (=)
-
-  interface quaternion
-    module procedure init__
-  end interface quaternion
-
-  interface abs
-    procedure abs__
-  end interface abs
-
-  interface dot_product
-    procedure dot_product__
-  end interface dot_product
-
-  interface conjg
-    module procedure conjg__
-  end interface conjg
-
-  interface exp
-    module procedure exp__
-  end interface exp
-
-  interface log
-    module procedure log__
-  end interface log
-
-  interface real
-    module procedure real__
-  end interface real
-
-  interface aimag
-    module procedure aimag__
-  end interface aimag
-
-  public :: &
-    quaternions_init, &
-    assignment(=), &
-    conjg, aimag, &
-    log, exp, &
-    abs, dot_product, &
-    inverse, &
-    real
-
-contains
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief Do self test.
-!--------------------------------------------------------------------------------------------------
-subroutine quaternions_init
-
-  print'(/,a)', ' <<<+-  quaternions init  -+>>>'; flush(6)
-
-  call selfTest
-
-end subroutine quaternions_init
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief construct a quaternion from a 4-vector
-!---------------------------------------------------------------------------------------------------
-type(quaternion) pure function init__(array)
-
-  real(pReal), intent(in), dimension(4) :: array
-
-  init__%w = array(1)
-  init__%x = array(2)
-  init__%y = array(3)
-  init__%z = array(4)
-
-end function init__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief assign a quaternion
-!---------------------------------------------------------------------------------------------------
-elemental pure subroutine assign_quat__(self,other)
-
-  type(quaternion), intent(out) :: self
-  type(quaternion), intent(in)  :: other
-
-  self = [other%w,other%x,other%y,other%z]
-
-end subroutine assign_quat__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief assign a 4-vector
-!---------------------------------------------------------------------------------------------------
-pure subroutine assign_vec__(self,other)
-
-  type(quaternion), intent(out)                :: self
-  real(pReal),       intent(in), dimension(4)  :: other
-
-  self%w = other(1)
-  self%x = other(2)
-  self%y = other(3)
-  self%z = other(4)
-
-end subroutine assign_vec__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief add a quaternion
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function add__(self,other)
-
-  class(quaternion), intent(in) :: self,other
-
-  add__ = [ self%w,  self%x,  self%y ,self%z] &
-        + [other%w, other%x, other%y,other%z]
-
-end function add__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief return (unary positive operator)
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function pos__(self)
-
-  class(quaternion), intent(in) :: self
-
-  pos__ = self * (+1.0_pReal)
-
-end function pos__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief subtract a quaternion
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function sub__(self,other)
-
-  class(quaternion), intent(in) :: self,other
-
-  sub__ = [ self%w,  self%x,  self%y ,self%z] &
-        - [other%w, other%x, other%y,other%z]
-
-end function sub__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief negate (unary negative operator)
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function neg__(self)
-
-  class(quaternion), intent(in) :: self
-
-  neg__ = self * (-1.0_pReal)
-
-end function neg__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief multiply with a quaternion
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function mul_quat__(self,other)
-
-  class(quaternion), intent(in) :: self, other
-
-  mul_quat__%w = self%w*other%w - self%x*other%x -      self%y*other%y - self%z*other%z
-  mul_quat__%x = self%w*other%x + self%x*other%w + P * (self%y*other%z - self%z*other%y)
-  mul_quat__%y = self%w*other%y + self%y*other%w + P * (self%z*other%x - self%x*other%z)
-  mul_quat__%z = self%w*other%z + self%z*other%w + P * (self%x*other%y - self%y*other%x)
-
-end function mul_quat__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief multiply with a scalar
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function mul_scal__(self,scal)
-
-  class(quaternion), intent(in) :: self
-  real(pReal),       intent(in) :: scal
-
-  mul_scal__ = [self%w,self%x,self%y,self%z]*scal
-
-end function mul_scal__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief divide by a quaternion
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function div_quat__(self,other)
-
-  class(quaternion), intent(in) :: self, other
-
-  div_quat__ = self * (conjg(other)/(abs(other)**2.0_pReal))
-
-end function div_quat__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief divide by a scalar
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function div_scal__(self,scal)
-
-  class(quaternion), intent(in) :: self
-  real(pReal),       intent(in) :: scal
-
-  div_scal__ = [self%w,self%x,self%y,self%z]/scal
-
-end function div_scal__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief test equality
-!---------------------------------------------------------------------------------------------------
-logical elemental pure function eq__(self,other)
-
-  class(quaternion), intent(in) :: self,other
-
-  eq__ = all(dEq([ self%w, self%x, self%y, self%z], &
-                 [other%w,other%x,other%y,other%z]))
-
-end function eq__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief test inequality
-!---------------------------------------------------------------------------------------------------
-logical elemental pure function neq__(self,other)
-
-  class(quaternion), intent(in) :: self,other
-
-  neq__ = .not. self%eq__(other)
-
-end function neq__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief raise to the power of a quaternion
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function pow_quat__(self,expon)
-
-  class(quaternion), intent(in) :: self
-  type(quaternion),  intent(in) :: expon
-
-  pow_quat__ = exp(log(self)*expon)
-
-end function pow_quat__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief raise to the power of a scalar
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function pow_scal__(self,expon)
-
-  class(quaternion), intent(in) :: self
-  real(pReal),       intent(in) :: expon
-
-  pow_scal__ = exp(log(self)*expon)
-
-end function pow_scal__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief take exponential
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function exp__(a)
-
-  class(quaternion), intent(in) :: a
-  real(pReal)                   :: absImag
-
-  absImag = norm2(aimag(a))
-
-  exp__ = merge(exp(a%w) * [               cos(absImag), &
-                             a%x/absImag * sin(absImag), &
-                             a%y/absImag * sin(absImag), &
-                             a%z/absImag * sin(absImag)], &
-                IEEE_value(1.0_pReal,IEEE_SIGNALING_NAN), &
-                dNeq0(absImag))
-
-end function exp__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief take logarithm
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function log__(a)
-
-  class(quaternion), intent(in) :: a
-  real(pReal)                   :: absImag
-
-  absImag = norm2(aimag(a))
-
-  log__ = merge([log(abs(a)), &
-                 a%x/absImag * acos(a%w/abs(a)), &
-                 a%y/absImag * acos(a%w/abs(a)), &
-                 a%z/absImag * acos(a%w/abs(a))], &
-                IEEE_value(1.0_pReal,IEEE_SIGNALING_NAN), &
-                dNeq0(absImag))
-
-end function log__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief return norm
-!---------------------------------------------------------------------------------------------------
-real(pReal) elemental pure function abs__(self)
-
-  class(quaternion), intent(in) :: self
-
-  abs__ = norm2([self%w,self%x,self%y,self%z])
-
-end function abs__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief calculate dot product
-!---------------------------------------------------------------------------------------------------
-real(pReal) elemental pure function dot_product__(a,b)
-
-  class(quaternion), intent(in) :: a,b
-
-  dot_product__ = a%w*b%w + a%x*b%x + a%y*b%y + a%z*b%z
-
-end function dot_product__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief take conjugate complex
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function conjg__(self)
-
-  class(quaternion), intent(in) :: self
-
-  conjg__ = [self%w,-self%x,-self%y,-self%z]
-
-end function conjg__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief homomorph
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function homomorphed(self)
-
-  class(quaternion), intent(in) :: self
-
-  homomorphed = - self
-
-end function homomorphed
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief return as plain array
-!---------------------------------------------------------------------------------------------------
-pure function asArray(self)
-
-  real(pReal), dimension(4)     :: asArray
-  class(quaternion), intent(in) :: self
-
-  asArray = [self%w,self%x,self%y,self%z]
-
-end function asArray
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief real part (scalar)
-!---------------------------------------------------------------------------------------------------
-pure function real__(self)
-
-  real(pReal)                   :: real__
-  class(quaternion), intent(in) :: self
-
-  real__ = self%w
-
-end function real__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief imaginary part (3-vector)
-!---------------------------------------------------------------------------------------------------
-pure function aimag__(self)
-
-  real(pReal), dimension(3)     :: aimag__
-  class(quaternion), intent(in) :: self
-
-  aimag__ = [self%x,self%y,self%z]
-
-end function aimag__
-
-
-!---------------------------------------------------------------------------------------------------
-!> @brief inverse
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function inverse(self)
-
-  class(quaternion), intent(in) :: self
-
-  inverse = conjg(self)/abs(self)**2.0_pReal
-
-end function inverse
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief check correctness of some quaternions functions
-!--------------------------------------------------------------------------------------------------
-subroutine selfTest
-
-  real(pReal), dimension(4) :: qu
-  type(quaternion)          :: q, q_2
-
-  if(dNeq(abs(P),1.0_pReal))                           error stop 'P not in {-1,+1}'
-
-  call random_number(qu)
-  qu = (qu-0.5_pReal) * 2.0_pReal
-  q  = quaternion(qu)
-
-  q_2= qu
-  if(any(dNeq(q%asArray(),q_2%asArray())))             error stop 'assign_vec__'
-
-  q_2 = q + q
-  if(any(dNeq(q_2%asArray(),2.0_pReal*qu)))            error stop 'add__'
-
-  q_2 = q - q
-  if(any(dNeq0(q_2%asArray())))                        error stop 'sub__'
-
-  q_2 = q * 5.0_pReal
-  if(any(dNeq(q_2%asArray(),5.0_pReal*qu)))            error stop 'mul__'
-
-  q_2 = q / 0.5_pReal
-  if(any(dNeq(q_2%asArray(),2.0_pReal*qu)))            error stop 'div__'
-
-  q_2 = q * 0.3_pReal
-  if(dNeq0(abs(q)) .and. q_2 == q)                     error stop 'eq__'
-
-  q_2 = q
-  if(q_2 /= q)                                         error stop 'neq__'
-
-  if(dNeq(abs(q),norm2(qu)))                           error stop 'abs__'
-  if(dNeq(abs(q)**2.0_pReal, real(q*q%conjg()),1.0e-14_pReal)) &
-                                                       error stop 'abs__/*conjg'
-
-  if(any(dNeq(q%asArray(),qu)))                        error stop 'eq__'
-  if(dNeq(q%real(),       qu(1)))                      error stop 'real()'
-  if(any(dNeq(q%aimag(),  qu(2:4))))                   error stop 'aimag()'
-
-  q_2 = q%homomorphed()
-  if(q                 /= q_2*    (-1.0_pReal))        error stop 'homomorphed'
-  if(dNeq(q_2%real(),     qu(1)*  (-1.0_pReal)))       error stop 'homomorphed/real'
-  if(any(dNeq(q_2%aimag(),qu(2:4)*(-1.0_pReal))))      error stop 'homomorphed/aimag'
-
-  q_2 = conjg(q)
-  if(dNeq(abs(q),abs(q_2)))                            error stop 'conjg/abs'
-  if(q /= conjg(q_2))                                  error stop 'conjg/involution'
-  if(dNeq(q_2%real(),     q%real()))                   error stop 'conjg/real'
-  if(any(dNeq(q_2%aimag(),q%aimag()*(-1.0_pReal))))    error stop 'conjg/aimag'
-
-  if(abs(q) > 0.0_pReal) then
-    q_2 = q * q%inverse()
-    if(     dNeq(real(q_2), 1.0_pReal,1.0e-15_pReal))  error stop 'inverse/real'
-    if(any(dNeq0(aimag(q_2),          1.0e-15_pReal))) error stop 'inverse/aimag'
-
-    q_2 = q/abs(q)
-    q_2 = conjg(q_2) - inverse(q_2)
-    if(any(dNeq0(q_2%asArray(),1.0e-15_pReal)))        error stop 'inverse/conjg'
-  endif
-  if(dNeq(dot_product(qu,qu),dot_product(q,q)))        error stop 'dot_product'
-
-#if !(defined(__GFORTRAN__) &&  __GNUC__ < 9)
-  if (norm2(aimag(q)) > 0.0_pReal) then
-    if (dNeq0(abs(q-exp(log(q))),1.0e-13_pReal))       error stop 'exp/log'
-    if (dNeq0(abs(q-log(exp(q))),1.0e-13_pReal))       error stop 'log/exp'
-  endif
-#endif
-
-end subroutine selfTest
-
-
-end module quaternions

src/results.f90

@@ -111,8 +111,6 @@ subroutine results_addIncrement(inc,time)
   call results_closeGroup(results_addGroup(trim('inc'//trim(adjustl(incChar)))))
   call results_setLink(trim('inc'//trim(adjustl(incChar))),'current')
   call results_addAttribute('time/s',time,trim('inc'//trim(adjustl(incChar))))
-  call results_closeGroup(results_addGroup('current/phase'))
-  call results_closeGroup(results_addGroup('current/homogenization'))
 
 end subroutine results_addIncrement
 

src/rotations.f90

@@ -47,16 +47,16 @@
 !---------------------------------------------------------------------------------------------------
 
 module rotations
-  use prec
   use IO
   use math
-  use quaternions
 
   implicit none
   private
 
+  real(pReal), parameter :: P = -1.0_pReal                                                          !< parameter for orientation conversion.
+
   type, public :: rotation
-    type(quaternion) :: q
+    real(pReal), dimension(4) :: q
     contains
       procedure, public :: asQuaternion
       procedure, public :: asEulers
@@ -103,7 +103,6 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine rotations_init
 
-  call quaternions_init
   print'(/,a)', ' <<<+-  rotations init  -+>>>'; flush(IO_STDOUT)
 
   print*, 'Rowenhorst et al., Modelling and Simulation in Materials Science and Engineering 23:083501, 2015'
@@ -122,7 +121,7 @@ pure function asQuaternion(self)
   class(rotation), intent(in) :: self
   real(pReal), dimension(4)   :: asQuaternion
 
-  asQuaternion = self%q%asArray()
+  asQuaternion = self%q
 
 end function asQuaternion
 !---------------------------------------------------------------------------------------------------
@@ -131,7 +130,7 @@ pure function asEulers(self)
   class(rotation), intent(in) :: self
   real(pReal), dimension(3)   :: asEulers
 
-  asEulers = qu2eu(self%q%asArray())
+  asEulers = qu2eu(self%q)
 
 end function asEulers
 !---------------------------------------------------------------------------------------------------
@@ -140,7 +139,7 @@ pure function asAxisAngle(self)
   class(rotation), intent(in) :: self
   real(pReal), dimension(4)   :: asAxisAngle
 
-  asAxisAngle = qu2ax(self%q%asArray())
+  asAxisAngle = qu2ax(self%q)
 
 end function asAxisAngle
 !---------------------------------------------------------------------------------------------------
@@ -149,7 +148,7 @@ pure function asMatrix(self)
   class(rotation), intent(in) :: self
   real(pReal), dimension(3,3) :: asMatrix
 
-  asMatrix = qu2om(self%q%asArray())
+  asMatrix = qu2om(self%q)
 
 end function asMatrix
 !---------------------------------------------------------------------------------------------------
@@ -158,7 +157,7 @@ pure function asRodrigues(self)
   class(rotation), intent(in) :: self
   real(pReal), dimension(4)   :: asRodrigues
 
-  asRodrigues = qu2ro(self%q%asArray())
+  asRodrigues = qu2ro(self%q)
 
 end function asRodrigues
 !---------------------------------------------------------------------------------------------------
@@ -167,7 +166,7 @@ pure function asHomochoric(self)
   class(rotation), intent(in) :: self
   real(pReal), dimension(3)   :: asHomochoric
 
-  asHomochoric = qu2ho(self%q%asArray())
+  asHomochoric = qu2ho(self%q)
 
 end function asHomochoric
 
@@ -259,7 +258,7 @@ pure elemental function rotRot__(self,R) result(rRot)
   type(rotation)              :: rRot
   class(rotation), intent(in) :: self,R
 
-  rRot = rotation(self%q*R%q)
+  rRot = rotation(multiply_quaternion(self%q,R%q))
   call rRot%standardize()
 
 end function rotRot__
@@ -272,14 +271,14 @@ pure elemental subroutine standardize(self)
 
   class(rotation), intent(inout) :: self
 
-  if (real(self%q) < 0.0_pReal) self%q = self%q%homomorphed()
+  if (self%q(1) < 0.0_pReal) self%q = - self%q
 
 end subroutine standardize
 
 
 !---------------------------------------------------------------------------------------------------
 !> @author Marc De Graef, Carnegie Mellon University
-!> @brief rotate a vector passively (default) or actively
+!> @brief Rotate a vector passively (default) or actively.
 !---------------------------------------------------------------------------------------------------
 pure function rotVector(self,v,active) result(vRot)
 
@@ -288,9 +287,8 @@ pure function rotVector(self,v,active) result(vRot)
   real(pReal),     intent(in), dimension(3) :: v
   logical,         intent(in), optional     :: active
 
-  real(pReal),    dimension(3) :: v_normed
+  real(pReal), dimension(4) :: v_normed, q
-  type(quaternion)             :: q
+  logical                   :: passive
-  logical                      :: passive
 
   if (present(active)) then
     passive = .not. active
@@ -301,13 +299,13 @@ pure function rotVector(self,v,active) result(vRot)
   if (dEq0(norm2(v))) then
     vRot = v
   else
-    v_normed = v/norm2(v)
+    v_normed = [0.0_pReal,v]/norm2(v)
     if (passive) then
-      q = self%q * (quaternion([0.0_pReal, v_normed(1), v_normed(2), v_normed(3)]) * conjg(self%q) )
+      q = multiply_quaternion(self%q, multiply_quaternion(v_normed, conjugate_quaternion(self%q)))
     else
-      q = conjg(self%q) * (quaternion([0.0_pReal, v_normed(1), v_normed(2), v_normed(3)]) * self%q )
+      q = multiply_quaternion(conjugate_quaternion(self%q), multiply_quaternion(v_normed, self%q))
     endif
-    vRot = q%aimag()*norm2(v)
+    vRot = q(2:4)*norm2(v)
   endif
 
 end function rotVector
@@ -315,8 +313,8 @@ end function rotVector
 
 !---------------------------------------------------------------------------------------------------
 !> @author Marc De Graef, Carnegie Mellon University
-!> @brief rotate a rank-2 tensor passively (default) or actively
+!> @brief Rotate a rank-2 tensor passively (default) or actively.
-!> @details: rotation is based on rotation matrix
+!> @details: Rotation is based on rotation matrix
 !---------------------------------------------------------------------------------------------------
 pure function rotTensor2(self,T,active) result(tRot)
 
@@ -403,7 +401,7 @@ pure elemental function misorientation(self,other)
   type(rotation)              :: misorientation
   class(rotation), intent(in) :: self, other
 
-  misorientation%q = other%q * conjg(self%q)
+  misorientation%q = multiply_quaternion(other%q, conjugate_quaternion(self%q))
 
 end function misorientation
 
@@ -1338,7 +1336,7 @@ end function cu2ho
 !--------------------------------------------------------------------------
 !> @author Marc De Graef, Carnegie Mellon University
 !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief determine to which pyramid a point in a cubic grid belongs
+!> @brief Determine to which pyramid a point in a cubic grid belongs.
 !--------------------------------------------------------------------------
 pure function GetPyramidOrder(xyz)
 
@@ -1362,7 +1360,39 @@ end function GetPyramidOrder
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief check correctness of some rotations functions
+!> @brief Multiply two quaternions.
+!--------------------------------------------------------------------------------------------------
+pure function multiply_quaternion(qu1,qu2)
+
+  real(pReal), dimension(4), intent(in) :: qu1, qu2
+  real(pReal), dimension(4) :: multiply_quaternion
+
+
+  multiply_quaternion(1) = qu1(1)*qu2(1) - qu1(2)*qu2(2) -      qu1(3)*qu2(3) - qu1(4)*qu2(4)
+  multiply_quaternion(2) = qu1(1)*qu2(2) + qu1(2)*qu2(1) + P * (qu1(3)*qu2(4) - qu1(4)*qu2(3))
+  multiply_quaternion(3) = qu1(1)*qu2(3) + qu1(3)*qu2(1) + P * (qu1(4)*qu2(2) - qu1(2)*qu2(4))
+  multiply_quaternion(4) = qu1(1)*qu2(4) + qu1(4)*qu2(1) + P * (qu1(2)*qu2(3) - qu1(3)*qu2(2))
+
+end function multiply_quaternion
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Calculate conjugate complex of a quaternion.
+!--------------------------------------------------------------------------------------------------
+pure function conjugate_quaternion(qu)
+
+  real(pReal), dimension(4), intent(in) :: qu
+  real(pReal), dimension(4) :: conjugate_quaternion
+
+
+  conjugate_quaternion = [qu(1), -qu(2), -qu(3), -qu(4)]
+
+
+end function conjugate_quaternion
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Check correctness of some rotations functions.
 !--------------------------------------------------------------------------------------------------
 subroutine selfTest
 
@@ -1374,7 +1404,8 @@ subroutine selfTest
   real    :: A,B
   integer :: i
 
-  do i=1,10
+
+  do i = 1, 10
 
 #if defined(__GFORTRAN__) && __GNUC__<9
     if(i<7) cycle

src/source_damage_anisoBrittle.f90

@@ -120,10 +120,10 @@ end function source_damage_anisoBrittle_init
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates derived quantities from state
 !--------------------------------------------------------------------------------------------------
-module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el)
+module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el)
 
   integer, intent(in) :: &
-    ipc, &                                                                                          !< component-ID of integration point
+    co, &                                                                                          !< component-ID of integration point
     ip, &                                                                                           !< integration point
     el                                                                                              !< element
   real(pReal),  intent(in), dimension(3,3) :: &
@@ -139,8 +139,8 @@ module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el)
   real(pReal) :: &
     traction_d, traction_t, traction_n, traction_crit
 
-  phase = material_phaseAt(ipc,el)
+  phase = material_phaseAt(co,el)
-  constituent = material_phasememberAt(ipc,ip,el)
+  constituent = material_phasememberAt(co,ip,el)
   sourceOffset = source_damage_anisoBrittle_offset(phase)
   homog = material_homogenizationAt(el)
   damageOffset = material_homogenizationMemberAt(ip,el)

src/source_damage_anisoDuctile.f90

@@ -107,10 +107,10 @@ end function source_damage_anisoDuctile_init
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates derived quantities from state
 !--------------------------------------------------------------------------------------------------
-module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el)
+module subroutine source_damage_anisoDuctile_dotState(co, ip, el)
 
   integer, intent(in) :: &
-    ipc, &                                                                                          !< component-ID of integration point
+    co, &                                                                                          !< component-ID of integration point
     ip, &                                                                                           !< integration point
     el                                                                                              !< element
 
@@ -121,8 +121,8 @@ module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el)
     damageOffset, &
     homog
 
-  phase = material_phaseAt(ipc,el)
+  phase = material_phaseAt(co,el)
-  constituent = material_phasememberAt(ipc,ip,el)
+  constituent = material_phasememberAt(co,ip,el)
   sourceOffset = source_damage_anisoDuctile_offset(phase)
   homog = material_homogenizationAt(el)
   damageOffset = material_homogenizationMemberAt(ip,el)

src/source_damage_isoBrittle.f90

@@ -94,10 +94,10 @@ end function source_damage_isoBrittle_init
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates derived quantities from state
 !--------------------------------------------------------------------------------------------------
-module subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
+module subroutine source_damage_isoBrittle_deltaState(C, Fe, co, ip, el)
 
   integer, intent(in) :: &
-    ipc, &                                                                                          !< component-ID of integration point
+    co, &                                                                                          !< component-ID of integration point
     ip, &                                                                                           !< integration point
     el                                                                                              !< element
   real(pReal),  intent(in), dimension(3,3) :: &
@@ -114,8 +114,8 @@ module subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
   real(pReal) :: &
     strainenergy
 
-  phase = material_phaseAt(ipc,el)                                                                  !< phase ID at ipc,ip,el
+  phase = material_phaseAt(co,el)                                                                  !< phase ID at co,ip,el
-  constituent = material_phasememberAt(ipc,ip,el)                                                   !< state array offset for phase ID at ipc,ip,el
+  constituent = material_phasememberAt(co,ip,el)                                                   !< state array offset for phase ID at co,ip,el
   sourceOffset = source_damage_isoBrittle_offset(phase)
 
   strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3)

src/source_damage_isoDuctile.f90

@@ -98,10 +98,10 @@ end function source_damage_isoDuctile_init
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates derived quantities from state
 !--------------------------------------------------------------------------------------------------
-module subroutine source_damage_isoDuctile_dotState(ipc, ip, el)
+module subroutine source_damage_isoDuctile_dotState(co, ip, el)
 
   integer, intent(in) :: &
-    ipc, &                                                                                          !< component-ID of integration point
+    co, &                                                                                          !< component-ID of integration point
     ip, &                                                                                           !< integration point
     el                                                                                              !< element
 
@@ -112,8 +112,8 @@ module subroutine source_damage_isoDuctile_dotState(ipc, ip, el)
     damageOffset, &
     homog
 
-  phase = material_phaseAt(ipc,el)
+  phase = material_phaseAt(co,el)
-  constituent = material_phasememberAt(ipc,ip,el)
+  constituent = material_phasememberAt(co,ip,el)
   sourceOffset = source_damage_isoDuctile_offset(phase)
   homog = material_homogenizationAt(el)
   damageOffset = material_homogenizationMemberAt(ip,el)