From 20755270a478c755d0b454513695594dc8e8d02b Mon Sep 17 00:00:00 2001 From: Alex Date: Mon, 24 Feb 2020 12:58:13 -0500 Subject: [PATCH] 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