Working element refinement

src/opt_group.f90

@@ -261,7 +261,7 @@ module opt_group
     subroutine remesh_group
         !This command is used to remesh the group to a desired element size
 
-        integer :: i, j, ix, iy, iz, inod, ibasis, ie, type_interp(max_basisnum*max_esize**3), add_atom_num, orig_atom_num,   &
+        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)
         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)
@@ -313,9 +313,6 @@ module opt_group
             end do 
             !Do a check to make sure the code is working and that lattice points aren't being written on top of each other. 
             !This is primarily a debugging statement
-            if((r_lat(1)==13).and.(r_lat(2)==13).and.(r_lat(3)==-13)) then 
-                print *, i, atom_index(i), r_atom(:,atom_index(i))
-            end if 
             if(lat_points(r_lat(1)-bd_in_lat(1)+5,r_lat(2)-bd_in_lat(3)+5,r_lat(3)-bd_in_lat(5)+5)) then 
                 stop "Multiple atoms share same position in lat point array, this shouldn't happen"
             else 
@@ -328,7 +325,11 @@ module opt_group
             ie = element_index(i)
             call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp)
             do j = 1, size_ele(ie)**3 * basisnum(lat_ele(ie))
-                r_lat = matmul(fcc_inv,matmul(ori_inv, r_interp(:,atom_index(i))/remesh_lat_pam)) + 1
+                r = r_interp(:,j)/remesh_lat_pam
+                r = matmul(fcc_inv,matmul(ori_inv,r))
+                do k = 1, 3
+                    r_lat(k) = nint(r(k))
+                end do 
                 !Do a check to make sure the code is working and that lattice points aren't being written on top of each other. 
                 !This is primarily a debugging statement
                 if(lat_points(r_lat(1)-bd_in_lat(1)+5,r_lat(2)-bd_in_lat(3)+5,r_lat(3)-bd_in_lat(5)+5)) then 
@@ -395,6 +396,7 @@ module opt_group
                         do i = 1, 3
                             vlat(i) = vlat(i) + bd_in_lat(2*i-1)-5
                         end do 
+                        lat_points(ix,iy,iz) = .false.
                        r = matmul(orient, matmul(fcc_mat, vlat))*remesh_lat_pam
                        call add_atom(remesh_type, sbox_atom(atom_index(1)), r)
                     end if