223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

fixed bug in nearest neighbor search, corrected error message for kdtree2.f90

This commit is contained in:
Taymor El Achkar 2012-05-08 13:16:59 +00:00
parent eb8265b914
commit a9a72cee97
4 changed files with 60 additions and 57 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
code/IO.f90
View File

@ -1280,6 +1280,8 @@ subroutine IO_error(error_ID,e,i,g,ext_msg)
msg = 'I_TO_HALTON-error: An input base BASE is <= 1' msg = 'I_TO_HALTON-error: An input base BASE is <= 1'
case (406_pInt) case (406_pInt)
msg = 'Prime-error: N must be between 0 and PRIME_MAX' msg = 'Prime-error: N must be between 0 and PRIME_MAX'
case (407_pInt)
msg = 'Dimension in nearest neigbor search wrong'
case (450_pInt) case (450_pInt)
msg = 'unknown symmetry type specified' msg = 'unknown symmetry type specified'
case (460_pInt) case (460_pInt)

17
code/damask.core.pyf
View File

@ -160,16 +160,19 @@ python module core ! in
real*8, dimension(res[0],res[1],res[2]),intent(out),depend(res[0],res[1],res[2]) :: vm real*8, dimension(res[0],res[1],res[2]),intent(out),depend(res[0],res[1],res[2]) :: vm
end subroutine math_equivStrain33_field end subroutine math_equivStrain33_field
subroutine math_nearestNeighborSearch(res,defgradAv,geomdim,,domainPoints,querySet,domainSet,indices) ! in :math:math.f90 subroutine math_nearestNeighborSearch(spatialDim,Favg,geomdim,queryPoints,domainPoints,querySet,domainSet,indices) ! in :math:math.f90
! input variables ! input variables
integer, dimension(3), intent(in) :: res integer, intent(in) :: spatialDim
integer, intent(in) :: domainPoints real*8, dimension(3,3), intent(in) :: Favg
real*8, dimension(3), intent(in) :: geomdim real*8, dimension(3), intent(in) :: geomdim
real*8, dimension(3,3), intent(in) :: defgradAv integer, intent(in) :: queryPoints
real*8, dimension(res[0],res[1],res[2],3), intent(in), depend(res[0],res[1],res[2]) :: querySet integer, intent(in) :: domainPoints
real*8, dimension(domainPoints,3), intent(in), depend(domainPoints) :: domainSet real*8, dimension(queryPoints,spatialDim), intent(in), depend(queryPoints,spatialDim) :: querySet
real*8, dimension(domainPoints,spatialDim), intent(in), depend(domainPoints,spatialDim) :: domainSet
! output variable ! output variable
integer, dimension(res[0]*res[1]*res[2]), intent(out), depend(res[0],res[1],res[2]) :: indices integer, dimension(queryPoints), intent(out), depend(queryPoints) :: indices
! depending on input
real*8, dimension(3,(3**spatialDim)*domainPoints), depend(domainPoints,spatialDim) :: domainSetLarge
end subroutine math_nearestNeighborSearch end subroutine math_nearestNeighborSearch
end module math end module math

12
code/kdtree2.f90
View File

