From 9c1cce828d28676d3f0dee539e2c03b258b530e0 Mon Sep 17 00:00:00 2001 From: achalhp Date: Wed, 24 Jan 2024 09:27:42 +0530 Subject: [PATCH] Added neighbor loop --- ...phase_mechanical_plastic_phenopowerlaw.f90 | 88 ++++++++++++++++++- 1 file changed, 86 insertions(+), 2 deletions(-) diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90 index ced0576a7..0ad9cd918 100644 --- a/src/phase_mechanical_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -6,6 +6,23 @@ !-------------------------------------------------------------------------------------------------- submodule(phase:plastic) phenopowerlaw + use geometry_plastic_nonlocal, only: & !< Achal Copied from nonlocal + 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, & + geometry_plastic_nonlocal_disable + + type :: tGeometry + real(pReal), dimension(:), allocatable :: V_0 + integer, dimension(:,:,:), allocatable :: IPneighborhood + real(pReal), dimension(:,:), allocatable :: IParea, IPcoordinates + real(pReal), dimension(:,:,:), allocatable :: IPareaNormal + end type tGeometry + + type(tGeometry), dimension(:), allocatable :: geom + type :: tParameters real(pReal) :: & dot_gamma_0_sl = 1.0_pReal, & !< reference shear strain rate for slip @@ -111,6 +128,9 @@ print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) phases => config_material%get('phase') + +allocate(geom(phases%length)) !< Achal + allocate(param(phases%length)) allocate(indexDotState(phases%length)) allocate(state(phases%length)) @@ -248,6 +268,13 @@ do ph = 1, phases%length call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState) deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely + allocate(geom(ph)%V_0(Nmembers)) !Achal + allocate(geom(ph)%IPneighborhood(3,nIPneighbors,Nmembers)) !Achal + allocate(geom(ph)%IPareaNormal(3,nIPneighbors,Nmembers)) + allocate(geom(ph)%IParea(nIPneighbors,Nmembers)) + allocate(geom(ph)%IPcoordinates(3,Nmembers)) + call storeGeometry(ph) !Achal + !-------------------------------------------------------------------------------------------------- ! state aliases and initialization startIndex = 1 @@ -392,7 +419,7 @@ associate(prm => param(ph), stt => state(ph), & !write(6,*)'characteristicShearTwin', prm%gamma_char !write(6,*)'Schmid_twin',prm%P_sl !if(en==1) write(6,*)'maxF',maxval(stt%gamma_tw(:,en)/prm%gamma_char) ! delete Achal - if(en==1) write(6,*)'f_twin',indexDotState(ph)%f_twin + !if(en==1) write(6,*)'f_twin',fdot_twin !if(en==1) write(6,*)'f_twin',f_twin sumF = sum(stt%gamma_tw(:,en)/prm%gamma_char) @@ -430,8 +457,17 @@ logical , intent(out) :: & real(pReal), dimension(3,3), intent(out) :: & deltaFp + integer :: & + n, & ! neighbor index + en, & + neighbor_e, & ! element index of my neighbor + neighbor_i, & ! integration point index of my neighbor + neighbor_me, & + neighbor_phase + real(pReal) :: & - random + random, & + nRealNeighbors !Achal number of really existing neighbors integer :: & twin_var real(pReal), dimension(param(ph)%sum_N_tw) :: & @@ -442,6 +478,18 @@ integer :: i twinJump = .false. deltaFp = math_I3 +!* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities !Achal + +nRealNeighbors = 0.0_pReal !Achal + +neighbors: do n = 1,nIPneighbors + neighbor_e = geom(ph)%IPneighborhood(1,n,en) + neighbor_i = geom(ph)%IPneighborhood(2,n,en) + neighbor_me = material_phaseEntry(1,(neighbor_e-1)*discretization_nIPs + neighbor_i) !Neighbour offset + neighbor_phase = material_phaseID(1,(neighbor_e-1)*discretization_nIPs + neighbor_i) + +enddo neighbors + tau_tw = [(math_tensordot(Mp,param(ph)%P_tw(1:3,1:3,i)),i=1,param(ph)%sum_N_tw)] twin_var = maxloc((0.05_pReal*(abs(tau_tw)/state(ph)%xi_tw(:,en))**param(ph)%n_tw)/param(ph)%gamma_char,dim=1) ! This prints values from 1 to 6, fdot0_twin is taken as 0.05 fdot_twin = (0.05_pReal*(abs(tau_tw)/state(ph)%xi_tw(:,en))**param(ph)%n_tw)/param(ph)%gamma_char ! This is sometimes >1 @@ -643,4 +691,40 @@ end associate end subroutine kinetics_tw +!-------------------------------------------------------------------------------------------------- +!> Achal: Added from nonlocal +!-------------------------------------------------------------------------------------------------- + +subroutine storeGeometry(ph) + + integer, intent(in) :: ph + + integer :: ce, co, nCell + real(pReal), dimension(:), allocatable :: V + integer, dimension(:,:,:), allocatable :: neighborhood + real(pReal), dimension(:,:), allocatable :: area, coords + real(pReal), dimension(:,:,:), allocatable :: areaNormal + + nCell = product(shape(IPvolume)) + + V = reshape(IPvolume,[nCell]) + neighborhood = reshape(IPneighborhood,[3,nIPneighbors,nCell]) + area = reshape(IParea,[nIPneighbors,nCell]) + areaNormal = reshape(IPareaNormal,[3,nIPneighbors,nCell]) + coords = reshape(discretization_IPcoords,[3,nCell]) + + do ce = 1, size(material_homogenizationEntry,1) + do co = 1, homogenization_maxNconstituents + if (material_phaseID(co,ce) == ph) then + geom(ph)%V_0(material_phaseEntry(co,ce)) = V(ce) + geom(ph)%IPneighborhood(:,:,material_phaseEntry(co,ce)) = neighborhood(:,:,ce) + geom(ph)%IParea(:,material_phaseEntry(co,ce)) = area(:,ce) + geom(ph)%IPareaNormal(:,:,material_phaseEntry(co,ce)) = areaNormal(:,:,ce) + geom(ph)%IPcoordinates(:,material_phaseEntry(co,ce)) = coords(:,ce) + end if + end do + end do + +end subroutine + end submodule phenopowerlaw