Added correspondance matrix

This commit is contained in:
achalhp 2023-11-02 22:15:15 +05:30
parent 9eeee3d39d
commit 2c13ca53bc
5 changed files with 127 additions and 20 deletions

BIN
src/io.mod Normal file

Binary file not shown.

View File

@ -400,7 +400,8 @@ module lattice
lattice_slip_direction, &
lattice_slip_transverse, &
lattice_labels_slip, &
lattice_labels_twin
lattice_labels_twin, &
lattice_CorrespondanceMatrix_twin !< Achal
contains
@ -2283,4 +2284,68 @@ subroutine selfTest
end subroutine selfTest
!--------------------------------------------------------------------------------------------------
!> @brief correspondance matrix for twinning
!> details only active twin systems are considered
!--------------------------------------------------------------------------------------------------
function lattice_CorrespondanceMatrix_twin(Ntwin,structure,cOverA) result(CorrespondanceMatrix)
use prec, only: &
tol_math_check
use IO, only: &
IO_error
use math, only: &
math_mul3333xx33, &
math_axisAngleToR, &
INRAD, &
math_I3
implicit none
integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family
character(len=*), intent(in) :: structure !< lattice structure
real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(3,3,sum(Ntwin)) :: CorrespondanceMatrix
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
if (len_trim(structure) /= 3) &
call IO_error(137,ext_msg='lattice_CorrespondanceMatrix_twin: '//trim(structure))
select case(structure(1:3))
case('fcc')
NtwinMax = CF_NTWINSYSTEM
twinSystems = CF_SYSTEMTWIN
case('bcc')
NtwinMax = CI_NTWINSYSTEM
twinSystems = CI_SYSTEMTWIN
case('hex')
NtwinMax = HP_NTWINSYSTEM
twinSystems = HP_SYSTEMTWIN !< the twin system matrix is different from V2.0
case default
call IO_error(137,ext_msg='lattice_CorrespondanceMatrix_twin: '//trim(structure))
end select
if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0)) &
call IO_error(145,ext_msg='Ntwin '//trim(structure))
if (any(Ntwin < 0)) &
call IO_error(144,ext_msg='Ntwin '//trim(structure))
coordinateSystem = buildCoordinateSystem(Ntwin,NtwinMax,twinSystems,structure,cOverA)
! characteristicShearTwin = 0.0_pReal*lattice_characteristicShear_Twin(Ntwin,structure,cOverA) ! for removing shear from CorrespondanceMatrix
characteristicShearTwin = lattice_characteristicShear_Twin(Ntwin,structure,cOverA)
SchmidMatrixTwin = lattice_SchmidMatrix_twin(Ntwin,structure,cOverA)
do i = 1, sum(Ntwin)
CorrespondanceMatrix(1:3,1:3,i) = math_mul3333xx33(math_axisAngleToR(coordinateSystem(1:3,2,i), &
180.0_pReal*INRAD), MATH_I3 + characteristicShearTwin(i)* &
SchmidMatrixTwin(1:3,1:3,i))
enddo
end function lattice_CorrespondanceMatrix_twin
end module lattice

View File

@ -1488,4 +1488,46 @@ subroutine selfTest()
end subroutine 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

View File

@ -280,15 +280,15 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
Nmembers = count(material_phaseID == ph)
sizeDotState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl &
+ size(['xi_tw ','gamma_tw']) * prm%sum_N_tw &
+ size(['xi_tw_nucl','xi_tw_grow']) * prm%sum_N_tw !& ! Why not size(['xi_tw_nucl','gamma_tw'])?
+ size(['xi_tw_nucl','xi_tw_grow']) * prm%sum_N_tw !& ! Achal Why not size(['xi_tw_nucl','gamma_tw'])?
!+ size(['f_tw_nucl','f_tw_grow']) * prm%sum_N_tw &
!+ 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(['xi_tw_nucl','xi_tw_grow']) * prm%sum_N_tw !& ! Why not size(['xi_tw_nucl','gamma_tw'])?
!+ size(['f_tw_nucl','f_tw_grow']) * prm%sum_N_tw &
!+ size(['fmc_tw_nucl','fmc_tw_grow']) * prm%sum_N_tw &
!+ size(['variant_twin','frozen']) * prm%sum_N_tw &
+ size(['xi_tw_nucl','xi_tw_grow']) * prm%sum_N_tw !& ! Achal Why not size(['xi_tw_nucl','gamma_tw'])?
!+ size(['f_tw_nucl','f_tw_grow']) * prm%sum_N_tw & !Achal
!+ size(['fmc_tw_nucl','fmc_tw_grow']) * prm%sum_N_tw & !Achal
!+ size(['variant_twin','frozen']) * prm%sum_N_tw & !Achal
call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely
@ -323,19 +323,19 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
stt%gamma_tw => plasticState(ph)%state(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw
idx_dot%xi_tw_nucl = [startIndex,endIndex]
stt%xi_tw_nucl => plasticState(ph)%state(startIndex:endIndex,:)
stt%xi_tw_nucl = spread(xi_0_tw, 2, Nmembers)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
startIndex = endIndex + 1 !Achal
endIndex = endIndex + prm%sum_N_tw !Achal
idx_dot%xi_tw_nucl = [startIndex,endIndex] !Achal
stt%xi_tw_nucl => plasticState(ph)%state(startIndex:endIndex,:) !Achal
stt%xi_tw_nucl = spread(xi_0_tw, 2, Nmembers) !Achal
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) !Achal
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw
idx_dot%xi_tw_grow = [startIndex,endIndex]
stt%xi_tw_grow => plasticState(ph)%state(startIndex:endIndex,:)
stt%xi_tw_grow = spread(xi_0_tw, 2, Nmembers)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
startIndex = endIndex + 1 !Achal
endIndex = endIndex + prm%sum_N_tw !Achal
idx_dot%xi_tw_grow = [startIndex,endIndex] !Achal
stt%xi_tw_grow => plasticState(ph)%state(startIndex:endIndex,:) !Achal
stt%xi_tw_grow = spread(xi_0_tw, 2, Nmembers) !Achal
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) !Achal
!startIndex = endIndex + 1
!endIndex = endIndex + prm%sum_N_tw
@ -601,7 +601,7 @@ module subroutine plastic_kinematic_deltaFp(twinJump,deltaFp,ipc, ip, el)
! ! enddo TwinLooptest
! !Saving the neighbor information in an array
!Saving the neighbor information in an array
! NeighborLoop1: do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el)))) ! only 4 neighbors for quasi 2D (1 element in z direction)
! neighbor_el = mesh_ipNeighborhood(1,n,ip,el)
! neighbor_ip = mesh_ipNeighborhood(2,n,ip,el)

BIN
src/prec.mod Normal file

Binary file not shown.