Debugging version of code
This commit is contained in:
Alex
parent
d0b6c595f0
commit
20755270a4
4 changed files with 167 additions and 34 deletions
|
|
@ -8,10 +8,10 @@ module opt_group
|
|||
use box
|
||||
implicit none
|
||||
|
||||
integer :: group_ele_num, group_atom_num, remesh_size
|
||||
integer :: group_ele_num, group_atom_num, remesh_size, remesh_type
|
||||
character(len=15) :: type, shape !Type indicates what element type is selected and shape is the group shape
|
||||
real(kind=dp) :: block_bd(6), disp_vec(3)
|
||||
logical :: displace, delete
|
||||
real(kind=dp) :: block_bd(6), disp_vec(3), remesh_lat_pam
|
||||
logical :: displace, delete, max_remesh, refine
|
||||
|
||||
integer, allocatable :: element_index(:), atom_index(:)
|
||||
|
||||
|
|
@ -29,6 +29,8 @@ module opt_group
|
|||
remesh_size=0
|
||||
displace=.false.
|
||||
delete=.false.
|
||||
! max_remesh=.false.
|
||||
refine = .false.
|
||||
|
||||
if(allocated(element_index)) deallocate(element_index)
|
||||
if(allocated(atom_index)) deallocate(atom_index)
|
||||
|
|
@ -43,6 +45,9 @@ module opt_group
|
|||
if(remesh_size > 0) call remesh_group
|
||||
|
||||
if(delete) call delete_group
|
||||
|
||||
if(refine) call refine_group
|
||||
|
||||
end subroutine group
|
||||
|
||||
subroutine parse_group(arg_pos)
|
||||
|
|
@ -101,13 +106,24 @@ module opt_group
|
|||
if (arglen==0) stop "Missing vector component for shift command"
|
||||
read(textholder, *) disp_vec(i)
|
||||
end do
|
||||
case('refine')
|
||||
refine=.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('delete')
|
||||
arg_pos = arg_pos + 1
|
||||
call get_command_argument(arg_pos, textholder, arglen)
|
||||
if (arglen==0) stop "Missing remesh lattice parameter in group command"
|
||||
read(textholder, *) remesh_lat_pam
|
||||
arg_pos = arg_pos + 1
|
||||
call get_command_argument(arg_pos, textholder, arglen)
|
||||
if (arglen==0) stop "Missing remesh type in group command"
|
||||
read(textholder, *) remesh_type
|
||||
case('max')
|
||||
max_remesh =.true.
|
||||
case('delete')
|
||||
delete=.true.
|
||||
case default
|
||||
!If it isn't an available option to opt_disl then we just exit
|
||||
|
|
@ -211,46 +227,152 @@ module opt_group
|
|||
|
||||
end subroutine displace_group
|
||||
|
||||
subroutine remesh_group
|
||||
subroutine refine_group
|
||||
!This command is used to remesh the group to a desired element size
|
||||
|
||||
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 and remeshing to elements are different processes so check which code we need to run
|
||||
select case(remesh_size)
|
||||
|
||||
!Refining to atoms
|
||||
case(2)
|
||||
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)
|
||||
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(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
|
||||
|
||||
!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(type_interp(j), sbox_ele(ie), r_interp(:,j))
|
||||
end do
|
||||
end subroutine
|
||||
|
||||
subroutine remesh_group
|
||||
!This command is used to remesh the group to a desired element size
|
||||
|
||||
integer :: i, j, 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)
|
||||
real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), ori_inv(3,3), &
|
||||
r_new_node(3,max_basisnum, max_ng_node), orient(3,3), group_in_lat(3,8)
|
||||
logical, allocatable :: lat_points(:,:,:)
|
||||
character(len=100) :: remesh_ele_type
|
||||
|
||||
|
||||
!Right now we just hardcode only remeshing to elements
|
||||
remesh_ele_type = 'fcc'
|
||||
|
||||
!Get the orientations, this assumes that the orientation of the subbox for the first atom is the
|
||||
!same as the rest of the atoms
|
||||
!If this assumption is false then the code will break and exit
|
||||
orient = sub_box_ori(:, :, sbox_atom(atom_index(1)))
|
||||
call matrix_inverse(orient,3,ori_inv)
|
||||
|
||||
!First calculate max position in lattice space to be able to allocate lat_points array, also sum the total numbers of
|
||||
!degrees of freedom which are added
|
||||
dof = 0
|
||||
|
||||
select case(trim(adjustl(shape)))
|
||||
case('block')
|
||||
group_in_lat = reshape((/ block_bd(1),block_bd(3),block_bd(5), &
|
||||
block_bd(2),block_bd(3),block_bd(5), &
|
||||
block_bd(2),block_bd(4),block_bd(5), &
|
||||
block_bd(1),block_bd(4),block_bd(5), &
|
||||
block_bd(1),block_bd(3),block_bd(6), &
|
||||
block_bd(2),block_bd(3),block_bd(6), &
|
||||
block_bd(2),block_bd(4),block_bd(6), &
|
||||
block_bd(1),block_bd(4),block_bd(6) /), [3,8])
|
||||
|
||||
group_in_lat = matmul(fcc_inv, matmul(ori_inv, group_in_lat/remesh_lat_pam))
|
||||
do i = 1, 3
|
||||
bd_in_lat(2*i-1) = minval(group_in_lat(i,:))
|
||||
bd_in_lat(2*i) = maxval(group_in_lat(i,:))
|
||||
end do
|
||||
|
||||
!Once all atoms are added we delete all of the elements
|
||||
call delete_elements(group_ele_num, element_index)
|
||||
end select
|
||||
|
||||
print *, group_ele_num, " elements of group are refined to ", atom_num -orig_atom_num, " atoms."
|
||||
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.
|
||||
|
||||
!Now place all group atoms and group interpolated atoms into lat_points
|
||||
do i = 1, group_atom_num
|
||||
r_lat = r_atom(:,atom_index(i))/remesh_lat_pam
|
||||
r_lat = matmul(fcc_inv,matmul(ori_inv,r_lat))
|
||||
!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((r_lat(1)==13).and.(r_lat(2)==13).and.(r_lat(3)==-13)) then
|
||||
print *, i, atom_index(i), r_atom(:,atom_index(i))
|
||||
end if
|
||||
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.
|
||||
end if
|
||||
end do
|
||||
|
||||
!Now place interpolated atoms within lat_points array
|
||||
do i =1, group_ele_num
|
||||
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_lat = matmul(fcc_inv,matmul(ori_inv, r_interp(:,atom_index(i))/remesh_lat_pam)) + 1
|
||||
!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.
|
||||
end if
|
||||
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)
|
||||
|
||||
!Now run remeshing algorithm, not the most optimized or efficient but gets the job done
|
||||
|
||||
ele = (remesh_size-1)*cubic_cell
|
||||
|
||||
!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
|
||||
|
||||
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
|
||||
if (all(lat_points(ix:ix+remesh_size-1,iy:iy+remesh_size-1,iz:iz+remesh_size-1))) then
|
||||
do inod = 1, 8
|
||||
vlat = ele(:,inod) + (/ix, iy, iz /)
|
||||
do i = 1, 3
|
||||
vlat = vlat + bd_in_lat(2*i-1)-5
|
||||
end do
|
||||
r_new_node(:,1,inod) = matmul(orient, matmul(fcc_mat, vlat))*remesh_lat_pam
|
||||
end do
|
||||
|
||||
lat_points(ix:ix+remesh_size-1,iy:iy+remesh_size-1,iz:iz+remesh_size-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.
|
||||
call add_element(remesh_ele_type, remesh_size, remesh_type, sbox_atom(atom_index(1)),r_new_node)
|
||||
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
end if
|
||||
!Remeshing to elements, currently not available
|
||||
case default
|
||||
print *, "Remeshing to elements is currently not available. Please refine to atoms by passing a remsh size of 2"
|
||||
end select
|
||||
|
||||
|
||||
end subroutine remesh_group
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue