From f63335708b5c38be318e8e5452e23b6e3ba7c102 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Sat, 24 Oct 2020 09:13:16 -0400 Subject: [PATCH] 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