From df5cdb6b017ddfce676f3880aee1e13df575a9d6 Mon Sep 17 00:00:00 2001
From: Su	Leen	Wong <s.wong@mpie.de>
Date: Wed, 8 Jul 2015 11:58:52 +0000
Subject: [PATCH] Elasticity matrices for transformed phase

---
 code/IO.f90                |  2 ++
 code/lattice.f90           | 35 +++++++++++++++++++++++++++++++++--
 code/plastic_dislotwin.f90 | 21 ++++++++++++---------
 3 files changed, 47 insertions(+), 11 deletions(-)

diff --git a/code/IO.f90 b/code/IO.f90
index a68c197d2..0ab8bfd90 100644
--- a/code/IO.f90
+++ b/code/IO.f90
@@ -1552,6 +1552,8 @@ subroutine IO_error(error_ID,el,ip,g,ext_msg)
    msg = 'transformed hex lattice structure with invalid c/a ratio'
  case (135_pInt)
    msg = 'zero entry on stiffness diagonal'
+ case (136_pInt)
+   msg = 'zero entry on stiffness diagonal for transformed phase'
 
 !--------------------------------------------------------------------------------------------------
 ! material error messages and related messages in mesh
diff --git a/code/lattice.f90 b/code/lattice.f90
index d5e8fe155..774ad6713 100644
--- a/code/lattice.f90
+++ b/code/lattice.f90
@@ -1642,6 +1642,37 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc)
    if (abs(lattice_C66(i,i,myPhase))<tol_math_check) &
      call IO_error(135_pInt,el=i,ip=myPhase,ext_msg='matrix diagonal "el"ement of phase "ip"')
  enddo
+
+ ! Elasticity matrices for transformed phase
+ select case(lattice_structure(myPhase))
+   case (LATTICE_fcc_ID)
+     select case(trans_lattice_structure(myPhase))
+       case (LATTICE_bcc_ID)
+         lattice_trans_C66(1:6,1:6,myPhase) = lattice_C66(1:6,1:6,myPhase)
+         lattice_trans_mu(myPhase) = lattice_mu(myPhase)
+         lattice_trans_nu(myPhase) = lattice_nu(myPhase)
+         lattice_trans_C3333(1:3,1:3,1:3,1:3,myPhase) = lattice_C3333(1:3,1:3,1:3,1:3,myPhase)
+       case (LATTICE_hex_ID)
+         lattice_trans_C66(1:6,1:6,myPhase) = lattice_symmetrizeC66(trans_lattice_structure(myPhase),&
+                                                                   lattice_trans_C66(1:6,1:6,myPhase))
+         lattice_trans_mu(myPhase) = 0.2_pReal *(  lattice_trans_C66(1,1,myPhase) &
+                                                 - lattice_trans_C66(1,2,myPhase) &
+                                                 + 3.0_pReal*lattice_trans_C66(4,4,myPhase))
+         lattice_trans_nu(myPhase) = (  lattice_trans_C66(1,1,myPhase) &
+                                      + 4.0_pReal*lattice_trans_C66(1,2,myPhase) &
+                                      - 2.0_pReal*lattice_trans_C66(4,4,myPhase)) &
+                                                     /(  4.0_pReal*lattice_trans_C66(1,1,myPhase) &
+                                                       + 6.0_pReal*lattice_trans_C66(1,2,myPhase) &
+                                                       + 2.0_pReal*lattice_trans_C66(4,4,myPhase))
+         lattice_trans_C3333(1:3,1:3,1:3,1:3,myPhase) = math_Voigt66to3333(lattice_trans_C66(1:6,1:6,myPhase))
+         lattice_trans_C66(1:6,1:6,myPhase) = math_Mandel3333to66(lattice_trans_C3333(1:3,1:3,1:3,1:3,myPhase))
+         do i = 1_pInt, 6_pInt
+           if (abs(lattice_trans_C66(i,i,myPhase))<tol_math_check) &
+           call IO_error(136_pInt,el=i,ip=myPhase,ext_msg='matrix diagonal "el"ement of phase "ip"')
+         enddo
+     end select
+ end select
+
  lattice_thermalConductivity33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),&
                                                                        lattice_thermalConductivity33(1:3,1:3,myPhase))
  lattice_thermalExpansion33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),&
