diff --git a/src/mode_create.f90 b/src/mode_create.f90 index e9d7fb2..4a69533 100644 --- a/src/mode_create.f90 +++ b/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)))