From 4b18278781dcb5a8ea961ba34ce37c621735e0d4 Mon Sep 17 00:00:00 2001 From: Krishna Komerla Date: Fri, 4 May 2012 13:07:37 +0000 Subject: [PATCH] reworked neares neigbor search subroutine, now also works for voronoi tesselation (general nn search in periodic cell) --- code/damask.core.pyf | 25 ++++++++-------- code/math.f90 | 51 ++++++++++++++++--------------- code/mesh.f90 | 71 ++++++++++++++++++++++++++++++++++++++++---- 3 files changed, 103 insertions(+), 44 deletions(-) diff --git a/code/damask.core.pyf b/code/damask.core.pyf index e7f97f50a..ec29502f8 100644 --- a/code/damask.core.pyf +++ b/code/damask.core.pyf @@ -137,7 +137,7 @@ python module core ! in subroutine logstrain_spat(res,defgrad,logstrain_field) ! in :math:math.f90 ! input variables - integer, dimension(3), intent(in) :: res + integer, dimension(3), intent(in) :: res real*8, dimension(res[0],res[1],res[2],3,3), intent(in), depend(res[0],res[1],res[2]) :: defgrad ! output variables real*8, dimension(res[0],res[1],res[2],3,3), intent(out), depend(res[0],res[1],res[2]) :: logstrain_field @@ -154,21 +154,22 @@ python module core ! in subroutine math_equivStrain33_field(res,tensor,vm) ! in :math:math.f90 ! input variables - integer, dimension(3), intent(in) :: res - real*8, dimension(res[0],res[1],res[2],3,3),intent(in),depend(res[0],res[1],res[2]) :: tensor + integer, dimension(3), intent(in) :: res + real*8, dimension(res[0],res[1],res[2],3,3), intent(in),depend(res[0],res[1],res[2]) :: tensor ! output variables 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_new,res_old,defgrad_av,geomdim,deformed_set,result_indices) ! in :math:math.f90 - ! input variables - integer, dimension(3), intent(in) :: res_new - integer, dimension(3), intent(in) :: res_old - real, dimension(3,3), intent(in) :: defgrad_av - real, dimension(3), intent(in) :: geomdim - real, dimension(res_old[0],res_old[1],res_old[2],3), intent(in), depend(res_old[0],res_old[1],res_old[2]) :: deformed_set - ! output variables - integer, dimension(res_new[0]*res_new[1]*res_new[2]), intent(out), depend(res_new[0],res_new[1],res_new[2]):: result_indices + subroutine math_nearestNeighborSearch(res,defgradAv,geomdim,,domainPoints,querySet,domainSet,indices) ! in :math:math.f90 + ! input variables + integer, dimension(3), intent(in) :: res + integer, intent(in) :: domainPoints + 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 + ! output variable + integer, dimension(res[0]*res[1]*res[2]), intent(out), depend(res[0],res[1],res[2]) :: indices end subroutine math_nearestNeighborSearch end module math diff --git a/code/math.f90 b/code/math.f90 index 75b8b9a93..b9af3cb6b 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -1,4 +1,4 @@ -! Copyright 2011 Max-Planck-Institut für Eisenforschung GmbH +! Copyright 2011 Max-Planck-Institut für Eisenforschung GmbH ! ! This file is part of DAMASK, ! the Düsseldorf Advanced MAterial Simulation Kit. @@ -3818,62 +3818,61 @@ subroutine calculate_cauchy(res,defgrad,p_stress,c_stress) end subroutine calculate_cauchy !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -subroutine math_nearestNeighborSearch(res_old, res_new, defgrad_av, geomdim, & - deformed_set, result_indices) +subroutine math_nearestNeighborSearch(res, defgradAv, geomdim, domainPoints, querySet, domainSet, indices) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !Obtain the nearest neighbour ! use kdtree2_module implicit none ! input variables - integer(pInt), dimension(3), intent(in) :: res_new, res_old + 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) :: defgrad_av - real(pReal), dimension( res_old(1), res_old(2),res_old(3),3), intent(in) :: deformed_set - ! other variables - integer(pInt), dimension(res_new(1)*res_new(2)*res_new(3)), intent(out) :: result_indices - ! other variables - real(pReal), dimension(:,:), allocatable :: deformed_set_large + 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 + ! output variable + integer(pInt), dimension(res(1)*res(2)*res(3)), intent(out) :: indices + ! other variables + real(pReal), dimension(:,:), allocatable :: querySetLarge integer(pInt) :: i,j,k, l,m,n, ielem_large, spatial_dim - real(pReal), dimension(3) :: shift, query_point + real(pReal), dimension(3) :: shift type(kdtree2), pointer :: tree type(kdtree2_result), dimension(1) :: Results - shift = math_mul33x3(defgrad_av,geomdim) + shift = math_mul33x3(defgradAv,geomdim) ielem_large = 0_pInt - if(res_old(3)==1_pInt) then + if(res(3) == 1_pInt) then spatial_dim = 2_pInt - allocate(deformed_set_large(2,(res_new(1)*res_new(2))*9_pInt)) - do j=1_pInt, res_old(2); do i=1_pInt, res_old(1) + 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 - deformed_set_large(1:2,ielem_large) = deformed_set(i,j,1,1:2) + real([l,m],pReal)* shift(1:2) + querySetLarge(1:2,ielem_large) = querySet(i,j,1,1:2) + real([l,m],pReal)* shift(1:2) enddo; enddo; enddo; enddo else - allocate(deformed_set_large(3,(res_new(1)*res_new(2)*res_new(3))*27)) + allocate(querySetLarge(3,(res(1)*res(2)*res(3))*27)) spatial_dim = 3_pInt - do k=1_pInt,res_old(3); do j=1_pInt, res_old(2); do i=1_pInt, res_old(1) + 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 - deformed_set_large(1:3,ielem_large) = deformed_set(i,j,k,1:3) + real([l,m,n],pReal)* shift + 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(deformed_set_large,sort=.true.,rearrange=.true.) + tree => kdtree2_create(querySetLarge,sort=.true.,rearrange=.true.) ielem_large = 0_pInt - do k=1_pInt,res_new(3); do j=1_pInt, res_new(2); do i=1_pInt, res_new(1) + do k=1_pInt,res(3); do j=1_pInt, res(2); do i=1_pInt, res(1) ielem_large = ielem_large + 1_pInt - query_point = math_mul33x3(defgrad_av, & - geomdim/real(res_new,pReal)*real([i,j,k],pReal) - geomdim/real(2_pInt*res_new,pReal)) - call kdtree2_n_nearest(tp=tree, qv=query_point(1:spatial_dim),nn=1_pInt, results = Results) - result_indices(ielem_large) = Results(1)%idx /3_pInt**spatial_dim +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(deformed_set_large) + deallocate(querySetLarge) end subroutine math_nearestNeighborSearch diff --git a/code/mesh.f90 b/code/mesh.f90 index 808223569..b1709dd40 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -3539,14 +3539,17 @@ subroutine mesh_regrid(resNew) !use new_res=[0.0,0.0,0.0] for automati IO_write_jobBinaryFile use math, only: & math_nearestNeighborSearch, & - deformed_FFT + deformed_FFT, & + math_mul33x3 integer(pInt), dimension(3), optional, intent(inout) :: resNew - integer(pInt):: maxsize, i, j, k, m, Npoints, NpointsNew + integer(pInt):: maxsize, i, j, k, ielem, Npoints, NpointsNew integer(pInt), dimension(3) :: res integer(pInt), dimension(:), allocatable :: indices, outputSize real(pReal), dimension(3) :: geomdim = 0.0_pReal real(pReal), dimension(3,3) :: Favg + real(pReal), dimension(:,:), allocatable :: & + coordinatesNew real(pReal), dimension(:,:,:,:,:), allocatable :: & F, FNew, & Fp, FpNew, & @@ -3559,7 +3562,7 @@ subroutine mesh_regrid(resNew) !use new_res=[0.0,0.0,0.0] for automati Tstar, TstarNew, & stateConst, stateConstNew, & stateHomog, stateHomogNew - + integer(pInt), dimension (:,:), allocatable :: SizeConst res = mesh_spectral_getResolution() geomdim = mesh_spectral_getDimension() @@ -3575,8 +3578,10 @@ subroutine mesh_regrid(resNew) !use new_res=[0.0,0.0,0.0] for automati print*, 'new resolution', resNew allocate(coordinates(res(1),res(2),res(3),3)) + allocate(coordinatesNew(resNew(1)*resNew(2)*resNew(3),3)) allocate(indices(Npoints)) - allocate(F(res(1),res(2),res(3),3,3)) + allocate(F(res(1),res(2),res(3),3,3)) + call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',trim(getSolverJobName()),size(F)) read (777,rec=1) F @@ -3589,10 +3594,21 @@ subroutine mesh_regrid(resNew) !use new_res=[0.0,0.0,0.0] for automati call deformed_fft(res,geomdim,Favg,1.0_pReal,F,coordinates) + ielem = 0_pInt + do k=1_pInt,resNew(3); do j=1_pInt, resNew(2); do i=1_pInt, resNew(1) + ielem = ielem + 1_pInt + coordinatesNew(ielem,1:3) = math_mul33x3(Favg, geomdim/real(resNew,pReal)*real([i,j,k],pReal) & + - geomdim/real(2_pInt*resNew,pReal)) + enddo; enddo; enddo + !----- Nearest neighbour search ----------------------- - call math_nearestNeighborSearch(res, res, Favg, geomdim, & - coordinates, indices) + call math_nearestNeighborSearch(res, Favg, geomdim, NpointsNew, & + coordinates, coordinatesNew, indices) + indices = indices /3_pInt**3_pInt +1_pInt ! adjust it for 2D case + deallocate(F) + deallocate(coordinates) + deallocate(coordinatesNew) ! --------------------------------------------------------------------------- allocate(F (3,3,1,1, Npoints)) @@ -3697,6 +3713,49 @@ subroutine mesh_regrid(resNew) !use new_res=[0.0,0.0,0.0] for automati close (777) deallocate(Tstar) deallocate(TstarNew) +!---------------------------------------------------------------------------- +!for homog and state const + ! allocate(SizeConst (1,Npoints)) + ! allocate(StateConst (1,Npoints)) + ! allocate(StateConstNew (1,Npoints)) + ! call IO_read_jobBinaryFile(777,'SizeConst',trim(getSolverJobName())) + ! read (777,rec=1) SizeConst + ! close(777) + ! maxsize = maxval(SizeConst) + ! call IO_read_jobBinaryFile(777,'convergedStateConst',trim(getSolverJobName())) + ! do i =1, Npoints + ! do l = 1,SizeConst(1,i) + ! read(777,rec=i) StateConst(1,1,i,l) + ! enddo + ! enddo + ! close(777) + + ! do i = 1, NpointsNew + ! do l = 1,SizeConst(1,i) + ! StateConstNew(1,1,i,l) = StateConst(1,1,indices(i),l) + ! enddo + ! enddo + + ! call IO_write_jobBinaryFile(777,'convergedStateConst',trim(getSolverJobName())) + ! m = 0_pInt + ! do i = 1,NpointsNew + ! do l = 1,SizeConst(1,i) + ! write(777,rec=m) StateConstNew(1,1,i,l) + ! enddo + ! enddo + ! close (777) + !---------------------------------------------------------------------------- + ! State homog + ! call IO_read_jobBinaryFile(777,'convergedStateHomog',FEmodelGeometry) + ! m = 0_pInt + ! do k = 1,Npoints; + ! do l = 1,homogenization_sizeState(j,k) + ! m = m+1_pInt + ! read(777,rec=m) stateHomog(1,k)%p(l) + ! enddo + ! enddo; enddo + ! close (777) + end subroutine mesh_regrid