diff --git a/src/io.mod b/src/io.mod new file mode 100644 index 000000000..c51696395 Binary files /dev/null and b/src/io.mod differ diff --git a/src/lattice.f90 b/src/lattice.f90 index 7bd846d0f..0fda143d1 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -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 diff --git a/src/math.f90 b/src/math.f90 index a875741b3..2902ed430 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -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 diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90 index 63f7aa37f..1005a9912 100644 --- a/src/phase_mechanical_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -280,18 +280,18 @@ 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 + deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely !-------------------------------------------------------------------------------------------------- ! state aliases and initialization @@ -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) diff --git a/src/prec.mod b/src/prec.mod new file mode 100644 index 000000000..9a27e158e Binary files /dev/null and b/src/prec.mod differ