From a9a72cee9799e98687d78cca3da0c6a7e6e23a90 Mon Sep 17 00:00:00 2001 From: Taymor El Achkar Date: Tue, 8 May 2012 13:16:59 +0000 Subject: [PATCH] fixed bug in nearest neighbor search, corrected error message for kdtree2.f90 --- code/IO.f90 | 2 ++ code/damask.core.pyf | 17 +++++---- code/kdtree2.f90 | 12 +++---- code/math.f90 | 86 ++++++++++++++++++++++---------------------- 4 files changed, 60 insertions(+), 57 deletions(-) diff --git a/code/IO.f90 b/code/IO.f90 index d0d822342..71b2a4de9 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -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' case (406_pInt) msg = 'Prime-error: N must be between 0 and PRIME_MAX' + case (407_pInt) + msg = 'Dimension in nearest neigbor search wrong' case (450_pInt) msg = 'unknown symmetry type specified' case (460_pInt) diff --git a/code/damask.core.pyf b/code/damask.core.pyf index ec29502f8..a439df716 100644 --- a/code/damask.core.pyf +++ b/code/damask.core.pyf @@ -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 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 - integer, dimension(3), intent(in) :: res - integer, intent(in) :: domainPoints + integer, intent(in) :: spatialDim + real*8, dimension(3,3), intent(in) :: Favg real*8, dimension(3), intent(in) :: geomdim - real*8, dimension(3,3), intent(in) :: defgradAv - real*8, dimension(res[0],res[1],res[2],3), intent(in), depend(res[0],res[1],res[2]) :: querySet - real*8, dimension(domainPoints,3), intent(in), depend(domainPoints) :: domainSet + integer, intent(in) :: queryPoints + integer, intent(in) :: domainPoints + real*8, dimension(queryPoints,spatialDim), intent(in), depend(queryPoints,spatialDim) :: querySet + real*8, dimension(domainPoints,spatialDim), intent(in), depend(domainPoints,spatialDim) :: domainSet ! 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 module math diff --git a/code/kdtree2.f90 b/code/kdtree2.f90 index 5fc769833..f13967f08 100644 --- a/code/kdtree2.f90 +++ b/code/kdtree2.f90 @@ -236,7 +236,7 @@ subroutine pq_max(a,e) if (a%heap_size .gt. 0) then e = a%elems(1) 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 return end subroutine pq_max @@ -250,7 +250,7 @@ real(pReal) function pq_maxpri(a) if (a%heap_size .gt. 0) then pq_maxpri = a%elems(1)%dis 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 return end function pq_maxpri @@ -281,7 +281,7 @@ subroutine pq_extract_max(a,e) call heapify(a,1_pInt) return 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 subroutine pq_extract_max @@ -456,7 +456,7 @@ subroutine pq_delete(a,i) integer(pInt) :: i 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 ! 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: with usually N >> D. If N =approx= D, then a k-d tree' write (*,*) 'KD_TREE_TRANS: is not an appropriate data structure.' - call IO_error (500_pInt) + call IO_error (460_pInt) end if call build_tree(mr) @@ -1356,7 +1356,7 @@ subroutine validate_query_storage(n) integer(pInt), intent(in) :: n 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 return diff --git a/code/math.f90 b/code/math.f90 index b9af3cb6b..ab7ff6f7b 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -3818,61 +3818,59 @@ subroutine calculate_cauchy(res,defgrad,p_stress,c_stress) 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 IO, only: & + IO_error implicit none ! input variables - integer(pInt), dimension(3), intent(in) :: res - integer(pInt), intent(in) :: domainPoints - real(pReal), dimension(3), intent(in) :: geomdim - real(pReal), dimension(3,3), intent(in) :: defgradAv - real(pReal), dimension(res(1),res(2),res(3),3), intent(in) :: querySet - real(pReal), dimension(domainPoints,3), intent(in) :: domainSet + integer(pInt), intent(in) :: spatialDim + real(pReal), dimension(3,3), intent(in) :: Favg + real(pReal), dimension(3), intent(in) :: geomdim + integer(pInt), intent(in) :: domainPoints + integer(pInt), intent(in) :: queryPoints + real(pReal), dimension(queryPoints,spatialDim), intent(in) :: querySet + real(pReal), dimension(domainPoints,spatialDim), intent(in) :: domainSet ! 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 - real(pReal), dimension(:,:), allocatable :: querySetLarge - integer(pInt) :: i,j,k, l,m,n, ielem_large, spatial_dim - real(pReal), dimension(3) :: shift + integer(pInt) :: i,j, l,m,n type(kdtree2), pointer :: tree type(kdtree2_result), dimension(1) :: Results + + 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') - shift = math_mul33x3(defgradAv,geomdim) - - ielem_large = 0_pInt - if(res(3) == 1_pInt) then - spatial_dim = 2_pInt - allocate(querySetLarge(2,(res(1)*res(2))*9_pInt)) - do j=1_pInt, res(2); do i=1_pInt, res(1) - do l = -1, 1; do m = -1, 1 - ielem_large = ielem_large + 1_pInt - querySetLarge(1:2,ielem_large) = querySet(i,j,1,1:2) + real([l,m],pReal)* shift(1:2) - enddo; enddo; - enddo; enddo - else - allocate(querySetLarge(3,(res(1)*res(2)*res(3))*27)) - spatial_dim = 3_pInt - do k=1_pInt,res(3); do j=1_pInt, res(2); do i=1_pInt, res(1) - do l = -1, 1; do m = -1, 1; do n = -1, 1 - 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 - endif - - tree => kdtree2_create(querySetLarge,sort=.true.,rearrange=.true.) - ielem_large = 0_pInt - do k=1_pInt,res(3); do j=1_pInt, res(2); do i=1_pInt, res(1) - ielem_large = ielem_large + 1_pInt - call kdtree2_n_nearest(tp=tree, qv=domainSet(ielem_large,1:spatial_dim),nn=1_pInt, results = Results) - indices(ielem_large) = Results(1)%idx !/3_pInt**spatial_dim +1_pInt - enddo; enddo; enddo - - deallocate(querySetLarge) + i = 0_pInt + if(spatialDim == 2_pInt) then + do j = 1_pInt, domainPoints + do l = -1_pInt, 1_pInt; do m = -1_pInt, 1_pInt + i = i + 1_pInt + domainSetLarge(1:3,i) = domainSet(j,1:3) + math_mul33x3(Favg,real([l,m,0_pInt],pReal)*geomdim) + enddo; enddo + enddo + else + do j = 1_pInt, domainPoints + do l = -1_pInt, 1_pInt; do m = -1_pInt, 1_pInt; do n = -1_pInt, 1_pInt + i = i + 1_pInt + domainSetLarge(1:3,i) = domainSet(j,1:3) + math_mul33x3(Favg,real([l,m,n],pReal)*geomdim) + enddo; enddo; enddo + enddo + endif + + tree => kdtree2_create(domainSetLarge,sort=.true.,rearrange=.true.) + + do j = 1_pInt, queryPoints + call kdtree2_n_nearest(tp=tree, qv=querySet(j,1:spatialDim),nn=1_pInt, results = Results) + indices(j) = Results(1)%idx + enddo + end subroutine math_nearestNeighborSearch