Merge branch 'ft--option-slip-plane' into development

src/Makefile

@@ -1,9 +1,15 @@
-FC=gfortran
+FC=ifort
-#FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin -heap-arrays
+
+#Ifort flags
+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 -O3 -g
+
+#gfortran flags
+#FFLAGS=-mcmodel=large -O3 -g
+#FFLAGS=-mcmodel=large -O0 -g -fbacktrace -fcheck=all
+
 MODES=mode_create.o mode_merge.o mode_convert.o
-OPTIONS=opt_disl.o opt_group.o opt_orient.o opt_delete.o opt_deform.o opt_redef_box.o
+OPTIONS=opt_disl.o opt_group.o opt_orient.o opt_delete.o opt_deform.o opt_redef_box.o opt_slip_plane.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:
@@ -17,7 +23,7 @@ cacmb: $(OBJECTS)
 
 .PHONY: clean
 clean: 
-	$(RM) cacmb  *.o 
+	$(RM) cacmb  *.o *.mod
 
 testfuncs: testfuncs.o functions.o subroutines.o
 	$(FC) testfuncs.o functions.o subroutines.o box.o elements.o -o $@

src/call_option.f90

@@ -6,6 +6,7 @@ subroutine call_option(option, arg_pos)
     use opt_deform
     use opt_delete
     use opt_redef_box
+    use opt_slip_plane
     use box
     implicit none
 
@@ -41,6 +42,8 @@ subroutine call_option(option, arg_pos)
         arg_pos=arg_pos +3
     case('-redef_box')
         call redef_box(arg_pos)
+    case('-slip_plane')
+        call run_slip_plane(arg_pos)
     case default
         print *, 'Option ', trim(adjustl(option)), ' is not currently accepted.'
         stop 3 

src/elements.f90

@@ -201,38 +201,39 @@ module elements
         
         !First check to make sure if it is allocated
         if (allocated(size_ele)) then 
+
             !Figure out the size of the atom and element arrays
             ele_size = size(size_ele)
 
             !Check if we need to grow the ele_size, if so grow all the variables
             if ( n+ele_num > size(size_ele)) then 
                 
-                allocate(temp_int(n+ele_num+buffer_size))
+                allocate(temp_int(n+ele_size+buffer_size))
-                temp_int(1:ele_size) = lat_ele
+                temp_int(1:ele_size) = lat_ele(1:ele_size)
                 temp_int(ele_size+1:) = 0
                 call move_alloc(temp_int, lat_ele)
 
-                allocate(temp_int(n+ele_num+buffer_size))
+                allocate(temp_int(n+ele_size+buffer_size))
-                temp_int(1:ele_size) = tag_ele
+                temp_int(1:ele_size) = tag_ele(1:ele_size)
                 temp_int(ele_size+1:) = 0
                 call move_alloc(temp_int, tag_ele)
 
-                allocate(temp_int(n+ele_num+buffer_size))
+                allocate(temp_int(n+ele_size+buffer_size))
-                temp_int(1:ele_size) = size_ele
+                temp_int(1:ele_size) = size_ele(1:ele_size)
                 temp_int(ele_size+1:) = 0
                 call move_alloc(temp_int, size_ele)
 
-                allocate(temp_int(n+ele_num+buffer_size))
+                allocate(temp_int(n+ele_size+buffer_size))
-                temp_int(1:ele_size) = lat_ele
+                temp_int(1:ele_size) = sbox_ele(1:ele_size)
                 temp_int(ele_size+1:) = 0
                 call move_alloc(temp_int, sbox_ele)
 
-                allocate(char_temp(n+ele_num+buffer_size))
+                allocate(char_temp(n+ele_size+buffer_size))
-                char_temp(1:ele_size) = type_ele
+                char_temp(1:ele_size) = type_ele(1:ele_size)
                 call move_alloc(char_temp, type_ele)
 
