Merge pull request #12 from aselimov/ft--opt-group

Ft  opt group

README.md

@@ -153,4 +153,23 @@ This option deletes vacancies on a plane which when minimized should result in a
 
 `bx by bz` - The burgers vector for the dislocation
 
-`poisson` - Poisson ratio for continuum solution
+`poisson` - Poisson ratio for continuum solution
+
+### Option Group
+
+`-group select_type group_shape shape_arguments additional keywords`
+
+This option selects a group of either elements, nodes, atoms and applies some transformation to them.
+
+`select_type` - Either `nodes`, `atoms`, `elements`, `nodes/atoms`, `all`. When using the option `nodes` only nodes which are within the group are selected, `elements` selects elements based on whether the element center is within the group, `nodes/atoms` selects both nodes and atoms for the group.  `all` selects elements based on the element center and atoms based on their position.
+
+`group_shape` - Specifies what shape the group takes and dictates which options must be passed. Each shape requires different arguments and these arguments are represented by the placeholder `shape_arguments`. The accepted group shapes and arguments are below:
+
+*Block:*
+
+`-group nodes block xlo xhi ylo yhi zlo zhi`
+
+This selects a group residing in a block with edges perpendicular to the simulation cell. The block boundaries are given by `xlo xhi ylo yhi zlo zhi`. 
+
+`additional keywords`- Represents the various transformations which can be performed on a group. These additional keywords are given below.
+

src/Makefile

@@ -2,8 +2,8 @@ FC=ifort
 FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin
 #FFLAGS=-mcmodel=large -Ofast -no-wrap-margin
 MODES=mode_create.o mode_merge.o mode_convert.o
-OPTIONS=opt_disl.o
+OPTIONS=opt_disl.o opt_group.o
-OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o call_option.o $(MODES) $(OPTIONS)
+OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o  $(MODES) $(OPTIONS) call_option.o
 
 .SUFFIXES:
 .SUFFIXES: .c .f .f90 .F90 .o
@@ -38,6 +38,6 @@ 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)
-$(MODES) io.o: atoms.o box.o
+$(MODES) $(OPTIONS) io.o : atoms.o box.o
 $(MODES) main.o : io.o
-testfuncs.o elements.o mode_create.o opt_disl.o: subroutines.o
+testfuncs.o elements.o mode_create.o $(OPTIONS): subroutines.o

src/call_option.f90

@@ -0,0 +1,18 @@
+subroutine call_option(option, arg_pos)
+    use parameters
+    use opt_disl
+    use opt_group
+    implicit none
+
+    integer, intent(inout) :: arg_pos
+    character(len=100), intent(in) :: option
+
+    select case(trim(adjustl(option)))
+    case('-dislgen', '-disloop')
+        call dislocation(option, arg_pos)
+    case('-group')
+        call group(arg_pos)
+    case default
+        print *, 'Option ', trim(adjustl(option)), ' is not currently accepted.'
+    end select
+end subroutine call_option

src/io.f90

@@ -399,6 +399,8 @@ module io
         select case(max_ng_node)
         case(8)
             interp_max = (max_esize)**3
+        case default
+            interp_max = 0
         end select
         write(11,20) interp_max
         write(11,3) node_num

src/opt_group.f90

