aselimov/cacmb
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Working changes

Browse source
This commit is contained in:
Alex Selimov 2020-10-24 09:13:16 -04:00
parent 22e250093b
commit f63335708b
5 changed files with 65 additions and 52 deletions
Download patch file Download diff file Expand all files Collapse all files

61
src/opt_slip_plane.f90
View file

@ -17,11 +17,11 @@ module opt_slip_plane
!Main calling function for the slip_plane option
integer, intent(inout) :: arg_pos
integer :: ie, ia, slip_enum, old_atom_num, esize, new_ele_num, n, m, o, ele(3,8), nump_ele, inod, vlat(3), &
rfill(3,max_basisnum, max_ng_node)
integer :: ie, ia, slip_enum, old_atom_num, esize, new_ele_num, n, m, o, ele(3,8), nump_ele, inod, vlat(3)
integer, allocatable :: slip_eles(:), temp_int(:)
real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3)
real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), rfill(3,max_basisnum, max_ng_node)
integer :: type_interp(max_basisnum*max_esize**3)
logical :: lat_points(max_esize,max_esize, max_esize)
@ -45,7 +45,7 @@ module opt_slip_plane
!Now loop over all elements, find which ones intersect
do ie = 1, ele_num
if( (spos < maxval(r_node(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie)),ie))).and. &
(spos > minval(r_node(sdim,1:basisnum(lat_ele(ie),1:ng_node(lat_ele(ie)),ie)))) then
(spos > minval(r_node(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie)),ie)))) then
slip_enum = slip_enum + 1
if (slip_enum > size(slip_eles)) then
allocate(temp_int(size(slip_eles)+1024))
@ -64,35 +64,36 @@ module opt_slip_plane
end do
!If we are efilling then the code is slightly more complex
else
!First populate the lat points array
lat_points(1:size_ele(ie),1:size_ele(ie), 1:size_ele(ie)) = .true.
!First populate the lat points array
lat_points(1:size_ele(ie),1:size_ele(ie), 1:size_ele(ie)) = .true.
!Now start trying to remesh the region, leaving the slip plane as a discontinuity
esize = size_ele(ie) - 2
nump_ele = size_ele(ie)**3
do while(esize > min_efillsize
if(nump_ele < esize**3) then
esize = esize - 2
else
ele = cubic_cell*(esize -1)
do o = 1, size_ele(ie) - esize
do n = 1, size_ele(ie) - esize
do m = 1, size_ele(ie) - esize
do inod = 1, ng_node(lat_ele(ie)
vlat = ele(:,inod) + (/ m, n, o /)
call get_interp_pos(vlat(1), vlat(2), vlat(3), ie, rfill(:,:,inod))
!Now start trying to remesh the region, leaving the slip plane as a discontinuity
esize = size_ele(ie) - 2
nump_ele = size_ele(ie)**3
do while(esize > min_efillsize)
if(nump_ele < esize**3) then
esize = esize - 2
else
ele = cubic_cell*(esize -1)
do o = 1, size_ele(ie) - esize
do n = 1, size_ele(ie) - esize
do m = 1, size_ele(ie) - esize
do inod = 1, ng_node(lat_ele(ie))
vlat = ele(:,inod) + (/ m, n, o /)
call get_interp_pos(vlat(1), vlat(2), vlat(3), ie, rfill(:,:,inod))
end do
!Check to see if the plane intersects this element if not then add it
if((spos < maxval(rfill(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie))))).and. &
(spos > minval(rfill(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie)))))) then
nump_ele = nump_ele - esize**3
lat_points(m:m+esize, n:n+esize, o:o+esize) = .false.
call add_element(0, type_ele(ie), esize, lat_ele(ie), sbox_ele(ie), rfill)
end if
end do
!Check to see if the plane intersects this element if not then add it
if((spos < maxval(rfill(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie)),ie))).and. &
(spos > minval(rfill(sdim,1:basisnum(lat_ele(ie),1:ng_node(lat_ele(ie)),ie)))) then
nump_ele = nump_ele - esize**3
end if
end do
end do
end do
end if
end if
end do
end if
end if
end do