aselimov
/
CACmb
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Restructure group code, add group shell, and add refinefill option

Browse Source
development
Alex Selimov 4 years ago
parent f92e80f9db
commit 2c42520102
2 changed files with 231 additions and 49 deletions
Show Stats Download Patch File Download Diff File

3
src/elements.f90
Unescape View File

@ -779,8 +779,7 @@ do i = 1, atom_num
end do end do
end do end do
end subroutine get_interp_pos
end subroutine
subroutine alloc_dat_arrays(n,m) subroutine alloc_dat_arrays(n,m)
!This subroutine used to provide initial allocation for the atom and element data arrays !This subroutine used to provide initial allocation for the atom and element data arrays

277
src/opt_group.f90
Unescape View File

@ -10,8 +10,8 @@ module opt_group
integer :: group_ele_num, group_atom_num, remesh_size,normal, dim1, dim2, random_num, group_type, notsize integer :: group_ele_num, group_atom_num, remesh_size,normal, dim1, dim2, random_num, group_type, notsize
character(len=15) :: type, shape !Type indicates what element type is selected and shape is the group shape character(len=15) :: type, shape !Type indicates what element type is selected and shape is the group shape
real(kind=dp) :: block_bd(6), centroid(3), vertices(3,3),disp_vec(3), radius, bwidth real(kind=dp) :: block_bd(6), centroid(3), vertices(3,3),disp_vec(3), radius, bwidth, shell_thickness
logical :: displace, delete, max_remesh, refine, group_nodes, flip logical :: displace, delete, max_remesh, refine, group_nodes, flip, efill, refinefill
integer, allocatable :: element_index(:), atom_index(:) integer, allocatable :: element_index(:), atom_index(:)
@ -22,7 +22,9 @@ module opt_group
!Main calling function for the group option !Main calling function for the group option
integer, intent(inout) :: arg_pos integer, intent(inout) :: arg_pos
print *, '-----------------------Option Group-------------------------' print *, '------------------------------------------------------------'
print *, 'Option Group'
print *, '------------------------------------------------------------'
group_ele_num = 0 group_ele_num = 0
group_atom_num = 0 group_atom_num = 0
@ -48,6 +50,11 @@ module opt_group
call refine_group call refine_group
end if end if
if(refinefill) then
call get_group
call refinefill_group
end if
if(displace)then if(displace)then
call get_group call get_group
call displace_group call displace_group
@ -335,6 +342,28 @@ module opt_group
if (arglen==0) STOP "Missing sphere radius in group command" if (arglen==0) STOP "Missing sphere radius in group command"
read(textholder, *) radius read(textholder, *) radius
case('shell')
!First extract the shell centroid
do i = 1, 3
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) STOP "Missing shell centroid in group command"
call parse_pos(i, textholder, centroid(i))
end do
!Now get the radius
arg_pos=arg_pos+1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) STOP "Missing shell radius in group command"
read(textholder, *) radius
!Now get the shell thickness
arg_pos=arg_pos+1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) STOP "Missing shell thickness in group command"
read(textholder, *) shell_thickness
case('all') case('all')
!Do nothing if the shape is all !Do nothing if the shape is all
continue continue
@ -364,6 +393,8 @@ module opt_group
end do end do
case('refine') case('refine')
refine=.true. refine=.true.
case('refinefill')
refinefill=.true.
case('remesh') case('remesh')
arg_pos = arg_pos + 1 arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen) call get_command_argument(arg_pos, textholder, arglen)
@ -382,6 +413,8 @@ module opt_group
read(textholder, *) random_num read(textholder, *) random_num
case('flip') case('flip')
flip=.true. flip=.true.
case('efill')
efill=.true.
case('type') case('type')
arg_pos = arg_pos + 1 arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen) call get_command_argument(arg_pos, textholder, arglen)
@ -403,7 +436,7 @@ module opt_group
subroutine get_group subroutine get_group
!This subroutine finds all elements and/or atoms within the group boundaries !This subroutine finds all elements and/or atoms within the group boundaries
!specified by the user. !specified by the user.
integer :: i, j, inod, ibasis, temp integer :: i, j, inod, ibasis, temp, node_in_out(max_ng_node)
integer, allocatable :: resize_array(:) integer, allocatable :: resize_array(:)
real(kind=dp) :: r_center(3), rand real(kind=dp) :: r_center(3), rand
@ -433,48 +466,19 @@ module opt_group
!Check the type to see whether we need to find the elements within the group !Check the type to see whether we need to find the elements within the group
select case(trim(adjustl(type))) select case(trim(adjustl(type)))
case('elements', 'both') case('elements', 'both')
if(.not.(group_nodes)) then do i = 1, ele_num
do i = 1, ele_num if(in_group_ele(size_ele(i), lat_ele(i), r_node(:,:,:,i))) then
r_center(:) = 0.0_dp group_ele_num = group_ele_num + 1
do inod = 1, ng_node(lat_ele(i)) if(group_ele_num > size(element_index)) then
do ibasis = 1, basisnum(lat_ele(i)) allocate(resize_array(size(element_index) + 1024))
r_center = r_center + r_node(:,ibasis,inod,i)/(basisnum(lat_ele(i))*ng_node(lat_ele(i))) resize_array(1:group_ele_num-1) = element_index
end do resize_array(group_ele_num:) = 0
end do call move_alloc(resize_array, element_index)
end if
if ((in_group(r_center).neqv.flip).and.(size_ele(i)/= notsize)) then
group_ele_num = group_ele_num + 1
if(group_ele_num > size(element_index)) then
allocate(resize_array(size(element_index) + 1024))
resize_array(1:group_ele_num-1) = element_index
resize_array(group_ele_num:) = 0
call move_alloc(resize_array, element_index)
end if
element_index(group_ele_num) = i
end if
end do
else if(group_nodes) then element_index(group_ele_num) = i
eleloop:do i = 1, ele_num end if
r_center(:) = 0.0_dp end do
do inod = 1, ng_node(lat_ele(i))
do ibasis = 1, basisnum(lat_ele(i))
if ((in_group(r_node(:,ibasis,inod,i)).neqv.flip).and.(size_ele(i)/=notsize)) then
group_ele_num = group_ele_num + 1
if(group_ele_num > size(element_index)) then
allocate(resize_array(size(element_index) + 1024))
resize_array(1:group_ele_num-1) = element_index
resize_array(group_ele_num:) = 0
call move_alloc(resize_array, element_index)
end if
element_index(group_ele_num) = i
cycle eleloop
end if
end do
end do
end do eleloop
end if
if(random_num > 0) then if(random_num > 0) then
!If we have the random option enabled then we select random_num number of elements from the group and overwrite !If we have the random option enabled then we select random_num number of elements from the group and overwrite
@ -569,7 +573,7 @@ module opt_group
end subroutine displace_group end subroutine displace_group
subroutine refine_group subroutine refine_group
!This command is used to remesh the group to a desired element size !This command is used to refine the group to full atomistics
integer :: i, j, ie, type_interp(max_basisnum*max_esize**3), add_atom_num, orig_atom_num integer :: i, j, ie, type_interp(max_basisnum*max_esize**3), add_atom_num, orig_atom_num
real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3) real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3)
@ -596,7 +600,104 @@ module opt_group
print *, group_ele_num, " elements of group are refined to ", atom_num -orig_atom_num, " atoms." print *, group_ele_num, " elements of group are refined to ", atom_num -orig_atom_num, " atoms."
end if end if
end subroutine end subroutine refine_group
subroutine refinefill_group
!This command is used to refine the group to full atomistics
integer :: i, j, ie, type_interp(max_basisnum*max_esize**3), add_atom_num, orig_atom_num, m, n, o, esize, &
ele(3,8), new_ele_num, ibasis, inod, vlat(3), nump_ele
real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), rfill(3,max_basisnum,max_ng_node), ravg(3), ratom(3,max_basisnum)
logical :: lat_points(max_esize, max_esize, max_esize)
!Refining to atoms
if(group_ele_num > 0) then
orig_atom_num = atom_num
new_ele_num = 0
!Estimate number of atoms we are adding, this doesn't have to be exact
add_atom_num = group_ele_num*basisnum(lat_ele(element_index(1)))*size_ele(element_index(1))**3
call grow_ele_arrays(0,add_atom_num)
do i = 1, group_ele_num
ie = element_index(i)
!Find all possible elements that we can make while making sure they aren't in the group
lat_points(1:size_ele(ie),1:size_ele(ie),1:size_ele(ie)) = .true.
!Now calculate the number of elemets which are available for remeshing
call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp)
nump_ele = size_ele(ie)**3
do o =1, size_ele(ie)
do n = 1, size_ele(ie)
do m =1, size_ele(ie)
call get_interp_pos(m,n,o,i,rfill(:,:,1))
ravg(:) = 0
do ibasis = 1, basisnum(lat_ele(ie))
ravg = ravg + rfill(:,ibasis, 1)/basisnum(lat_ele(ie))
end do
if( in_group(ravg)) then
nump_ele = nump_ele - 1
end if
end do
end do
end do
!Now start the remeshing loop for the element
esize = size_ele(ie) - 2
do while(esize > min_efillsize)
if(nump_ele < min_efillsize**3) then
exit
else 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
if (.not.in_group_ele(esize, lat_ele(ie), rfill)) 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
esize=esize-2
end if
end do
!Now add the leftover lattice points as atoms
do o = 1, size_ele(i)
do n = 1, size_ele(i)
do m = 1, size_ele(i)
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 do
!Once all atoms are added we delete all of the elements
call delete_elements(group_ele_num, element_index)
print *, group_ele_num, " elements of group are refined to ", atom_num -orig_atom_num, " atoms and ", new_ele_num, &
" elements."
end if
end subroutine refinefill_group
subroutine remesh_group subroutine remesh_group
!This command is used to remesh the group to a desired element size !This command is used to remesh the group to a desired element size
@ -900,6 +1001,11 @@ module opt_group
end subroutine change_group_type end subroutine change_group_type
subroutine split_group_elements
!
end subroutine split_group_elements
function in_group(r) function in_group(r)
!This subroutine determines if a point is within the group boundaries !This subroutine determines if a point is within the group boundaries
real(kind=dp), intent(in) :: r(3) real(kind=dp), intent(in) :: r(3)
@ -933,8 +1039,85 @@ module opt_group
else else
in_group = .false. in_group = .false.
end if end if
case('shell')
rnorm = norm2(r(:) - centroid(:))
if ((rnorm >= radius).and.(rnorm<=(radius + shell_thickness))) then
in_group = .true.
else
in_group = .false.
end if
case('all') case('all')
in_group = .true. in_group = .true.
end select end select
end function in_group end function in_group
function in_group_ele(esize, elat, rn)
!This figures out if the elements are in the group boundaries
real(kind=dp), intent(in) :: rn(3,max_basisnum, max_ng_node)
integer, intent(in) :: esize, elat
logical :: in_group_ele
integer :: i, inod, ibasis, node_in_out(max_ng_node)
real(kind=dp) :: r_center(3)
in_group_ele=.false.
if(trim(adjustl(shape)) == 'shell') then
node_in_out(:) = -1
!First calculate whether each element node is within the shell region, inside the shell sphere, or outside the
!shell region
nodeloop:do inod = 1, ng_node(elat)
r_center(:)=0.0_dp
do ibasis = 1, basisnum(elat)
r_center(:)= r_center(:) + rn(:,ibasis,inod)/basisnum(elat)
end do
if((in_group(rn(:, ibasis, inod)).neqv.flip).and.(size_ele(i)/=notsize)) then
node_in_out(inod) = 2
exit nodeloop
end if
shape ='sphere'
if((in_group(rn(:, ibasis, inod)).neqv.flip).and.(esize/=notsize)) then
node_in_out(inod) = 1
else
node_in_out(inod) = 0
end if
shape='shell'
end do nodeloop
!If any nodes are within the shell region, or if the shell region interescts an element then add it to the group
if(any(node_in_out == 2).or.(any(node_in_out==1).and.(any(node_in_out==0)))) then
in_group_ele=.true.
return
end if
else if(.not.(group_nodes)) then
r_center(:) = 0.0_dp
do inod = 1, ng_node(elat)
do ibasis = 1, basisnum(elat)
r_center = r_center + rn(:,ibasis,inod)/(basisnum(elat)*ng_node(elat))
end do
end do
if ((in_group(r_center).neqv.flip).and.(esize/= notsize)) then
in_group_ele=.true.
return
end if
else if(group_nodes) then
r_center(:) = 0.0_dp
do inod = 1, ng_node(elat)
do ibasis = 1, basisnum(elat)
if ((in_group(rn(:,ibasis,inod)).neqv.flip).and.(esize/=notsize)) then
in_group_ele=.true.
return
end if
end do
end do
end if
end function in_group_ele
end module opt_group end module opt_group

Write Preview
Loading…
Cancel
Save