@ -236,7 +236,7 @@ subroutine pq_max(a,e)
if (a%heap_size .gt. 0) then if (a%heap_size .gt. 0) then
e = a%elems(1) e = a%elems(1)
else else
call IO_error (500_pInt, ext_msg='PQ_MAX: heap_size < 1') call IO_error (460_pInt, ext_msg='PQ_MAX: heap_size < 1')
endif endif
return return
end subroutine pq_max end subroutine pq_max
@ -250,7 +250,7 @@ real(pReal) function pq_maxpri(a)
if (a%heap_size .gt. 0) then if (a%heap_size .gt. 0) then
pq_maxpri = a%elems(1)%dis pq_maxpri = a%elems(1)%dis
else else
call IO_error (500_pInt,ext_msg='PPQ_MAX_PRI: heap_size < 1') call IO_error (460_pInt,ext_msg='PPQ_MAX_PRI: heap_size < 1')
endif endif
return return
end function pq_maxpri end function pq_maxpri
@ -281,7 +281,7 @@ subroutine pq_extract_max(a,e)
call heapify(a,1_pInt) call heapify(a,1_pInt)
return return
else else
call IO_error (500_pInt,ext_msg='PQ_EXTRACT_MAX: attempted to pop non-positive PQ') call IO_error (460_pInt,ext_msg='PQ_EXTRACT_MAX: attempted to pop non-positive PQ')
end if end if
end subroutine pq_extract_max end subroutine pq_extract_max
@ -456,7 +456,7 @@ subroutine pq_delete(a,i)
integer(pInt) :: i integer(pInt) :: i
if ((i .lt. 1) .or. (i .gt. a%heap_size)) then if ((i .lt. 1) .or. (i .gt. a%heap_size)) then
call IO_error (500_pInt,ext_msg='PQ_DELETE: attempt to remove out of bounds element') call IO_error (460_pInt,ext_msg='PQ_DELETE: attempt to remove out of bounds element')
endif endif
! swap the item to be deleted with the last element ! swap the item to be deleted with the last element
@ -659,7 +659,7 @@ function kdtree2_create(input_data,myDim,sort,rearrange) result (mr)
write (*,*) 'KD_TREE_TRANS: note, that new format is myData(1:D,1:N)' write (*,*) 'KD_TREE_TRANS: note, that new format is myData(1:D,1:N)'
write (*,*) 'KD_TREE_TRANS: with usually N >> D. If N =approx= D, then a k-d tree' write (*,*) 'KD_TREE_TRANS: with usually N >> D. If N =approx= D, then a k-d tree'
write (*,*) 'KD_TREE_TRANS: is not an appropriate data structure.' write (*,*) 'KD_TREE_TRANS: is not an appropriate data structure.'
call IO_error (500_pInt) call IO_error (460_pInt)
end if end if
call build_tree(mr) call build_tree(mr)
@ -1356,7 +1356,7 @@ subroutine validate_query_storage(n)
integer(pInt), intent(in) :: n integer(pInt), intent(in) :: n
if (int(size(sr%results,1),pInt) .lt. n) then if (int(size(sr%results,1),pInt) .lt. n) then
call IO_error (500_pInt,ext_msg='KD_TREE_TRANS: not enough storage for results(1:n)') call IO_error (460_pInt,ext_msg='KD_TREE_TRANS: not enough storage for results(1:n)')
endif endif
return return

70
code/math.f90
View File

