diff --git a/code/damask.core.pyf b/code/damask.core.pyf index b6be470ff..e7f97f50a 100644 --- a/code/damask.core.pyf +++ b/code/damask.core.pyf @@ -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 - real, dimension(3,3), intent(in) :: defgrad_av - real, dimension(spatial_dim,Npoints_old), intent(in),depend(spatial_dim,Npoints_old) :: deformed_set + 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 + 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 diff --git a/code/math.f90 b/code/math.f90 index ef9b22d41..75b8b9a93 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -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 + ! 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 - 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 - 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 = 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 - + 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: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.) - - 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 - enddo; enddo; enddo - + + 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) + 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 diff --git a/code/mesh.f90 b/code/mesh.f90 index ef09fd3c0..808223569 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -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 - - 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 + 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 + 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)