@@ -1686,9 +1717,9 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc)
      select case(trans_lattice_structure(myPhase))
        case (LATTICE_bcc_ID)                                                                         ! fcc to bcc transformation                         
          do i = 1_pInt,myNtrans
-           Rtr(1:3,1:3,i) = math_axisAngleToR(lattice_fccTobcc_systemTrans(1:3,i), &                      ! Pitsch rotation
+           Rtr(1:3,1:3,i) = math_axisAngleToR(lattice_fccTobcc_systemTrans(1:3,i), &                 ! Pitsch rotation
                               lattice_fccTobcc_systemTrans(4,i)*INRAD)
-           Btr(1:3,1:3,i) = math_axisAngleToR(lattice_fccTobcc_bainRot(1:3,i), &                          ! Rotation of fcc to Bain coordinate system
+           Btr(1:3,1:3,i) = math_axisAngleToR(lattice_fccTobcc_bainRot(1:3,i), &                     ! Rotation of fcc to Bain coordinate system
                               lattice_fccTobcc_bainRot(4,i)*INRAD)
            xtr(1:3,i) = real(LATTICE_fccTobcc_bainVariant(1:3,i),pReal)
            ytr(1:3,i) = real(LATTICE_fccTobcc_bainVariant(4:6,i),pReal)
diff --git a/code/plastic_dislotwin.f90 b/code/plastic_dislotwin.f90
index 0ed40590c..844a884a2 100644
--- a/code/plastic_dislotwin.f90
+++ b/code/plastic_dislotwin.f90
@@ -982,7 +982,7 @@ subroutine plastic_dislotwin_init(fileUnit)
            do p = 1_pInt,3_pInt; do q = 1_pInt,3_pInt; do r = 1_pInt,3_pInt; do s = 1_pInt,3_pInt
              plastic_dislotwin_Ctrans3333(l,m,n,o,index_myFamily+j,instance) = &
              plastic_dislotwin_Ctrans3333(l,m,n,o,index_myFamily+j,instance) + &
-               lattice_C3333(p,q,r,s,instance) * &
+               lattice_trans_C3333(p,q,r,s,instance) * &
                lattice_Qtrans(l,p,index_otherFamily+j,phase) * &
                lattice_Qtrans(m,q,index_otherFamily+j,phase) * &
                lattice_Qtrans(n,r,index_otherFamily+j,phase) * &
@@ -993,14 +993,17 @@ subroutine plastic_dislotwin_init(fileUnit)
            math_Mandel3333to66(plastic_dislotwin_Ctrans3333(1:3,1:3,1:3,1:3,index_myFamily+j,instance))
 
          !* Projection matrices for shear from slip systems to fault-band (twin) systems for strain-induced martensite nucleation
-         do o = 1_pInt,lattice_maxNtransFamily
-           index_otherFamily = sum(plastic_dislotwin_Nslip(1:o-1_pInt,instance))
-           do k = 1_pInt,plastic_dislotwin_Nslip(o,instance)                                   ! loop over (active) systems in other family (trans)
-             plastic_dislotwin_projectionMatrix_Trans(index_myFamily+j,index_otherFamily+k,instance) = &
-                   lattice_projectionTrans( sum(lattice_NtransSystem(1:f-1,phase))+j, &
-                                            sum(lattice_NslipSystem(1:o-1,phase))+k, phase)
-         enddo; enddo
- 
+         select case(trans_lattice_structure(phase))
+           case (LATTICE_bcc_ID)
+             do o = 1_pInt,lattice_maxNtransFamily
+               index_otherFamily = sum(plastic_dislotwin_Nslip(1:o-1_pInt,instance))
+               do k = 1_pInt,plastic_dislotwin_Nslip(o,instance)                                   ! loop over (active) systems in other family (trans)
+                 plastic_dislotwin_projectionMatrix_Trans(index_myFamily+j,index_otherFamily+k,instance) = &
+                       lattice_projectionTrans( sum(lattice_NtransSystem(1:f-1,phase))+j, &
+                                                sum(lattice_NslipSystem(1:o-1,phase))+k, phase)
+             enddo; enddo
+         end select 
+
        enddo transSystemsLoop
      enddo transFamiliesLoop