Current working changes to option-slip-plane

src/Makefile

@@ -1,8 +1,8 @@
 FC=ifort
 
 #Ifort flags
-#FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin -heap-arrays
-FFLAGS=-mcmodel=large -Ofast -no-wrap-margin -heap-arrays
+FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin -heap-arrays
+#FFLAGS=-mcmodel=large -Ofast -no-wrap-margin -heap-arrays
 
 #gfortran flags
 #FFLAGS=-mcmodel=large -O3 -g

src/elements.f90

@@ -201,24 +201,25 @@ module elements
         
         !First check to make sure if it is allocated
         if (allocated(size_ele)) then 
+
             !Figure out the size of the atom and element arrays
             ele_size = size(size_ele)
 
             !Check if we need to grow the ele_size, if so grow all the variables
-            if ( n+ele_size > size(size_ele)) then 
+            if ( n+ele_num > size(size_ele)) then 
                 
                 allocate(temp_int(n+ele_size+buffer_size))
-                temp_int(1:ele_size) = lat_ele
+                temp_int(1:ele_size) = lat_ele(1:ele_size)
                 temp_int(ele_size+1:) = 0
                 call move_alloc(temp_int, lat_ele)
 
                 allocate(temp_int(n+ele_size+buffer_size))
-                temp_int(1:ele_size) = tag_ele
+                temp_int(1:ele_size) = tag_ele(1:ele_size)
                 temp_int(ele_size+1:) = 0
                 call move_alloc(temp_int, tag_ele)
 
                 allocate(temp_int(n+ele_size+buffer_size))
-                temp_int(1:ele_size) = size_ele
+                temp_int(1:ele_size) = size_ele(1:ele_size)
                 temp_int(ele_size+1:) = 0
                 call move_alloc(temp_int, size_ele)
 
@@ -278,6 +279,7 @@ module elements
         integer :: newtag
 
         ele_num = ele_num + 1
+        node_num = node_num + ng_node(lat)
 
         if (tag==0) then
             newtag = ele_num !If we don't assign a tag then pass the tag as the ele_num
@@ -293,7 +295,6 @@ module elements
         lat_ele(ele_num) = lat
         sbox_ele(ele_num) = sbox
         r_node(:,:,:,ele_num) = r(:,:,:)
-        node_num = node_num + ng_node(lat)
 
 
     end subroutine add_element

src/opt_slip_plane.f90

@@ -17,10 +17,10 @@ module opt_slip_plane
         !Main calling function for the slip_plane option
         integer, intent(inout) :: arg_pos 
 
-        integer :: ie, ia, slip_enum, old_atom_num, esize, new_ele_num, n, m, o, ele(3,8), nump_ele, inod, vlat(3)
+        integer :: ie, ia, slip_enum, old_atom_num, esize, new_ele_num, n, m, o, ele(3,8), nump_ele, inod, vlat(3), ibasis
                 
         integer, allocatable :: slip_eles(:), temp_int(:)
-        real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), rfill(3,max_basisnum, max_ng_node)
+        real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), rfill(3,max_basisnum, max_ng_node), ratom(3,max_basisnum)
 
         integer :: type_interp(max_basisnum*max_esize**3)
         logical :: lat_points(max_esize,max_esize, max_esize)
@@ -86,14 +86,28 @@ module opt_slip_plane
                                         if((spos < maxval(rfill(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie))))).and. &
                                             (spos > minval(rfill(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie)))))) then 
                                             nump_ele = nump_ele - esize**3 
-                                            lat_points(m:m+esize, n:n+esize, o:o+esize) = .false.
+                                            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
                                 end do
                             end do
                         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 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 if
             end if
         end do 
@@ -104,6 +118,8 @@ module opt_slip_plane
         !Output data
         if(.not.efill) then 
             print *, "We refine ", slip_enum, " elements into ", atom_num - old_atom_num , " atoms"
+        else 
+            print *, "We refine ", slip_enum, " elements into ", atom_num - old_atom_num , " atoms and ", new_ele_num, " elements"
         end if
 
     end subroutine run_slip_plane