diff --git a/src/mode_create.f90 b/src/mode_create.f90 index ed60daf..546e345 100644 --- a/src/mode_create.f90 +++ b/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,60 +439,66 @@ module mode_create !Now build the finite element region lat_ele_num = 0 lat_atom_num = 0 - allocate(r_lat(3,8,numlatpoints/esize)) - - !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 - end do - - !Now start the element at rzero - do inod=1, 8 - ele(:,inod) = ele(:,inod) + rzero - end do - do iz = -bd_in_array(6), bd_in_array(6) - do iy = -bd_in_array(5), bd_in_array(5) - do ix = -bd_in_array(4), bd_in_array(4) - node_in_bd(:) = .false. - temp_nodes(:,:,:) = 0.0_dp - temp_lat(:,:) = 0 - do inod = 1, 8 - vlat= ele(:,inod) + (/ ix*(esize), iy*(esize), iz*(esize) /) - !Transform point back to real space for easier checking - ! v = matmul(orient, matmul(transform_matrix,v)) - do i = 1,3 - v(i) = real(vlat(i) + bd_in_lat(2*i-1) - 5) - end do - temp_nodes(:,1, inod) = matmul(orient, matmul(transform_matrix, v)) - temp_lat(:,inod) = vlat + curr_esize=esize/(2**(esize_nums-1)) + 1 + allocate(r_lat(3,8,numlatpoints/curr_esize)) - !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 - !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. - end if - end do + 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)/curr_esize + end do - if(all(node_in_bd)) then - lat_ele_num = lat_ele_num+1 - r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:) - - !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 + !Now start the element at rzero + do inod=1, 8 + ele(:,inod) = ele(:,inod) + rzero + end do + do iz = -bd_in_array(6), bd_in_array(6) + do iy = -bd_in_array(5), bd_in_array(5) + do ix = -bd_in_array(4), bd_in_array(4) + node_in_bd(:) = .false. + temp_nodes(:,:,:) = 0.0_dp + temp_lat(:,:) = 0 + do inod = 1, 8 + 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 + v(i) = real(vlat(i) + bd_in_lat(2*i-1) - 5) end do + temp_nodes(:,1, inod) = matmul(orient, matmul(transform_matrix, v)) + temp_lat(:,inod) = vlat + + !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 + !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. + end if end do - end if + if(all(node_in_bd)) then + lat_ele_num = lat_ele_num+1 + r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:) + + !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 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