diff --git a/src/Makefile b/src/Makefile index 18881ae..4d6971d 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,9 +1,15 @@ -FC=gfortran -#FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin -heap-arrays +FC=ifort + +#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 $@ diff --git a/src/call_option.f90 b/src/call_option.f90 index d55e1ee..af63333 100644 --- a/src/call_option.f90 +++ b/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 diff --git a/src/elements.f90 b/src/elements.f90 index 4b3d5e0..7e658c2 100644 --- a/src/elements.f90 +++ b/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)) - temp_int(1:ele_size) = lat_ele + allocate(temp_int(n+ele_size+buffer_size)) + 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)) - temp_int(1:ele_size) = tag_ele + allocate(temp_int(n+ele_size+buffer_size)) + 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)) - temp_int(1:ele_size) = size_ele + allocate(temp_int(n+ele_size+buffer_size)) + 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)) - temp_int(1:ele_size) = lat_ele + allocate(temp_int(n+ele_size+buffer_size)) + 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)) - char_temp(1:ele_size) = type_ele + allocate(char_temp(n+ele_size+buffer_size)) + 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)) - temp_ele_real(:,:,:,1:ele_size) = r_node + allocate(temp_ele_real(3, max_basisnum, max_ng_node, n+ele_size+buffer_size)) + 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 @@ -292,8 +294,7 @@ module elements size_ele(ele_num) = size lat_ele(ele_num) = lat sbox_ele(ele_num) = sbox - r_node(:,:,:,ele_num) = r(:,:,:) - node_num = node_num + ng_node(lat) + r_node(:,:,:,ele_num) = r(:,:,:) 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 diff --git a/src/mode_create.f90 b/src/mode_create.f90 index 6f1ca28..7994f55 100644 --- a/src/mode_create.f90 +++ b/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 diff --git a/src/opt_slip_plane.f90 b/src/opt_slip_plane.f90 new file mode 100644 index 0000000..80da930 --- /dev/null +++ b/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 diff --git a/src/parameters.f90 b/src/parameters.f90 index f261552..8677381 100644 --- a/src/parameters.f90 +++ b/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), &