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

reworked neares neigbor search subroutine, now also works for voronoi tesselation (general nn search in periodic cell)

This commit is contained in:
Krishna Komerla 2012-05-04 13:07:37 +00:00
parent 863f0c766e
commit 4b18278781
3 changed files with 103 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
code/damask.core.pyf
View File

@ -160,15 +160,16 @@ 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_new,res_old,defgrad_av,geomdim,deformed_set,result_indices) ! in :math:math.f90 subroutine math_nearestNeighborSearch(res,defgradAv,geomdim,,domainPoints,querySet,domainSet,indices) ! in :math:math.f90
! input variables ! input variables
integer, dimension(3), intent(in) :: res_new integer, dimension(3), intent(in) :: res
integer, dimension(3), intent(in) :: res_old integer, intent(in) :: domainPoints
real, dimension(3,3), intent(in) :: defgrad_av real*8, dimension(3), intent(in) :: geomdim
real, dimension(3), intent(in) :: geomdim real*8, dimension(3,3), intent(in) :: defgradAv
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 real*8, dimension(res[0],res[1],res[2],3), intent(in), depend(res[0],res[1],res[2]) :: querySet
! output variables real*8, dimension(domainPoints,3), intent(in), depend(domainPoints) :: domainSet
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 ! 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 subroutine math_nearestNeighborSearch
end module math end module math

47
code/math.f90
View File

@ -3818,62 +3818,61 @@ subroutine calculate_cauchy(res,defgrad,p_stress,c_stress)
end subroutine calculate_cauchy end subroutine calculate_cauchy
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine math_nearestNeighborSearch(res_old, res_new, defgrad_av, geomdim, & subroutine math_nearestNeighborSearch(res, defgradAv, geomdim, domainPoints, querySet, domainSet, indices)
deformed_set, result_indices)
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!Obtain the nearest neighbour !Obtain the nearest neighbour
! !
use kdtree2_module use kdtree2_module
implicit none implicit none
! input variables ! 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), intent(in) :: geomdim
real(pReal), dimension(3,3), intent(in) :: defgrad_av real(pReal), dimension(3,3), intent(in) :: defgradAv
real(pReal), dimension( res_old(1), res_old(2),res_old(3),3), intent(in) :: deformed_set 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 ! other variables
integer(pInt), dimension(res_new(1)*res_new(2)*res_new(3)), intent(out) :: result_indices real(pReal), dimension(:,:), allocatable :: querySetLarge
! other variables
real(pReal), dimension(:,:), allocatable :: deformed_set_large
integer(pInt) :: i,j,k, l,m,n, ielem_large, spatial_dim 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), pointer :: tree
type(kdtree2_result), dimension(1) :: Results type(kdtree2_result), dimension(1) :: Results
shift = math_mul33x3(defgrad_av,geomdim) shift = math_mul33x3(defgradAv,geomdim)
ielem_large = 0_pInt ielem_large = 0_pInt
if(res_old(3)==1_pInt) then if(res(3) == 1_pInt) then
spatial_dim = 2_pInt spatial_dim = 2_pInt
allocate(deformed_set_large(2,(res_new(1)*res_new(2))*9_pInt)) allocate(querySetLarge(2,(res(1)*res(2))*9_pInt))
do j=1_pInt, res_old(2); do i=1_pInt, res_old(1) do j=1_pInt, res(2); do i=1_pInt, res(1)
do l = -1, 1; do m = -1, 1 do l = -1, 1; do m = -1, 1
ielem_large = ielem_large + 1_pInt 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;
enddo; enddo enddo; enddo
else 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 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 do l = -1, 1; do m = -1, 1; do n = -1, 1
ielem_large = ielem_large + 1_pInt 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;
enddo; enddo; enddo enddo; enddo; enddo
endif endif
tree => kdtree2_create(deformed_set_large,sort=.true.,rearrange=.true.) tree => kdtree2_create(querySetLarge,sort=.true.,rearrange=.true.)
ielem_large = 0_pInt 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 ielem_large = ielem_large + 1_pInt
query_point = math_mul33x3(defgrad_av, & call kdtree2_n_nearest(tp=tree, qv=domainSet(ielem_large,1:spatial_dim),nn=1_pInt, results = Results)
geomdim/real(res_new,pReal)*real([i,j,k],pReal) - geomdim/real(2_pInt*res_new,pReal)) indices(ielem_large) = Results(1)%idx !/3_pInt**spatial_dim +1_pInt
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
enddo; enddo; enddo enddo; enddo; enddo
deallocate(deformed_set_large) deallocate(querySetLarge)
end subroutine math_nearestNeighborSearch end subroutine math_nearestNeighborSearch

