From 7196c697e2b8c1ccb0d308178fed498c3719f341 Mon Sep 17 00:00:00 2001
From: achalhp <achal_hp@hotmail.com>
Date: Sun, 4 Feb 2024 19:23:46 +0530
Subject: [PATCH] Edit1 for Fp multiply and adjust plasticstate indexing

---
 src/phase_mechanical.f90                       | 13 ++++++++++++-
 src/phase_mechanical_plastic_phenopowerlaw.f90 | 10 ++++++----
 2 files changed, 18 insertions(+), 5 deletions(-)

diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90
index c143adbfe..bfbaa1133 100644
--- a/src/phase_mechanical.f90
+++ b/src/phase_mechanical.f90
@@ -471,10 +471,21 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
                                         Fe, Fi_new, ph, en)
       ! achal call Kinematic DeltaFp here
       !** Starting to implement changes for accommodating large shear and reorientation caused by twinning**
-      if(.not. FpJumped .and. NiterationStressLp>1) then
+      if(.not. FpJumped .and. NiterationStressLp>1) then                !Achal: Reason for this if statement?
         call plastic_KinematicJump(ph, en, FpJumped,deltaFp)
+        if(en==15) write(6,*)'deltaFp',deltaFp                          !Achal Delete
+        
+        !converged = .true. means no more iteration
 
+        if(FpJumped) then
+          !crystallite_converged(ipc,ip,el) = .true.  !> See "phase_mechanical_constitutive" and "homogenization_mechanical_response"
+          !crystallite_todo(ipc,ip,el) = .false. !> Can't find this
+          ! _converged = .not. broken
+          subFp0 = matmul(deltaFp,subFp0) 
+! subFp0 is input need to change "phase_mechanical_Fp(ph)%data(1:3,1:3,en) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal)"
 
+          !plasticState(ph)%state()
+        endif
 
       endif
                                         
diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90
index 57f6ea55b..75390e6be 100644
--- a/src/phase_mechanical_plastic_phenopowerlaw.f90
+++ b/src/phase_mechanical_plastic_phenopowerlaw.f90
@@ -439,9 +439,9 @@ associate(prm => param(ph), stt => state(ph), &
   call kinetics_tw(Mp,ph,en,dot_gamma_tw,fdot_twin)
   
   
-  if(en==1) call plastic_kinematic_deltaFp(ph,en,twinJump,deltaFp)                                        ! delete this  
+  !if(en==1) call plastic_kinematic_deltaFp(ph,en,twinJump,deltaFp)                                        ! delete this  
   
-  if(en==1) write(6,*)'f_twin',dotState(indexDotState(ph)%f_twin(1):indexDotState(ph)%f_twin(2))                  !Achal delete
+  !if(en==1) write(6,*)'f_twin',dotState(indexDotState(ph)%f_twin(1):indexDotState(ph)%f_twin(2))                  !Achal delete
 
   sumF = sum(stt%gamma_tw(:,en)/prm%gamma_char)
   xi_sl_sat_offset = prm%f_sat_sl_tw*sqrt(sumF)
@@ -531,7 +531,8 @@ associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dlt => deltas
     Success_Nucleation: if (random*0.00000000000000000000001_pReal <= stt%f_twin(twin_var,en)) then          ! Instead of sum take max
       twinJump = .true.
       deltaFp  = prm%CorrespondanceMatrix(:,:,twin_var)
-      dlt%f_twin(:,en) = 0.0_pReal - stt%f_twin(:,en)    
+      dlt%f_twin(:,en)     = 0.0_pReal - stt%f_twin(:,en)
+      dlt%fmc_twin(:,en)   = 0.0_pReal - stt%fmc_twin(:,en) 
     end if Success_Nucleation
 
   endif Ability_Nucleation
@@ -551,7 +552,8 @@ integer, intent(in)::&
   ph, &
   en
 
-deltastate(ph)%f_twin=0.0_pReal
+deltastate(ph)%f_twin   = 0.0_pReal
+deltastate(ph)%fmc_twin = 0.0_pReal
   
 
 end subroutine plastic_phenopowerlaw_deltaState