Working remesh command and remesh max command

src/io.f90

@@ -310,6 +310,9 @@ module io
 20      format('SCALARS lattice_type float', /&
                'LOOKUP_TABLE default')
 
+21      format('SCALARS esize float', /&
+               'LOOKUP_TABLE default')
+
         !First we write the vtk file containing the atoms
         open(unit=11, file='atoms_'//trim(adjustl(file)), action='write', status='replace',position='rewind')
 
@@ -358,6 +361,10 @@ module io
         do i = 1, ele_num
             write(11, '(i16)') lat_ele(i)
         end do
+        write(11,21)
+        do i = 1, ele_num
+            write(11, '(i16)') size_ele(i)
+        end do
         close(11)
     end subroutine
 
@@ -624,8 +631,13 @@ module io
         !Read in the box boundary and grow the current active box bd 
         read(11, *) temp_box_bd(:)
 
+        print *, "displace", displace
         do i = 1, 3
-            newdisplace(i) = displace(i) - temp_box_bd(2*i-1)
+            if (displace(i) > lim_zero) then 
+                newdisplace(i) = displace(i) - temp_box_bd(2*i-1)
+            else
+                newdisplace=displace(i)
+            end if 
             temp_box_bd(2*i-1) = temp_box_bd(2*i-1) + newdisplace(i)
             temp_box_bd(2*i) = temp_box_bd(2*i) + newdisplace(i)
         end do 
@@ -717,7 +729,8 @@ module io
         !Read the atoms 
         do i = 1, in_atoms
             read(11,*) j, type, sbox, r(:)
-            call add_atom(new_type_to_type(type), sbox, r+newdisplace )
+            r = r+newdisplace
+            call add_atom(new_type_to_type(type), sbox, r)
         end do
 
         !Read the elements

src/opt_group.f90

@@ -29,7 +29,7 @@ module opt_group
         remesh_size=0
         displace=.false.
         delete=.false.
-        ! max_remesh=.false.
+        max_remesh=.false.
         refine = .false.
 
         if(allocated(element_index)) deallocate(element_index)
@@ -262,7 +262,9 @@ module opt_group
         !This command is used to remesh the group to a desired element size
 
         integer :: i, j, k, ix, iy, iz, inod, ibasis, ie, type_interp(max_basisnum*max_esize**3), add_atom_num, orig_atom_num,   &
-                   current_esize, dof, max_lat(3), r_lat(3), ele(3,8), vlat(3), bd_in_lat(6), bd_in_array(3), old_ele, old_atom
+                   current_esize, dof, max_lat(3), r_lat(3), ele(3,8), vlat(3), bd_in_lat(6), bd_in_array(3), old_ele, old_atom, &
+                   max_loops, working_esize
+
         real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3),  ori_inv(3,3),  r(3), &
                          r_new_node(3,max_basisnum, max_ng_node), orient(3,3), group_in_lat(3,8)
         logical, allocatable :: lat_points(:,:,:)
@@ -352,48 +354,58 @@ module opt_group
         old_atom = atom_num
         old_ele = ele_num
 
-        !Now run remeshing algorithm, not the most optimized or efficient but gets the job done
 
-        ele = (remesh_size-1)*cubic_cell
 
+        !Now run remeshing algorithm, not the most optimized or efficient but gets the job done
         !Figure out new looping boundaries
         bd_in_array(1) = bd_in_lat(2) - bd_in_lat(1) + 10
         bd_in_array(2) = bd_in_lat(4) - bd_in_lat(3) + 10
         bd_in_array(3) = bd_in_lat(6) - bd_in_lat(5) + 10
 
-        zloop: do iz = 1, bd_in_array(3)
-            yloop: do iy = 1, bd_in_array(2)
-                xloop: do ix = 1, bd_in_array(1)
-                    if (lat_points(ix, iy,iz)) then 
-                        
-                        !Check to see if the element overshoots the bound
-                        if (iz+remesh_size-1 > bd_in_array(3)) then
-                            exit zloop
-                        else if (iy+remesh_size-1 > bd_in_array(2)) then
-                            cycle zloop
-                        else if (ix+remesh_size-1 > bd_in_array(1)) then 
-                            cycle yloop
-                        end if  
 
-                        if (all(lat_points(ix:ix+remesh_size-1,iy:iy+remesh_size-1,iz:iz+remesh_size-1))) then 
+        if (max_remesh) then 
-                            do inod = 1, 8
+            max_loops = (remesh_size-2)/2
-                                vlat = ele(:,inod) + (/ix, iy, iz /)
+        else 
-                                do i = 1, 3
+            max_loops = 1
-                                    vlat(i) = vlat(i) + bd_in_lat(2*i-1)-5
+        end if
+        do j = 1, max_loops
+            working_esize = remesh_size - 2*(j-1)
+            ele = (working_esize-1)*cubic_cell
+            zloop: do iz = 1, bd_in_array(3)
+                yloop: do iy = 1, bd_in_array(2)
+                    xloop: do ix = 1, bd_in_array(1)
+                        if (lat_points(ix, iy,iz)) then 
+                            r_new_node(:,:,:) = 0.0_dp
+                            
+                            !Check to see if the element overshoots the bound
+                            if (iz+working_esize-1 > bd_in_array(3)) then
+                                exit zloop
+                            else if (iy+working_esize-1 > bd_in_array(2)) then
+                                cycle zloop
+                            else if (ix+working_esize-1 > bd_in_array(1)) then 
+                                cycle yloop
+                            end if  
+
+                            if (all(lat_points(ix:ix+working_esize-1,iy:iy+working_esize-1,iz:iz+working_esize-1))) then 
+                                do inod = 1, ng_node(remesh_type)
+                                    vlat = ele(:,inod) + (/ix, iy, iz /)
+                                    do i = 1, 3
+                                        vlat(i) = vlat(i) + bd_in_lat(2*i-1)-5
+                                    end do 
+                                   r_new_node(:,1,inod) = matmul(orient, matmul(fcc_mat, vlat))*remesh_lat_pam
                                 end do 
-                               r_new_node(:,1,inod) = matmul(orient, matmul(fcc_mat, vlat))*remesh_lat_pam
-                            end do 
 
-                            lat_points(ix:ix+remesh_size-1,iy:iy+remesh_size-1,iz:iz+remesh_size-1) = .false.
+                                lat_points(ix:ix+working_esize-1,iy:iy+working_esize-1,iz:iz+working_esize-1) = .false.
-                            !Add the element, for the sbox we just set it to the same sbox that we get the orientation from.
+                                !Add the element, for the sbox we just set it to the same sbox that we get the orientation from.
-                            !In this case it is from the sbox of the first atom in the group.
+                                !In this case it is from the sbox of the first atom in the group.
-                            call add_element(remesh_ele_type, remesh_size, remesh_type, sbox_atom(atom_index(1)),r_new_node)
+                                call add_element(remesh_ele_type, working_esize, remesh_type, sbox_atom(atom_index(1)),r_new_node)
 
+                            end if
                         end if
-                    end if
+                    end do xloop
-                end do xloop
+                end do yloop
-            end do yloop
+            end do zloop
-        end do zloop
+        end do 
 
         !Now we have to add any leftover lattice points as atoms
         do iz = 1, bd_in_array(3)