From 20755270a478c755d0b454513695594dc8e8d02b Mon Sep 17 00:00:00 2001 From: Alex Date: Mon, 24 Feb 2020 12:58:13 -0500 Subject: [PATCH 1/6] Debugging version of code --- README.md | 10 ++- src/io.f90 | 3 + src/mode_create.f90 | 2 +- src/opt_group.f90 | 190 ++++++++++++++++++++++++++++++++++++-------- 4 files changed, 169 insertions(+), 36 deletions(-) diff --git a/README.md b/README.md index 8ad643f..e56e857 100644 --- a/README.md +++ b/README.md @@ -217,7 +217,15 @@ This command wraps atoms back into the simulation cell as though periodic bounda remesh esize ``` -This command remeshes the atoms/elements within the group to the new element size `esize`. Currently only accepts an `esize` of 2 which refines it to full atomistics. +This command remeshes the atoms/elements within the group to the new element size `esize`. Currently only accepts an `esize` of 2 which refines it to full atomistics. When remeshing to atomistics the group can contain any orientations of elements but when remeshing to different finite elements, the group must contain all atoms/elements with the same orientation. + +**Max** + +``` +max +``` + +This command attempts to reduce the degrees of freedom in the model by replacing them with graded elements. This code works by starting at elements with size `esize` and then checks all degrees of freedom to see which ones can be replaced by inserting the element. It then iterates over elements of `esize-2` to `esize=2` which is full atomic resolution. **Delete** diff --git a/src/io.f90 b/src/io.f90 index cb33508..3d03416 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -87,6 +87,9 @@ module io call set_max_esize do i = 1, outfilenum + + print *, "Writing data out to ", trim(adjustl(outfiles(i))) + !Pull out the extension of the file and call the correct write subroutine select case(trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:)))) case('xyz') diff --git a/src/mode_create.f90 b/src/mode_create.f90 index db6fdd0..3cbafa2 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -393,7 +393,7 @@ module mode_create outerloop: do iz = 1, bd_in_array(3) do iy = 1, bd_in_array(2) do ix = 1, bd_in_array(1) - node_in_bd(8) = .false. + node_in_bd(:) = .false. do inod = 1, 8 vlat = ele(:,inod) + (/ ix, iy, iz /) diff --git a/src/opt_group.f90 b/src/opt_group.f90 index 48bd2ed..cfd6499 100644 --- a/src/opt_group.f90 +++ b/src/opt_group.f90 @@ -8,10 +8,10 @@ module opt_group use box implicit none - integer :: group_ele_num, group_atom_num, remesh_size + integer :: group_ele_num, group_atom_num, remesh_size, remesh_type character(len=15) :: type, shape !Type indicates what element type is selected and shape is the group shape - real(kind=dp) :: block_bd(6), disp_vec(3) - logical :: displace, delete + real(kind=dp) :: block_bd(6), disp_vec(3), remesh_lat_pam + logical :: displace, delete, max_remesh, refine integer, allocatable :: element_index(:), atom_index(:) @@ -29,6 +29,8 @@ module opt_group remesh_size=0 displace=.false. delete=.false. + ! max_remesh=.false. + refine = .false. if(allocated(element_index)) deallocate(element_index) if(allocated(atom_index)) deallocate(atom_index) @@ -43,6 +45,9 @@ module opt_group if(remesh_size > 0) call remesh_group if(delete) call delete_group + + if(refine) call refine_group + end subroutine group subroutine parse_group(arg_pos) @@ -101,13 +106,24 @@ module opt_group if (arglen==0) stop "Missing vector component for shift command" read(textholder, *) disp_vec(i) end do + case('refine') + refine=.true. case('remesh') arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) if (arglen==0) stop "Missing remesh element size in group command" read(textholder, *) remesh_size - case('delete') arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) stop "Missing remesh lattice parameter in group command" + read(textholder, *) remesh_lat_pam + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) stop "Missing remesh type in group command" + read(textholder, *) remesh_type + case('max') + max_remesh =.true. + case('delete') delete=.true. case default !If it isn't an available option to opt_disl then we just exit @@ -211,46 +227,152 @@ module opt_group end subroutine displace_group - subroutine remesh_group + subroutine refine_group !This command is used to remesh the group to a desired element size integer :: i, j, ie, type_interp(max_basisnum*max_esize**3), add_atom_num, orig_atom_num real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3) - !Refining to atoms and remeshing to elements are different processes so check which code we need to run - select case(remesh_size) - !Refining to atoms - case(2) - if(group_ele_num > 0) then - orig_atom_num = atom_num - !Estimate number of atoms we are adding, this doesn't have to be exact - add_atom_num = group_ele_num*basisnum(lat_ele(element_index(1)))*size_ele(element_index(1))**3 - call grow_ele_arrays(0,add_atom_num) - - do i = 1, group_ele_num - ie = element_index(i) - !Get the interpolated atom positions - call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp) - - !Loop over all interpolated atoms and add them to the system, we apply periodic boundaries - !here as well to make sure they are in the box - do j = 1, basisnum(lat_ele(ie))*size_ele(ie)**3 - call apply_periodic(r_interp(:,j)) - call add_atom(type_interp(j), sbox_ele(ie), r_interp(:,j)) - end do + + if(group_ele_num > 0) then + orig_atom_num = atom_num + !Estimate number of atoms we are adding, this doesn't have to be exact + add_atom_num = group_ele_num*basisnum(lat_ele(element_index(1)))*size_ele(element_index(1))**3 + call grow_ele_arrays(0,add_atom_num) + do i = 1, group_ele_num + ie = element_index(i) + !Get the interpolated atom positions + call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp) + !Loop over all interpolated atoms and add them to the system, we apply periodic boundaries + !here as well to make sure they are in the box + do j = 1, basisnum(lat_ele(ie))*size_ele(ie)**3 + call apply_periodic(r_interp(:,j)) + call add_atom(type_interp(j), sbox_ele(ie), r_interp(:,j)) + end do + end do + !Once all atoms are added we delete all of the elements + call delete_elements(group_ele_num, element_index) + print *, group_ele_num, " elements of group are refined to ", atom_num -orig_atom_num, " atoms." + end if + + end subroutine + + 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, & + 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_new_node(3,max_basisnum, max_ng_node), orient(3,3), group_in_lat(3,8) + logical, allocatable :: lat_points(:,:,:) + character(len=100) :: remesh_ele_type + + + !Right now we just hardcode only remeshing to elements + remesh_ele_type = 'fcc' + + !Get the orientations, this assumes that the orientation of the subbox for the first atom is the + !same as the rest of the atoms + !If this assumption is false then the code will break and exit + orient = sub_box_ori(:, :, sbox_atom(atom_index(1))) + call matrix_inverse(orient,3,ori_inv) + + !First calculate max position in lattice space to be able to allocate lat_points array, also sum the total numbers of + !degrees of freedom which are added + dof = 0 + + select case(trim(adjustl(shape))) + case('block') + group_in_lat = reshape((/ block_bd(1),block_bd(3),block_bd(5), & + block_bd(2),block_bd(3),block_bd(5), & + block_bd(2),block_bd(4),block_bd(5), & + block_bd(1),block_bd(4),block_bd(5), & + block_bd(1),block_bd(3),block_bd(6), & + block_bd(2),block_bd(3),block_bd(6), & + block_bd(2),block_bd(4),block_bd(6), & + block_bd(1),block_bd(4),block_bd(6) /), [3,8]) + + group_in_lat = matmul(fcc_inv, matmul(ori_inv, group_in_lat/remesh_lat_pam)) + do i = 1, 3 + bd_in_lat(2*i-1) = minval(group_in_lat(i,:)) + bd_in_lat(2*i) = maxval(group_in_lat(i,:)) end do - !Once all atoms are added we delete all of the elements - call delete_elements(group_ele_num, element_index) + end select - print *, group_ele_num, " elements of group are refined to ", atom_num -orig_atom_num, " atoms." + allocate(lat_points(bd_in_lat(2)-bd_in_lat(1)+10, bd_in_lat(4)-bd_in_lat(3)+10, bd_in_lat(6)-bd_in_lat(5)+10)) + lat_points(:,:,:) = .false. + + !Now place all group atoms and group interpolated atoms into lat_points + do i = 1, group_atom_num + r_lat = r_atom(:,atom_index(i))/remesh_lat_pam + r_lat = matmul(fcc_inv,matmul(ori_inv,r_lat)) + !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 + 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) = .true. + end if + end do + + !Now place interpolated atoms within lat_points array + do i =1, group_ele_num + 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 + !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 + stop "Multiple atoms/interpolated atoms share same position in lat point array, this shouldn't happen" + else + 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) = .true. + end if + end do + end do + + !Delete all elements and atoms to make space for new elements and atoms + call delete_atoms(group_atom_num, atom_index) + call delete_elements(group_ele_num, element_index) + + !Now run remeshing algorithm, not the most optimized or efficient but gets the job done + + ele = (remesh_size-1)*cubic_cell + + !Figure out new looping boundaries + bd_in_array(1) = bd_in_lat(2) - bd_in_lat(1) + 10 + bd_in_array(2) = bd_in_lat(4) - bd_in_lat(3) + 10 + bd_in_array(3) = bd_in_lat(6) - bd_in_lat(5) + 10 + + do iz = 1, bd_in_array(3) + do iy = 1, bd_in_array(2) + do ix = 1, bd_in_array(1) + if (lat_points(ix, iy,iz)) then + if (all(lat_points(ix:ix+remesh_size-1,iy:iy+remesh_size-1,iz:iz+remesh_size-1))) then + do inod = 1, 8 + vlat = ele(:,inod) + (/ix, iy, iz /) + do i = 1, 3 + vlat = vlat + bd_in_lat(2*i-1)-5 + end do + r_new_node(:,1,inod) = matmul(orient, matmul(fcc_mat, vlat))*remesh_lat_pam + end do + + lat_points(ix:ix+remesh_size-1,iy:iy+remesh_size-1,iz:iz+remesh_size-1) = .false. + !Add the element, for the sbox we just set it to the same sbox that we get the orientation from. + !In this case it is from the sbox of the first atom in the group. + call add_element(remesh_ele_type, remesh_size, remesh_type, sbox_atom(atom_index(1)),r_new_node) + + end if + end if + end do + end do + end do - end if - !Remeshing to elements, currently not available - case default - print *, "Remeshing to elements is currently not available. Please refine to atoms by passing a remsh size of 2" - end select end subroutine remesh_group From c698c31ede5a5d90ce9ccf62665a4c503a1a439e Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Mon, 24 Feb 2020 14:23:32 -0500 Subject: [PATCH 2/6] First version that doesn't crash, incorrectly places elements and atoms --- src/opt_group.f90 | 51 ++++++++++++++++++++++++++++++++++++----------- 1 file changed, 39 insertions(+), 12 deletions(-) diff --git a/src/opt_group.f90 b/src/opt_group.f90 index cfd6499..8ff789a 100644 --- a/src/opt_group.f90 +++ b/src/opt_group.f90 @@ -263,7 +263,7 @@ module opt_group integer :: i, j, 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), & + 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) logical, allocatable :: lat_points(:,:,:) character(len=100) :: remesh_ele_type @@ -295,8 +295,8 @@ module opt_group group_in_lat = matmul(fcc_inv, matmul(ori_inv, group_in_lat/remesh_lat_pam)) do i = 1, 3 - bd_in_lat(2*i-1) = minval(group_in_lat(i,:)) - bd_in_lat(2*i) = maxval(group_in_lat(i,:)) + bd_in_lat(2*i-1) = nint(minval(group_in_lat(i,:))) + bd_in_lat(2*i) = nint(maxval(group_in_lat(i,:))) end do end select @@ -306,8 +306,11 @@ module opt_group !Now place all group atoms and group interpolated atoms into lat_points do i = 1, group_atom_num - r_lat = r_atom(:,atom_index(i))/remesh_lat_pam - r_lat = matmul(fcc_inv,matmul(ori_inv,r_lat)) + r = r_atom(:,atom_index(i))/remesh_lat_pam + r = matmul(fcc_inv,matmul(ori_inv,r)) + do j = 1, 3 + r_lat(j) = nint(r(j)) + 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 @@ -349,10 +352,20 @@ module opt_group bd_in_array(2) = bd_in_lat(4) - bd_in_lat(3) + 10 bd_in_array(3) = bd_in_lat(6) - bd_in_lat(5) + 10 - do iz = 1, bd_in_array(3) - do iy = 1, bd_in_array(2) - do ix = 1, bd_in_array(1) + zloop: do iz = 1, bd_in_array(3) + yloop: do iy = 1, bd_in_array(2) + xloop: do ix = 1, bd_in_array(1) if (lat_points(ix, iy,iz)) then + + !Check to see if the element overshoots the bound + if (iz+remesh_size-1 > bd_in_array(3)) then + exit zloop + else if (iy+remesh_size-1 > bd_in_array(2)) then + cycle zloop + else if (ix+remesh_size-1 > bd_in_array(1)) then + cycle yloop + end if + if (all(lat_points(ix:ix+remesh_size-1,iy:iy+remesh_size-1,iz:iz+remesh_size-1))) then do inod = 1, 8 vlat = ele(:,inod) + (/ix, iy, iz /) @@ -369,11 +382,25 @@ module opt_group end if end if - end do - end do - end do - + end do xloop + end do yloop + end do zloop + !Now we have to add any leftover lattice points as atoms + do iz = 1, bd_in_array(3) + do iy=1, bd_in_array(2) + do ix = 1, bd_in_array(1) + if(lat_points(ix,iy,iz)) then + vlat = (/ ix, iy, iz /) + do i = 1, 3 + vlat = vlat + bd_in_lat(2*i-1)-5 + end do + r = matmul(orient, matmul(fcc_mat, vlat))*remesh_lat_pam + call add_atom(remesh_type, sbox_atom(atom_index(1)), r) + end if + end do + end do + end do end subroutine remesh_group From 56994393a07c4c58d53354e0f646e9cafcf58b24 Mon Sep 17 00:00:00 2001 From: Alex Date: Tue, 25 Feb 2020 10:11:12 -0500 Subject: [PATCH 3/6] Fix to bug causing incorrrect positions. Remeshing algorithm works --- src/opt_group.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/opt_group.f90 b/src/opt_group.f90 index 8ff789a..2176a88 100644 --- a/src/opt_group.f90 +++ b/src/opt_group.f90 @@ -370,7 +370,7 @@ module opt_group do inod = 1, 8 vlat = ele(:,inod) + (/ix, iy, iz /) do i = 1, 3 - vlat = vlat + bd_in_lat(2*i-1)-5 + vlat(i) = vlat(i) + bd_in_lat(2*i-1)-5 end do r_new_node(:,1,inod) = matmul(orient, matmul(fcc_mat, vlat))*remesh_lat_pam end do @@ -393,7 +393,7 @@ module opt_group if(lat_points(ix,iy,iz)) then vlat = (/ ix, iy, iz /) do i = 1, 3 - vlat = vlat + bd_in_lat(2*i-1)-5 + vlat(i) = vlat(i) + bd_in_lat(2*i-1)-5 end do r = matmul(orient, matmul(fcc_mat, vlat))*remesh_lat_pam call add_atom(remesh_type, sbox_atom(atom_index(1)), r) From 09c2e6315500256d7b32a21eb9b221590910562f Mon Sep 17 00:00:00 2001 From: Alex Date: Tue, 25 Feb 2020 10:34:09 -0500 Subject: [PATCH 4/6] Working element refinement --- src/opt_group.f90 | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/opt_group.f90 b/src/opt_group.f90 index 2176a88..3b40f14 100644 --- a/src/opt_group.f90 +++ b/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 From 5b925122df056f81cc89bec839bff845516ad297 Mon Sep 17 00:00:00 2001 From: Alex Date: Tue, 25 Feb 2020 10:47:07 -0500 Subject: [PATCH 5/6] Added verbosity to remesh command --- src/opt_group.f90 | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/src/opt_group.f90 b/src/opt_group.f90 index 3b40f14..59270b5 100644 --- a/src/opt_group.f90 +++ b/src/opt_group.f90 @@ -262,7 +262,7 @@ module opt_group !This command is used to remesh the group to a desired element size 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) + current_esize, dof, max_lat(3), r_lat(3), ele(3,8), vlat(3), bd_in_lat(6), bd_in_array(3), old_ele, old_atom 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) logical, allocatable :: lat_points(:,:,:) @@ -303,6 +303,7 @@ module opt_group allocate(lat_points(bd_in_lat(2)-bd_in_lat(1)+10, bd_in_lat(4)-bd_in_lat(3)+10, bd_in_lat(6)-bd_in_lat(5)+10)) lat_points(:,:,:) = .false. + dof = 0 !Now place all group atoms and group interpolated atoms into lat_points do i = 1, group_atom_num @@ -317,6 +318,7 @@ module opt_group stop "Multiple atoms share same position in lat point array, this shouldn't happen" else 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) = .true. + dof = dof + 1 end if end do @@ -336,14 +338,20 @@ module opt_group stop "Multiple atoms/interpolated atoms share same position in lat point array, this shouldn't happen" else 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) = .true. + dof = dof + 1 end if end do end do + print *, "Group has ", dof, " degrees of freedom to remesh" + !Delete all elements and atoms to make space for new elements and atoms call delete_atoms(group_atom_num, atom_index) call delete_elements(group_ele_num, element_index) + old_atom = atom_num + old_ele = ele_num + !Now run remeshing algorithm, not the most optimized or efficient but gets the job done ele = (remesh_size-1)*cubic_cell @@ -404,6 +412,8 @@ module opt_group end do end do + + print *, "Remeshing has created ", ele_num-old_ele, " elements and ", atom_num-old_atom, " atoms." end subroutine remesh_group subroutine delete_group From d670bb50830edd9fe1a296de8a96df2105567773 Mon Sep 17 00:00:00 2001 From: Alex Date: Fri, 28 Feb 2020 09:36:33 -0500 Subject: [PATCH 6/6] Working remesh command and remesh max command --- src/io.f90 | 17 +++++++++-- src/opt_group.f90 | 78 +++++++++++++++++++++++++++-------------------- 2 files changed, 60 insertions(+), 35 deletions(-) diff --git a/src/io.f90 b/src/io.f90 index 3d03416..b1bc084 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -310,6 +310,9 @@ module io 20 format('SCALARS lattice_type float', /& 'LOOKUP_TABLE default') +21 format('SCALARS esize float', /& + 'LOOKUP_TABLE default') + !First we write the vtk file containing the atoms open(unit=11, file='atoms_'//trim(adjustl(file)), action='write', status='replace',position='rewind') @@ -358,6 +361,10 @@ module io do i = 1, ele_num write(11, '(i16)') lat_ele(i) end do + write(11,21) + do i = 1, ele_num + write(11, '(i16)') size_ele(i) + end do close(11) end subroutine @@ -624,8 +631,13 @@ module io !Read in the box boundary and grow the current active box bd read(11, *) temp_box_bd(:) + print *, "displace", displace do i = 1, 3 - newdisplace(i) = displace(i) - temp_box_bd(2*i-1) + if (displace(i) > lim_zero) then + newdisplace(i) = displace(i) - temp_box_bd(2*i-1) + else + newdisplace=displace(i) + end if temp_box_bd(2*i-1) = temp_box_bd(2*i-1) + newdisplace(i) temp_box_bd(2*i) = temp_box_bd(2*i) + newdisplace(i) end do @@ -717,7 +729,8 @@ module io !Read the atoms do i = 1, in_atoms read(11,*) j, type, sbox, r(:) - call add_atom(new_type_to_type(type), sbox, r+newdisplace ) + r = r+newdisplace + call add_atom(new_type_to_type(type), sbox, r) end do !Read the elements diff --git a/src/opt_group.f90 b/src/opt_group.f90 index 59270b5..03b3063 100644 --- a/src/opt_group.f90 +++ b/src/opt_group.f90 @@ -29,7 +29,7 @@ module opt_group remesh_size=0 displace=.false. delete=.false. - ! max_remesh=.false. + max_remesh=.false. refine = .false. if(allocated(element_index)) deallocate(element_index) @@ -262,7 +262,9 @@ module opt_group !This command is used to remesh the group to a desired element size 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), old_ele, old_atom + current_esize, dof, max_lat(3), r_lat(3), ele(3,8), vlat(3), bd_in_lat(6), bd_in_array(3), old_ele, old_atom, & + max_loops, working_esize + 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) logical, allocatable :: lat_points(:,:,:) @@ -352,48 +354,58 @@ module opt_group old_atom = atom_num old_ele = ele_num - !Now run remeshing algorithm, not the most optimized or efficient but gets the job done - ele = (remesh_size-1)*cubic_cell + !Now run remeshing algorithm, not the most optimized or efficient but gets the job done !Figure out new looping boundaries bd_in_array(1) = bd_in_lat(2) - bd_in_lat(1) + 10 bd_in_array(2) = bd_in_lat(4) - bd_in_lat(3) + 10 bd_in_array(3) = bd_in_lat(6) - bd_in_lat(5) + 10 - zloop: do iz = 1, bd_in_array(3) - yloop: do iy = 1, bd_in_array(2) - xloop: do ix = 1, bd_in_array(1) - if (lat_points(ix, iy,iz)) then - - !Check to see if the element overshoots the bound - if (iz+remesh_size-1 > bd_in_array(3)) then - exit zloop - else if (iy+remesh_size-1 > bd_in_array(2)) then - cycle zloop - else if (ix+remesh_size-1 > bd_in_array(1)) then - cycle yloop - end if - - if (all(lat_points(ix:ix+remesh_size-1,iy:iy+remesh_size-1,iz:iz+remesh_size-1))) then - do inod = 1, 8 - vlat = ele(:,inod) + (/ix, iy, iz /) - do i = 1, 3 - vlat(i) = vlat(i) + bd_in_lat(2*i-1)-5 + + if (max_remesh) then + max_loops = (remesh_size-2)/2 + else + max_loops = 1 + end if + do j = 1, max_loops + working_esize = remesh_size - 2*(j-1) + ele = (working_esize-1)*cubic_cell + zloop: do iz = 1, bd_in_array(3) + yloop: do iy = 1, bd_in_array(2) + xloop: do ix = 1, bd_in_array(1) + if (lat_points(ix, iy,iz)) then + r_new_node(:,:,:) = 0.0_dp + + !Check to see if the element overshoots the bound + if (iz+working_esize-1 > bd_in_array(3)) then + exit zloop + else if (iy+working_esize-1 > bd_in_array(2)) then + cycle zloop + else if (ix+working_esize-1 > bd_in_array(1)) then + cycle yloop + end if + + if (all(lat_points(ix:ix+working_esize-1,iy:iy+working_esize-1,iz:iz+working_esize-1))) then + do inod = 1, ng_node(remesh_type) + vlat = ele(:,inod) + (/ix, iy, iz /) + do i = 1, 3 + vlat(i) = vlat(i) + bd_in_lat(2*i-1)-5 + end do + r_new_node(:,1,inod) = matmul(orient, matmul(fcc_mat, vlat))*remesh_lat_pam end do - r_new_node(:,1,inod) = matmul(orient, matmul(fcc_mat, vlat))*remesh_lat_pam - end do - lat_points(ix:ix+remesh_size-1,iy:iy+remesh_size-1,iz:iz+remesh_size-1) = .false. - !Add the element, for the sbox we just set it to the same sbox that we get the orientation from. - !In this case it is from the sbox of the first atom in the group. - call add_element(remesh_ele_type, remesh_size, remesh_type, sbox_atom(atom_index(1)),r_new_node) + lat_points(ix:ix+working_esize-1,iy:iy+working_esize-1,iz:iz+working_esize-1) = .false. + !Add the element, for the sbox we just set it to the same sbox that we get the orientation from. + !In this case it is from the sbox of the first atom in the group. + call add_element(remesh_ele_type, working_esize, remesh_type, sbox_atom(atom_index(1)),r_new_node) + end if end if - end if - end do xloop - end do yloop - end do zloop + end do xloop + end do yloop + end do zloop + end do !Now we have to add any leftover lattice points as atoms do iz = 1, bd_in_array(3)