-                allocate(temp_ele_real(3, max_basisnum, max_ng_node, n+ele_num+buffer_size))
+                allocate(temp_ele_real(3, max_basisnum, max_ng_node, n+ele_size+buffer_size))
-                temp_ele_real(:,:,:,1:ele_size) = r_node
+                temp_ele_real(:,:,:,1:ele_size) = r_node(:,:,:,1:ele_size)
                 temp_ele_real(:,:,:,ele_size+1:) = 0.0_dp
                 call move_alloc(temp_ele_real, r_node)
             end if
@@ -244,22 +245,22 @@ module elements
         if (allocated(type_atom)) then 
             atom_size = size(type_atom)
             if (m+atom_num > atom_size) then 
-                allocate(temp_int(m+atom_num+buffer_size))
+                allocate(temp_int(m+atom_size+buffer_size))
                 temp_int(1:atom_size) = type_atom
                 temp_int(atom_size+1:) = 0
                 call move_alloc(temp_int, type_atom)
 
-                allocate(temp_int(m+atom_num+buffer_size))
+                allocate(temp_int(m+atom_size+buffer_size))
                 temp_int(1:atom_size) = tag_atom
                 temp_int(atom_size+1:) = 0
                 call move_alloc(temp_int, tag_atom)
 
-                allocate(temp_int(m+atom_num+buffer_size))
+                allocate(temp_int(m+atom_size+buffer_size))
                 temp_int(1:atom_size) = sbox_atom
                 temp_int(atom_size+1:) = 0
                 call move_alloc(temp_int, sbox_atom)
 
-                allocate(temp_real(3,m+atom_num+buffer_size))
+                allocate(temp_real(3,m+atom_size+buffer_size))
                 temp_real(:,1:atom_size) = r_atom
                 temp_real(:, atom_size+1:) = 0.0_dp
                 call move_alloc(temp_real, r_atom)
@@ -278,6 +279,7 @@ module elements
         integer :: newtag
 
         ele_num = ele_num + 1
+        node_num = node_num + ng_node(lat)
 
         if (tag==0) then
             newtag = ele_num !If we don't assign a tag then pass the tag as the ele_num
@@ -293,7 +295,6 @@ module elements
         lat_ele(ele_num) = lat
         sbox_ele(ele_num) = sbox
         r_node(:,:,:,ele_num) = r(:,:,:)
-       node_num = node_num + ng_node(lat)
 
 
     end subroutine add_element
@@ -750,4 +751,28 @@ module elements
 
     end subroutine lattice_map
 
+    subroutine get_interp_pos(i,j,k, ie, rout)
+        !This returns the position of an interpolated basis from an element ie.
+        !i, j, k should be in natural coordinates
+
+        integer, intent(in) :: i, j, k
+        real(kind=dp), dimension(3,max_basisnum), intent(out) :: rout
+
+        integer :: ie, ibasis, inod
+        real(kind=dp) :: a_shape(8), r, s, t
+
+        r = (1.0_dp*(i-1)-(size_ele(ie)-1)/2)/(1.0_dp*(size_ele(ie)-1)/2)
+        s = (1.0_dp*(j-1)-(size_ele(ie)-1)/2)/(1.0_dp*(size_ele(ie)-1)/2)
+        t = (1.0_dp*(k-1)-(size_ele(ie)-1)/2)/(1.0_dp*(size_ele(ie)-1)/2)
+        rout(:,:) = 0
+        do ibasis = 1, basisnum(lat_ele(ie))
+            do inod = 1, ng_node(lat_ele(ie)) 
+                call rhombshape(r,s,t,a_shape)
+                rout(:,ibasis) = rout(:,ibasis) + a_shape(inod) * r_node(:,ibasis,inod,ie)
+            end do
+        end do
+
+
+    end subroutine 
+
 end module elements

src/mode_create.f90

@@ -528,7 +528,7 @@ module mode_create
                             do i = 1, 3
                                 filzero(i) = bd_ele_lat(2*i-1) -1
                             end do
