Correspondence Matrix applied

src/crystal.f90

@@ -404,6 +404,7 @@ module crystal
     crystal_slip_direction, &
     crystal_slip_transverse, &
     crystal_labels_slip, &
+    crystal_CorrespondenceMatrix_twin, &
     crystal_labels_twin
 
 contains
@@ -2111,6 +2112,60 @@ function getlabels(active,potential,system) result(labels)
 
 end function getlabels
 
+!--------------------------------------------------------------------------------------------------
+!> @brief correspondence matrix for twinning
+!> details only active twin systems are considered
+!--------------------------------------------------------------------------------------------------
+function crystal_CorrespondenceMatrix_twin(Ntwin,lattice,cOverA) result(CorrespondenceMatrix)
+
+  integer,    dimension(:),            intent(in) :: Ntwin                                     !< number of active twin systems per family
+  character(len=*),                          intent(in) :: lattice                                 !< lattice structure
+  real(pReal),                               intent(in) :: cOverA                                    !< c/a ratio
+  real(pReal),     dimension(3,3,sum(Ntwin))            :: CorrespondenceMatrix
+
+  real(pReal),     dimension(3,3,sum(Ntwin))            :: coordinateSystem
+  real(pReal),     dimension(sum(Ntwin))                :: characteristicShearTwin
+  real(pReal),     dimension(3,3,sum(Ntwin))            :: SchmidMatrixTwin
+  real(pReal),     dimension(:,:), allocatable          :: twinSystems
+  integer,   dimension(:),   allocatable          :: NtwinMax
+  integer :: i
+
+  select case(lattice)
+    case('cF')
+      NtwinMax    = CF_NTWINSYSTEM
+      twinSystems = CF_SYSTEMTWIN
+    case('cI')
+      NtwinMax    = CI_NTWINSYSTEM
+      twinSystems = CI_SYSTEMTWIN
+    case('hP')
+      NtwinMax    = HP_NTWINSYSTEM
+      twinSystems = HP_SYSTEMTWIN                                                              !< the twin system matrix is different from V2.0
+    case default
+      allocate(NtwinMax(0))
+      call IO_error(137,ext_msg='crystal_CorrespondenceMatrix_twin: '//trim(lattice))
+  end select
+
+  if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0)) &
+    call IO_error(145,ext_msg='Ntwin '//trim(lattice))
+  if (any(Ntwin < 0)) &
+    call IO_error(144,ext_msg='Ntwin '//trim(lattice))
+
+  coordinateSystem        =  buildCoordinateSystem(Ntwin,NtwinMax,twinSystems,lattice,cOverA)
+  !  characteristicShearTwin =  0.0_pReal*lattice_characteristicShear_Twin(Ntwin,lattice,cOverA)          ! for removing shear from CorrespondenceMatrix
+  characteristicShearTwin =  crystal_characteristicShear_Twin(Ntwin,lattice,cOverA)
+  SchmidMatrixTwin        =  crystal_SchmidMatrix_twin(Ntwin,lattice,cOverA)
+
+  !write(6,*)'coordinate system', coordinateSystem(1:3,2,1)
+  
+  !CorrespondenceMatrix(1:3,1:3,1) = math_axisAngleToR(coordinateSystem(1:3,2,6), 180.0_pReal*INRAD)                           ! delete this
+
+  do i = 1, sum(Ntwin)
+    CorrespondenceMatrix(1:3,1:3,i) = matmul(math_axisAngleToR(coordinateSystem(1:3,2,i), &
+                                      180.0_pReal*INRAD), MATH_I3 + characteristicShearTwin(i)* &
+                                      SchmidMatrixTwin(1:3,1:3,i))
+  enddo
+
+end function crystal_CorrespondenceMatrix_twin
 
 !--------------------------------------------------------------------------------------------------
 !> @brief Equivalent Poisson's ratio (ν)
@@ -2307,4 +2362,7 @@ subroutine crystal_selfTest()
 
 end subroutine crystal_selfTest
 
+
+
+
 end module crystal

src/math.f90

@@ -1470,4 +1470,46 @@ subroutine math_selfTest()
 
 end subroutine math_selfTest
 
+!--------------------------------------------------------------------------------------------------
+!> @brief rotation matrix from axis and angle (in radians)
+!> @details rotation matrix is meant to represent a ACTIVE rotation
+!> @details (see http://en.wikipedia.org/wiki/Euler_angles for definitions)
+!> @details formula for active rotation taken from http://mathworld.wolfram.com/RodriguesRotationFormula.html
+!> @details equivalent to eu2om (P=-1) from "D Rowenhorst et al. Consistent representations of and
+!> @details conversions between 3D rotations, Model. Simul. Mater. Sci. Eng. 23-8 (2015)"
+!--------------------------------------------------------------------------------------------------
+pure function math_axisAngleToR(axis,omega)
+
+ implicit none
+ real(pReal), dimension(3,3) :: math_axisAngleToR
+ real(pReal), dimension(3), intent(in) :: axis
+ real(pReal), intent(in) :: omega
+ real(pReal), dimension(3) :: n
+ real(pReal) :: norm,s,c,c1
+
+ norm = norm2(axis)
+ wellDefined: if (norm > 1.0e-8_pReal) then
+   n = axis/norm                                                                                    ! normalize axis to be sure
+
+   s = sin(omega)
+   c = cos(omega)
+   c1 = 1.0_pReal - c
+
+   math_axisAngleToR(1,1) =  c + c1*n(1)**2.0_pReal
+   math_axisAngleToR(1,2) =  c1*n(1)*n(2) - s*n(3)
+   math_axisAngleToR(1,3) =  c1*n(1)*n(3) + s*n(2)
+                             
+   math_axisAngleToR(2,1) =  c1*n(1)*n(2) + s*n(3)
+   math_axisAngleToR(2,2) =  c + c1*n(2)**2.0_pReal
+   math_axisAngleToR(2,3) =  c1*n(2)*n(3) - s*n(1)
+                             
+   math_axisAngleToR(3,1) =  c1*n(1)*n(3) - s*n(2)
+   math_axisAngleToR(3,2) =  c1*n(2)*n(3) + s*n(1)
+   math_axisAngleToR(3,3) =  c + c1*n(3)**2.0_pReal
+ else wellDefined
+   math_axisAngleToR = math_I3
+ endif wellDefined
+ 
+end function math_axisAngleToR
+
 end module math

src/phase_mechanical_plastic_phenopowerlaw.f90

@@ -32,7 +32,8 @@ submodule(phase:plastic) phenopowerlaw
       P_sl, &
       P_tw, &
       P_nS_pos, &
-      P_nS_neg
+      P_nS_neg, &
+      CorrespondenceMatrix
     integer :: &
       sum_N_sl, &                                                                                   !< total number of active slip system
       sum_N_tw                                                                                      !< total number of active twin systems
@@ -203,6 +204,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
       prm%c_4            = math_expand(pl%get_as1dReal('c_4',            requiredSize=size(N_tw), &
                                                                          defaultVal=misc_zeros(size(N_tw))), N_tw)
 
+      prm%CorrespondenceMatrix = crystal_CorrespondenceMatrix_twin(N_tw,phase_lattice(ph),phase_cOverA(ph))
       prm%gamma_char = crystal_characteristicShear_twin(N_tw,phase_lattice(ph),phase_cOverA(ph))
       prm%h_tw_tw    = crystal_interaction_TwinByTwin(N_tw,pl%get_as1dReal('h_tw-tw'),phase_lattice(ph))
 
@@ -225,6 +227,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
                prm%h_0_tw_tw, &
                source=emptyRealArray)
       allocate(prm%h_tw_tw(0,0))
+      allocate(prm%CorrespondenceMatrix(0,0,0))
     end if twinActive
 
 !--------------------------------------------------------------------------------------------------
@@ -511,7 +514,7 @@ pure subroutine kinetics_tw(Mp,ph,en,&
 
   associate(prm => param(ph), stt => state(ph))
 
-    tau_tw = [(math_tensordot(Mp,prm%P_tw(1:3,1:3,i)),i=1,prm%sum_N_tw)]
+    tau_tw = [(math_tensordot(Mp,prm%CorrespondenceMatrix(1:3,1:3,i)),i=1,prm%sum_N_tw)]
 
     where(tau_tw > 0.0_pREAL)
       dot_gamma_tw = (1.0_pREAL-sum(stt%gamma_tw(:,en)/prm%gamma_char)) &                           ! only twin in untwinned volume fraction