diff --git a/src/lattice.f90 b/src/lattice.f90 index 7bd846d0f..8d8ca6be1 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -2283,4 +2283,67 @@ 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 f0dc04869..6904f4790 100644 --- a/src/phase_mechanical_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -377,6 +377,169 @@ module function phenopowerlaw_dotState(Mp,ph,en) result(dotState) end function phenopowerlaw_dotState +!-------------------------------------------------------------------------------------------------- +!> @brief calculates instantaneous incremental change of kinematics and associated jump state +!-------------------------------------------------------------------------------------------------- +module subroutine plastic_kinematic_deltaFp(twinJump,deltaFp,ipc, ip, el) + use prec, only: & + dNeq, & + dEq0 +! #ifdef DEBUG +! use debug, only: & +! debug_level, & +! debug_constitutive,& +! debug_levelExtensive, & +! debug_levelSelective +! #endif + + use geometry_plastic_nonlocal, only: & + nIPneighbors => geometry_plastic_nonlocal_nIPneighbors, & + IPneighborhood => geometry_plastic_nonlocal_IPneighborhood, & + IPvolume => geometry_plastic_nonlocal_IPvolume0, & + IParea => geometry_plastic_nonlocal_IParea0, & + IPareaNormal => geometry_plastic_nonlocal_IPareaNormal0 + + + !use mesh, only: & + ! mesh_element, & !name changed + ! mesh_ipNeighborhood, & + ! mesh_ipCoordinates, & + ! mesh_ipVolume, & + ! mesh_ipAreaNormal, & + ! mesh_ipArea, & + ! FE_NipNeighbors, & + ! mesh_maxNipNeighbors, & + ! FE_geomtype, & + ! FE_celltype + + use lattice + use math, only: & + math_I3 + +! use material, only: & +! phaseAt, phasememberAt, & !name changed +! phase_plasticityInstance + + implicit none + integer :: & + ph, of, instance, & + neighbor_el, & !< element number of neighboring material point + neighbor_ip, & !< integration point of neighboring material point + np, & !< neighbor phase + no, n !< nieghbor offset and index for loop at neighbor + + logical , intent(out) :: & + twinJump + + real(pReal), dimension(3,3), intent(out) :: & + deltaFp + + integer, intent(in) :: & + ipc, & !< element index + ip, & !< integration point index + el !< grain index +! ! real(pReal), dimension(3,3,param(instance)%totalNslip) :: & +! ! CorrespondanceMatrix + integer, dimension(52) :: & + twin_el_incl + real(pReal), dimension(6) :: & + neighbor_stt + real(pReal) :: & + random, random1 + integer :: & + i,j,var_growth,var_nucl + var_growth = 0 + var_nucl = 0 + !ph = phaseAt(ipc, ip, el) + !of = phasememberAt(ipc, ip, el) + !instance = phase_plasticityInstance(ph) + +! associate(prm => param(instance), stt => state(instance), dlt => deltaState(instance)) + +! twinJump = .false. +! deltaFp = math_I3 + +! ! for eshelby circular inclusion +! twin_el_incl = (/ 10913,10914,10915,10916,10917,10918,10919,10920,10921,10922,10923,10924,10925, & +! 10993,10994,10995,10996,10997,10998,10999,11000,11001,11002,11074,11075,11076, & +! 11077,11078,11079,10751,10752,10753,10754,10755,10756,10757,10758,10759,10760, & +! 10670,10671,10672,10673,10674,10675,10676,10677,10678,10679,10680,10681,10682 /) +! ! TwinLooptest: do i=1_pInt, prm%totalNtwin +! ! write(6,*)'CorrespondenceMatrix for system',i, prm%CorrespondanceMatrix(:,:,i) +! ! enddo TwinLooptest + + +!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) ! Integer +! neighbor_ip = mesh_ipNeighborhood(2,n,ip,el) ! Integer +! np = phaseAt(1,neighbor_ip,neighbor_el) ! Integer +! no = phasememberAt(1,neighbor_ip,neighbor_el) ! Integer +! neighbor_stt(n) = state(phase_plasticityInstance(np))%variant_twin(no) ! Integer +! enddo NeighborLoop1 + +! !checking if any of my neighbor is twinned if yes recognize the variant and exit +! ! NeighborLoop2: 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) +! ! np = phaseAt(1,neighbor_ip,neighbor_el) +! ! ! if(of == 1) write(6,*)'phaseAt neighbor_ip of neighbor_el', np +! ! no = phasememberAt(1,neighbor_ip,neighbor_el) +! ! ! if(of == 1) write(6,*)'phasememberAt at neighbor_ip of neighbor_el', no +! ! if (state(phase_plasticityInstance(np))%variant_twin(no) > 0_pInt) then +! ! var_growth = state(phase_plasticityInstance(np))%variant_twin(no) +! ! +! ! exit NeighborLoop2 +! ! endif +! ! enddo NeighborLoop2 + +! call RANDOM_NUMBER(random) +! call RANDOM_NUMBER(random1) +! ! Sampling: if (var_growth > 0_pInt) then +! ! ! write(6,*)'I am sampling for growth with variant',var_growth +! ! Ability_Growth: if (stt%f_twin_grow(var_growth,of) > stt%fmc_twin_grow(var_growth,of) & +! ! + prm%checkstep_grow) then +! ! stt%fmc_twin_grow(var_growth,of) = stt%fmc_twin_grow(var_growth,of) & +! ! + prm%checkstep_grow +! ! Success_Growth: if (random <= stt%f_twin_grow(var_growth,of) .or. ALL(neighbor_stt > 0_pReal)) then +! ! write(6,*)'growth sampling is successful for elem',el +! ! twinJump = .true. +! ! deltaFp = prm%CorrespondanceMatrix(:,:,var_growth) +! ! dlt%f_twin_grow(:,of) = 0.0_pReal - stt%f_twin_grow(:,of) +! ! dlt%f_twin_nucl(:,of) = 0.0_pReal - stt%f_twin_nucl(:,of) +! ! dlt%fmc_twin_grow(:,of) = 0.0_pReal - stt%fmc_twin_grow(:,of) +! ! dlt%fmc_twin_nucl(:,of) = 0.0_pReal - stt%fmc_twin_nucl(:,of) +! ! dlt%frozen(of) = 1.0_pReal - stt%frozen(of) +! ! dlt%variant_twin(of) = var_growth - stt%variant_twin(of) +! ! endif Success_Growth +! ! endif Ability_Growth + +! ! elseif (var_growth == 0_pInt .and. prm%checkgrowth_twin > 0_pReal ) then +! if (var_growth == 0_pInt .and. prm%checkgrowth_twin > 0_pReal ) then +! var_nucl = maxloc(stt%f_twin_nucl(:,of), dim=1) +! ! write(6,*)'I am sampling for nucleation with variant',var_nucl,stt%f_twin_nucl(var_nucl,of) +! Ability_Nucleation: if (stt%f_twin_nucl(var_nucl,of) > stt%fmc_twin_nucl(var_nucl,of) & +! + prm%checkstep_nucl) then +! stt%fmc_twin_nucl(var_nucl,of) = stt%fmc_twin_nucl(var_nucl,of) & +! + prm%checkstep_nucl +! Success_Nucleation: if (random <= stt%f_twin_nucl(var_nucl,of) & +! .and. random1 <= 0.20) then +! write(6,*)'nucleation sampling is successful for elem',el +! twinJump = .true. +! deltaFp = prm%CorrespondanceMatrix(:,:,var_nucl) +! dlt%f_twin_nucl(:,of) = 0.0_pReal - stt%f_twin_nucl(:,of) +! dlt%f_twin_grow(:,of) = 0.0_pReal - stt%f_twin_grow(:,of) +! dlt%fmc_twin_nucl(:,of) = 0.0_pReal - stt%fmc_twin_nucl(:,of) +! dlt%fmc_twin_grow(:,of) = 0.0_pReal - stt%fmc_twin_grow(:,of) +! dlt%frozen(of) = 1.0_pReal - stt%frozen(of) +! dlt%variant_twin(of) = var_nucl - stt%variant_twin(of) +! endif Success_Nucleation +! endif Ability_Nucleation +! endif +! ! endif Sampling +! end associate + +end subroutine plastic_kinematic_deltaFp !-------------------------------------------------------------------------------------------------- !> @brief Write results to HDF5 output file.