From 2c425201028e967be1b261562607dae6508db1fe Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Wed, 2 Dec 2020 16:28:10 -0500 Subject: [PATCH 1/2] Restructure group code, add group shell, and add refinefill option --- src/elements.f90 | 3 +- src/opt_group.f90 | 277 ++++++++++++++++++++++++++++++++++++++-------- 2 files changed, 231 insertions(+), 49 deletions(-) diff --git a/src/elements.f90 b/src/elements.f90 index 5d3b3c6..fffa9cb 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -779,8 +779,7 @@ do i = 1, atom_num end do end do - - end subroutine + end subroutine get_interp_pos subroutine alloc_dat_arrays(n,m) !This subroutine used to provide initial allocation for the atom and element data arrays diff --git a/src/opt_group.f90 b/src/opt_group.f90 index 6fe2896..384b369 100644 --- a/src/opt_group.f90 +++ b/src/opt_group.f90 @@ -10,8 +10,8 @@ module opt_group 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 - logical :: displace, delete, max_remesh, refine, group_nodes, flip + 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(:) @@ -22,7 +22,9 @@ module opt_group !Main calling function for the group option integer, intent(inout) :: arg_pos - print *, '-----------------------Option Group-------------------------' + print *, '------------------------------------------------------------' + print *, 'Option Group' + print *, '------------------------------------------------------------' group_ele_num = 0 group_atom_num = 0 @@ -48,6 +50,11 @@ module opt_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 @@ -335,6 +342,28 @@ module opt_group 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 @@ -364,6 +393,8 @@ module opt_group 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) @@ -382,6 +413,8 @@ module opt_group 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) @@ -403,7 +436,7 @@ module opt_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 + integer :: i, j, inod, ibasis, temp, node_in_out(max_ng_node) integer, allocatable :: resize_array(:) 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 select case(trim(adjustl(type))) case('elements', 'both') - if(.not.(group_nodes)) then - do i = 1, ele_num - r_center(:) = 0.0_dp - do inod = 1, ng_node(lat_ele(i)) - do ibasis = 1, basisnum(lat_ele(i)) - r_center = r_center + r_node(:,ibasis,inod,i)/(basisnum(lat_ele(i))*ng_node(lat_ele(i))) - end do - end do - - 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 + 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 - else if(group_nodes) then - eleloop:do i = 1, ele_num - r_center(:) = 0.0_dp - 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 + 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 @@ -569,7 +573,7 @@ module opt_group end subroutine displace_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 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." 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 !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 + + 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) @@ -933,8 +1039,85 @@ module opt_group 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 From 29df0f9eb2ae79d8444e49d33e58df975c099bd7 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Thu, 3 Dec 2020 09:21:31 -0500 Subject: [PATCH 2/2] Fix to opt_group which got the code working --- src/opt_group.f90 | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/src/opt_group.f90 b/src/opt_group.f90 index 384b369..603bf6e 100644 --- a/src/opt_group.f90 +++ b/src/opt_group.f90 @@ -606,7 +606,7 @@ module opt_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 + 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) @@ -623,7 +623,6 @@ module opt_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) @@ -644,6 +643,7 @@ module opt_group !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 @@ -667,6 +667,7 @@ module opt_group 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 @@ -676,20 +677,25 @@ module opt_group 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) + 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 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)