You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
313 lines
12 KiB
313 lines
12 KiB
module mode_da
|
|
!This mode is used to calculate the dislocation analysis algorithm for both CG and atomistic regions
|
|
|
|
use parameters
|
|
use io
|
|
use elements
|
|
use neighbors
|
|
|
|
implicit none
|
|
|
|
integer :: num_d_points
|
|
integer, allocatable :: is_slipped(:,:,:)
|
|
real(kind=dp), allocatable :: disl_line(:,:), disl_data(:)
|
|
|
|
private :: write_xyz
|
|
public
|
|
contains
|
|
|
|
subroutine da(arg_pos)
|
|
!This is the main calling subroutine for the dislocation analysis code
|
|
integer, intent(out) :: arg_pos
|
|
real(kind=dp), dimension(6) :: temp_box_bd
|
|
|
|
integer :: ppos
|
|
character(len=100) :: outfile
|
|
|
|
!Now call the dislocation analysis code for the coarse-grained elements
|
|
call parse_command(arg_pos)
|
|
|
|
call read_in(1, (/0.0_dp,0.0_dp,0.0_dp/), temp_box_bd)
|
|
|
|
!Now allocate the necessary variables
|
|
if(allocated(is_slipped)) deallocate(is_slipped)
|
|
allocate(is_slipped(max_basisnum, max_ng_node, ele_num))
|
|
|
|
if(allocated(disl_line)) deallocate(disl_line)
|
|
num_d_points=0
|
|
allocate(disl_line(3,1024))
|
|
allocate(disl_data(1024))
|
|
|
|
call cga_da
|
|
|
|
!Now create the output file num and write out to xyz format
|
|
ppos = scan(trim(infiles(1)),".", BACK= .true.)
|
|
if ( ppos > 0 ) then
|
|
outfile = 'da_'//infiles(1)(1:ppos)//'vtk'
|
|
else
|
|
outfile = 'da_'//infiles(1)//'.vtk'
|
|
end if
|
|
call write_da_vtk(outfile)
|
|
|
|
ppos = scan(trim(infiles(1)),".", BACK= .true.)
|
|
if ( ppos > 0 ) then
|
|
outfile = 'da_'//infiles(1)(1:ppos)//'xyz'
|
|
else
|
|
outfile = 'da_'//infiles(1)//'.xyz'
|
|
end if
|
|
call write_da_xyz(outfile)
|
|
|
|
return
|
|
end subroutine da
|
|
|
|
|
|
subroutine parse_command(arg_pos)
|
|
!This subroutine parses the arguments for mode command
|
|
integer, intent(out) :: arg_pos
|
|
integer :: i, arglen
|
|
character(len = 100) :: textholder
|
|
logical :: file_exists
|
|
|
|
!Read the input file
|
|
call get_command_argument(2,textholder, arglen)
|
|
if (arglen == 0) stop "Missing file for dislocation analysis"
|
|
call get_in_file(textholder)
|
|
|
|
arg_pos = 3
|
|
return
|
|
end subroutine parse_command
|
|
|
|
subroutine cga_da
|
|
|
|
integer :: i, j, k, l, ibasis, nei, close_neighbors(2,4,6), far_neighbor(4,6), minindex, minindices(2), linevec(3)
|
|
|
|
integer :: face_types(ele_num*6), face_ele(6*ele_num), diff_pair(2,6)
|
|
real(kind = dp) :: face_centroids(3, ele_num*6), r(3), rc, vnode(3, max_basisnum, 4), vnorm, &
|
|
max_node_dist, ndiff, rmax(3), rmin(3), rnorm, v1(3), v2(3), v1norm, v2norm, theta1, theta2, &
|
|
diff_mag(6), diffvec(3)
|
|
real(kind=dp), allocatable :: disl_line_temp(:,:), disl_data_temp(:)
|
|
logical :: is_edge
|
|
|
|
!Initialize variables
|
|
l = 0
|
|
max_node_dist = 0
|
|
|
|
!Now save the close and far neighbors of the nodes. This is done to attempt to figure out the character of the dislocation
|
|
!by comparing how the burgers vector is distributed over the face.
|
|
do j = 1, 6
|
|
close_neighbors(1, 1, j) = cubic_faces(4, j)
|
|
close_neighbors(2, 1, j) = cubic_faces(2, j)
|
|
far_neighbor(1,j) = cubic_faces(3,j)
|
|
|
|
close_neighbors(1, 2, j) = cubic_faces(1, j)
|
|
close_neighbors(2, 2, j) = cubic_faces(3, j)
|
|
far_neighbor(2,j) = cubic_faces(4,j)
|
|
|
|
close_neighbors(1, 3, j) = cubic_faces(2, j)
|
|
close_neighbors(2, 3, j) = cubic_faces(4, j)
|
|
far_neighbor(3,j) = cubic_faces(1,j)
|
|
|
|
close_neighbors(1, 4, j) = cubic_faces(3, j)
|
|
close_neighbors(2, 4, j) = cubic_faces(1, j)
|
|
far_neighbor(4,j) = cubic_faces(2,j)
|
|
end do
|
|
|
|
!First calculate all of the face centroids
|
|
do i = 1, ele_num
|
|
do j = 1, 6
|
|
r(:) = 0.0_dp
|
|
rmax= -Huge(1.0)
|
|
rmin= Huge(1.0)
|
|
do k = 1, 4
|
|
do ibasis = 1, basisnum(lat_ele(i))
|
|
r = r + r_node(:, ibasis, cubic_faces(k,j), i)
|
|
rmax= max(rmax, r_node(:, ibasis, cubic_faces(k,j), i))
|
|
rmin= min(rmin, r_node(:, ibasis, cubic_faces(k,j), i))
|
|
end do
|
|
end do
|
|
|
|
rnorm=norm2(rmax-rmin)
|
|
max_node_dist = max(max_node_dist, rnorm)
|
|
r = r/(basisnum(lat_ele(i))*4)
|
|
|
|
!add the face centroids, the type, and map the elements faces to the face arrays
|
|
l = l + 1
|
|
face_centroids(:, l) = r
|
|
face_types(l) = j
|
|
face_ele(l) = i
|
|
|
|
|
|
end do
|
|
end do
|
|
|
|
!Now set the cut off distance as half the distance from opposite nodes in the largest element
|
|
rc = 0.5*max_node_dist
|
|
call calc_NN(l, face_centroids(:,1:l), rc)
|
|
|
|
!Now loop overall the faces and make sure the nearest neighbor is the opposite type. If it isn't than we dscard it
|
|
is_slipped = 0
|
|
do i = 1, l
|
|
nei = nn(i)
|
|
|
|
!Skip if it's 0 or the closest face belongs to the same element
|
|
if ((nei == 0).or.(face_ele(i) == face_ele(nei))) cycle
|
|
|
|
|
|
!Check the face types, the way that the faces are laid out in the cubic_faces array face 1's opposite is 6 and
|
|
! face 2's opposite is 5 and etc
|
|
vnode = 0
|
|
if(face_types(i) == (7-face_types(nei))) then
|
|
vnorm = 0
|
|
do j = 1, 4
|
|
do ibasis = 1, basisnum(lat_ele(face_ele(i)))
|
|
!Compute the vectors between all nodes at the face.
|
|
vnode(:,ibasis,j) = r_node(:,ibasis, cubic_faces(j,face_types(i)), face_ele(i)) - &
|
|
r_node(:,ibasis, cubic_faces(j,face_types(nei)), face_ele(nei))
|
|
end do
|
|
end do
|
|
|
|
!Now calculate the maximum difference between nodes
|
|
vnorm = 0
|
|
do j = 1, 4
|
|
do k=j,4
|
|
v1=vnode(:,1,j) - vnode(:,1,k)
|
|
v1norm=norm2(v1)
|
|
if (vnorm < v1norm) then
|
|
vnorm = v1norm
|
|
diffvec= v1
|
|
end if
|
|
end do
|
|
end do
|
|
!Now calculate the difference between the largest norm and the smallest norm, if it's larger than 0.5 then mark it
|
|
!as slipped. This value probably can be converted to a variable value that depends on the current lattice parameter
|
|
!I think 0.5 works ok though.
|
|
if (vnorm>0.5_dp) then
|
|
print *, "Element number ", face_ele(i), " is dislocated along face ", face_types(i), &
|
|
" with neighbor ", face_ele(nei), " with max displacement of ", vnorm
|
|
l=0
|
|
do j = 1, 4
|
|
is_slipped(:, cubic_faces(j,face_types(i)), face_ele(i)) = 1
|
|
|
|
!This portion of the code is used to determine what the character of the dislocation most likely is
|
|
!effectively how this works is we look for the two nodes with the most similar burgers vector.
|
|
!The vector between these two nodes is close to the line direction and as a result we can probably estimate
|
|
!if it's close to edge, screw, or if it's pretty fairly mixed.
|
|
do k = j+1, 4
|
|
l=l+1
|
|
diff_mag(l) = norm2(vnode(:, 1, j)-vnode(:,1,k))
|
|
diff_pair(1,l)=cubic_faces(j, face_types(i))
|
|
diff_pair(2,l)=cubic_faces(k, face_types(i))
|
|
end do
|
|
end do
|
|
minindex = minloc(diff_mag,1)
|
|
!Now figure out if the min difference between nodes is associated with a 112 direction
|
|
is_edge = .false.
|
|
do j = 1,6
|
|
if((diff_pair(1,minindex) == oneonetwopairs(1,j)).and.(diff_pair(2,minindex)==oneonetwopairs(2,j))) then
|
|
is_edge=.true.
|
|
exit
|
|
else if((diff_pair(2,minindex)==oneonetwopairs(1,j)).and.(diff_pair(1,minindex)==oneonetwopairs(2,j)))then
|
|
is_edge=.true.
|
|
exit
|
|
end if
|
|
end do
|
|
|
|
if(is_edge) then
|
|
print *, 'Dislocation has primarily edge character'
|
|
else
|
|
print *, 'Dislocation has primarily screw character'
|
|
end if
|
|
|
|
if( i < nei) then
|
|
num_d_points=num_d_points+2
|
|
if(num_d_points > size(disl_line,2)) then
|
|
allocate(disl_line_temp(3,size(disl_line,2)+1024))
|
|
disl_line_temp=0.0_dp
|
|
disl_line_temp(:, 1:size(disl_line,2))=disl_line
|
|
call move_alloc(disl_line_temp, disl_line)
|
|
end if
|
|
if(num_d_points/2 > size(disl_data)) then
|
|
allocate(disl_data_temp(size(disl_data)+1024))
|
|
disl_data_temp=0
|
|
disl_data_temp(1:size(disl_data)) = disl_data
|
|
call move_alloc(disl_data_temp, disl_data)
|
|
end if
|
|
print *, num_d_points/2, vnorm
|
|
disl_data(num_d_points/2)=vnorm
|
|
disl_line(:,num_d_points-1) =r_node(:, 1, diff_pair(1,minindex), face_ele(i))
|
|
disl_line(:,num_d_points) =r_node(:, 1, diff_pair(2,minindex), face_ele(i))
|
|
end if
|
|
end if
|
|
end if
|
|
end do
|
|
|
|
end subroutine cga_da
|
|
|
|
subroutine write_da_vtk(outfile)
|
|
character (len=*), intent(in) :: outfile
|
|
|
|
integer :: i
|
|
open(unit=11, file=outfile, action='write', status='replace', position='rewind')
|
|
|
|
write(11, '(a)') '# vtk DataFile Version 4.0.1'
|
|
write(11, '(a)') 'Dislocation line analyzed via cacmb dislocation analysis'
|
|
write(11, '(a)') 'ASCII'
|
|
write(11, '(a)') 'DATASET UNSTRUCTURED_GRID'
|
|
write(11, *) 'POINTS', num_d_points, 'float'
|
|
|
|
!Write the points to the file
|
|
do i = 1, num_d_points
|
|
write(11, '(3f23.15)') disl_line(:, i)
|
|
end do
|
|
|
|
write(11,*) 'CELLS', num_d_points/2, num_d_points/2*3
|
|
do i=1, num_d_points, 2
|
|
write(11,'(3i16)') 2, i-1, i
|
|
end do
|
|
|
|
write(11,*) 'CELL_TYPES', num_d_points/2
|
|
do i=1, num_d_points, 2
|
|
write(11,'(i16)') 3
|
|
end do
|
|
|
|
write(11,*) 'CELL_DATA', num_d_points/2
|
|
write(11, '(a)') 'SCALARS burgers_vec, float'
|
|
write(11, '(a)') 'LOOKUP_TABLE default'
|
|
do i=1, num_d_points, 2
|
|
write(11, '(f23.15)') disl_data((i+1)/2)
|
|
end do
|
|
|
|
close(11)
|
|
|
|
end subroutine write_da_vtk
|
|
|
|
subroutine write_da_xyz(outfile)
|
|
!This subroutine write the element positions to a .xyz file and marks whether they are slipped or not
|
|
character(len=*), intent(in) :: outfile
|
|
integer :: i, ibasis, inod, outn
|
|
|
|
open(unit=11, file=outfile, action='write', status='replace', position='rewind')
|
|
|
|
!Write number of node_atoms
|
|
write(11, '(i16)') node_atoms
|
|
!Write comment
|
|
write(11, '(a)') "is_slipped x y z"
|
|
!Write nodal positions
|
|
outn = 0
|
|
do i = 1, ele_num
|
|
do inod = 1, ng_node(lat_ele(i))
|
|
do ibasis = 1, basisnum(lat_ele(i))
|
|
write(11, '(1i16, 3f23.15)') is_slipped(ibasis,inod,i), r_node(:,ibasis,inod,i)
|
|
outn = outn + 1
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
if(outn /= node_atoms) then
|
|
print *, "outn", outn, " doesn't equal node_atoms ", node_atoms
|
|
end if
|
|
|
|
close(11)
|
|
end subroutine write_da_xyz
|
|
end module mode_da
|