From d670bb50830edd9fe1a296de8a96df2105567773 Mon Sep 17 00:00:00 2001 From: Alex Date: Fri, 28 Feb 2020 09:36:33 -0500 Subject: [PATCH] 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)