From b5629b1563940375f69b770a9af9a7b996e5cead Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Fri, 23 Oct 2020 14:13:35 -0400 Subject: [PATCH 1/6] Working changes to slip_plane code --- src/Makefile | 2 +- src/call_option.f90 | 3 +++ src/elements.f90 | 20 ++++++++++++++++++++ src/mode_create.f90 | 2 +- src/parameters.f90 | 3 ++- 5 files changed, 27 insertions(+), 3 deletions(-) diff --git a/src/Makefile b/src/Makefile index feb933f..3b170e7 100644 --- a/src/Makefile +++ b/src/Makefile @@ -2,7 +2,7 @@ FC=ifort 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 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: 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 6df859b..a1ab9bc 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -736,4 +736,24 @@ module elements end subroutine lattice_map + subroutine get_interp_pos(i,j,k, ie, r) + !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, r, s, t, ie, inod -= + real(kind=dp), dimension(3,max_basisnum), intent(out) :: r + + 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) + r(:) = 0 + do ibasis = 1, bnum + do inod = 1, 8 + r(:,ibasis) = r(:,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 bcec2d5..00d7d8b 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/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), & From 22e250093b1a9fc1c7f9c5fcf60e7162aef13fad Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Fri, 23 Oct 2020 17:17:40 -0400 Subject: [PATCH 2/6] Working changes. Working on the efill code --- src/elements.f90 | 4 +- src/opt_slip_plane.f90 | 151 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 153 insertions(+), 2 deletions(-) create mode 100644 src/opt_slip_plane.f90 diff --git a/src/elements.f90 b/src/elements.f90 index a1ab9bc..42eec2d 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -747,8 +747,8 @@ module elements 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) r(:) = 0 - do ibasis = 1, bnum - do inod = 1, 8 + do ibasis = 1, basisnum(lat_ele(ie)) + do inod = 1, ng_node(lat_ele(ie)) r(:,ibasis) = r(:,ibasis) + a_shape(inod) * r_node(:,ibasis,inod,ie) end do end do diff --git a/src/opt_slip_plane.f90 b/src/opt_slip_plane.f90 new file mode 100644 index 0000000..36a3d4d --- /dev/null +++ b/src/opt_slip_plane.f90 @@ -0,0 +1,151 @@ +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), & + rfill(3,max_basisnum, max_ng_node) + + integer, allocatable :: slip_eles(:), temp_int(:) + real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3) + 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 + do m = 1, size_ele(ie) - esize + do inod = 1, ng_node(lat_ele(ie) + vlat = ele(:,inod) + (/ m, n, o /) + call get_interp_pos(vlat(1), vlat(2), vlat(3), ie, rfill(:,:,inod)) + + end do + !Check to see if the plane intersects this element if not then add it + if((spos < maxval(rfill(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie)),ie))).and. & + (spos > minval(rfill(sdim,1:basisnum(lat_ele(ie),1:ng_node(lat_ele(ie)),ie)))) then + nump_ele = nump_ele - esize**3 + + end if + end do + end do + end do + end if + 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" + 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 From f63335708b5c38be318e8e5452e23b6e3ba7c102 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Sat, 24 Oct 2020 09:13:16 -0400 Subject: [PATCH 3/6] Working changes --- src/Makefile | 13 +++++++-- src/elements.f90 | 36 +++++++++++++----------- src/mode_create.f90 | 2 +- src/opt_group.f90 | 5 ++-- src/opt_slip_plane.f90 | 63 +++++++++++++++++++++--------------------- 5 files changed, 66 insertions(+), 53 deletions(-) diff --git a/src/Makefile b/src/Makefile index 3b170e7..b7b0c5e 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,6 +1,13 @@ -FC=ifort -FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin -heap-arrays +FC=gfortran + +#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 + +#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 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 @@ -16,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/elements.f90 b/src/elements.f90 index 42eec2d..47bc22c 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -211,17 +211,17 @@ module elements call move_alloc(temp_int, lat_ele) allocate(temp_int(n+ele_num+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)) - 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)) - 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, sbox_ele) @@ -282,8 +282,8 @@ 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(:,:,:) + node_num = node_num + ng_node(lat) end subroutine add_element @@ -669,17 +669,17 @@ module elements esize = size_ele(ie) select case(iface) case(1) - pos = (/ real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp, -10.0_dp**-2.0_dp /) + pos = (/ real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp, -10.0_dp**(-2.0_dp) /) case(2) - pos = (/ real(esize-1,dp)/2.0_dp, -10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp /) + pos = (/ real(esize-1,dp)/2.0_dp, -10.0_dp**(-2.0_dp), real(esize-1,dp)/2.0_dp /) case(3) - pos = (/ (esize-1)+10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp /) + pos = (/ (esize-1)+10.0_dp**(-2.0_dp), real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp /) case(4) - pos = (/ real(esize-1,dp)/2.0_dp, (esize-1)+10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp /) + pos = (/ real(esize-1,dp)/2.0_dp, (esize-1)+10.0_dp**(-2.0_dp), real(esize-1,dp)/2.0_dp /) case(5) - pos = (/ -10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp /) + pos = (/ -10.0_dp**(-2.0_dp), real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp /) case(6) - pos = (/ real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp, (esize-1)+10.0_dp**-2.0_dp /) + pos = (/ real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp, (esize-1)+10.0_dp**(-2.0_dp) /) end select !Now transform it to real space and adjust it to the position of the element in the first node. @@ -736,20 +736,24 @@ module elements end subroutine lattice_map - subroutine get_interp_pos(i,j,k, ie, r) + 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, r, s, t, ie, inod -= - real(kind=dp), dimension(3,max_basisnum), intent(out) :: r + 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) - r(:) = 0 + rout(:,:) = 0 do ibasis = 1, basisnum(lat_ele(ie)) do inod = 1, ng_node(lat_ele(ie)) - r(:,ibasis) = r(:,ibasis) + a_shape(inod) * r_node(:,ibasis,inod,ie) + call rhombshape(r,s,t,a_shape) + rout(:,ibasis) = rout(:,ibasis) + a_shape(inod) * r_node(:,ibasis,inod,ie) end do end do diff --git a/src/mode_create.f90 b/src/mode_create.f90 index 00d7d8b..7994f55 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -135,7 +135,7 @@ module mode_create case('bcc') call build_with_rhomb(box_lat_vert, bcc_mat) case default - print *, "Element type ", trim(adjustl(element_type)), " not accepted in mode create, please specify a supported ", & + print *, "Element type ", trim(adjustl(element_type)), " not accepted in mode create, please specify a supported ",& "element type" stop 3 end select diff --git a/src/opt_group.f90 b/src/opt_group.f90 index 855e04e..0956677 100644 --- a/src/opt_group.f90 +++ b/src/opt_group.f90 @@ -300,7 +300,7 @@ module opt_group case('elements','element') if (group_ele_num > 0) then - print *, "Elements specifier used more than once in group id command with type both, either use type ", & + print *, "Elements specifier used more than once in group id command with type both, either use type ",& "elements or include atoms identifier" stop 3 @@ -824,7 +824,8 @@ module opt_group !Add the element, for the sbox we just set it to the same sbox that we get the orientation !from. In this case it is from the sbox of the first atom in the group. new_ele = new_ele+1 - call add_element(0,remesh_ele_type, working_esize, ilat, sbox_atom(atom_index(1)),r_new_node) + call add_element(0,remesh_ele_type, working_esize, ilat, & + sbox_atom(atom_index(1)),r_new_node) end if end if diff --git a/src/opt_slip_plane.f90 b/src/opt_slip_plane.f90 index 36a3d4d..4e4f1c1 100644 --- a/src/opt_slip_plane.f90 +++ b/src/opt_slip_plane.f90 @@ -17,11 +17,11 @@ module opt_slip_plane !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), & - rfill(3,max_basisnum, max_ng_node) - + integer :: ie, ia, slip_enum, old_atom_num, esize, new_ele_num, n, m, o, ele(3,8), nump_ele, inod, vlat(3) + integer, allocatable :: slip_eles(:), temp_int(:) - real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3) + real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), rfill(3,max_basisnum, max_ng_node) + integer :: type_interp(max_basisnum*max_esize**3) logical :: lat_points(max_esize,max_esize, max_esize) @@ -45,7 +45,7 @@ module opt_slip_plane !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 + (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)) @@ -64,35 +64,36 @@ module opt_slip_plane 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 - do m = 1, size_ele(ie) - esize - do inod = 1, ng_node(lat_ele(ie) - vlat = ele(:,inod) + (/ m, n, o /) - call get_interp_pos(vlat(1), vlat(2), vlat(3), ie, rfill(:,:,inod)) - + !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 + do m = 1, size_ele(ie) - esize + do inod = 1, ng_node(lat_ele(ie)) + vlat = ele(:,inod) + (/ m, n, o /) + call get_interp_pos(vlat(1), vlat(2), vlat(3), ie, rfill(:,:,inod)) + end do + !Check to see if the plane intersects this element if not then add it + if((spos < maxval(rfill(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie))))).and. & + (spos > minval(rfill(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie)))))) then + nump_ele = nump_ele - esize**3 + lat_points(m:m+esize, n:n+esize, o:o+esize) = .false. + call add_element(0, type_ele(ie), esize, lat_ele(ie), sbox_ele(ie), rfill) + end if end do - !Check to see if the plane intersects this element if not then add it - if((spos < maxval(rfill(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie)),ie))).and. & - (spos > minval(rfill(sdim,1:basisnum(lat_ele(ie),1:ng_node(lat_ele(ie)),ie)))) then - nump_ele = nump_ele - esize**3 - - end if end do end do - end do - end if + end if + end do end if end if end do From b9ce916e42ecc864259d3a6c44b43a4b267cf65e Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Tue, 27 Oct 2020 13:23:00 -0400 Subject: [PATCH 4/6] Current working changes to option-slip-plane --- src/Makefile | 4 ++-- src/elements.f90 | 11 ++++++----- src/opt_slip_plane.f90 | 22 +++++++++++++++++++--- 3 files changed, 27 insertions(+), 10 deletions(-) diff --git a/src/Makefile b/src/Makefile index 95f074a..4d6971d 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,8 +1,8 @@ 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 -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 #gfortran flags #FFLAGS=-mcmodel=large -O3 -g diff --git a/src/elements.f90 b/src/elements.f90 index 429a07b..8a5f54e 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -201,24 +201,25 @@ 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_size > size(size_ele)) then + if ( n+ele_num > size(size_ele)) then 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_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_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) @@ -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 diff --git a/src/opt_slip_plane.f90 b/src/opt_slip_plane.f90 index 4e4f1c1..c6c8468 100644 --- a/src/opt_slip_plane.f90 +++ b/src/opt_slip_plane.f90 @@ -17,10 +17,10 @@ module opt_slip_plane !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) + 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) + real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), rfill(3,max_basisnum, max_ng_node), ratom(3,max_basisnum) integer :: type_interp(max_basisnum*max_esize**3) logical :: lat_points(max_esize,max_esize, max_esize) @@ -86,14 +86,28 @@ module opt_slip_plane if((spos < maxval(rfill(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie))))).and. & (spos > minval(rfill(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie)))))) then nump_ele = nump_ele - esize**3 - lat_points(m:m+esize, n:n+esize, o:o+esize) = .false. + 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 end do end do end if 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 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 @@ -104,6 +118,8 @@ module opt_slip_plane !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 From 84a84e578a28f7bdbf652cf51b747d04b9c55a72 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Tue, 27 Oct 2020 16:44:03 -0400 Subject: [PATCH 5/6] Working version of slip_plane ooption with effill --- src/opt_slip_plane.f90 | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/src/opt_slip_plane.f90 b/src/opt_slip_plane.f90 index c6c8468..e607ae8 100644 --- a/src/opt_slip_plane.f90 +++ b/src/opt_slip_plane.f90 @@ -20,7 +20,8 @@ module opt_slip_plane 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) + 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) @@ -77,23 +78,29 @@ module opt_slip_plane ele = cubic_cell*(esize -1) do o = 1, size_ele(ie) - esize do n = 1, size_ele(ie) - esize - do m = 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(lat_points(m:m+esize-1, n:n+esize-1, o:o+esize-1) == 0)) cycle latloop !Check to see if the plane intersects this element if not then add it - if((spos < maxval(rfill(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie))))).and. & - (spos > minval(rfill(sdim,1:basisnum(lat_ele(ie)),1:ng_node(lat_ele(ie)))))) then + 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 + 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) @@ -102,6 +109,7 @@ module opt_slip_plane 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 From a528dc4c5245d4a6a4d1c033547b920dd1eed4f8 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Wed, 28 Oct 2020 17:53:12 -0400 Subject: [PATCH 6/6] Final fix to get working opt_slip_plane --- src/opt_slip_plane.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/opt_slip_plane.f90 b/src/opt_slip_plane.f90 index e607ae8..80da930 100644 --- a/src/opt_slip_plane.f90 +++ b/src/opt_slip_plane.f90 @@ -86,7 +86,7 @@ module opt_slip_plane end do !Check to make sure all lattice points exist for the current element - if(any(lat_points(m:m+esize-1, n:n+esize-1, o:o+esize-1) == 0)) cycle latloop + 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))))