added HDF5 outputs

src/phase_mechanical.f90

@@ -473,15 +473,15 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
       !** Starting to implement changes for accommodating large shear and reorientation caused by twinning**
       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
+        !if(en==15) write(6,*)'deltaFp',deltaFp                          !Achal Delete
-        
+        if(en==15) write(6,*)'FpJumped',FpJumped
         !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) 
+          Fp_new = 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()
@@ -1031,7 +1031,7 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_)
     formerSubStep
   integer :: &
     ph, en, sizeDotState
-  logical :: todo
+  logical :: todo, FpJumped                                                                  !Achal
   real(pReal) :: subFrac,subStep
   real(pReal), dimension(3,3) :: &
     subFp0, &
@@ -1039,7 +1039,8 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_)
     subLp0, &
     subLi0, &
     subF0, &
-    subF
+    subF, &
+    deltaFp                                                                                 !Achal
   real(pReal), dimension(plasticState(material_phaseID(co,ce))%sizeState) :: subState0
 
 
@@ -1060,6 +1061,18 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_)
   todo = .true.
   cutbackLooping: do while (todo)
 
+  !matmul(deltaFp,subFp0)
+
+  ! achal calling Kinematic DeltaFp here
+  !** Starting to implement changes for accommodating large shear and reorientation caused by twinning**
+  !if(.not. FpJumped .and. NiterationStressLp>1) then                !Achal: Reason for this if statement?
+    call plastic_KinematicJump(ph, en, FpJumped,deltaFp)
+    if(FpJumped) then
+      subFp0 = matmul(deltaFp,phase_mechanical_Fp0(ph)%data(1:3,1:3,en))  
+    endif
+
+  !endif
+
     if (converged_) then
       formerSubStep = subStep
       subFrac = subFrac + subStep

src/phase_mechanical_plastic_phenopowerlaw.f90

@@ -82,8 +82,11 @@ type :: tPhenopowerlawState
     xi_tw, &
     gamma_sl, &
     gamma_tw, &
-    f_twin, &
+    f_twin, &                                                                                   !< Twin volume fraction
-    fmc_twin                                                                                       !< Achal, To control sampling frequency
+    fmc_twin                                                                                    !< Achal, To control sampling frequency
+  real(pReal), pointer, dimension(:) :: &
+    variant_twin, &
+    frozen
 end type tPhenopowerlawState
 
 !--------------------------------------------------------------------------------------------------
@@ -261,11 +264,13 @@ do ph = 1, phases%length
   sizeDotState = size(['xi_sl   ','gamma_sl']) * prm%sum_N_sl &
                + size(['xi_tw   ','gamma_tw','f_twin  ']) * prm%sum_N_tw      !Achal
   
-  sizeDeltaState = size(['f_twin  ','fmc_twin']) * prm%sum_N_tw                         !Achal
+  sizeDeltaState = size(['f_twin  ','fmc_twin']) * prm%sum_N_tw  &                       !Achal
+                   + size(['variant_twin','frozen      ']) * prm%sum_N_tw
                
   sizeState = size(['xi_sl   ','gamma_sl']) * prm%sum_N_sl &
               + size(['xi_tw   ','gamma_tw']) * prm%sum_N_tw &
-              + size(['f_twin  ','fmc_twin']) * prm%sum_N_tw                !Achal
+              + size(['f_twin  ','fmc_twin']) * prm%sum_N_tw &
+              + size(['variant_twin','frozen      ']) * prm%sum_N_tw              !Achal
 
 
 
@@ -326,12 +331,27 @@ do ph = 1, phases%length
   dlt%f_twin => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:)         ! Achal
   plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
 
+  startIndex =  endIndex + 1
+  endIndex   =  endIndex + prm%sum_N_tw
+  stt%frozen => plasticState(ph)%state(startIndex,:)
+  stt%frozen = 0.0_pReal-1.0_pReal 
+  dlt%frozen => plasticState(ph)%deltaState(startIndex-o,:)
+  plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
+
   startIndex =  endIndex + 1
   endIndex   =  endIndex + prm%sum_N_tw
   stt%fmc_twin => plasticState(ph)%state(startIndex:endIndex,:)
   dlt%fmc_twin => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:)      
   plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
 
+  startIndex       =  endIndex + 1
+  endIndex         =  endIndex + prm%sum_N_tw
+  stt%variant_twin => plasticState(ph)%state(startIndex,:)
+  dlt%variant_twin => plasticState(ph)%deltaState(startIndex-o,:)      
+  plasticState(ph)%atol(startIndex:endIndex) = 0.0_pReal
+
+
+  
   write(6,*)"offset", o                          ! Achal Delete
   
 
@@ -532,6 +552,8 @@ associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dlt => deltas
       deltaFp  = prm%CorrespondanceMatrix(:,:,twin_var)
       dlt%f_twin(:,en)     = 0.0_pReal - stt%f_twin(:,en)
       dlt%fmc_twin(:,en)   = 0.0_pReal - stt%fmc_twin(:,en)
+      dlt%frozen(en)       = 1.0_pReal - stt%frozen(en)
+      dlt%variant_twin(en) = twin_var - stt%variant_twin(en) 
     end if Success_Nucleation
 
   endif Ability_Nucleation
@@ -553,6 +575,8 @@ integer, intent(in)::&
 
 deltastate(ph)%f_twin   = 0.0_pReal
 deltastate(ph)%fmc_twin = 0.0_pReal
+deltastate(ph)%variant_twin = 0
+deltastate(ph)%frozen = 0.0_pReal
   
 
 end subroutine plastic_phenopowerlaw_deltaState
@@ -589,6 +613,17 @@ associate(prm => param(ph), stt => state(ph))
         call results_writeDataset(stt%gamma_tw,group,trim(prm%output(ou)), &
                                   'twinning shear','1',prm%systems_tw)
 
+      case('f_twin')
+        call results_writeDataset(stt%f_twin,group,trim(prm%output(ou)), &
+                                  'volume fraction','1',prm%systems_tw)
+
+      case('variant_twin')
+        call results_writeDataset(stt%variant_twin,group,trim(prm%output(ou)), &
+                                  'twin variant','1',prm%systems_tw)
+
+      case('fbinary_twin')
+        call results_writeDataset(stt%frozen,group,trim(prm%output(ou)), &
+                                  'binary twin flag','1',prm%systems_tw)                                  
     end select
 
   end do