From 22e250093b1a9fc1c7f9c5fcf60e7162aef13fad Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Fri, 23 Oct 2020 17:17:40 -0400 Subject: [PATCH] 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