|
|
|
|
module opt_slip_plane
|
|
|
|
|
use parameters
|
|
|
|
|
use elements
|
|
|
|
|
use functions
|
|
|
|
|
use subroutines
|
|
|
|
|
|
|
|
|
|
implicit none
|
|
|
|
|
|
|
|
|
|
integer :: sdim
|
|
|
|
|
real(kind=dp) :: spos
|
|
|
|
|
logical :: efill
|
|
|
|
|
|
|
|
|
|
public
|
|
|
|
|
contains
|
|
|
|
|
|
|
|
|
|
subroutine run_slip_plane(arg_pos)
|
|
|
|
|
!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), ibasis
|
|
|
|
|
|
|
|
|
|
integer, allocatable :: slip_eles(:), temp_int(:)
|
|
|
|
|
real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), rfill(3,max_basisnum, max_ng_node), ratom(3,max_basisnum), &
|
|
|
|
|
maxp, minp
|
|
|
|
|
|
|
|
|
|
integer :: type_interp(max_basisnum*max_esize**3)
|
|
|
|
|
logical :: lat_points(max_esize,max_esize, max_esize)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
print *, '---------------------Option Slip_Plane----------------------'
|
|
|
|
|
|
|
|
|
|
!Initialize variables
|
|
|
|
|
efill = .false.
|
|
|
|
|
slip_enum = 0
|
|
|
|
|
old_atom_num = atom_num
|
|
|
|
|
|
|
|
|
|
!!Parse the argument
|
|
|
|
|
call parse(arg_pos)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
!If we are running the efill code then we have to initiate some variables
|
|
|
|
|
if(efill) then
|
|
|
|
|
new_ele_num = 0
|
|
|
|
|
end if
|
|
|
|
|
allocate(slip_eles(1024))
|
|
|
|
|
!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
|
|
|
|
|
slip_enum = slip_enum + 1
|
|
|
|
|
if (slip_enum > size(slip_eles)) then
|
|
|
|
|
allocate(temp_int(size(slip_eles)+1024))
|
|
|
|
|
temp_int(1:size(slip_eles)) = slip_eles
|
|
|
|
|
temp_int(size(slip_eles)+1:) = 0
|
|
|
|
|
call move_alloc(temp_int, slip_eles)
|
|
|
|
|
end if
|
|
|
|
|
slip_eles(slip_enum) = ie
|
|
|
|
|
|
|
|
|
|
!If we aren't efilling then just refine the element
|
|
|
|
|
if(.not.efill) then
|
|
|
|
|
call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp)
|
|
|
|
|
do ia = 1, basisnum(lat_ele(ie)) * size_ele(ie)**3
|
|
|
|
|
call apply_periodic(r_interp(:,ia))
|
|
|
|
|
call add_atom(0, type_interp(ia), sbox_ele(ie), r_interp(:,ia))
|
|
|
|
|
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.
|
|
|
|
|
|
|
|
|
|
!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
|
|
|
|
|
latloop:do m = 1, size_ele(ie) - esize
|
|
|
|
|
do inod = 1, ng_node(lat_ele(ie))
|
|
|
|
|
vlat = ele(:,inod) + (/ m, n, o /)
|
|
|
|
|
if (.not.lat_points(vlat(1), vlat(2),vlat(3))) cycle latloop
|
|
|
|
|
call get_interp_pos(vlat(1), vlat(2), vlat(3), ie, rfill(:,:,inod))
|
|
|
|
|
end do
|
|
|
|
|
|
|
|
|
|
!Check to make sure all lattice points exist for the current element
|
|
|
|
|
if(any(.not.lat_points(m:m+esize-1, n:n+esize-1, o:o+esize-1))) cycle latloop
|
|
|
|
|
!Check to see if the plane intersects this element if not then add it
|
|
|
|
|
maxp = maxval(rfill(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie))))
|
|
|
|
|
minp = minval(rfill(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie))))
|
|
|
|
|
if(.not.(spos < maxp).and.(spos > minp))then
|
|
|
|
|
nump_ele = nump_ele - esize**3
|
|
|
|
|
lat_points(m:m+esize-1, n:n+esize-1, o:o+esize-1) = .false.
|
|
|
|
|
call add_element(0, type_ele(ie), esize, lat_ele(ie), sbox_ele(ie), rfill)
|
|
|
|
|
new_ele_num = new_ele_num + 1
|
|
|
|
|
end if
|
|
|
|
|
end do latloop
|
|
|
|
|
end do
|
|
|
|
|
end do
|
|
|
|
|
end if
|
|
|
|
|
esize= esize-2
|
|
|
|
|
end do
|
|
|
|
|
! Now add the leftover lattice points as atoms
|
|
|
|
|
do o = 1, size_ele(ie)
|
|
|
|
|
do n = 1, size_ele(ie)
|
|
|
|
|
do m = 1, size_ele(ie)
|
|
|
|
|
if(lat_points(m,n,o)) then
|
|
|
|
|
call get_interp_pos(m,n,o, ie, ratom(:,:))
|
|
|
|
|
do ibasis = 1, basisnum(lat_ele(ie))
|
|
|
|
|
call apply_periodic(r_atom(:,ibasis))
|
|
|
|
|
call add_atom(0, basis_type(ibasis,lat_ele(ie)), sbox_ele(ie), ratom(:,ibasis))
|
|
|
|
|
end do
|
|
|
|
|
end if
|
|
|
|
|
end do
|
|
|
|
|
end do
|
|
|
|
|
end do
|
|
|
|
|
end if
|
|
|
|
|
end if
|
|
|
|
|
end do
|
|
|
|
|
|
|
|
|
|
!Once we finish adding elements delete the old ones
|
|
|
|
|
call delete_elements(slip_enum, slip_eles(1:slip_enum))
|
|
|
|
|
|
|
|
|
|
!Output data
|
|
|
|
|
if(.not.efill) then
|
|
|
|
|
print *, "We refine ", slip_enum, " elements into ", atom_num - old_atom_num , " atoms"
|
|
|
|
|
else
|
|
|
|
|
print *, "We refine ", slip_enum, " elements into ", atom_num - old_atom_num , " atoms and ", new_ele_num, " elements"
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
end subroutine run_slip_plane
|
|
|
|
|
|
|
|
|
|
subroutine parse(arg_pos)
|
|
|
|
|
!This subroutine parses the input arguments to the mode
|
|
|
|
|
integer, intent(inout) :: arg_pos
|
|
|
|
|
|
|
|
|
|
integer :: arglen
|
|
|
|
|
character(len = 100) :: textholder
|
|
|
|
|
|
|
|
|
|
!First read the dimension
|
|
|
|
|
arg_pos = arg_pos +1
|
|
|
|
|
call get_command_argument(arg_pos,textholder, arglen)
|
|
|
|
|
if(arglen == 0) stop "Incorrect slip_plane command. Please check documentation for correct format"
|
|
|
|
|
|
|
|
|
|
!Check to make sure that the dimension is correct
|
|
|
|
|
select case(trim(adjustl(textholder)))
|
|
|
|
|
case('x','X')
|
|
|
|
|
sdim = 1
|
|
|
|
|
case('y','Y')
|
|
|
|
|
sdim = 2
|
|
|
|
|
case('z','Z')
|
|
|
|
|
sdim = 3
|
|
|
|
|
case default
|
|
|
|
|
print *, "Error: dimension ", trim(adjustl(textholder)), " is not accepted. Please select from x, y, or z"
|
|
|
|
|
end select
|
|
|
|
|
|
|
|
|
|
!Now parse the position of the slip plane
|
|
|
|
|
arg_pos = arg_pos + 1
|
|
|
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
|
|
|
if(arglen == 0) stop "Incorrect slip_plane command. Please check documentation for correct format"
|
|
|
|
|
call parse_pos(sdim, textholder, spos)
|
|
|
|
|
|
|
|
|
|
!Now check to see if efill was passed
|
|
|
|
|
arg_pos = arg_pos + 1
|
|
|
|
|
if(.not.(arg_pos > command_argument_count())) then
|
|
|
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
|
|
|
if(arglen == 0) stop "Incorrect slip_plane command. Please check documentation for correct format"
|
|
|
|
|
if(trim(adjustl(textholder)) == "efill") then
|
|
|
|
|
arg_pos = arg_pos +1
|
|
|
|
|
efill = .true.
|
|
|
|
|
end if
|
|
|
|
|
end if
|
|
|
|
|
end subroutine parse
|
|
|
|
|
end module opt_slip_plane
|
