Added efill zigzag code which adjusts the efill option

src/mode_create.f90

@@ -10,12 +10,12 @@ module mode_create
 
     implicit none
 
-    character(len=100) :: name, element_type
+    character(len=100) :: name, element_type    
     real(kind = dp) :: lattice_parameter, orient(3,3), cell_mat(3,8), box_len(3), basis(3,3), origin(3), maxlen(3), &
                        orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3), duplicate(3)
     integer :: esize, ix, iy, iz, box_lat_vert(3,8), lat_ele_num, lat_atom_num, bd_in_lat(6), &
                basis_pos(3,10), esize_nums, esize_index(10)
-    logical :: dup_flag, dim_flag, efill
+    logical :: dup_flag, dim_flag, efill(3)
 
     real(kind=dp), allocatable :: r_lat(:,:,:), r_atom_lat(:,:)
     integer, allocatable :: elat(:)
@@ -46,6 +46,7 @@ module mode_create
         basisnum = 0
         lat_ele_num = 0
         lat_atom_num = 0
+        efill(:) = .false.
 
         !First we parse the command
         call parse_command(arg_pos)
@@ -264,7 +265,30 @@ module mode_create
                 end do
 
             case('efill')
-                efill = .true.
+                call get_command_argument(arg_pos, textholder)
+                select case(trim(adjustl(textholder)))
+                case('x')
+                    efill(1) = .true.
+                case('y')
+                    efill(2) = .true.
+                case('z')
+                    efill(3) = .true.
+                case('xy','yx')
+                    efill(1) = .true.
+                    efill(2) = .true.
+                case('yz','zy')
+                    efill(2) = .true.
+                    efill(3) = .true.
+                case('xz','zx')
+                    efill(1) = .true.
+                    efill(3) = .true.
+                case('xyz','xzy','yxz','yzx','zxy','zyx')
+                    efill(:) = .true.
+                case default
+                    print *, "Error: ", trim(adjustl(textholder)), " is not an acceptable argument for the efill argument"
+                    stop 3
+                end select
+                arg_pos = arg_pos + 1
             case default
                 !If it isn't an option then you have to exit
                 arg_pos = arg_pos -1
@@ -339,16 +363,20 @@ module mode_create
         !Internal variables
         integer :: i,  inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz,  numlatpoints, ele(3,8), rzero(3), &
                    vlat(3), temp_lat(3,8), m, n, o, curr_esize, ei
-        real(kind=dp) :: v(3), temp_nodes(3,1,8)
+        real(kind=dp) :: v(3), temp_nodes(3,1,8), r(3), centroid_bd(6)
         logical, allocatable :: lat_points(:,:,:)
-        logical :: node_in_bd(8)
+        logical :: node_in_bd(8), add
 
         !Do some value initialization
         max_esize = esize
+        do i = 1,3
+            centroid_bd(2*i) = -huge(1.0_dp)
+            centroid_bd(2*i-1) = huge(1.0_dp)
+        end do
 
         !Now initialize the code if we are doing efill. This means calculate the number of times we can divide the esize in 2 with
         !the value still being > 7
-        if(efill) then 
+        if(any(efill)) then 
             curr_esize=esize
             esize_nums=0
             do while (curr_esize >= 7) 
@@ -485,14 +513,18 @@ module mode_create
 
                                 !Check to see if the lattice point values are greater then the array limits
                                 if(any(vlat > shape(lat_points)).or.any(vlat < 1)) then 
-                                    continue
+                                    exit
                                 !If within array boundaries check to see if it is a lattice point
                                 else if(lat_points(vlat(1),vlat(2),vlat(3))) then 
                                     node_in_bd(inod) = .true.
+                                else
+                                    exit
                                 end if
                             end do
 
-                            if(all(node_in_bd)) then 
+                            !If we are on the first round of element building then we can just add the element if all(node_in_bd) is
+                            !true
+                            if(all(node_in_bd).and.(ei==1)) then 
                                 lat_ele_num = lat_ele_num+1
                                 r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:)
                                 elat(lat_ele_num) = curr_esize
@@ -505,11 +537,60 @@ module mode_create
                                     end do
                                 end do
 
+                            !Otherwise we have to also do a box boundary check
+                            else if(all(node_in_bd)) then
+                                r(:) = 0.0_dp
+                                add = .false.
+                                do inod = 1,8
+                                    r = r+ temp_nodes(:,1,inod)/8.0_dp
+                                end do
+
+                                !Here we check to make sure the centroid of the element we are adding is outside of the bounds set
+                                !by the centroids of the elements of the initial iteration
+
+                                do i = 1,3
+                                    if(efill(i)) then 
+                                        if((r(i) > centroid_bd(2*i)).or.(r(i) < centroid_bd(2*i-1)))then 
+                                            add = .true.
+                                            exit
+                                        end if
+                                    end if
+                                end do
+                                if(add) then 
+                                    lat_ele_num = lat_ele_num+1
+                                    r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:)
+                                    elat(lat_ele_num) = curr_esize
+                                    !Now set all the lattice points contained within an element to false
+                                    do o = minval(temp_lat(3,:)), maxval(temp_lat(3,:))
+                                        do n = minval(temp_lat(2,:)), maxval(temp_lat(2,:))
+                                            do m = minval(temp_lat(1,:)), maxval(temp_lat(1,:))
+                                                lat_points(m,n,o) = .false.
+                                            end do 
+                                        end do
+                                    end do
+                                end if
                             end if
                         end do
                     end do
                 end do
                 curr_esize=curr_esize-2
+                !If we are running efill code, after the first iteration we have to calculate the min and max element centroids in
+                !each dimension
+                if((ei == 1).and.(any(efill))) then 
+                    do i = 1, lat_ele_num
+                        !Calculate the current element centroid
+                        r(:) = 0.0_dp
+                        do inod = 1,8
+                            r = r + r_lat(:,inod,i)/8.0_dp
+                        end do
+                        
+                        !Check to see if it's a min or max
+                        do o = 1,3
+                            if(r(o) > centroid_bd(2*o)) centroid_bd(2*o) = r(o)
+                            if(r(o) < centroid_bd(2*o-1)) centroid_bd(2*o-1) = r(o)
+                        end do
+                    end do
+                end if
             end do
             !Now figure out how many lattice points could not be contained in elements
             allocate(r_atom_lat(3,count(lat_points)))