From 66670671b56dd38acf5a3943b136792f11081fe6 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Mon, 18 May 2020 13:50:14 -0400 Subject: [PATCH] Notch group shape --- README.md | 15 +++++++- src/Makefile | 7 ++-- src/opt_group.f90 | 97 ++++++++++++++++++++++++++++++++++++++++++++--- 3 files changed, 108 insertions(+), 11 deletions(-) diff --git a/README.md b/README.md index 3a4630e..2bb3db2 100644 --- a/README.md +++ b/README.md @@ -205,7 +205,8 @@ This selects a group residing in a block with edges perpendicular to the simulat `additional keywords`- Represents the various transformations which can be performed on a group. These additional keywords are given below. -*Wedge* +*Wedge:* + `-group nodes wedge dim1 dim2 bx by bz bw` This selects a group which are within a wedge shape. The options are given as follows: `dim1` - The dimension containing the plane normal of the wedge base. @@ -213,7 +214,17 @@ This selects a group which are within a wedge shape. The options are given as fo `bx by bz` - Centroid of the center of the base `bw` - Base width -**Displace** +*Notch:* + +`-group nodes notch dim1 dim2 bx by bz bw tr` +This shape is similar to a wedge shape except instead of becoming atomically sharp, it finishes in a rounded tip with tip radius `tr`. Options are as follows. +`dim1` - The dimension containing the plane normal of the wedge base. +`dim2` - The thickness dimension. Wedge groups are currently required to span the entire cell thickness in one dimensions which is normal to the triangular face. This through thickness dimension is dim2. +`bx by bz` - Centroid of the center of the base +`bw` - Base width +`tr` - Tip radius + +**Displace:** ``` displace x y z diff --git a/src/Makefile b/src/Makefile index 9a3fced..642fce5 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,9 +1,9 @@ FC=ifort -FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin -heap-arrays -#FFLAGS=-mcmodel=large -Ofast -no-wrap-margin -heap-arrays +#FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin -heap-arrays +FFLAGS=-mcmodel=large -Ofast -no-wrap-margin -heap-arrays MODES=mode_create.o mode_merge.o mode_convert.o OPTIONS=opt_disl.o opt_group.o opt_orient.o opt_delete.o -OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o $(MODES) $(OPTIONS) call_option.o +OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o $(MODES) $(OPTIONS) call_option.o sorts.o .SUFFIXES: .SUFFIXES: .c .f .f90 .F90 .o @@ -38,6 +38,7 @@ atoms.o subroutines.o testfuncs.o box.o : functions.o main.o io.o $(MODES) $(OPTIONS) : elements.o call_mode.o : $(MODES) call_option.o : $(OPTIONS) +elements.o : sorts.o $(MODES) $(OPTIONS) subroutines.o io.o : atoms.o box.o $(MODES) main.o : io.o testfuncs.o elements.o mode_create.o $(OPTIONS) $(MODES): subroutines.o diff --git a/src/opt_group.f90 b/src/opt_group.f90 index dbc984a..dd5a1fb 100644 --- a/src/opt_group.f90 +++ b/src/opt_group.f90 @@ -10,7 +10,7 @@ module opt_group integer :: group_ele_num, group_atom_num, remesh_size,normal, dim1, dim2 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) + real(kind=dp) :: block_bd(6), centroid(3), vertices(3,3),disp_vec(3), tip_radius, bwidth logical :: displace, delete, max_remesh, refine integer, allocatable :: element_index(:), atom_index(:) @@ -56,7 +56,7 @@ module opt_group integer :: i, j, arglen, in_num character(len=100) :: textholder, type_spec - real(kind=dp) bwidth, wheight + real(kind=dp) H !Parse type and shape command arg_pos = arg_pos + 1 @@ -131,6 +131,65 @@ module opt_group end if end do + + case('notch') + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) STOP "Missing normal dim in group notch command" + read(textholder,*) dim1 + + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) STOP "Missing normal dim in group notch command" + read(textholder,*) dim2 + + do i = 1, 3 + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) STOP "Missing centroid in group notch command" + call parse_pos(i, textholder, centroid(i)) + end do + + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) STOP "Missing base width in group notch command" + read(textholder,*) bwidth + + 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 + + !Calculate the vertex positions + vertices(:,1) = centroid + vertices(dim2,1) = 0.0_dp + do i = 1, 3 + if (i == dim1) then + if (bwidth > 0) then + vertices(i,2) = box_bd(2*i) + vertices(i,3) = box_bd(2*i) + H= (box_bd(2*i) - centroid(i)) !Calculate the height of the wedge + else if (bwidth < 0) then + vertices(i,2) = box_bd(2*i-1) + vertices(i,3) = box_bd(2*i-1) + H = (centroid(i) - box_bd(2*i-1)) + else + print *, "bwidth cannot be 0 in wedge shaped group" + stop 3 + end if + else if (i == dim2) then + vertices(i,2) = 0.0_dp + vertices(i,3) = 0.0_dp + else + vertices(i,2) = centroid(i) + 0.5_dp*bwidth + vertices(i,3) = centroid(i) - 0.5_dp*bwidth + end if + end do + + !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 + end if !Read the ID type shape for group case('id') arg_pos = arg_pos + 1 @@ -285,15 +344,15 @@ module opt_group subroutine get_group !This subroutine finds all elements and/or atoms within the group boundaries !specified by the user. - integer :: i, inod, ibasis + integer :: i, j, inod, ibasis integer, allocatable :: resize_array(:) real(kind=dp) :: r_center(3) select case(trim(adjustl(shape))) case('block') print *, "Group has block shape with boundaries: ", block_bd - case ('crack') - print *, "Group has crack shape with dim1", dim1, "and dim2", dim2, "and vertices ", vertices + case ('wedge') + print *, "Group has wedge shape with dim1", dim1, "and dim2", dim2, "and vertices ", vertices case('id') print *, 'Group contains ', group_ele_num, " elements and ", group_atom_num, " atoms." return @@ -343,6 +402,15 @@ module opt_group end if end do end select + + j = 0 + do i = 1, group_atom_num + if (atom_index(i) == 23318348) then + j = j + 1 + end if + end do + + if (j > 1) stop "Code broken" print *, 'Group contains ', group_ele_num, " elements and ", group_atom_num, " atoms." end subroutine get_group @@ -691,12 +759,29 @@ 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 + + integer :: dim3, i logical :: in_group select case(trim(adjustl(shape))) case('block') in_group=in_block_bd(r,block_bd) case('wedge') in_group = in_wedge_bd(r,vertices) - end select + case('notch') + do i = 1, 3 + if (.not.((dim1==i).or.(dim2==i))) dim3 = i + end do + + in_group = in_wedge_bd(r,vertices) + !Do a check to make sure the wedge isn't used if it should be the tip radius + if (bwidth>0) then + 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) + + end select end function in_group end module opt_group