diff --git a/README.md b/README.md index fe71ed0..3455000 100644 --- a/README.md +++ b/README.md @@ -203,6 +203,14 @@ 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* +`-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. +`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 + **Displace** ``` @@ -323,4 +331,4 @@ Face 6: [5,6,7,8] Image for vertex numbers is: -![](/home/alexselimov/Documents/CACmb/Numbered_element.png) \ No newline at end of file +![](/home/alexselimov/Documents/CACmb/Numbered_element.png) diff --git a/src/functions.f90 b/src/functions.f90 index 82bcbcd..6fed2b7 100644 --- a/src/functions.f90 +++ b/src/functions.f90 @@ -155,8 +155,7 @@ END FUNCTION StrDnCase end function in_bd_lat function in_block_bd(v, box_bd) - !This function returns whether a point is within a block in 3d - + !This function determines whether a point is within a block in 3d !Input/output real(kind=dp), dimension(3), intent(in) :: v real(kind=dp), dimension(6), intent(in) :: box_bd @@ -180,6 +179,28 @@ END FUNCTION StrDnCase end do end function in_block_bd + function in_wedge_bd(r,vertex) + !This code determines whether the 2dimensional projection of a point lies within the 2 dimensional projection of a wedge. + real(kind=dp), intent(in) :: r(3) !This is the point position + real(kind=dp), intent(in) :: vertex(3,3) !These are the relevant vertex positions for the wedge + real(kind=dp) :: v1(3), v2(3), v3(3), c1(3), c2(3) !Vertex vector to point and cross products + integer :: i + logical :: in_wedge_bd + + in_wedge_bd = .true. + do i = 1, 3 + v1 = vertex(:,mod(i,3)+1) - vertex(:,i) + v2 = r - vertex(:,i) + v3 = vertex(:,mod(i+1,3)+1) - vertex(:,i) + c1 = cross_product(v1,v2) + c2 = cross_product(v1,v3) + if(dot_product(c1,c2) < 0) then + in_wedge_bd=.false. + exit + end if + end do + end function in_wedge_bd + function lcm(a,b) !This function returns the smallest least common multiple of two numbers diff --git a/src/opt_group.f90 b/src/opt_group.f90 index be87611..80bdc26 100644 --- a/src/opt_group.f90 +++ b/src/opt_group.f90 @@ -8,9 +8,9 @@ module opt_group use box implicit none - integer :: group_ele_num, group_atom_num, remesh_size, remesh_type + integer :: group_ele_num, group_atom_num, remesh_size,normal, remesh_type, 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), disp_vec(3), remesh_lat_pam + real(kind=dp) :: block_bd(6), centroid(3), vertices(3,3),disp_vec(3), remesh_lat_pam logical :: displace, delete, max_remesh, refine integer, allocatable :: element_index(:), atom_index(:) @@ -56,6 +56,7 @@ module opt_group integer :: i,arglen character(len=100) :: textholder + real(kind=dp) bwidth, wheight !Parse type and shape command arg_pos = arg_pos + 1 @@ -79,10 +80,56 @@ module opt_group do i= 1, 6 arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) - if (arglen==0) STOP "Missing block boundary in dislgen command" + if (arglen==0) STOP "Missing block boundary in group command" call parse_pos(int((i+1)/2), textholder, block_bd(i)) end do + case('wedge') + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) STOP "Missing normal dim in group wedge 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 wedge 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 wedge 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 wedge command" + read(textholder,*) bwidth + + !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) + else if (bwidth < 0) then + vertices(i,2) = box_bd(2*i-1) + vertices(i,3) = 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) + bwidth + vertices(i,3) = centroid(i) - bwidth + end if + end do case default print *, "Group shape ", trim(adjustl(shape)), " is not currently accepted. Please check documentation ", & "for accepted group shapes." @@ -141,54 +188,55 @@ module opt_group 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 + end select - !Allocate variables to arbitrary size - allocate(element_index(1024), atom_index(1024)) - !Check the type to see whether we need to find the elements within the group - select case(trim(adjustl(type))) - case('elements', 'both') - do i = 1, ele_num - r_center(:) = 0.0_dp - do inod = 1, ng_node(lat_ele(i)) - do ibasis = 1, basisnum(lat_ele(i)) - r_center = r_center + r_node(:,ibasis,inod,i)/(basisnum(lat_ele(i))*ng_node(lat_ele(i))) - end do + !Allocate variables to arbitrary size + allocate(element_index(1024), atom_index(1024)) + !Check the type to see whether we need to find the elements within the group + select case(trim(adjustl(type))) + case('elements', 'both') + do i = 1, ele_num + r_center(:) = 0.0_dp + do inod = 1, ng_node(lat_ele(i)) + do ibasis = 1, basisnum(lat_ele(i)) + r_center = r_center + r_node(:,ibasis,inod,i)/(basisnum(lat_ele(i))*ng_node(lat_ele(i))) end do + end do - if (in_block_bd(r_center, block_bd)) then - group_ele_num = group_ele_num + 1 - if(group_ele_num > size(element_index)) then - allocate(resize_array(size(element_index) + 1024)) - resize_array(1:group_ele_num-1) = element_index - resize_array(group_ele_num:) = 0 - call move_alloc(resize_array, element_index) - end if - - element_index(group_ele_num) = i - end if - end do - end select + if (in_group(r_center)) then + group_ele_num = group_ele_num + 1 + if(group_ele_num > size(element_index)) then + allocate(resize_array(size(element_index) + 1024)) + resize_array(1:group_ele_num-1) = element_index + resize_array(group_ele_num:) = 0 + call move_alloc(resize_array, element_index) + end if - !Check the type to see if we need to find the atoms within the group - select case(trim(adjustl(type))) - case('atoms','both') - do i = 1, atom_num - if(in_block_bd(r_atom(:,i),block_bd)) then - group_atom_num = group_atom_num + 1 - if (group_atom_num > size(atom_index)) then - allocate(resize_array(size(atom_index) + 1024)) - resize_array(1:group_atom_num -1) = atom_index - resize_array(group_atom_num:) = 0 - call move_alloc(resize_array, atom_index) - end if + element_index(group_ele_num) = i + end if + end do + end select - atom_index(group_atom_num) = i + !Check the type to see if we need to find the atoms within the group + select case(trim(adjustl(type))) + case('atoms','both') + do i = 1, atom_num + if(in_group(r_atom(:,i))) then + group_atom_num = group_atom_num + 1 + if (group_atom_num > size(atom_index)) then + allocate(resize_array(size(atom_index) + 1024)) + resize_array(1:group_atom_num -1) = atom_index + resize_array(group_atom_num:) = 0 + call move_alloc(resize_array, atom_index) end if - end do - end select - end select + + atom_index(group_atom_num) = i + end if + end do + end select print *, 'Group contains ', group_ele_num, " elements and ", group_atom_num, " atoms." end subroutine get_group @@ -440,4 +488,16 @@ module opt_group call delete_elements(group_ele_num, element_index) end subroutine delete_group -end module opt_group \ No newline at end of file + + function in_group(r) + !This subroutine determines if a point is within the group boundaries + real(kind=dp), intent(in) :: r(3) + 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 + end function in_group +end module opt_group