Merge branch 'ft--group-shell' into development

src/elements.f90

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

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
+    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(:) 
 
@@ -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
-                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 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,110 @@ 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, 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
@@ -900,6 +1007,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 +1045,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