You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
1130 lines
47 KiB
1130 lines
47 KiB
module opt_group
|
|
|
|
!This module contains all code associated with dislocations
|
|
|
|
use parameters
|
|
use elements
|
|
use subroutines
|
|
use box
|
|
implicit none
|
|
|
|
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
|
|
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, efill, refinefill
|
|
|
|
integer, allocatable :: element_index(:), atom_index(:)
|
|
|
|
public
|
|
contains
|
|
|
|
subroutine group(arg_pos)
|
|
!Main calling function for the group option
|
|
integer, intent(inout) :: arg_pos
|
|
|
|
print *, '------------------------------------------------------------'
|
|
print *, 'Option Group'
|
|
print *, '------------------------------------------------------------'
|
|
|
|
group_ele_num = 0
|
|
group_atom_num = 0
|
|
remesh_size=0
|
|
random_num=0
|
|
group_type=0
|
|
notsize=0
|
|
displace=.false.
|
|
delete=.false.
|
|
max_remesh=.false.
|
|
refine = .false.
|
|
flip = .false.
|
|
|
|
if(allocated(element_index)) deallocate(element_index)
|
|
if(allocated(atom_index)) deallocate(atom_index)
|
|
|
|
call parse_group(arg_pos)
|
|
|
|
|
|
!Now call the transformation functions for the group
|
|
if(refine) then
|
|
call get_group
|
|
call refine_group
|
|
end if
|
|
|
|
if(refinefill) then
|
|
call get_group
|
|
call refinefill_group
|
|
end if
|
|
|
|
if(displace)then
|
|
call get_group
|
|
call displace_group
|
|
end if
|
|
|
|
if(delete)then
|
|
call get_group
|
|
call delete_group
|
|
end if
|
|
|
|
if(remesh_size > 0)then
|
|
call get_group
|
|
call remesh_group
|
|
end if
|
|
|
|
if(group_type > 0) then
|
|
call get_group
|
|
call change_group_type
|
|
end if
|
|
|
|
end subroutine group
|
|
|
|
subroutine parse_group(arg_pos)
|
|
!Parse the group command
|
|
integer, intent(inout) :: arg_pos
|
|
|
|
integer :: i, j, arglen, in_num
|
|
character(len=100) :: textholder, type_spec
|
|
real(kind=dp) H
|
|
|
|
!Parse type and shape command
|
|
arg_pos = arg_pos + 1
|
|
call get_command_argument(arg_pos, type, arglen)
|
|
if (arglen==0) STOP "Missing select_type in group command"
|
|
select case(trim(adjustl(type)))
|
|
case('atoms', 'elements', 'both')
|
|
continue
|
|
case default
|
|
print *, "Select_type ", trim(adjustl(type)), " is not an accept group selection criteria. ", &
|
|
"Please select from atoms, elements, or both."
|
|
end select
|
|
|
|
arg_pos = arg_pos + 1
|
|
call get_command_argument(arg_pos, shape, arglen)
|
|
if (arglen==0) STOP "Missing group_shape in group command"
|
|
|
|
!Now parse the arguments required by the user selected shape
|
|
select case(trim(adjustl(shape)))
|
|
case('block')
|
|
do i= 1, 6
|
|
arg_pos = arg_pos + 1
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
if (arglen==0) STOP "Missing block boundary in group command"
|
|
call parse_pos(int((i+1)/2), textholder, block_bd(i))
|
|
end do
|
|
|
|
case('wedge')
|
|
arg_pos = arg_pos + 1
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
if (arglen==0) STOP "Missing normal dim in group wedge command"
|
|
read(textholder,*) dim1
|
|
|
|
arg_pos = arg_pos + 1
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
if (arglen==0) STOP "Missing normal dim in group wedge command"
|
|
read(textholder,*) dim2
|
|
|
|
do i = 1, 3
|
|
arg_pos = arg_pos + 1
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
if (arglen==0) STOP "Missing centroid in group wedge command"
|
|
call parse_pos(i, textholder, centroid(i))
|
|
end do
|
|
|
|
arg_pos = arg_pos + 1
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
if (arglen==0) STOP "Missing base width in group wedge command"
|
|
read(textholder,*) bwidth
|
|
|
|
!Calculate the vertex positions
|
|
vertices(:,1) = centroid
|
|
vertices(dim2,1) = 0.0_dp
|
|
do i = 1, 3
|
|
if (i == dim1) then
|
|
if (bwidth > 0) then
|
|
vertices(i,2) = box_bd(2*i)
|
|
vertices(i,3) = box_bd(2*i)
|
|
else if (bwidth < 0) then
|
|
vertices(i,2) = box_bd(2*i-1)
|
|
vertices(i,3) = box_bd(2*i-1)
|
|
else
|
|
print *, "bwidth cannot be 0 in wedge shaped group"
|
|
stop 3
|
|
end if
|
|
else if (i == dim2) then
|
|
vertices(i,2) = 0.0_dp
|
|
vertices(i,3) = 0.0_dp
|
|
else
|
|
vertices(i,2) = centroid(i) + bwidth
|
|
vertices(i,3) = centroid(i) - bwidth
|
|
end if
|
|
end do
|
|
|
|
|
|
case('notch')
|
|
arg_pos = arg_pos + 1
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
if (arglen==0) STOP "Missing normal dim in group notch command"
|
|
read(textholder,*) dim1
|
|
|
|
arg_pos = arg_pos + 1
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
if (arglen==0) STOP "Missing normal dim in group notch command"
|
|
read(textholder,*) dim2
|
|
|
|
do i = 1, 3
|
|
arg_pos = arg_pos + 1
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
if (arglen==0) STOP "Missing centroid in group notch command"
|
|
call parse_pos(i, textholder, centroid(i))
|
|
end do
|
|
|
|
arg_pos = arg_pos + 1
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
if (arglen==0) STOP "Missing base width in group notch command"
|
|
read(textholder,*) bwidth
|
|
|
|
arg_pos = arg_pos + 1
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
if (arglen==0) STOP "Missing tip radius in group notch command"
|
|
read(textholder,*) radius
|
|
|
|
!Calculate the vertex positions
|
|
vertices(:,1) = centroid
|
|
vertices(dim2,1) = 0.0_dp
|
|
do i = 1, 3
|
|
if (i == dim1) then
|
|
if (bwidth > 0) then
|
|
vertices(i,2) = box_bd(2*i)
|
|
vertices(i,3) = box_bd(2*i)
|
|
H= (box_bd(2*i) - centroid(i)) !Calculate the height of the wedge
|
|
else if (bwidth < 0) then
|
|
vertices(i,2) = box_bd(2*i-1)
|
|
vertices(i,3) = box_bd(2*i-1)
|
|
H = (centroid(i) - box_bd(2*i-1))
|
|
else
|
|
print *, "bwidth cannot be 0 in wedge shaped group"
|
|
stop 3
|
|
end if
|
|
else if (i == dim2) then
|
|
vertices(i,2) = 0.0_dp
|
|
vertices(i,3) = 0.0_dp
|
|
else
|
|
vertices(i,2) = centroid(i) + 0.5_dp*bwidth
|
|
vertices(i,3) = centroid(i) - 0.5_dp*bwidth
|
|
end if
|
|
end do
|
|
|
|
!Now update the centroid so that the desire tip diameter matches the wedge with
|
|
if (bwidth > 0) then
|
|
centroid(dim1) = centroid(dim1) + 2*radius*(H)/bwidth
|
|
end if
|
|
!Read the ID type shape for group
|
|
case('id')
|
|
arg_pos = arg_pos + 1
|
|
!For this type we have to call different options if type is atoms, elements, or both. Both is the most complex.
|
|
select case(trim(adjustl(type)))
|
|
case('atoms')
|
|
!Read number of ids
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
if(arglen == 0) then
|
|
print *, "Missing number of input atoms for group id"
|
|
end if
|
|
read(textholder,*) in_num
|
|
|
|
!allocate arrays
|
|
allocate(atom_index(in_num))
|
|
!Read ids
|
|
do i = 1, in_num
|
|
arg_pos = arg_pos + 1
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
if(arglen == 0) then
|
|
print *, "Missing atom id in group atom id"
|
|
stop 3
|
|
end if
|
|
read(textholder,*) atom_index(i)
|
|
end do
|
|
group_atom_num = in_num
|
|
|
|
case('elements')
|
|
!Read number of ids
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
if(arglen == 0) then
|
|
print *, "Missing number of input elements for group id"
|
|
end if
|
|
read(textholder, *) in_num
|
|
|
|
!allocate arrays
|
|
allocate(element_index(in_num))
|
|
!Read ids
|
|
do i = 1, in_num
|
|
arg_pos = arg_pos + 1
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
if(arglen == 0) then
|
|
print *, "Missing element id in group element id"
|
|
stop 3
|
|
end if
|
|
read(textholder,*) element_index(i)
|
|
end do
|
|
group_ele_num = in_num
|
|
|
|
case('both')
|
|
!We repeat this code twice, once for the atoms and once for the elements
|
|
allocate(element_index(1024),atom_index(1024))
|
|
do i = 1, 2
|
|
!Read the first id type (either atoms or elements)
|
|
call get_command_argument(arg_pos, type_spec, arglen)
|
|
if (arglen == 0) then
|
|
print *, "Missing type specifier in group id command"
|
|
stop 3
|
|
end if
|
|
|
|
!Now read the first in_num
|
|
arg_pos = arg_pos + 1
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
if(arglen == 0) then
|
|
print *, "Missing input number in group_id"
|
|
stop 3
|
|
end if
|
|
read(textholder, *) in_num
|
|
|
|
select case(trim(adjustl(type_spec)))
|
|
case('atoms','atom')
|
|
if (group_atom_num > 0) then
|
|
print *, "Atoms specifier used more than once in group id command with type both, either use type ", &
|
|
"atoms or include elements identifier"
|
|
stop 3
|
|
|
|
do j = 1, in_num
|
|
arg_pos = arg_pos + 1
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
if (arglen == 0) then
|
|
print *, "Missing atom id in group atom id"
|
|
stop 3
|
|
end if
|
|
read(textholder, *) atom_index(j)
|
|
end do
|
|
group_atom_num = in_num
|
|
end if
|
|
|
|
case('elements','element')
|
|
if (group_ele_num > 0) then
|
|
print *, "Elements specifier used more than once in group id command with type both, either use type ",&
|
|
"elements or include atoms identifier"
|
|
stop 3
|
|
|
|
do j = 1, in_num
|
|
arg_pos = arg_pos + 1
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
if (arglen == 0) then
|
|
print *, "Missing element id in group element id"
|
|
stop 3
|
|
end if
|
|
read(textholder, *) element_index(j)
|
|
end do
|
|
group_ele_num = in_num
|
|
end if
|
|
|
|
end select
|
|
if(i ==1) arg_pos = arg_pos + 1
|
|
end do
|
|
end select
|
|
|
|
case('sphere')
|
|
!First extract the sphere centroid
|
|
do i = 1, 3
|
|
arg_pos = arg_pos + 1
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
if (arglen==0) STOP "Missing sphere centroid in group command"
|
|
call parse_pos(i, textholder, centroid(i))
|
|
end do
|
|
!Now get the tip radius
|
|
arg_pos=arg_pos+1
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
if (arglen==0) STOP "Missing sphere radius in group command"
|
|
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')
|
|
!Do nothing if the shape is all
|
|
continue
|
|
|
|
case default
|
|
print *, "Group shape ", trim(adjustl(shape)), " is not currently accepted. Please check documentation ", &
|
|
"for accepted group shapes."
|
|
end select
|
|
|
|
!Now parse the additional options which may be present
|
|
do while(.true.)
|
|
if(arg_pos > command_argument_count()) exit
|
|
!Pull out the next argument which should either be a keyword or an option
|
|
arg_pos=arg_pos+1
|
|
call get_command_argument(arg_pos, textholder)
|
|
textholder=adjustl(textholder)
|
|
|
|
!Choose what to based on what the option string is
|
|
select case(trim(textholder))
|
|
case('displace')
|
|
displace = .true.
|
|
do i = 1,3
|
|
arg_pos = arg_pos + 1
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
if (arglen==0) stop "Missing vector component for shift command"
|
|
read(textholder, *) disp_vec(i)
|
|
end do
|
|
case('refine')
|
|
refine=.true.
|
|
case('refinefill')
|
|
refinefill=.true.
|
|
case('remesh')
|
|
arg_pos = arg_pos + 1
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
if (arglen==0) stop "Missing remesh element size in group command"
|
|
read(textholder, *) remesh_size
|
|
case('max')
|
|
max_remesh =.true.
|
|
case('delete')
|
|
delete=.true.
|
|
case('nodes')
|
|
group_nodes=.true.
|
|
case('random')
|
|
arg_pos = arg_pos + 1
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
if (arglen==0) stop "Missing number of random atoms in group command"
|
|
read(textholder, *) random_num
|
|
case('flip')
|
|
flip=.true.
|
|
case('efill')
|
|
efill=.true.
|
|
case('type')
|
|
arg_pos = arg_pos + 1
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
if (arglen==0) stop "Missing atom type for group"
|
|
call add_atom_type(textholder, group_type)
|
|
case('notsize')
|
|
arg_pos = arg_pos + 1
|
|
call get_command_argument(arg_pos, textholder, arglen)
|
|
if(arglen ==0) stop "Missing notsize size"
|
|
read(textholder, *) notsize
|
|
print *, "Ignoring elements with size ", notsize
|
|
case default
|
|
!If it isn't an available option to opt_disl then we just exit
|
|
exit
|
|
end select
|
|
end do
|
|
end subroutine parse_group
|
|
|
|
subroutine get_group
|
|
!This subroutine finds all elements and/or atoms within the group boundaries
|
|
!specified by the user.
|
|
integer :: i, j, inod, ibasis, temp, node_in_out(max_ng_node)
|
|
integer, allocatable :: resize_array(:)
|
|
real(kind=dp) :: r_center(3), rand
|
|
|
|
select case(trim(adjustl(shape)))
|
|
case('block')
|
|
print *, "Group has block shape with boundaries: ", block_bd
|
|
case ('wedge')
|
|
print *, "Group has wedge shape with dim1", dim1, "and dim2", dim2, "and vertices ", vertices
|
|
case ('notch')
|
|
print *, "Group has notch shape with dim1", dim1, "and dim2", dim2, " tip radius ", radius, "and vertices ", &
|
|
vertices
|
|
case('id')
|
|
print *, 'Group contains ', group_ele_num, " elements and ", group_atom_num, " atoms."
|
|
return
|
|
case('sphere')
|
|
print *, "Group has sphere shape with centroid ", centroid(:), " and radius ", radius
|
|
case('all')
|
|
print *, "Group has all of type ", type
|
|
end select
|
|
|
|
!Reset group if needed
|
|
if(allocated(element_index)) deallocate(element_index,atom_index)
|
|
|
|
group_ele_num = 0
|
|
group_atom_num = 0
|
|
allocate(element_index(1024), atom_index(1024))
|
|
!Check the type to see whether we need to find the elements within the group
|
|
select case(trim(adjustl(type)))
|
|
case('elements', 'both')
|
|
do i = 1, ele_num
|
|
if(in_group_ele(size_ele(i), lat_ele(i), r_node(:,:,:,i))) 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
|
|
|
|
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
|
|
!the group with those elements
|
|
do i = 1, random_num
|
|
call random_number(rand)
|
|
j = i + floor((group_ele_num+1-i)*rand)
|
|
temp = element_index(j)
|
|
element_index(j) = element_index(i)
|
|
element_index(i) = temp
|
|
end do
|
|
|
|
group_ele_num = random_num
|
|
end if
|
|
|
|
end select
|
|
!Check the type to see if we need to find the atoms within the group
|
|
select case(trim(adjustl(type)))
|
|
case('atoms','both')
|
|
do i = 1, atom_num
|
|
if(in_group(r_atom(:,i)).neqv.flip) then
|
|
group_atom_num = group_atom_num + 1
|
|
if (group_atom_num > size(atom_index)) then
|
|
allocate(resize_array(size(atom_index) + 1024))
|
|
resize_array(1:group_atom_num -1) = atom_index
|
|
resize_array(group_atom_num:) = 0
|
|
call move_alloc(resize_array, atom_index)
|
|
end if
|
|
|
|
atom_index(group_atom_num) = i
|
|
end if
|
|
end do
|
|
if(random_num > 0) then
|
|
!If we have the random option enabled then we select random_num number of atom from the group and overwrite
|
|
!the group with those atom
|
|
do i = 1, random_num
|
|
call random_number(rand)
|
|
j = i + floor((group_atom_num+1-i)*rand)
|
|
temp = atom_index(j)
|
|
atom_index(j) = atom_index(i)
|
|
atom_index(i) = temp
|
|
end do
|
|
|
|
group_atom_num = random_num
|
|
end if
|
|
end select
|
|
|
|
j = 0
|
|
do i = 1, group_atom_num
|
|
if (atom_index(i) == 23318348) then
|
|
j = j + 1
|
|
end if
|
|
end do
|
|
|
|
if (j > 1) stop "Code broken"
|
|
|
|
print *, 'Group contains ', group_ele_num, " elements and ", group_atom_num, " atoms."
|
|
end subroutine get_group
|
|
|
|
subroutine displace_group
|
|
!This subroutine applies a displacement to elements/atoms in the groups
|
|
|
|
integer :: i, inod, ibasis
|
|
|
|
print *, "Elements/atoms in group displaced by ", disp_vec
|
|
|
|
!Displace atoms
|
|
do i = 1, group_atom_num
|
|
r_atom(:,atom_index(i)) = r_atom(:,atom_index(i)) + disp_vec
|
|
end do
|
|
|
|
!Displace elements
|
|
do i = 1, group_ele_num
|
|
do inod = 1, ng_node(lat_ele(element_index(i)))
|
|
do ibasis = 1, basisnum(lat_ele(element_index(i)))
|
|
r_node(:,ibasis,inod,element_index(i)) = r_node(:,ibasis,inod,element_index(i)) + disp_vec
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
!If we don't include the wrap command then we have to increase the size of the box
|
|
if (.not.(wrap_flag)) then
|
|
do i = 1,3
|
|
if (disp_vec(i) < -lim_zero) then
|
|
box_bd(2*i-1) = box_bd(2*i-1) - disp_vec(i)
|
|
else if (disp_vec(i) > lim_zero) then
|
|
box_bd(2*i) = box_bd(2*i) + disp_vec(i)
|
|
end if
|
|
end do
|
|
end if
|
|
|
|
end subroutine displace_group
|
|
|
|
subroutine refine_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
|
|
real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3)
|
|
|
|
!Refining to atoms
|
|
if(group_ele_num > 0) then
|
|
orig_atom_num = atom_num
|
|
!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)
|
|
!Get the interpolated atom positions
|
|
call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp)
|
|
!Loop over all interpolated atoms and add them to the system, we apply periodic boundaries
|
|
!here as well to make sure they are in the box
|
|
do j = 1, basisnum(lat_ele(ie))*size_ele(ie)**3
|
|
call apply_periodic(r_interp(:,j))
|
|
call add_atom(0,type_interp(j), sbox_ele(ie), r_interp(:,j))
|
|
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."
|
|
end if
|
|
|
|
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, added_points
|
|
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
|
|
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
|
|
added_points=0
|
|
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
|
|
added_points = added_points + esize**3
|
|
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(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(ratom(:,ibasis))
|
|
call add_atom(0, basis_type(ibasis,lat_ele(ie)), sbox_ele(ie), ratom(:,ibasis))
|
|
added_points=added_points + 1
|
|
end do
|
|
end if
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
if (added_points /= (size_ele(ie)**3)) then
|
|
|
|
print *, "Element ", ie, " is refined incorrectly in refinefill"
|
|
end if
|
|
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
|
|
!This command is used to remesh the group to a desired element size
|
|
|
|
integer :: i, j, k, ix, iy, iz, inod, ibasis, ie, type_interp(max_basisnum*max_esize**3), add_atom_num, orig_atom_num, &
|
|
current_esize, dof, max_lat(3), r_lat(3), ele(3,8), vlat(3), bd_in_lat(6), bd_in_array(3), new_ele, new_atom, &
|
|
max_loops, working_esize, group_lat_num, lat_list(10), group_sbox_num, sbox_list(100), is, ilat, tot_dof
|
|
|
|
real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), ori_inv(3,3), r(3), &
|
|
r_new_node(3,max_basisnum, max_ng_node), orient(3,3), group_in_lat(3,8), group_bd(6)
|
|
logical, allocatable :: lat_points(:,:,:), new_sbox, new_lat
|
|
character(len=100) :: remesh_ele_type
|
|
|
|
|
|
!Right now we just hardcode only remeshing to elements
|
|
remesh_ele_type = 'fcc'
|
|
|
|
! Determine which sub_boxes and lattices types are within in the group
|
|
group_sbox_num = 0
|
|
group_lat_num = 0
|
|
do i = 1, group_atom_num
|
|
new_sbox=.true.
|
|
new_lat=.true.
|
|
do j = 1, group_sbox_num
|
|
if (sbox_list(j) == sbox_atom(atom_index(i))) then
|
|
new_sbox=.false.
|
|
exit
|
|
end if
|
|
end do
|
|
|
|
if(new_sbox) then
|
|
group_sbox_num = group_sbox_num + 1
|
|
sbox_list(group_sbox_num) = sbox_atom(atom_index(i))
|
|
end if
|
|
|
|
do j = 1, group_lat_num
|
|
if (basis_type(1,lat_list(j)) == type_atom(atom_index(i))) then
|
|
new_lat = .false.
|
|
exit
|
|
end if
|
|
end do
|
|
|
|
if (new_lat) then
|
|
group_lat_num = group_lat_num + 1
|
|
do k = 1, lattice_types
|
|
if (basis_type(1,k) == type_atom(atom_index(i))) lat_list(group_lat_num) = k
|
|
end do
|
|
end if
|
|
|
|
end do
|
|
|
|
do i = 1, group_ele_num
|
|
new_sbox=.true.
|
|
new_lat = .true.
|
|
do j = 1, group_sbox_num
|
|
if (sbox_list(j) == sbox_ele(element_index(i))) then
|
|
new_sbox = .false.
|
|
exit
|
|
end if
|
|
end do
|
|
|
|
if (new_sbox) then
|
|
group_sbox_num = group_sbox_num + 1
|
|
sbox_list(group_sbox_num) = sbox_ele(element_index(i))
|
|
end if
|
|
|
|
do j = 1, group_lat_num
|
|
if (lat_list(group_lat_num) == lat_ele(element_index(i))) then
|
|
new_lat=.false.
|
|
exit
|
|
end if
|
|
end do
|
|
|
|
if (new_lat) then
|
|
group_lat_num = group_lat_num + 1
|
|
lat_list(group_lat_num) = lat_ele(element_index(i))
|
|
end if
|
|
end do
|
|
|
|
new_atom = 0
|
|
new_ele=0
|
|
tot_dof=0
|
|
do is = 1, group_sbox_num
|
|
|
|
orient = sub_box_ori(:, :, sbox_list(is))
|
|
call matrix_inverse(orient,3,ori_inv)
|
|
|
|
do ilat = 1, group_lat_num
|
|
|
|
!First calculate max position in lattice space to be able to allocate lat_points array, also sum the total number
|
|
!of degrees of freedom which are added
|
|
dof = 0
|
|
do j = 1, 3
|
|
group_bd(2*j) = -huge(1.0_dp)
|
|
group_bd(2*j-1) = huge(1.0_dp)
|
|
end do
|
|
do i = 1, group_atom_num
|
|
if ((type_atom(atom_index(i)) == basis_type(1,ilat)).and.(sbox_atom(atom_index(i)) == is)) then
|
|
do j =1 ,3
|
|
if (r_atom(j,atom_index(i)) > group_bd(2*j)) group_bd(2*j) = r_atom(j,atom_index(i))
|
|
if (r_atom(j,atom_index(i)) < group_bd(2*j-1)) group_bd(2*j-1) = r_atom(j,atom_index(i))
|
|
end do
|
|
dof = dof + 1
|
|
end if
|
|
end do
|
|
|
|
do i = 1, group_ele_num
|
|
if ((lat_ele(element_index(i)) == ilat).and.(sbox_ele(element_index(i)) == is)) then
|
|
do inod =1, ng_node(ilat)
|
|
do ibasis = 1, basisnum(ilat)
|
|
do j = 1, 3
|
|
r =r_node(j,ibasis,inod,element_index(i))
|
|
if (r(j) > group_bd(2*j)) group_bd(2*j) = r(j)
|
|
if (r(j) < group_bd(2*j-1)) group_bd(2*j-1) = r(j)
|
|
end do
|
|
end do
|
|
end do
|
|
dof = dof + size_ele(element_index(i))**3
|
|
end if
|
|
end do
|
|
|
|
!If for some reason there are no dof in this loop then cycle out
|
|
if(dof == 0) cycle
|
|
tot_dof = tot_dof+dof
|
|
|
|
group_in_lat = reshape((/ group_bd(1),group_bd(3),group_bd(5), &
|
|
group_bd(2),group_bd(3),group_bd(5), &
|
|
group_bd(2),group_bd(4),group_bd(5), &
|
|
group_bd(1),group_bd(4),group_bd(5), &
|
|
group_bd(1),group_bd(3),group_bd(6), &
|
|
group_bd(2),group_bd(3),group_bd(6), &
|
|
group_bd(2),group_bd(4),group_bd(6), &
|
|
group_bd(1),group_bd(4),group_bd(6) /), [3,8])
|
|
|
|
group_in_lat = matmul(fcc_inv, matmul(ori_inv, group_in_lat/lapa(ilat)))
|
|
do i = 1, 3
|
|
bd_in_lat(2*i-1) = nint(minval(group_in_lat(i,:)))
|
|
bd_in_lat(2*i) = nint(maxval(group_in_lat(i,:)))
|
|
end do
|
|
|
|
allocate(lat_points(bd_in_lat(2)-bd_in_lat(1)+10, bd_in_lat(4)-bd_in_lat(3)+10, bd_in_lat(6)-bd_in_lat(5)+10))
|
|
lat_points(:,:,:) = .false.
|
|
dof = 0
|
|
|
|
!Now place all group atoms and group interpolated atoms into lat_points
|
|
do i = 1, group_atom_num
|
|
if (.not.((sbox_atom(atom_index(i)) == is).and.(basis_type(1,ilat) == type_atom(atom_index(i))))) cycle
|
|
r = r_atom(:,atom_index(i))/lapa(ilat)
|
|
r = matmul(fcc_inv,matmul(ori_inv,r))
|
|
do j = 1, 3
|
|
r_lat(j) = nint(r(j))
|
|
end do
|
|
!Do a check to make sure the code is working and that lattice points aren't being written on top of each other.
|
|
!This is primarily a debugging statement
|
|
if(lat_points(r_lat(1)-bd_in_lat(1)+5,r_lat(2)-bd_in_lat(3)+5,r_lat(3)-bd_in_lat(5)+5)) then
|
|
stop "Multiple atoms share same position in lat point array, this shouldn't happen"
|
|
else
|
|
lat_points(r_lat(1)-bd_in_lat(1)+5, r_lat(2)-bd_in_lat(3)+5, r_lat(3)-bd_in_lat(5)+5) = .true.
|
|
dof = dof + 1
|
|
end if
|
|
end do
|
|
|
|
!Now place interpolated atoms within lat_points array
|
|
do i =1, group_ele_num
|
|
if (.not.((sbox_ele(element_index(i)) == is).and.( lat_ele(element_index(i)) == ilat))) cycle
|
|
ie = element_index(i)
|
|
call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp)
|
|
do j = 1, size_ele(ie)**3 * basisnum(lat_ele(ie))
|
|
r = r_interp(:,j)/lapa(ilat)
|
|
r = matmul(fcc_inv,matmul(ori_inv,r))
|
|
do k = 1, 3
|
|
r_lat(k) = nint(r(k))
|
|
end do
|
|
!Do a check to make sure the code is working and that lattice points aren't being written on top of each
|
|
!other. This is primarily a debugging statement
|
|
if(lat_points(r_lat(1)-bd_in_lat(1)+5,r_lat(2)-bd_in_lat(3)+5,r_lat(3)-bd_in_lat(5)+5)) then
|
|
stop "Multiple atoms/interpolated atoms share same position in lat point array, this shouldn't happen"
|
|
else
|
|
lat_points(r_lat(1)-bd_in_lat(1)+5, r_lat(2)-bd_in_lat(3)+5, r_lat(3)-bd_in_lat(5)+5) = .true.
|
|
dof = dof + 1
|
|
end if
|
|
end do
|
|
end do
|
|
|
|
!Now run remeshing algorithm, not the most optimized or efficient but gets the job done
|
|
!Figure out new looping boundaries
|
|
bd_in_array(1) = bd_in_lat(2) - bd_in_lat(1) + 10
|
|
bd_in_array(2) = bd_in_lat(4) - bd_in_lat(3) + 10
|
|
bd_in_array(3) = bd_in_lat(6) - bd_in_lat(5) + 10
|
|
|
|
|
|
if (max_remesh) then
|
|
max_loops = (remesh_size-3)/2
|
|
else
|
|
max_loops = 1
|
|
end if
|
|
do j = 1, max_loops
|
|
working_esize = remesh_size - 2*(j-1)
|
|
ele = (working_esize-1)*cubic_cell
|
|
zloop: do iz = 1, bd_in_array(3)
|
|
yloop: do iy = 1, bd_in_array(2)
|
|
xloop: do ix = 1, bd_in_array(1)
|
|
if (lat_points(ix, iy,iz)) then
|
|
r_new_node(:,:,:) = 0.0_dp
|
|
|
|
!Check to see if the element overshoots the bound
|
|
if (iz+working_esize-1 > bd_in_array(3)) then
|
|
exit zloop
|
|
else if (iy+working_esize-1 > bd_in_array(2)) then
|
|
cycle zloop
|
|
else if (ix+working_esize-1 > bd_in_array(1)) then
|
|
cycle yloop
|
|
end if
|
|
|
|
if (all(lat_points(ix:ix+working_esize-1,iy:iy+working_esize-1,iz:iz+working_esize-1))) then
|
|
do inod = 1, ng_node(ilat)
|
|
vlat = ele(:,inod) + (/ix, iy, iz /)
|
|
do i = 1, 3
|
|
vlat(i) = vlat(i) + bd_in_lat(2*i-1)-5
|
|
end do
|
|
r_new_node(:,1,inod) = matmul(orient, matmul(fcc_mat, vlat))*lapa(ilat)
|
|
end do
|
|
|
|
lat_points(ix:ix+working_esize-1,iy:iy+working_esize-1,iz:iz+working_esize-1) = .false.
|
|
!Add the element, for the sbox we just set it to the same sbox that we get the orientation
|
|
!from. In this case it is from the sbox of the first atom in the group.
|
|
new_ele = new_ele+1
|
|
call add_element(0,remesh_ele_type, working_esize, ilat, &
|
|
sbox_atom(atom_index(1)),r_new_node)
|
|
|
|
end if
|
|
end if
|
|
end do xloop
|
|
end do yloop
|
|
end do zloop
|
|
end do
|
|
|
|
!Now we have to add any leftover lattice points as atoms
|
|
do iz = 1, bd_in_array(3)
|
|
do iy=1, bd_in_array(2)
|
|
do ix = 1, bd_in_array(1)
|
|
if(lat_points(ix,iy,iz)) then
|
|
vlat = (/ ix, iy, iz /)
|
|
do i = 1, 3
|
|
vlat(i) = vlat(i) + bd_in_lat(2*i-1)-5
|
|
end do
|
|
lat_points(ix,iy,iz) = .false.
|
|
r = matmul(orient, matmul(fcc_mat, vlat))*lapa(ilat)
|
|
new_atom=new_atom+1
|
|
call add_atom(0,basis_type(1,ilat), is, r)
|
|
end if
|
|
end do
|
|
end do
|
|
end do
|
|
deallocate(lat_points)
|
|
end do
|
|
end do
|
|
|
|
!Delete all elements and atoms to make space for new elements and atoms
|
|
call delete_atoms(group_atom_num, atom_index)
|
|
call delete_elements(group_ele_num, element_index)
|
|
|
|
print *, tot_dof, " degrees of freedom in group"
|
|
print *, "remesh_group creates ", new_atom, " atoms and ", new_ele, " elements."
|
|
end subroutine remesh_group
|
|
|
|
subroutine delete_group
|
|
!This subroutine deletes all atoms/elements within a group
|
|
|
|
print *, "Deleting group containing ", group_atom_num, " atoms and ", group_ele_num, " elements."
|
|
|
|
!Delete atoms
|
|
call delete_atoms(group_atom_num, atom_index)
|
|
|
|
!Delete elements
|
|
call delete_elements(group_ele_num, element_index)
|
|
|
|
end subroutine delete_group
|
|
|
|
subroutine change_group_type
|
|
!This subroutine changes all atoms and nodes at atoms within a group to a specific type
|
|
integer :: i, j, ltype,ibasis, inod, basis_type(10)
|
|
|
|
print *, "Changing ", group_atom_num, " atoms and ", group_ele_num, " elements to atom type ", group_type
|
|
|
|
!Change all atom group types
|
|
do i = 1, group_atom_num
|
|
j = atom_index(i)
|
|
type_atom(j) = group_type
|
|
end do
|
|
|
|
!Map to a new lattice type for all element
|
|
do i =1, group_ele_num
|
|
j = element_index(i)
|
|
ltype = lat_ele(j)
|
|
do ibasis = 1, basisnum(ltype)
|
|
basis_type(ibasis) = group_type
|
|
end do
|
|
call lattice_map(basisnum(ltype), basis_type, ng_node(ltype), lapa(ltype), lat_ele(j))
|
|
end do
|
|
|
|
end subroutine change_group_type
|
|
|
|
|
|
subroutine split_group_elements
|
|
!
|
|
end subroutine split_group_elements
|
|
|
|
function in_group(r)
|
|
!This subroutine determines if a point is within the group boundaries
|
|
real(kind=dp), intent(in) :: r(3)
|
|
real(kind=dp) :: rnorm
|
|
|
|
integer :: dim3, i
|
|
logical :: in_group
|
|
select case(trim(adjustl(shape)))
|
|
case('block')
|
|
in_group=in_block_bd(r,block_bd)
|
|
case('wedge')
|
|
in_group = in_wedge_bd(r,vertices)
|
|
case('notch')
|
|
do i = 1, 3
|
|
if (.not.((dim1==i).or.(dim2==i))) dim3 = i
|
|
end do
|
|
|
|
in_group = in_wedge_bd(r,vertices)
|
|
!Do a check to make sure the wedge isn't used if it should be the tip radius
|
|
if (bwidth>0) then
|
|
if (r(dim1) < centroid(dim1)) in_group = .false.
|
|
end if
|
|
|
|
rnorm = sqrt((r(dim1) - centroid(dim1))**2 + (r(dim3)-centroid(dim3))**2)
|
|
in_group = in_group.or.(rnorm < radius)
|
|
|
|
case('sphere')
|
|
rnorm = norm2(r(:) - centroid(:))
|
|
if (rnorm <= radius) then
|
|
in_group = .true.
|
|
else
|
|
in_group = .false.
|
|
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')
|
|
in_group = .true.
|
|
end select
|
|
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
|