69
code/mesh.f90
View File

@ -3539,14 +3539,17 @@ subroutine mesh_regrid(resNew) !use new_res=[0.0,0.0,0.0] for automati
IO_write_jobBinaryFile IO_write_jobBinaryFile
use math, only: & use math, only: &
math_nearestNeighborSearch, & math_nearestNeighborSearch, &
deformed_FFT deformed_FFT, &
math_mul33x3
integer(pInt), dimension(3), optional, intent(inout) :: resNew 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(3) :: res
integer(pInt), dimension(:), allocatable :: indices, outputSize integer(pInt), dimension(:), allocatable :: indices, outputSize
real(pReal), dimension(3) :: geomdim = 0.0_pReal real(pReal), dimension(3) :: geomdim = 0.0_pReal
real(pReal), dimension(3,3) :: Favg real(pReal), dimension(3,3) :: Favg
real(pReal), dimension(:,:), allocatable :: &
coordinatesNew
real(pReal), dimension(:,:,:,:,:), allocatable :: & real(pReal), dimension(:,:,:,:,:), allocatable :: &
F, FNew, & F, FNew, &
Fp, FpNew, & Fp, FpNew, &
@ -3559,7 +3562,7 @@ subroutine mesh_regrid(resNew) !use new_res=[0.0,0.0,0.0] for automati
Tstar, TstarNew, & Tstar, TstarNew, &
stateConst, stateConstNew, & stateConst, stateConstNew, &
stateHomog, stateHomogNew stateHomog, stateHomogNew
integer(pInt), dimension (:,:), allocatable :: SizeConst
res = mesh_spectral_getResolution() res = mesh_spectral_getResolution()
geomdim = mesh_spectral_getDimension() geomdim = mesh_spectral_getDimension()
@ -3575,9 +3578,11 @@ subroutine mesh_regrid(resNew) !use new_res=[0.0,0.0,0.0] for automati
print*, 'new resolution', resNew print*, 'new resolution', resNew
allocate(coordinates(res(1),res(2),res(3),3)) allocate(coordinates(res(1),res(2),res(3),3))
allocate(coordinatesNew(resNew(1)*resNew(2)*resNew(3),3))
allocate(indices(Npoints)) 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)) call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',trim(getSolverJobName()),size(F))
read (777,rec=1) F read (777,rec=1) F
close (777) close (777)
@ -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) 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 ----------------------- !----- Nearest neighbour search -----------------------
call math_nearestNeighborSearch(res, res, Favg, geomdim, & call math_nearestNeighborSearch(res, Favg, geomdim, NpointsNew, &
coordinates, indices) coordinates, coordinatesNew, indices)
indices = indices /3_pInt**3_pInt +1_pInt ! adjust it for 2D case
deallocate(F) deallocate(F)
deallocate(coordinates)
deallocate(coordinatesNew)
! --------------------------------------------------------------------------- ! ---------------------------------------------------------------------------
allocate(F (3,3,1,1, Npoints)) 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) close (777)
deallocate(Tstar) deallocate(Tstar)
deallocate(TstarNew) 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 end subroutine mesh_regrid