Initial testing of efill code

src/mode_create.f90

@@ -14,8 +14,8 @@ module mode_create
     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)
-    logical :: dup_flag, dim_flag
+               basis_pos(3,10), esize_nums, esize_index(10)
+    logical :: dup_flag, dim_flag, efill
 
     real(kind=dp), allocatable :: r_lat(:,:,:), r_atom_lat(:,:)
     public 
@@ -26,7 +26,7 @@ module mode_create
 
         integer, intent(out) :: arg_pos
 
-        integer :: i, ibasis, inod
+        integer :: i, ibasis, inod, ei, curr_esize
         real(kind=dp), allocatable :: r_node_temp(:,:,:)
 
         print *, '-----------------------Mode Create---------------------------'
@@ -148,7 +148,15 @@ module mode_create
                             r_node_temp(:,ibasis,inod) = (r_lat(:,inod,i)*lattice_parameter)+basis_pos(:,ibasis)
                         end do
                     end do 
-                    call add_element(element_type, esize, 1, 1, r_node_temp)
+
+                    curr_esize=esize
+                    do ei = 1, esize_nums
+                        if(i < esize_index(ei)) then 
+                            call add_element(element_type, curr_esize, 1, 1, r_node_temp)
+                            exit
+                        end if
+                        curr_esize=curr_esize/2 + 1
+                    end do
                 end do 
             end if
         end if
@@ -248,6 +256,9 @@ module mode_create
                     end do
                 end do
 
+            case('efill')
+                arg_pos=arg_pos+1
+                efill = .true.
             case default
                 !If it isn't an option then you have to exit
                 arg_pos = arg_pos -1
@@ -314,7 +325,7 @@ module mode_create
         real(kind=dp), dimension(3,3), intent(in) :: transform_matrix !The transformation matrix from lattice_space to real space
         !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
+                   vlat(3), temp_lat(3,8), m, n, o, curr_esize, ei
         real(kind=dp) :: v(3), temp_nodes(3,1,8)
         logical, allocatable :: lat_points(:,:,:)
         logical :: node_in_bd(8)
@@ -322,6 +333,19 @@ module mode_create
         !Do some value initialization
         max_esize = esize
 
+        !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 
+            curr_esize=esize
+            esize_nums=0
+            do while (curr_esize >= 7) 
+                esize_nums=esize_nums+1
+                curr_esize = curr_esize/2 + 1 
+            end do
+        else
+            esize_nums=1
+        end if
+
         !First find the bounding lattice points (min and max points for the box in each dimension)
         numlatpoints = 1
         do i = 1, 3
@@ -415,11 +439,15 @@ module mode_create
             !Now build the finite element region
             lat_ele_num = 0
             lat_atom_num = 0
-            allocate(r_lat(3,8,numlatpoints/esize))
+            curr_esize=esize/(2**(esize_nums-1)) + 1
+            allocate(r_lat(3,8,numlatpoints/curr_esize))
 
+            curr_esize=esize 
+            do ei = 1, esize_nums
+                ele(:,:) = (curr_esize-1) * cubic_cell(:,:)
                 !Redefined the second 3 indices as the number of elements that fit in the bounds
                 do i = 1, 3
-                bd_in_array(3+i) = bd_in_array(i)/esize
+                    bd_in_array(3+i) = bd_in_array(i)/curr_esize
                 end do
 
                 !Now start the element at rzero
@@ -433,7 +461,7 @@ module mode_create
                             temp_nodes(:,:,:) = 0.0_dp
                             temp_lat(:,:) = 0
                             do inod = 1, 8
-                            vlat= ele(:,inod) + (/ ix*(esize), iy*(esize), iz*(esize) /)
+                                vlat= ele(:,inod) + (/ ix*(curr_esize), iy*(curr_esize), iz*(curr_esize) /)
                                 !Transform point back to real space for easier checking
                                 ! v = matmul(orient, matmul(transform_matrix,v))
                                 do i = 1,3
@@ -468,7 +496,9 @@ module mode_create
                         end do
                     end do
                 end do
-
+                esize_index(ei) = lat_ele_num
+                curr_esize=curr_esize/2 + 1
+            end do
             !Now figure out how many lattice points could not be contained in elements
             allocate(r_atom_lat(3,count(lat_points)))
             lat_atom_num = 0