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

Added sphere as a possible group shape

Browse Source
This commit is contained in:
Alex Selimov 2020-06-04 12:13:42 -04:00
parent d01d0b0c57
commit 338587b7b7
2 changed files with 37 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
README.md
View File

@ -218,6 +218,11 @@ This shape is similar to a wedge shape except instead of becoming atomically sha
`bw` - Base width
`tr` - Tip radius
*Sphere*
`-group atoms sphere x y z r`
This shape selects all atoms within a sphere centered at `(x,y,z)` with radius `r`.
**Displace:**
```

40
src/opt_group.f90
View File

@ -10,7 +10,7 @@ module opt_group
integer :: group_ele_num, group_atom_num, remesh_size,normal, dim1, dim2, random_num
character(len=15) :: type, shape !Type indicates what element type is selected and shape is the group shape
real(kind=dp) :: block_bd(6), centroid(3), vertices(3,3),disp_vec(3), tip_radius, bwidth
real(kind=dp) :: block_bd(6), centroid(3), vertices(3,3),disp_vec(3), radius, bwidth
logical :: displace, delete, max_remesh, refine, group_nodes, flip
integer, allocatable :: element_index(:), atom_index(:)
@ -171,7 +171,7 @@ module opt_group
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) STOP "Missing tip radius in group notch command"
read(textholder,*) tip_radius
read(textholder,*) radius
!Calculate the vertex positions
vertices(:,1) = centroid
@ -201,7 +201,7 @@ module opt_group
!Now update the centroid so that the desire tip diameter matches the wedge with
if (bwidth > 0) then
centroid(dim1) = centroid(dim1) + 2*tip_radius*(H)/bwidth
centroid(dim1) = centroid(dim1) + 2*radius*(H)/bwidth
end if
!Read the ID type shape for group
case('id')
@ -313,6 +313,21 @@ module opt_group
if(i ==1) arg_pos = arg_pos + 1
end do
end select
case('sphere')
!First extract the sphere centroid
do i = 1, 3
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) STOP "Missing sphere centroid in group command"
call parse_pos(i, textholder, centroid(i))
end do
!Now get the tip radius
arg_pos=arg_pos+1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) STOP "Missing sphere radius in group command"
read(textholder, *) radius
case default
print *, "Group shape ", trim(adjustl(shape)), " is not currently accepted. Please check documentation ", &
"for accepted group shapes."
@ -376,11 +391,13 @@ module opt_group
case ('wedge')
print *, "Group has wedge shape with dim1", dim1, "and dim2", dim2, "and vertices ", vertices
case ('notch')
print *, "Group has notch shape with dim1", dim1, "and dim2", dim2, " tip radius ", tip_radius, "and vertices ", &
print *, "Group has notch shape with dim1", dim1, "and dim2", dim2, " tip radius ", radius, "and vertices ", &
vertices
case('id')
print *, 'Group contains ', group_ele_num, " elements and ", group_atom_num, " atoms."
return
case('sphere')
print *, "Group has sphere shape with centroid ", centroid(:), " and radius ", radius
end select
!Reset group if needed
@ -837,7 +854,7 @@ module opt_group
function in_group(r)
!This subroutine determines if a point is within the group boundaries
real(kind=dp), intent(in) :: r(3)
real(kind=dp) :: r_notch
real(kind=dp) :: rnorm
integer :: dim3, i
logical :: in_group
@ -857,9 +874,16 @@ module opt_group
if (r(dim1) < centroid(dim1)) in_group = .false.
end if
r_notch = sqrt((r(dim1) - centroid(dim1))**2 + (r(dim3)-centroid(dim3))**2)
in_group = in_group.or.(r_notch < tip_radius)
rnorm = sqrt((r(dim1) - centroid(dim1))**2 + (r(dim3)-centroid(dim3))**2)
in_group = in_group.or.(rnorm < radius)
end select
case('sphere')
rnorm = norm2(r(:) - centroid(:))
if (rnorm <= radius) then
in_group = .true.
else
in_group = .false.
end if
end select
end function in_group
end module opt_group