@ -3818,61 +3818,59 @@ subroutine calculate_cauchy(res,defgrad,p_stress,c_stress)
end subroutine calculate_cauchy end subroutine calculate_cauchy
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine math_nearestNeighborSearch(res, defgradAv, geomdim, domainPoints, querySet, domainSet, indices) subroutine math_nearestNeighborSearch(spatialDim, Favg, geomdim, queryPoints, domainPoints, querySet, domainSet, indices)
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!Obtain the nearest neighbour ! Obtain the nearest neighbor in domain set for all points in querySet
! !
use kdtree2_module use kdtree2_module
use IO, only: &
IO_error
implicit none implicit none
! input variables ! input variables
integer(pInt), dimension(3), intent(in) :: res integer(pInt), intent(in) :: spatialDim
integer(pInt), intent(in) :: domainPoints real(pReal), dimension(3,3), intent(in) :: Favg
real(pReal), dimension(3), intent(in) :: geomdim real(pReal), dimension(3), intent(in) :: geomdim
real(pReal), dimension(3,3), intent(in) :: defgradAv integer(pInt), intent(in) :: domainPoints
real(pReal), dimension(res(1),res(2),res(3),3), intent(in) :: querySet integer(pInt), intent(in) :: queryPoints
real(pReal), dimension(domainPoints,3), intent(in) :: domainSet real(pReal), dimension(queryPoints,spatialDim), intent(in) :: querySet
real(pReal), dimension(domainPoints,spatialDim), intent(in) :: domainSet
! output variable ! output variable
integer(pInt), dimension(res(1)*res(2)*res(3)), intent(out) :: indices integer(pInt), dimension(queryPoints), intent(out) :: indices
! other variables depending on input
real(pReal), dimension(3,(3_pInt**spatialDim)*domainPoints) :: domainSetLarge
! other variables ! other variables
real(pReal), dimension(:,:), allocatable :: querySetLarge integer(pInt) :: i,j, l,m,n
integer(pInt) :: i,j,k, l,m,n, ielem_large, spatial_dim
real(pReal), dimension(3) :: shift
type(kdtree2), pointer :: tree type(kdtree2), pointer :: tree
type(kdtree2_result), dimension(1) :: Results type(kdtree2_result), dimension(1) :: Results
shift = math_mul33x3(defgradAv,geomdim) if (size(querySet(1,:)) /= spatialDim) call IO_error(407_pInt,ext_msg='query set')
if (size(domainSet(1,:)) /= spatialDim) call IO_error(407_pInt,ext_msg='domain set')
ielem_large = 0_pInt
if(res(3) == 1_pInt) then i = 0_pInt
spatial_dim = 2_pInt if(spatialDim == 2_pInt) then
allocate(querySetLarge(2,(res(1)*res(2))*9_pInt)) do j = 1_pInt, domainPoints
do j=1_pInt, res(2); do i=1_pInt, res(1) do l = -1_pInt, 1_pInt; do m = -1_pInt, 1_pInt
do l = -1, 1; do m = -1, 1 i = i + 1_pInt
ielem_large = ielem_large + 1_pInt domainSetLarge(1:3,i) = domainSet(j,1:3) + math_mul33x3(Favg,real([l,m,0_pInt],pReal)*geomdim)
querySetLarge(1:2,ielem_large) = querySet(i,j,1,1:2) + real([l,m],pReal)* shift(1:2)
enddo; enddo;
enddo; enddo enddo; enddo
enddo
else else
allocate(querySetLarge(3,(res(1)*res(2)*res(3))*27)) do j = 1_pInt, domainPoints
spatial_dim = 3_pInt do l = -1_pInt, 1_pInt; do m = -1_pInt, 1_pInt; do n = -1_pInt, 1_pInt
do k=1_pInt,res(3); do j=1_pInt, res(2); do i=1_pInt, res(1) i = i + 1_pInt
do l = -1, 1; do m = -1, 1; do n = -1, 1 domainSetLarge(1:3,i) = domainSet(j,1:3) + math_mul33x3(Favg,real([l,m,n],pReal)*geomdim)
ielem_large = ielem_large + 1_pInt
querySetLarge(1:3,ielem_large) = querySet(i,j,k,1:3) + real([l,m,n],pReal)* shift
enddo; enddo; enddo;
enddo; enddo; enddo enddo; enddo; enddo
enddo
endif endif
tree => kdtree2_create(querySetLarge,sort=.true.,rearrange=.true.) tree => kdtree2_create(domainSetLarge,sort=.true.,rearrange=.true.)
ielem_large = 0_pInt do j = 1_pInt, queryPoints
do k=1_pInt,res(3); do j=1_pInt, res(2); do i=1_pInt, res(1) call kdtree2_n_nearest(tp=tree, qv=querySet(j,1:spatialDim),nn=1_pInt, results = Results)
ielem_large = ielem_large + 1_pInt indices(j) = Results(1)%idx
call kdtree2_n_nearest(tp=tree, qv=domainSet(ielem_large,1:spatial_dim),nn=1_pInt, results = Results) enddo
indices(ielem_large) = Results(1)%idx !/3_pInt**spatial_dim +1_pInt
enddo; enddo; enddo
deallocate(querySetLarge)
end subroutine math_nearestNeighborSearch end subroutine math_nearestNeighborSearch