improved regridding, especially interfacing to python (also for nearest neighbor search)

This commit is contained in:
Krishna Komerla 2012-04-24 17:02:27 +00:00
parent e33c34d86b
commit 83e89fba3c
3 changed files with 220 additions and 64 deletions

View File

@ -1,7 +1,7 @@
! $Id$
! -*- f90 -*-
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Note: the context of this file is case sensitive.
! Note: the syntax of this file is case sensitive.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This file was auto-generated with f2py (version:2_5972).
! See http://cens.ioc.ee/projects/f2py2e/
@ -160,27 +160,22 @@ 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_new,Npoints_old,defgrad_av,geomdim,spatial_dim,deformed_set,result_indices) ! in :math:math.f90
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, intent(in) :: Npoints_old
integer, intent(in) :: spatial_dim
real, dimension(3), intent(in) :: geomdim
integer, dimension(3), intent(in) :: res_old
real, dimension(3,3), intent(in) :: defgrad_av
real, dimension(spatial_dim,Npoints_old), intent(in),depend(spatial_dim,Npoints_old) :: deformed_set
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
end subroutine math_nearestNeighborSearch
end module math
module mesh ! in :mesh:mesh.f90
subroutine mesh_regrid(res,resNew) ! in :mesh:mesh.f90
integer, dimension(3), intent(in) :: res
integer, dimension(3), intent(in,out) :: resNew
real*8, dimension(res[0],res[1],res[2],3,3), depend(res[0],res[1],res[2]) :: F
subroutine mesh_regrid(resNew) ! in :mesh:mesh.f90
integer, dimension(3), intent(in,out), optional :: resNew
end subroutine mesh_regrid
end module mesh

View File

@ -3817,49 +3817,64 @@ subroutine calculate_cauchy(res,defgrad,p_stress,c_stress)
end subroutine calculate_cauchy
!############################################################################
! subroutine to find nearest_neighbor.
!############################################################################
subroutine math_nearestNeighborSearch(res_new, Npoints_old, defgrad_av, geomdim, &
spatial_dim, deformed_set, result_indices)
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine math_nearestNeighborSearch(res_old, res_new, defgrad_av, geomdim, &
deformed_set, result_indices)
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!Obtain the nearest neighbour
!
use kdtree2_module
implicit none
integer(pInt), dimension(3), intent(in) :: res_new
integer(pInt), intent(in):: spatial_dim, Npoints_old
! input variables
integer(pInt), dimension(3), intent(in) :: res_new, res_old
real(pReal), dimension(3), intent(in) :: geomdim
real(pReal), dimension(3,3), intent(in) :: defgrad_av
real(pReal), dimension(spatial_dim,Npoints_old), intent(in) :: deformed_set
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
real(pReal), dimension(spatial_dim,Npoints_old*3_pInt**spatial_dim) :: deformed_set_large
integer(pInt):: i, j, k, ielem_small, ielem_large
! other variables
real(pReal), dimension(:,:), allocatable :: deformed_set_large
integer(pInt) :: i,j,k, l,m,n, ielem_large, spatial_dim
real(pReal), dimension(3) :: shift, query_point
type(kdtree2), pointer :: tree
type(kdtree2_result), dimension(1) :: Results
shift = math_mul33x3(defgrad_av,geomdim)
ielem_large = 0_pInt
do ielem_small=1_pInt, Npoints_old ! making copies (27 for 3D, 9 for 2D)
do k = -1, 1
do j = -1, 1
do i = -1, 1
ielem_large = 0_pInt
if(res_old(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)
do l = -1, 1; do m = -1, 1
ielem_large = ielem_large + 1_pInt
deformed_set_large(1:spatial_dim,ielem_large) = &
deformed_set(1:spatial_dim,ielem_small) + real([i,j,k],pReal)* shift
enddo; enddo; enddo; enddo
deformed_set_large(1:2,ielem_large) = deformed_set(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))
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 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
enddo; enddo; enddo;
enddo; enddo; enddo
endif
tree => kdtree2_create(deformed_set_large,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)
query_point = math_mul33x3(defgrad_av,(real([i,j,k],pReal)-0.5_pReal)/geomdim*real(res_new,pReal))
call kdtree2_n_nearest(tp=tree, qv=query_point(1_pInt:spatial_dim),nn=1_pInt, results = Results)
result_indices(i) = Results(1)%idx
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
enddo; enddo; enddo
deallocate(deformed_set_large)
end subroutine math_nearestNeighborSearch
end module math

View File