@@ -0,0 +1,200 @@
+module opt_group
+
+    !This module contains all code associated with dislocations
+
+    use parameters
+    use elements
+    use subroutines
+    use box 
+    implicit none
+
+    integer :: group_ele_num, group_atom_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), disp_vec(3)
+    logical :: displace, wrap
+
+    integer, allocatable :: element_index(:), atom_index(:) 
+
+    public 
+    contains
+
+    subroutine group(arg_pos)
+        !Main calling function for the group option
+        integer, intent(inout) :: arg_pos
+
+        group_ele_num = 0
+        group_atom_num = 0
+        if(allocated(element_index)) deallocate(element_index)
+        if(allocated(atom_index)) deallocate(atom_index)
+
+        call parse_group(arg_pos)
+
+        call get_group
+
+        !Now call the transformation functions for the group
+        if(displace) call displace_group
+
+    end subroutine group
+
+    subroutine parse_group(arg_pos)
+        !Parse the group command
+        integer, intent(inout) :: arg_pos
+
+        integer :: i,arglen
+        character(len=100) :: textholder
+
+        !Parse type and shape command
+        arg_pos = arg_pos + 1
+        call get_command_argument(arg_pos, type, arglen)
+        if (arglen==0) STOP "Missing select_type in group command"
+        select case(trim(adjustl(type)))
+        case('atoms', 'elements', 'both')
+            continue
+        case default
+            print *, "Select_type ", trim(adjustl(type)), " is not an accept group selection criteria. ", &
+                     "Please select from atoms, nodes, or both."
+        end select
+
+        arg_pos = arg_pos + 1
+        call get_command_argument(arg_pos, shape, arglen)
+        if (arglen==0) STOP "Missing group_shape in group command"
+
+        !Now parse the arguments required by the user selected shape
+        select case(trim(adjustl(shape)))
+        case('block')
+            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"
+                call parse_pos(int((i+1)/2), textholder, block_bd(i))   
+            end do 
+
+        case default
+            print *, "Group shape ", trim(adjustl(shape)), " is not currently accepted. Please check documentation ", &
+                     "for accepted group shapes."
+        end select
+
+        !Now parse the additional options which may be present
+        do while(.true.)
+            if(arg_pos > command_argument_count()) exit
+            !Pull out the next argument which should either be a keyword or an option
+            arg_pos=arg_pos+1
+            call get_command_argument(arg_pos, textholder)
+            textholder=adjustl(textholder)
+
+            !Choose what to based on what the option string is 
+            select case(trim(textholder))
+            case('displace')
+                displace = .true.
+                do i = 1,3 
+                    arg_pos = arg_pos + 1
+                    call get_command_argument(arg_pos, textholder, arglen)
+                    if (arglen==0) stop "Missing vector component for shift command"
+                    read(textholder, *) disp_vec(i)
+                end do
+            case('wrap')
+                wrap = .true.
+            case default
+                !If it isn't an available option to opt_disl then we just exit
+                exit 
+            end select
+        end do
+    end subroutine parse_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, allocatable :: resize_array(:) 
+        real(kind=dp) :: r_center(3)
+
+        select case(trim(adjustl(shape)))
+        case('block')
+            !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
+
+            !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
+
+                        atom_index(group_atom_num) = i
+                    end if
+                end do
+            end select
+        end select 
+    end subroutine get_group 
+
+    subroutine displace_group
+        !This subroutine applies a displacement to elements/atoms in the groups
+
+        integer :: i, inod, ibasis
+
+        !Displace atoms
+        do i = 1, group_atom_num
+            r_atom(:,atom_index(i)) = r_atom(:,atom_index(i)) + disp_vec
+        end do  
+
+        !Displace elements
+        do i = 1, group_ele_num
+            do inod = 1, ng_node(lat_ele(element_index(i)))
+                do ibasis = 1, basisnum(lat_ele(element_index(i)))
+                    r_node(:,ibasis,inod,element_index(i)) = r_node(:,ibasis,inod,element_index(i)) + disp_vec
+                end do 
+            end do 
+        end do
+
+        !Now either apply periodic boundaries if wrap command was passed or adjust box dimensions
+                !Now we check if we have to wrap the atoms, nodes are not wrapped. For elements the periodic 
+        !boundary conditions are applied in the actual CAC codes
+        if(wrap) then   
+            do i = 1, atom_num
+                call apply_periodic(r_atom(:,i))
+            end do
+
+        !If we don't include the wrap command then we have to increase the size of the box
+        else
+            do i = 1,3
+                if (disp_vec(i) < -lim_zero) then 
+                    box_bd(2*i-1) = box_bd(2*i-1) - disp_vec(i)
+                else if (disp_vec(i) > lim_zero) then 
+                    box_bd(2*i) = box_bd(2*i) + disp_vec(i)
+                end if
+            end do
+        end if
+
+    end subroutine displace_group
+
+end module opt_group

src/subroutines.f90

@@ -178,6 +178,8 @@ module subroutines
         real(kind=dp), intent(out) :: pos !The output parsed position value
 
         integer :: iospara
+
+        iospara = 0
         if(trim(adjustl(pos_string)) == 'inf') then 
             pos=box_bd(2*i)
         else if(trim(adjustl(pos_string)) == '-inf') then
@@ -258,12 +260,13 @@ module subroutines
         real(kind=dp), dimension(3), intent(inout) :: r
 
         integer :: j 
-
+        real(kind=dp) ::box_len
         do j = 1, 3
+            box_len = box_bd(2*j) - box_bd(2*j-1)
             if (r(j) > box_bd(2*j)) then 
-                r(j) = r(j) - box_bd(2*j)
+                r(j) = r(j) - box_len
             else if (r(j) < box_bd(2*j-1)) then 
-                r(j) = r(j) + box_bd(2*j-1)
+                r(j) = r(j) + box_len
             end if
         end do
     end subroutine