-                            do while(efill_size>9)
+                            do while(efill_size>min_efillsize)
                                 !First check whether there are enough lattice points to house the current element size
                                 efill_ele=cubic_cell*(efill_size-1)
                                 if  (nump_ele < efill_size**3) then 

src/opt_slip_plane.f90

@@ -0,0 +1,176 @@
+module opt_slip_plane
+    use parameters 
+    use elements
+    use functions 
+    use subroutines
+
+    implicit none
+
+    integer :: sdim
+    real(kind=dp) :: spos
+    logical :: efill
+    
+    public 
+    contains
+
+    subroutine run_slip_plane(arg_pos)
+        !Main calling function for the slip_plane option
+        integer, intent(inout) :: arg_pos 
+
+        integer :: ie, ia, slip_enum, old_atom_num, esize, new_ele_num, n, m, o, ele(3,8), nump_ele, inod, vlat(3), ibasis
+                
+        integer, allocatable :: slip_eles(:), temp_int(:)
+        real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), rfill(3,max_basisnum, max_ng_node), ratom(3,max_basisnum), &
+                         maxp, minp
+
+        integer :: type_interp(max_basisnum*max_esize**3)
+        logical :: lat_points(max_esize,max_esize, max_esize)
+    
+
+        print *, '---------------------Option Slip_Plane----------------------'
+
+        !Initialize variables
+        efill = .false.
+        slip_enum = 0
+        old_atom_num = atom_num
+
+        !!Parse the argument
+        call parse(arg_pos)
+
+
+        !If we are running the efill code then we have to initiate some variables
+        if(efill) then 
+            new_ele_num = 0
+        end if
+        allocate(slip_eles(1024))
+        !Now loop over all elements, find which ones intersect 
+        do ie = 1, ele_num
+            if( (spos < maxval(r_node(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie)),ie))).and. &
+                (spos > minval(r_node(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie)),ie)))) then 
+                slip_enum = slip_enum + 1
+                if (slip_enum > size(slip_eles)) then 
+                    allocate(temp_int(size(slip_eles)+1024))
+                    temp_int(1:size(slip_eles)) = slip_eles
+                    temp_int(size(slip_eles)+1:)  = 0
+                    call move_alloc(temp_int, slip_eles)
+                end if
+                slip_eles(slip_enum) = ie
+                
+                !If we aren't efilling then just refine the element
+                if(.not.efill) then 
+                    call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp)
+                    do ia = 1, basisnum(lat_ele(ie)) * size_ele(ie)**3           
+                        call apply_periodic(r_interp(:,ia))
+                        call add_atom(0, type_interp(ia), sbox_ele(ie), r_interp(:,ia))
+                    end do
+                !If we are efilling then the code is slightly more complex
+                else
+                    !First populate the lat points array 
+                    lat_points(1:size_ele(ie),1:size_ele(ie), 1:size_ele(ie)) = .true.
+
+                    !Now start trying to remesh the region, leaving the slip plane as a discontinuity
+                    esize = size_ele(ie) - 2
+                    nump_ele = size_ele(ie)**3
+                    do while(esize > min_efillsize)
+                        if(nump_ele < esize**3) then 
+                            esize = esize - 2
+                        else
+                            ele = cubic_cell*(esize -1) 
+                            do o = 1, size_ele(ie) - esize
+                                do n = 1, size_ele(ie) - esize
+                                    latloop:do m = 1, size_ele(ie) - esize
+                                        do inod = 1, ng_node(lat_ele(ie))
+                                            vlat = ele(:,inod) + (/ m, n, o /)
+                                            if (.not.lat_points(vlat(1), vlat(2),vlat(3))) cycle latloop
+                                            call get_interp_pos(vlat(1), vlat(2), vlat(3), ie, rfill(:,:,inod))
+                                        end do
+
+                                        !Check to make sure all lattice points exist for the current element
+                                        if(any(.not.lat_points(m:m+esize-1, n:n+esize-1, o:o+esize-1))) cycle latloop 
+                                        !Check to see if the plane intersects this element if not then add it 
+                                        maxp = maxval(rfill(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie))))
+                                        minp = minval(rfill(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie))))
+                                        if(.not.(spos < maxp).and.(spos > minp))then
+                                            nump_ele = nump_ele - esize**3 
+                                            lat_points(m:m+esize-1, n:n+esize-1, o:o+esize-1) = .false.
+                                            call add_element(0, type_ele(ie), esize, lat_ele(ie), sbox_ele(ie), rfill) 
+                                            new_ele_num = new_ele_num + 1
+                                        end if
+                                    end do latloop
+                                end do
+                            end do
+                        end if
+                        esize= esize-2
+                    end do
+                    ! Now add the leftover lattice points as atoms
+                    do o = 1, size_ele(ie)
+                        do n = 1, size_ele(ie)
+                            do m = 1, size_ele(ie)
+                                if(lat_points(m,n,o)) then 
+                                    call get_interp_pos(m,n,o, ie, ratom(:,:))
+                                    do ibasis = 1, basisnum(lat_ele(ie))
+                                        call apply_periodic(r_atom(:,ibasis))
+                                        call add_atom(0, basis_type(ibasis,lat_ele(ie)), sbox_ele(ie), ratom(:,ibasis))
+                                    end do
+                                end if
+                            end do
+                        end do
+                    end do
+                end if
+            end if
+        end do 
+
+        !Once we finish adding elements delete the old ones
+        call delete_elements(slip_enum, slip_eles(1:slip_enum))
+
+        !Output data
+        if(.not.efill) then 
+            print *, "We refine ", slip_enum, " elements into ", atom_num - old_atom_num , " atoms"
+        else 
+            print *, "We refine ", slip_enum, " elements into ", atom_num - old_atom_num , " atoms and ", new_ele_num, " elements"
+        end if
+
+    end subroutine run_slip_plane 
+
+    subroutine parse(arg_pos)
+        !This subroutine parses the input arguments to the mode
+        integer, intent(inout) :: arg_pos 
+
+        integer :: arglen
+        character(len = 100) :: textholder
+
+        !First read the dimension
+        arg_pos = arg_pos +1
+        call get_command_argument(arg_pos,textholder, arglen) 
+        if(arglen == 0) stop "Incorrect slip_plane command. Please check documentation for correct format"
+        
+        !Check to make sure that the dimension is correct
+        select case(trim(adjustl(textholder))) 
+        case('x','X')
+            sdim = 1
+        case('y','Y')
+            sdim = 2
+        case('z','Z')
+            sdim = 3
+        case default
+            print *, "Error: dimension ", trim(adjustl(textholder)), " is not accepted. Please select from x, y, or z"
+        end select 
+
+        !Now parse the position of the slip plane
+        arg_pos = arg_pos + 1
+        call get_command_argument(arg_pos, textholder, arglen)
+        if(arglen == 0) stop "Incorrect slip_plane command. Please check documentation for correct format"
+        call parse_pos(sdim, textholder, spos)
+
+        !Now check to see if efill was passed 
+        arg_pos = arg_pos + 1
+        if(.not.(arg_pos > command_argument_count())) then 
+            call get_command_argument(arg_pos, textholder, arglen)
+            if(arglen == 0) stop "Incorrect slip_plane command. Please check documentation for correct format"
+            if(trim(adjustl(textholder)) == "efill") then 
+                arg_pos = arg_pos +1 
+                efill = .true.
+            end if
+        end if
+    end subroutine parse
+end module opt_slip_plane

src/parameters.f90

@@ -3,7 +3,8 @@ module parameters
     implicit none
 
     !Default precision
-    integer, parameter :: dp= selected_real_kind(15,307)
+    integer, parameter :: dp= selected_real_kind(15,307), &
+                          min_efillsize = 11
     !Parameters for floating point tolerance
     real(kind=dp), parameter :: lim_zero = epsilon(1.0_dp), &
                                 lim_large = huge(1.0_dp), &