From ff9fa1d4f7eb9b318b5b35cc804ba2d0add5f4ff Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Fri, 19 Nov 2021 07:29:40 +0100
Subject: [PATCH] using Voigt notation instead of proprietary scaled 6x6
 notation

Note: This results in a change of behavior for the transformation
systems of dislotwin. I assume that this fixes a bug, but still need to
confirm where the equations in lattice_C66_trans come from
---
 src/lattice.f90                            |  8 ++++----
 src/phase_mechanical_elastic.f90           |  5 +++--
 src/phase_mechanical_plastic_dislotwin.f90 | 18 ++++++++++--------
 3 files changed, 17 insertions(+), 14 deletions(-)

diff --git a/src/lattice.f90 b/src/lattice.f90
index 7fa6b8ab5..4b87fccbf 100644
--- a/src/lattice.f90
+++ b/src/lattice.f90
@@ -522,7 +522,7 @@ function lattice_C66_twin(Ntwin,C66,lattice,CoverA)
 
   do i = 1, sum(Ntwin)
     call R%fromAxisAngle([coordinateSystem(1:3,2,i),PI],P=1)                                        ! ToDo: Why always 180 deg?
-    lattice_C66_twin(1:6,1:6,i) = math_sym3333to66(R%rotTensor4(math_66toSym3333(C66)))
+    lattice_C66_twin(1:6,1:6,i) = math_3333toVoigt66(R%rotTensor4(math_Voigt66to3333(C66)))
   enddo
 
 end function lattice_C66_twin
@@ -572,16 +572,16 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
     call IO_error(137,ext_msg='lattice_C66_trans : '//trim(lattice_target))
   endif
 
-  do i = 1, 6
+  do i = 1,6
     if (abs(C_target_unrotated66(i,i))<tol_math_check) &
     call IO_error(135,el=i,ext_msg='matrix diagonal "el"ement in transformation')
   enddo
 
   call buildTransformationSystem(Q,S,Ntrans,cOverA_trans,a_fcc,a_bcc)
 
-  do i = 1, sum(Ntrans)
+  do i = 1,sum(Ntrans)
     call R%fromMatrix(Q(1:3,1:3,i))
-    lattice_C66_trans(1:6,1:6,i) = math_sym3333to66(R%rotTensor4(math_66toSym3333(C_target_unrotated66)))
+    lattice_C66_trans(1:6,1:6,i) = math_3333toVoigt66(R%rotTensor4(math_Voigt66to3333(C_target_unrotated66)))
   enddo
 
  end function lattice_C66_trans
diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90
index 38bca44b4..5792ea3cd 100644
--- a/src/phase_mechanical_elastic.f90
+++ b/src/phase_mechanical_elastic.f90
@@ -130,6 +130,7 @@ end function elastic_nu
 !--------------------------------------------------------------------------------------------------
 !> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
 !> the elastic and intermediate deformation gradients using Hooke's law
+! ToDo: Use Voigt matrix directly
 !--------------------------------------------------------------------------------------------------
 module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
                                               Fe, Fi, ph, en)
@@ -166,7 +167,7 @@ end subroutine phase_hooke_SandItsTangents
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief returns the homogenized elasticity matrix
+!> @brief Return the homogenized elasticity matrix.
 !--------------------------------------------------------------------------------------------------
 module function phase_homogenizedC66(ph,en) result(C)
 
@@ -176,7 +177,7 @@ module function phase_homogenizedC66(ph,en) result(C)
 
   plasticType: select case (phase_plasticity(ph))
     case (PLASTICITY_DISLOTWIN_ID) plasticType
-     C = math_3333toVoigt66(math_66toSym3333(plastic_dislotwin_homogenizedC(ph,en)))
+     C = plastic_dislotwin_homogenizedC(ph,en)
     case default plasticType
      C = elastic_C66(ph,en)
   end select plasticType
diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90
index 934773ad7..de027ae44 100644
--- a/src/phase_mechanical_plastic_dislotwin.f90
+++ b/src/phase_mechanical_plastic_dislotwin.f90
@@ -473,27 +473,29 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC)
   integer,     intent(in) :: &
     ph, en
   real(pReal), dimension(6,6) :: &
-    homogenizedC, &
-    C66
+    homogenizedC
+
+  real(pReal), dimension(6,6) :: &
+    C
   real(pReal), dimension(:,:,:), allocatable :: &
     C66_tw, &
     C66_tr
-
   integer :: i
   real(pReal) :: f_unrotated
 
-  associate(prm => param(ph), stt => state(ph))
 
-    C66 = math_sym3333to66(math_Voigt66to3333(elastic_C66(ph,en)))
+  C = elastic_C66(ph,en)
+
+  associate(prm => param(ph), stt => state(ph))
 
     f_unrotated = 1.0_pReal &
                 - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
                 - sum(stt%f_tr(1:prm%sum_N_tr,en))
 
-    homogenizedC = f_unrotated * C66
+    homogenizedC = f_unrotated * C
 
     twinActive: if (prm%sum_N_tw > 0) then
-      C66_tw    = lattice_C66_twin(prm%N_tw,C66,phase_lattice(ph),phase_cOverA(ph))
+      C66_tw    = lattice_C66_twin(prm%N_tw,C,phase_lattice(ph),phase_cOverA(ph))
       do i=1,prm%sum_N_tw
         homogenizedC = homogenizedC &
                      + stt%f_tw(i,en)*C66_tw(1:6,1:6,i)
@@ -501,7 +503,7 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC)
      end if twinActive 
     
     transActive: if (prm%sum_N_tr > 0) then
-      C66_tr    = lattice_C66_trans(prm%N_tr,C66,prm%lattice_tr,0.0_pReal,prm%a_cI,prm%a_cF)
+      C66_tr    = lattice_C66_trans(prm%N_tr,C,prm%lattice_tr,0.0_pReal,prm%a_cI,prm%a_cF)
       do i=1,prm%sum_N_tr
         homogenizedC = homogenizedC &
                      + stt%f_tr(i,en)*C66_tr(1:6,1:6,i)