@ -3526,32 +3526,178 @@ deallocate(mesh_HomogMicro)
end subroutine mesh_tell_statistics
subroutine mesh_regrid(res,resNew) !use new_res=0.0 for automatic determination of new grid
subroutine mesh_regrid(resNew) !use new_res=[0.0,0.0,0.0] for automatic determination of new grid
use prec, only: &
pInt, &
pReal
use DAMASK_interface, only: &
getSolverJobName
use IO, only: &
IO_read_jobBinaryFile
IO_read_jobBinaryFile ,&
IO_write_jobBinaryFile
use math, only: &
math_nearestNeighborSearch, &
deformed_FFT
integer(pInt), dimension(3), intent(in) :: res
integer(pInt), dimension(3), intent(inout) :: resNew
real(pReal), dimension(res(1),res(2),res(3),3,3) :: F
integer(pInt), dimension(3), optional, intent(inout) :: resNew
integer(pInt):: maxsize, i, j, k, m, 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 :: &
F, FNew, &
Fp, FpNew, &
Lp, LpNew, &
dcsdE, dcsdENew
real(pReal), dimension (:,:,:,:,:,:,:), allocatable :: &
dPdF, dPdFNew
real(pReal), dimension (:,:,:,:), allocatable :: &
coordinates, &
Tstar, TstarNew, &
stateConst, stateConstNew, &
stateHomog, stateHomogNew
real(pReal), dimension(:,:,:,:,:), allocatable :: crystallite_F0, &
CPFEM_dcsdE, &
crystallite_Fp0, &
crystallite_Lp0
real(pReal), dimension (:,:,:,:,:,:,:), allocatable :: crystallite_dPdF0
real(pReal), dimension (:,:,:,:), allocatable :: crystallite_Tstar0_v, &
convergedStateConst
integer(pInt), dimension (:,:), allocatable :: convergedSizeConst
res = mesh_spectral_getResolution()
geomdim = mesh_spectral_getDimension()
if (.not. present(resNew)) then
resNew=res
endif
Npoints = res(1)*res(2)*res(3)
NpointsNew = resNew(1)*resNew(2)*resNew(3)
print*, 'resolution ', res
print*, 'new resolution', resNew
allocate(coordinates(res(1),res(2),res(3),3))
allocate(indices(Npoints))
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
close (777)
! ----Calculate average deformation for remesh--------
do i= 1_pInt,3_pInt; do j = 1_pInt,3_pInt
Favg(i,j) = sum(F(1:res(1),1:res(2),1:res(3),i,j)) / Npoints
enddo; enddo
call deformed_fft(res,geomdim,Favg,1.0_pReal,F,coordinates)
!----- Nearest neighbour search -----------------------
call math_nearestNeighborSearch(res, res, Favg, geomdim, &
coordinates, indices)
deallocate(F)
! ---------------------------------------------------------------------------
allocate(F (3,3,1,1, Npoints))
allocate(FNew(3,3,1,1, NpointsNew))
call IO_read_jobBinaryFile(777,'convergedF',trim(getSolverJobName()),size(F))
read (777,rec=1) F
close (777)
call flush(6)
do i = 1, NpointsNew
FNew(1:3,1:3,1,1,i) = F(1:3,1:3,1,1,indices(i))
enddo
call flush(6)
call IO_write_jobBinaryFile(777,'convergedF',size(FNew))
write (777,rec=1) FNew
close (777)
deallocate(F)
deallocate(FNew)
!---------------------------------------------------------------------
allocate(Fp (3,3,1,1,Npoints))
allocate(FpNew(3,3,1,1,NpointsNew))
call IO_read_jobBinaryFile(777,'convergedFp',trim(getSolverJobName()),size(Fp))
read (777,rec=1) Fp
close (777)
do i = 1, NpointsNew
FpNew(1:3,1:3,1,1,i) = Fp(1:3,1:3,1,1,indices(i))
enddo
call IO_write_jobBinaryFile(777,'convergedFp',size(FpNew))
write (777,rec=1) FpNew
close (777)
deallocate(Fp)
deallocate(FpNew)
!------------------------------------------------------------------------
allocate(Lp (3,3,1,1,Npoints))
allocate(LpNew(3,3,1,1,NpointsNew))
call IO_read_jobBinaryFile(777,'convergedLp',trim(getSolverJobName()),size(Lp))
read (777,rec=1) Lp
close (777)
do i = 1, NpointsNew
LpNew(1:3,1:3,1,1,i) = Lp(1:3,1:3,1,1,indices(i))
enddo
call IO_write_jobBinaryFile(777,'convergedLp',size(LpNew))
write (777,rec=1) LpNew
close (777)
deallocate(Lp)
deallocate(LpNew)
!----------------------------------------------------------------------------
allocate(dcsdE (6,6,1,1,Npoints))
allocate(dcsdENew(6,6,1,1,NpointsNew))
call IO_read_jobBinaryFile(777,'convergeddcsdE',trim(getSolverJobName()),size(dcsdE))
read (777,rec=1) dcsdE
close (777)
do i = 1, NpointsNew
dcsdENew(1:6,1:6,1,1,i) = dcsdE(1:6,1:6,1,1,indices(i))
enddo
call IO_write_jobBinaryFile(777,'convergeddcsdE',size(dcsdENew))
write (777,rec=1) dcsdENew
close (777)
deallocate(dcsdE)
deallocate(dcsdENew)
! ---------------------------------------------------------------------------
allocate(dPdF (3,3,3,3,1,1,Npoints))
allocate(dPdFNew(3,3,3,3,1,1,NpointsNew))
call IO_read_jobBinaryFile(777,'convergeddPdF',trim(getSolverJobName()),size(dPdF))
read (777,rec=1) dPdF
close (777)
do i = 1, NpointsNew
dPdFNew(1:3,1:3,1:3,1:3,1,1,i) = dPdF(1:3,1:3,1:3,1:3,1,1,indices(i))
enddo
call IO_write_jobBinaryFile(777,'convergeddPdF',size(dPdFNew))
write (777,rec=1) dPdFNew
close (777)
deallocate(dPdF)
deallocate(dPdFNew)
!---------------------------------------------------------------------------
allocate(Tstar (6,1,1,Npoints))
allocate(TstarNew(6,1,1,NpointsNew))
call IO_read_jobBinaryFile(777,'convergedTstar',trim(getSolverJobName()),size(Tstar))
read (777,rec=1) Tstar
close (777)
do i = 1, NpointsNew
TstarNew(1:6,1,1,i) = Tstar(1:6,1,1,indices(i))
enddo
call IO_write_jobBinaryFile(777,'convergedTstar',size(TstarNew))
write (777,rec=1) TstarNew
close (777)
deallocate(Tstar)
deallocate(TstarNew)
end subroutine mesh_regrid
function mesh_spectral_getDimension(fileUnit)