Merge pull request #16 from aselimov/ft--group-refine

Added remesh option to group

src/elements.f90

@@ -223,6 +223,8 @@ module elements
         real(kind=dp), intent(in) :: r(3, max_basisnum, max_ng_node)
 
         ele_num = ele_num + 1
+        !Check to see if we need to grow the arrays
+        call grow_ele_arrays(1,0)
         type_ele(ele_num) = type
         size_ele(ele_num) = size
         lat_ele(ele_num) = lat
@@ -239,6 +241,8 @@ module elements
         real(kind=dp), intent(in), dimension(3) :: r
 
         atom_num = atom_num+1
+        !Check to see if we need to grow the arrays
+        call grow_ele_arrays(0,1)
         type_atom(atom_num) = type
         r_atom(:,atom_num) = r(:)
 
@@ -356,7 +360,7 @@ module elements
             end do 
         end select
 
-        if (ia /= esize**3) then 
+        if (ia /= bnum*esize**3) then 
             print *, "Incorrect interpolation"
             stop 3
         end if
@@ -403,4 +407,34 @@ module elements
             atom_num = atom_num - 1
         end do 
     end subroutine
+
+    subroutine delete_elements(num, index)
+        !This subroutine deletes elements from the element array 
+        integer, intent(in) :: num
+        integer, intent(inout), dimension(num) :: index
+
+        integer :: i, j
+
+        call heapsort(index)
+
+        !We go from largest index to smallest index just to make sure that we don't miss 
+        !accidentally overwrite values which need to be deleted
+        do i = num, 1, -1
+            if(index(i) == ele_num) then 
+                r_node(:,:,:,index(i)) = 0.0_dp
+                type_ele(index(i)) =''
+                size_ele(index(i)) = 0
+                lat_ele(index(i)) = 0
+                sbox_ele(index(i)) = 0
+            else
+                node_num = node_num - ng_node(lat_ele(index(i)))
+                r_node(:,:,:,index(i)) = r_node(:,:,:,ele_num)
+                type_ele(index(i)) = type_ele(ele_num)
+                size_ele(index(i)) = size_ele(ele_num)
+                lat_ele(index(i)) = lat_ele(ele_num)
+                sbox_ele(index(i)) = sbox_ele(ele_num)
+            end if            
+            ele_num = ele_num - 1
+        end do 
+    end subroutine delete_elements
 end module elements

src/opt_group.f90

@@ -8,7 +8,7 @@ module opt_group
     use box 
     implicit none
 
-    integer :: group_ele_num, group_atom_num
+    integer :: group_ele_num, group_atom_num, remesh_size
     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, wrap
@@ -24,6 +24,7 @@ module opt_group
 
         group_ele_num = 0
         group_atom_num = 0
+        remesh_size=0
         if(allocated(element_index)) deallocate(element_index)
         if(allocated(atom_index)) deallocate(atom_index)
 
@@ -34,6 +35,7 @@ module opt_group
         !Now call the transformation functions for the group
         if(displace) call displace_group
 
+        if(remesh_size > 0) call remesh_group
     end subroutine group
 
     subroutine parse_group(arg_pos)
@@ -94,6 +96,11 @@ module opt_group
                 end do
             case('wrap')
                 wrap = .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 default
                 !If it isn't an available option to opt_disl then we just exit
                 exit 
@@ -197,4 +204,44 @@ module opt_group
 
     end subroutine displace_group
 
+    subroutine remesh_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
+        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 
+                !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), 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)
+
+            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
 end module opt_group