From d8da2f60d870430fa63aceb5d4a37244f055f384 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 19 Feb 2013 14:56:26 +0000 Subject: [PATCH] added core module function math_periodicNearestNeighborDistances --- code/damask.core.pyf | 14 ++++++++- code/math.f90 | 66 ++++++++++++++++++++++++++++++++++++++---- lib/damask/__init__.py | 9 +++--- lib/kdtree2.f90 | 4 +-- 4 files changed, 81 insertions(+), 12 deletions(-) diff --git a/code/damask.core.pyf b/code/damask.core.pyf index 347d3b344..70634b485 100644 --- a/code/damask.core.pyf +++ b/code/damask.core.pyf @@ -85,6 +85,18 @@ python module core ! in real*8, dimension(size(domainSet,1),(3_pInt**size(domainSet,1))*size(domainSet,2)), depend(domainSet) :: domainSetLarge end function math_periodicNearestNeighbor + function math_periodicNearestNeighborDistances(geomdim,Favg,querySet,domainSet,Ndist) ! in :math:math.f90 + ! input variables + real*8, dimension(3), intent(in) :: geomdim + real*8, dimension(3,3), intent(in) :: Favg + integer, intent(in) :: Ndist + real*8, dimension(:,:), intent(in) :: querySet + real*8, dimension(:,:), intent(in) :: domainSet + real*8, dimension(Ndist,size(querySet,2)), depend(Ndist,querySet) :: math_periodicNearestNeighborDistances + ! depending on input + real*8, dimension(size(domainSet,1),(3_pInt**size(domainSet,1))*size(domainSet,2)), depend(domainSet) :: domainSetLarge + end function math_periodicNearestNeighborDistances + function math_tensorAvg(field) ! in :math:math.f90 ! input variables real*8 dimension(:,:,:,:,:), intent(in), :: field @@ -155,7 +167,7 @@ python module core ! in real*8, dimension(3), intent(in) :: gDim real*8, intent(in), optional, :: scalingIn = -1.0 real*8, dimension(3,3), intent(in), optional :: FavgIn = -1.0 - real*8, dimension(3,size(F,3),size(F,4),size(F,5)), depend(F) :: mesh_deformedCoordsFFT + real*8, dimension(3,size(F,3),size(F,4),size(F,5)), depend(F) :: mesh_deformedCoordsFFT end function mesh_deformedCoordsFFT function mesh_volumeMismatch(gDim,F,nodes) ! in :mesh:mesh.f90 diff --git a/code/math.f90 b/code/math.f90 index 994aaea27..c970adbce 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -232,7 +232,8 @@ real(pReal), dimension(4,36), parameter, private :: & math_logstrainSpat, & math_logstrainMat, & math_cauchy, & - math_periodicNearestNeighbor + math_periodicNearestNeighbor, & + math_periodicNearestNeighborDistances #endif private :: & math_partition, & @@ -1965,7 +1966,7 @@ end function math_sampleGaussVar !-------------------------------------------------------------------------------------------------- -!> @brief symmetric Euler angles for given symmetry string 'triclinic' or '', 'monoclinic', 'orthotropic' +!> @brief symmetric Euler angles for given symmetry 1:triclinic, 2:monoclinic, 4:orthotropic !-------------------------------------------------------------------------------------------------- pure function math_symmetricEulers(sym,Euler) @@ -3120,7 +3121,7 @@ function math_divergenceFDM(geomdim,order,field) end function math_divergenceFDM !-------------------------------------------------------------------------------------------------- -!> @brief Obtain the nearest neighbor in domain set for all points in querySet +!> @brief Obtain the nearest neighbor from domainSet at points in querySet !-------------------------------------------------------------------------------------------------- function math_periodicNearestNeighbor(geomdim, Favg, querySet, domainSet) use kdtree2_module @@ -3145,7 +3146,6 @@ function math_periodicNearestNeighbor(geomdim, Favg, querySet, domainSet) if (size(querySet,1) /= size(domainSet,1)) call IO_error(407_pInt,ext_msg='query set') spatialDim = size(querySet,1) - print*, geomdim i = 0_pInt if(spatialDim == 2_pInt) then do j = 1_pInt, size(domainSet,2) @@ -3170,9 +3170,65 @@ function math_periodicNearestNeighbor(geomdim, Favg, querySet, domainSet) math_periodicNearestNeighbor(j) = Results(1)%idx enddo math_periodicNearestNeighbor = math_periodicNearestNeighbor -1_pInt ! let them run from 0 to domainPoints -1 - print*, math_periodicNearestNeighbor end function math_periodicNearestNeighbor + + +!-------------------------------------------------------------------------------------------------- +!> @brief Obtain the distances to the next N nearest neighbors from domainSet at points in querySet +!-------------------------------------------------------------------------------------------------- +function math_periodicNearestNeighborDistances(geomdim, Favg, querySet, domainSet, Ndist) result(distances) + use kdtree2_module + use IO, only: & + IO_error + implicit none + ! input variables + real(pReal), dimension(3), intent(in) :: geomdim + real(pReal), dimension(3,3), intent(in) :: Favg + integer(pInt), intent(in) :: Ndist + real(pReal), dimension(:,:), intent(in) :: querySet + real(pReal), dimension(:,:), intent(in) :: domainSet + ! output variable + real(pReal), dimension(Ndist,size(querySet,2)) :: distances + + real(pReal), dimension(size(domainSet,1),(3_pInt**size(domainSet,1))*size(domainSet,2)) & + :: domainSetLarge + + integer(pInt) :: i,j, l,m,n, spatialDim + type(kdtree2), pointer :: tree + type(kdtree2_result), dimension(:), allocatable :: Results + + allocate(Results(Ndist)) + if (size(querySet,1) /= size(domainSet,1)) call IO_error(407_pInt,ext_msg='query set') + spatialDim = size(querySet,1) + + i = 0_pInt + if(spatialDim == 2_pInt) then + do j = 1_pInt, size(domainSet,2) + do l = -1_pInt, 1_pInt; do m = -1_pInt, 1_pInt + i = i + 1_pInt + domainSetLarge(1:2,i) = domainSet(1:2,j) +matmul(Favg(1:2,1:2),real([l,m],pReal)*geomdim(1:2)) + enddo; enddo + enddo + else + do j = 1_pInt, size(domainSet,2) + 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(1:3,j) + 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, size(querySet,2) + call kdtree2_n_nearest(tp=tree, qv=querySet(1:spatialDim,j),nn=Ndist, results = Results) + distances(1:Ndist,j) = sqrt(Results(1:Ndist)%dis) + enddo + + deallocate(Results) + +end function math_periodicNearestNeighborDistances #endif diff --git a/lib/damask/__init__.py b/lib/damask/__init__.py index a997a4193..8506533fe 100644 --- a/lib/damask/__init__.py +++ b/lib/damask/__init__.py @@ -19,7 +19,7 @@ try: core.IO = core.io core.FEsolving = core.fesolving core.DAMASK_interface = core.damask_interface -# remove XXX_ +# remove modulePrefix_ core.prec.init = core.prec.prec_init core.DAMASK_interface.init = core.DAMASK_interface.DAMASK_interface_init core.IO.init = core.IO.IO_init @@ -30,6 +30,7 @@ try: core.math.divergenceFFT = core.math.math_divergenceFFT core.math.divergenceFDM = core.math.math_divergenceFDM core.math.periodicNearestNeighbor = core.math.math_periodicNearestNeighbor + core.math.periodicNearestNeighborDistances = core.math.math_periodicNearestNeighborDistances core.math.tensorAvg = core.math.math_tensorAvg core.math.logstrainSpat = core.math.math_logstrainSpat core.math.logstrainMat = core.math.math_logstrainMat @@ -45,8 +46,8 @@ try: except (ImportError,AttributeError) as e: core = None # from http://www.python.org/dev/peps/pep-0008/ - if('setup_processing' not in sys.argv[0]): + if os.path.split(sys.argv[0])[1] not in ('symLink_Processing.py', + 'compile_CoreModule.py', + ): sys.stderr.write('\nWARNING: Core module (Fortran code) not available, '\ 'try to run setup_processing.py\nError Message when importing core.so: %s\n\n'%e) - - diff --git a/lib/kdtree2.f90 b/lib/kdtree2.f90 index f13967f08..2baae0bb0 100644 --- a/lib/kdtree2.f90 +++ b/lib/kdtree2.f90 @@ -1355,8 +1355,8 @@ subroutine validate_query_storage(n) implicit none integer(pInt), intent(in) :: n - if (int(size(sr%results,1),pInt) .lt. n) then - call IO_error (460_pInt,ext_msg='KD_TREE_TRANS: not enough storage for results(1:n)') + if (int(size(sr%results,1),pInt) < n) then + call IO_error (460_pInt,e=n,i=int(size(sr%results,1),pInt),ext_msg='KD_TREE_TRANS: not enough storage for results(1:n)') endif return