From 6b11d438424ccc74241719434de678d1d156758d Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Wed, 1 Apr 2020 08:02:08 +0200
Subject: [PATCH] handle error checking centrally

---
 src/constitutive.f90 | 22 +++++++++++++++-------
 src/crystallite.f90  | 13 ++-----------
 2 files changed, 17 insertions(+), 18 deletions(-)

diff --git a/src/constitutive.f90 b/src/constitutive.f90
index 4402820e8..582222dae 100644
--- a/src/constitutive.f90
+++ b/src/constitutive.f90
@@ -327,7 +327,7 @@ module constitutive
     constitutive_initialFi, &
     constitutive_SandItsTangents, &
     constitutive_collectDotState, &
-    constitutive_collectDeltaState, &
+    constitutive_deltaState, &
     constitutive_results
 
 contains
@@ -794,7 +794,7 @@ end function constitutive_collectDotState
 !> @brief for constitutive models having an instantaneous change of state
 !> will return false if delta state is not needed/supported by the constitutive model
 !--------------------------------------------------------------------------------------------------
-subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
+function constitutive_deltaState(S, Fe, Fi, ipc, ip, el) result(broken)
 
   integer, intent(in) :: &
     ipc, &                                                                                          !< component-ID of integration point
@@ -808,13 +808,18 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
     Mp
   integer :: &
     i, &
-    instance, of
+    instance, of, &
+    phase
+  logical :: &
+    broken
+
 
   Mp  = matmul(matmul(transpose(Fi),Fi),S)
   of = material_phasememberAt(ipc,ip,el)
+  phase = material_phaseAt(ipc,el)
   instance = phase_plasticityInstance(material_phaseAt(ipc,el))
 
-  plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el)))
+  plasticityType: select case (phase_plasticity(phase))
 
     case (PLASTICITY_KINEHARDENING_ID) plasticityType
       call plastic_kinehardening_deltaState(Mp,instance,of)
@@ -823,10 +828,11 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
       call plastic_nonlocal_deltaState(Mp,instance,of,ip,el)
 
   end select plasticityType
+  broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of)))
 
-  sourceLoop: do i = 1, phase_Nsources(material_phaseAt(ipc,el))
+  sourceLoop: do i = 1, phase_Nsources(phase)
 
-     sourceType: select case (phase_source(i,material_phaseAt(ipc,el)))
+     sourceType: select case (phase_source(i,phase))
 
       case (SOURCE_damage_isoBrittle_ID) sourceType
         call source_damage_isoBrittle_deltaState  (constitutive_homogenizedC(ipc,ip,el), Fe, &
@@ -834,9 +840,11 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
 
     end select sourceType
 
+    broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%deltaState(:,of)))
+
   enddo SourceLoop
 
-end subroutine constitutive_collectDeltaState
+end function constitutive_deltaState
 
 
 !--------------------------------------------------------------------------------------------------
diff --git a/src/crystallite.f90 b/src/crystallite.f90
index 03de4dea3..81119f2e2 100644
--- a/src/crystallite.f90
+++ b/src/crystallite.f90
@@ -1654,18 +1654,15 @@ logical function stateJump(ipc,ip,el)
   c = material_phaseMemberAt(ipc,ip,el)
   p = material_phaseAt(ipc,el)
 
-  call constitutive_collectDeltaState(crystallite_S(1:3,1:3,ipc,ip,el), &
+  stateJump = .not. constitutive_deltaState(crystallite_S(1:3,1:3,ipc,ip,el), &
                                       crystallite_Fe(1:3,1:3,ipc,ip,el), &
                                       crystallite_Fi(1:3,1:3,ipc,ip,el), &
                                       ipc,ip,el)
+  if(.not. stateJump) return
 
   myOffset = plasticState(p)%offsetDeltaState
   mySize   = plasticState(p)%sizeDeltaState
 
-  if( any(IEEE_is_NaN(plasticState(p)%deltaState(1:mySize,c)))) then
-    stateJump = .false.
-    return
-  endif
 
   plasticState(p)%state(myOffset + 1:myOffset + mySize,c) = &
   plasticState(p)%state(myOffset + 1:myOffset + mySize,c) + plasticState(p)%deltaState(1:mySize,c)
@@ -1673,16 +1670,10 @@ logical function stateJump(ipc,ip,el)
   do mySource = 1, phase_Nsources(p)
     myOffset = sourceState(p)%p(mySource)%offsetDeltaState
     mySize   = sourceState(p)%p(mySource)%sizeDeltaState
-    if (any(IEEE_is_NaN(sourceState(p)%p(mySource)%deltaState(1:mySize,c)))) then
-      stateJump = .false.
-      return
-    endif
     sourceState(p)%p(mySource)%state(myOffset + 1: myOffset + mySize,c) = &
     sourceState(p)%p(mySource)%state(myOffset + 1: myOffset + mySize,c) + sourceState(p)%p(mySource)%deltaState(1:mySize,c)
   enddo
 
-  stateJump = .true.
-
 end function stateJump