aselimov/CACmb
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Fixed the efill option to preserve interelement discontinuities

Browse Source
This commit is contained in:
Alex Selimov 2020-10-15 20:41:15 -04:00
parent c42db27f57
commit dadc1f7a4a
1 changed files with 127 additions and 146 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

273
src/mode_create.f90
View File

@ -15,7 +15,7 @@ module mode_create
orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3), duplicate(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), & 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) 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(:,:) real(kind=dp), allocatable :: r_lat(:,:,:), r_atom_lat(:,:)
integer, allocatable :: elat(:) integer, allocatable :: elat(:)
@ -46,7 +46,7 @@ module mode_create
basisnum = 0 basisnum = 0
lat_ele_num = 0 lat_ele_num = 0
lat_atom_num = 0 lat_atom_num = 0
efill(:) = .false. efill = .false.
!First we parse the command !First we parse the command
call parse_command(arg_pos) call parse_command(arg_pos)
@ -119,8 +119,8 @@ module mode_create
end do end do
end do end do
do i = 1,3 do i = 1,3
box_bd(2*i) = maxval(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 box_bd(2*i-1) = minval(r_node_temp(i,:,:)) - 10.0_dp**(-6.0_dp)
end do end do
call add_element(0,element_type, esize, 1, 1, r_node_temp) call add_element(0,element_type, esize, 1, 1, r_node_temp)
end if end if
@ -265,30 +265,7 @@ module mode_create
end do end do
case('efill') case('efill')
call get_command_argument(arg_pos, textholder) efill=.true.
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 case default
!If it isn't an option then you have to exit !If it isn't an option then you have to exit
arg_pos = arg_pos -1 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 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 real(kind=dp), dimension(3,3), intent(in) :: transform_matrix !The transformation matrix from lattice_space to real space
!Internal variables !Internal variables
integer :: i, inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz, numlatpoints, ele(3,8), rzero(3), & 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, curr_esize, ei 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) real(kind=dp) :: v(3), temp_nodes(3,1,8), r(3), centroid_bd(6)
logical, allocatable :: lat_points(:,:,:) 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 !Do some value initialization
max_esize = esize max_esize = esize
@ -374,19 +352,6 @@ module mode_create
centroid_bd(2*i-1) = huge(1.0_dp) centroid_bd(2*i-1) = huge(1.0_dp)
end do 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) !First find the bounding lattice points (min and max points for the box in each dimension)
numlatpoints = 1 numlatpoints = 1
do i = 1, 3 do i = 1, 3
@ -480,117 +445,134 @@ module mode_create
!Now build the finite element region !Now build the finite element region
lat_ele_num = 0 lat_ele_num = 0
lat_atom_num = 0 lat_atom_num = 0
curr_esize=esize - 2*(esize_nums-1) allocate(r_lat(3,8,numlatpoints/esize), elat(numlatpoints/esize))
allocate(r_lat(3,8,numlatpoints/curr_esize), elat(numlatpoints/curr_esize))
curr_esize=esize ele(:,:) = (esize-1) * cubic_cell(:,:)
do ei = 1, esize_nums !Redefined the second 3 indices as the number of elements that fit in the bounds
ele(:,:) = (curr_esize-1) * cubic_cell(:,:) do i = 1, 3
!Redefined the second 3 indices as the number of elements that fit in the bounds bd_in_array(3+i) = bd_in_array(i)/esize
do i = 1, 3 end do
bd_in_array(3+i) = bd_in_array(i)/curr_esize
end do
!Now start the element at rzero !Now start the element at rzero
do inod=1, 8 do inod=1, 8
ele(:,inod) = ele(:,inod) + rzero ele(:,inod) = ele(:,inod) + rzero
end do end do
do iz = -bd_in_array(6), bd_in_array(6) do iz = -bd_in_array(6), bd_in_array(6)
do iy = -bd_in_array(5), bd_in_array(5) do iy = -bd_in_array(5), bd_in_array(5)
do ix = -bd_in_array(4), bd_in_array(4) do ix = -bd_in_array(4), bd_in_array(4)
node_in_bd(:) = .false. node_in_bd(:) = .false.
temp_nodes(:,:,:) = 0.0_dp temp_nodes(:,:,:) = 0.0_dp
temp_lat(:,:) = 0 temp_lat(:,:) = 0
do inod = 1, 8 do inod = 1, 8
vlat= ele(:,inod) + (/ ix*(curr_esize), iy*(curr_esize), iz*(curr_esize) /) vlat= ele(:,inod) + (/ ix*(esize), iy*(esize), iz*(esize) /)
!Transform point back to real space for easier checking !Transform point back to real space for easier checking
! v = matmul(orient, matmul(transform_matrix,v)) ! v = matmul(orient, matmul(transform_matrix,v))
do i = 1,3 do i = 1,3
v(i) = real(vlat(i) + bd_in_lat(2*i-1) - 5) 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
!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
temp_nodes(:,1, inod) = matmul(orient, matmul(transform_matrix, v)) end do
temp_lat(:,inod) = vlat
!Check to see if the lattice point values are greater then the array limits !If any nodes are in the boundary and we want to efill then run the efill code
if(any(vlat > shape(lat_points)).or.any(vlat < 1)) then else if(any(node_in_bd).and.efill) then
exit
!If within array boundaries check to see if it is a lattice point !Pull out the section of the lat points array
else if(lat_points(vlat(1),vlat(2),vlat(3))) then lat_points_ele(:,:,:)=.false.
node_in_bd(inod) = .true. do i = 1,3
else if (minval(temp_lat(i,:)) <lim_zero) then
exit bd_ele_lat(2*i-1) = 1
else
bd_ele_lat(2*i-1) = minval(temp_lat(i,:))
end if
if(maxval(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 if
end do end do
!If we are on the first round of element building then we can just add the element if all(node_in_bd) is lat_points_ele(1:(bd_ele_lat(2)-bd_ele_lat(1)),1:(bd_ele_lat(4)-bd_ele_lat(3)),&
!true 1:(bd_ele_lat(6)-bd_ele_lat(5)))= lat_points(bd_ele_lat(1):bd_ele_lat(2), &
if(all(node_in_bd).and.(ei==1)) then bd_ele_lat(3):bd_ele_lat(4), &
lat_ele_num = lat_ele_num+1 bd_ele_lat(5):bd_ele_lat(6))
r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:) !Now start looping through elements and try to fit as many as you can
elat(lat_ele_num) = curr_esize efill_size = esize-2
!Now set all the lattice points contained within an element to false i=1
do o = minval(temp_lat(3,:)), maxval(temp_lat(3,:)) j=1
do n = minval(temp_lat(2,:)), maxval(temp_lat(2,:)) k=1
do m = minval(temp_lat(1,:)), maxval(temp_lat(1,:)) nump_ele = count(lat_points_ele)
lat_points(m,n,o) = .false. do i = 1, 3
end do filzero(i) = bd_ele_lat(2*i-1) -1
end do end do
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
!Otherwise we have to also do a box boundary check do o = 1,3
else if(all(node_in_bd)) then ebd(2*o-1) = minval(efill_temp_lat(o,:))
r(:) = 0.0_dp ebd(2*o) = maxval(efill_temp_lat(o,:))
add = .false. end do
do inod = 1,8 lat_ele_num = lat_ele_num+1
r = r+ temp_nodes(:,1,inod)/8.0_dp r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:)
end do elat(lat_ele_num) = efill_size
nump_ele = nump_ele - efill_size**3
!Here we check to make sure the centroid of the element we are adding is outside of the bounds set !Now set all the lattice points contained within an element to false
!by the centroids of the elements of the initial iteration do o = ebd(5), ebd(6)
do n = ebd(3), ebd(4)
do i = 1,3 do m = ebd(1), ebd(2)
if(efill(i)) then lat_points(m+filzero(1),n+filzero(2),o+filzero(3)) = .false.
if((r(i) > centroid_bd(2*i)).or.(r(i) < centroid_bd(2*i-1)))then lat_points_ele(m,n,o) = .false.
add = .true. end do
exit end do
end if end do
end if end do xe
end do end do ye
if(add) then end do ze
lat_ele_num = lat_ele_num+1 efill_size = efill_size-2
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 if end do
end do end if
end do 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 end do
!Now figure out how many lattice points could not be contained in elements !Now figure out how many lattice points could not be contained in elements
allocate(r_atom_lat(3,count(lat_points))) allocate(r_atom_lat(3,count(lat_points)))
@ -641,6 +623,5 @@ module mode_create
STOP 3 STOP 3
end subroutine error_message end subroutine error_message
end module mode_create end module mode_create