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, remesh_size
    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, delete

    integer, allocatable :: element_index(:), atom_index(:) 

    public 
    contains

    subroutine group(arg_pos)
        !Main calling function for the group option
        integer, intent(inout) :: arg_pos

        print *, '-----------------------Option Group-------------------------'

        group_ele_num = 0
        group_atom_num = 0
        remesh_size=0
        displace=.false.
        delete=.false.

        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

        if(remesh_size > 0) call remesh_group

        if(delete) call delete_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('remesh')
                arg_pos = arg_pos + 1
                call get_command_argument(arg_pos, textholder, arglen)
                if (arglen==0) stop "Missing remesh element size in group command"
                read(textholder, *) remesh_size
            case('delete')
                arg_pos = arg_pos + 1
                delete=.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')

            print *, "Group has block shape with boundaries: ", block_bd

            !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 

        print *, 'Group contains ', group_ele_num, " elements and ", group_atom_num, " atoms."
    end subroutine get_group 

    subroutine displace_group
        !This subroutine applies a displacement to elements/atoms in the groups

        integer :: i, inod, ibasis

        print *, "Elements/atoms in group displaced by ", disp_vec

        !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

        !If we don't include the wrap command then we have to increase the size of the box
        if (.not.(wrap_flag)) then 
            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

    subroutine remesh_group
        !This command is used to remesh the group to a desired element size

        integer :: i, j, ie, type_interp(max_basisnum*max_esize**3), add_atom_num, orig_atom_num
        real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3)

        !Refining to atoms and remeshing to elements are different processes so check which code we need to run
        select case(remesh_size)

        !Refining to atoms
        case(2)
            if(group_ele_num > 0) then 
                orig_atom_num = atom_num
                !Estimate number of atoms we are adding, this doesn't have to be exact
                add_atom_num = group_ele_num*basisnum(lat_ele(element_index(1)))*size_ele(element_index(1))**3
                call grow_ele_arrays(0,add_atom_num)

                do i = 1, group_ele_num
                    ie = element_index(i)
                    !Get the interpolated atom positions
                    call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp)

                    !Loop over all interpolated atoms and add them to the system, we apply periodic boundaries 
                    !here as well to make sure they are in the box
                    do j = 1, basisnum(lat_ele(ie))*size_ele(ie)**3
                        call apply_periodic(r_interp(:,j))
                        call add_atom(type_interp(j), sbox_ele(ie), r_interp(:,j))
                    end do 
                end do

                !Once all atoms are added we delete all of the elements
                call delete_elements(group_ele_num, element_index)

                print *, group_ele_num, " elements of group are refined to ", atom_num -orig_atom_num, " atoms."

            end if                
        !Remeshing to elements, currently not available
        case default
            print *, "Remeshing to elements is currently not available. Please refine to atoms by passing a remsh size of 2"
        end select

        
    end subroutine remesh_group

    subroutine delete_group
        !This subroutine deletes all atoms/elements within a group

        print *, "Deleting group containing ", group_atom_num, " atoms and ", group_ele_num, " elements."

        !Delete atoms
        call delete_atoms(group_atom_num, atom_index)

        !Delete elements
        call delete_elements(group_ele_num, element_index)

    end subroutine delete_group
end module opt_group