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..b1bc084 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') @@ -307,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') @@ -355,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 @@ -621,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 @@ -714,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/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..03b3063 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,48 +227,205 @@ 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, 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, & + 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(:,:,:) + 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) = nint(minval(group_in_lat(i,:))) + bd_in_lat(2*i) = nint(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. + dof = 0 - 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 + !Now place all group atoms and group interpolated atoms into lat_points + do i = 1, group_atom_num + 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(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. + dof = dof + 1 + 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 = 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 + 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 + !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 + + + 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 + + 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 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) + 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(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 + end do + 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