From dadc1f7a4a332462c1f17acada7c9090370ad62e Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Thu, 15 Oct 2020 20:41:15 -0400 Subject: [PATCH] Fixed the efill option to preserve interelement discontinuities --- src/mode_create.f90 | 273 +++++++++++++++++++++----------------------- 1 file changed, 127 insertions(+), 146 deletions(-) diff --git a/src/mode_create.f90 b/src/mode_create.f90 index 4a69533..bcec2d5 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -15,7 +15,7 @@ module mode_create 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(3) + logical :: dup_flag, dim_flag, efill real(kind=dp), allocatable :: r_lat(:,:,:), r_atom_lat(:,:) integer, allocatable :: elat(:) @@ -46,7 +46,7 @@ module mode_create basisnum = 0 lat_ele_num = 0 lat_atom_num = 0 - efill(:) = .false. + efill = .false. !First we parse the command call parse_command(arg_pos) @@ -119,8 +119,8 @@ module mode_create end do end do do i = 1,3 - box_bd(2*i) = maxval(r_node_temp(i,:,:))+10.0_dp**-6.0_dp - box_bd(2*i-1) = minval(r_node_temp(i,:,:)) - 10.0_dp**-6.0_dp + box_bd(2*i) = maxval(r_node_temp(i,:,:))+10.0_dp**(-6.0_dp) + box_bd(2*i-1) = minval(r_node_temp(i,:,:)) - 10.0_dp**(-6.0_dp) end do call add_element(0,element_type, esize, 1, 1, r_node_temp) end if @@ -265,30 +265,7 @@ module mode_create end do case('efill') - 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 + efill=.true. case default !If it isn't an option then you have to exit arg_pos = arg_pos -1 @@ -361,11 +338,12 @@ module mode_create integer, dimension(3,8), intent(in) :: box_in_lat !The box vertices transformed to lattice space 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, curr_esize, ei + integer :: i, inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz, numlatpoints, ele(3,8), rzero(3), efill_size, & + vlat(3), temp_lat(3,8), m, n, o, j, k, nump_ele, efill_temp_lat(3,8), filzero(3), bd_ele_lat(6),& + efill_ele(3,8), ebd(6) real(kind=dp) :: v(3), temp_nodes(3,1,8), r(3), centroid_bd(6) logical, allocatable :: lat_points(:,:,:) - logical :: node_in_bd(8), add + logical :: node_in_bd(8), add, lat_points_ele(esize,esize,esize), intersect_bd(3) !Do some value initialization max_esize = esize @@ -374,19 +352,6 @@ module mode_create 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(any(efill)) then - curr_esize=esize - esize_nums=0 - do while (curr_esize >= 7) - esize_nums=esize_nums+1 - curr_esize = curr_esize -2 - 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 @@ -480,117 +445,134 @@ module mode_create !Now build the finite element region lat_ele_num = 0 lat_atom_num = 0 - curr_esize=esize - 2*(esize_nums-1) - allocate(r_lat(3,8,numlatpoints/curr_esize), elat(numlatpoints/curr_esize)) + allocate(r_lat(3,8,numlatpoints/esize), elat(numlatpoints/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)/curr_esize - end do + ele(:,:) = (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 + 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 - 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 + !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 - !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 - !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 + !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 - !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 + !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)) then + lat_ele_num = lat_ele_num+1 + r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:) + elat(lat_ele_num) = 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 - !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 + !If any nodes are in the boundary and we want to efill then run the efill code + else if(any(node_in_bd).and.efill) then - 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 + !Pull out the section of the lat points array + lat_points_ele(:,:,:)=.false. + do i = 1,3 + if (minval(temp_lat(i,:)) size(lat_points,i)) then + bd_ele_lat(2*i) = size(temp_lat(i,:)) + else + bd_ele_lat(2*i) = maxval(temp_lat(i,:)) + end if + end do + + lat_points_ele(1:(bd_ele_lat(2)-bd_ele_lat(1)),1:(bd_ele_lat(4)-bd_ele_lat(3)),& + 1:(bd_ele_lat(6)-bd_ele_lat(5)))= lat_points(bd_ele_lat(1):bd_ele_lat(2), & + bd_ele_lat(3):bd_ele_lat(4), & + bd_ele_lat(5):bd_ele_lat(6)) + !Now start looping through elements and try to fit as many as you can + efill_size = esize-2 + i=1 + j=1 + k=1 + nump_ele = count(lat_points_ele) + do i = 1, 3 + filzero(i) = bd_ele_lat(2*i-1) -1 + end do + do while(efill_size>9) + !First check whether there are enough lattice points to house the current element size + efill_ele=cubic_cell*(efill_size-1) + if (nump_ele < efill_size**3) then + efill_size = efill_size - 2 + else + ze: do k = 1, (esize-efill_size) + ye: do j = 1, (esize-efill_size) + xe: do i = 1, (esize-efill_size) + do inod = 1,8 + vlat = efill_ele(:,inod) + (/ i, j, k /) + if (.not.lat_points_ele(vlat(1),vlat(2),vlat(3))) cycle xe + do o = 1,3 + v(o) = real(vlat(o) + filzero(o) + bd_in_lat(2*o-1) -5) + end do + temp_nodes(:,1, inod) = matmul(orient, matmul(transform_matrix, v)) + efill_temp_lat(:,inod) = vlat + end do + + do o = 1,3 + ebd(2*o-1) = minval(efill_temp_lat(o,:)) + ebd(2*o) = maxval(efill_temp_lat(o,:)) + end do + lat_ele_num = lat_ele_num+1 + r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:) + elat(lat_ele_num) = efill_size + nump_ele = nump_ele - efill_size**3 + !Now set all the lattice points contained within an element to false + do o = ebd(5), ebd(6) + do n = ebd(3), ebd(4) + do m = ebd(1), ebd(2) + lat_points(m+filzero(1),n+filzero(2),o+filzero(3)) = .false. + lat_points_ele(m,n,o) = .false. + end do + end do + end do + end do xe + end do ye + end do ze + efill_size = efill_size-2 + end if + end do + end if 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))) @@ -641,6 +623,5 @@ module mode_create STOP 3 end subroutine error_message - - + end module mode_create