From 96ba5ecae4ff40790a28dccfe53d5544ba20feef Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 27 Aug 2012 08:04:47 +0000 Subject: [PATCH] moved some more 'mesh related' functions for post processing from math.f90 to mesh.f90 f2py functions remaining in math.f90 now uses assumed size arrays in order to have simpler interfaces. This is only working with python 2.7! changed python pre- and postprocessing scripts. If you encounter any problems whith core modules, try to remove the old core.so in the lib/damask --- code/DAMASK_spectral.f90 | 6 +- code/DAMASK_spectral_Driver.f90 | 1 - code/DAMASK_spectral_SolverAL.f90 | 14 +- code/DAMASK_spectral_SolverBasic.f90 | 12 +- code/DAMASK_spectral_SolverBasicPETSC.f90 | 13 +- code/DAMASK_spectral_Utilities.f90 | 5 +- code/damask.core.pyf | 234 +++-- code/math.f90 | 997 ++++++-------------- code/mesh.f90 | 810 ++++++++++++---- lib/damask/__init__.py | 22 +- processing/post/3Dvisualize.py | 15 +- processing/post/addCompatibilityMismatch.py | 12 +- processing/post/addCurl.py | 2 +- processing/post/addDeformedConfiguration.py | 4 +- processing/post/addDivergence.py | 4 +- processing/post/spectral_buildElements.py | 6 +- 16 files changed, 1072 insertions(+), 1085 deletions(-) diff --git a/code/DAMASK_spectral.f90 b/code/DAMASK_spectral.f90 index cafcde85b..ea7429d56 100644 --- a/code/DAMASK_spectral.f90 +++ b/code/DAMASK_spectral.f90 @@ -83,7 +83,8 @@ program DAMASK_spectral mesh_NcpElems, & wgt, & geomdim, & - virt_dim + virt_dim, & + deformed_FFT use CPFEM, only: & CPFEM_general, & @@ -801,8 +802,7 @@ program DAMASK_spectral !-------------------------------------------------------------------------------------------------- ! call function to calculate divergence from math (for post processing) to check results if (debugDivergence) & - call divergence_fft(res,virt_dim,3_pInt,& - P_real(1:res(1),1:res(2),1:res(3),1:3,1:3),divergence_post) ! padding + divergence_post = math_divergenceFFT(virt_dim,P_real(1:res(1),1:res(2),1:res(3),1:3,1:3)) ! padding !-------------------------------------------------------------------------------------------------- ! doing the FT because it simplifies calculation of average stress in real space also diff --git a/code/DAMASK_spectral_Driver.f90 b/code/DAMASK_spectral_Driver.f90 index deea4e441..63c1bc92c 100644 --- a/code/DAMASK_spectral_Driver.f90 +++ b/code/DAMASK_spectral_Driver.f90 @@ -8,7 +8,6 @@ !-------------------------------------------------------------------------------------------------- program DAMASK_spectral_Driver use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) - use DAMASK_interface, only: & DAMASK_interface_init, & loadCaseFile, & diff --git a/code/DAMASK_spectral_SolverAL.f90 b/code/DAMASK_spectral_SolverAL.f90 index ff25ff2b6..0999a361c 100644 --- a/code/DAMASK_spectral_SolverAL.f90 +++ b/code/DAMASK_spectral_SolverAL.f90 @@ -7,9 +7,6 @@ !> @brief AL scheme solver !-------------------------------------------------------------------------------------------------- module DAMASK_spectral_SolverAL - - use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) - use prec, only: & pInt, & pReal @@ -81,8 +78,9 @@ module DAMASK_spectral_SolverAL !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields and fills them with data, potentially from restart info !-------------------------------------------------------------------------------------------------- - subroutine AL_init() - +subroutine AL_init() + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) + use IO, only: & IO_read_JobBinaryFile, & IO_write_JobBinaryFile @@ -224,11 +222,11 @@ module DAMASK_spectral_SolverAL update_gamma use math, only: & math_mul33x33 ,& - math_rotate_backward33, & - deformed_fft + math_rotate_backward33 use mesh, only: & res,& - geomdim + geomdim,& + deformed_fft use IO, only: & IO_write_JobBinaryFile diff --git a/code/DAMASK_spectral_SolverBasic.f90 b/code/DAMASK_spectral_SolverBasic.f90 index 15eefd1f6..a5f81585b 100644 --- a/code/DAMASK_spectral_SolverBasic.f90 +++ b/code/DAMASK_spectral_SolverBasic.f90 @@ -7,9 +7,6 @@ !> @brief Basic scheme solver !-------------------------------------------------------------------------------------------------- module DAMASK_spectral_SolverBasic - - use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) - use prec, only: & pInt, & pReal @@ -44,7 +41,8 @@ module DAMASK_spectral_SolverBasic !> @brief allocates all neccessary fields and fills them with data, potentially from restart info !-------------------------------------------------------------------------------------------------- subroutine basic_init() - + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) + use IO, only: & IO_read_JobBinaryFile, & IO_write_JobBinaryFile @@ -148,11 +146,11 @@ type(solutionState) function basic_solution(guessmode,timeinc,timeinc_old,P_BC,F math_mul33x33 ,& math_rotate_backward33, & math_transpose33, & - math_mul3333xx33, & - deformed_fft + math_mul3333xx33 use mesh, only: & res,& - geomdim + geomdim, & + deformed_fft use IO, only: & IO_write_JobBinaryFile diff --git a/code/DAMASK_spectral_SolverBasicPETSC.f90 b/code/DAMASK_spectral_SolverBasicPETSC.f90 index 2f835ea8f..ff1c2c79d 100644 --- a/code/DAMASK_spectral_SolverBasicPETSC.f90 +++ b/code/DAMASK_spectral_SolverBasicPETSC.f90 @@ -7,9 +7,6 @@ !> @brief Basic scheme PETSc solver !-------------------------------------------------------------------------------------------------- module DAMASK_spectral_SolverBasicPETSC - - use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) - use prec, only: & pInt, & pReal @@ -79,8 +76,10 @@ module DAMASK_spectral_SolverBasicPETSC !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields and fills them with data, potentially from restart info !-------------------------------------------------------------------------------------------------- - subroutine BasicPETSC_init() +subroutine BasicPETSC_init() + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) + use IO, only: & IO_read_JobBinaryFile, & IO_write_JobBinaryFile @@ -208,11 +207,11 @@ module DAMASK_spectral_SolverBasicPETSC update_gamma use math, only: & math_mul33x33 ,& - math_rotate_backward33, & - deformed_fft + math_rotate_backward33 use mesh, only: & res,& - geomdim + geomdim,& + deformed_fft use IO, only: & IO_write_JobBinaryFile diff --git a/code/DAMASK_spectral_Utilities.f90 b/code/DAMASK_spectral_Utilities.f90 index b0146ecfa..0d05f0359 100644 --- a/code/DAMASK_spectral_Utilities.f90 +++ b/code/DAMASK_spectral_Utilities.f90 @@ -7,7 +7,6 @@ !> @brief Utilities used by the different spectral solver variants !-------------------------------------------------------------------------------------------------- module DAMASK_spectral_Utilities - use prec, only: & pReal, & pInt @@ -78,7 +77,7 @@ contains !> Initializes FFTW. !-------------------------------------------------------------------------------------------------- subroutine Utilities_init() - + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) use numerics, only: & DAMASK_NumThreadsInt, & fftw_planner_flag, & @@ -258,7 +257,7 @@ subroutine Utilities_forwardFFT(row,column) !-------------------------------------------------------------------------------------------------- ! call function to calculate divergence from math (for post processing) to check results if (debugDivergence) & - call divergence_fft(res,virt_dim,3_pInt,field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),divergence_post) + divergence_post = math_divergenceFFT(virt_dim,field_real(1:res(1),1:res(2),1:res(3),1:3,1:3)) ! padding !-------------------------------------------------------------------------------------------------- ! doing the FT diff --git a/code/damask.core.pyf b/code/damask.core.pyf index 98ad5f63e..81d34d2e0 100644 --- a/code/damask.core.pyf +++ b/code/damask.core.pyf @@ -45,27 +45,99 @@ python module core ! in subroutine math_init end subroutine math_init - subroutine volume_compare(res,geomdim,defgrad,nodes,volume_mismatch) ! in :math:math.f90 + function math_curlFFT(geomdim,field) ! in :math:math.f90 ! input variables - integer, dimension(3), intent(in) :: res - real*8, dimension(3), intent(in) :: geomdim - real*8, dimension(res[0], res[1], res[2], 3,3),intent(in), depend(res[0],res[1],res[2]) :: defgrad - real*8, dimension(res[0]+1,res[1]+1,res[2]+1,3), intent(in), depend(res[0],res[1],res[2]) :: nodes - ! output variables - real*8, dimension(res[0], res[1], res[2]), intent(out), depend(res[0],res[1],res[2])) :: volume_mismatch - end subroutine volume_compare + real*8, dimension(3), intent(in) :: geomdim + real*8, dimension(:,:,:,:,:), intent(in), :: field + ! function definition + real*8, dimension(size(field,1),size(field,2),size(field,3),size(field,4),size(field,5)), depend(field) :: math_curlFFT + ! variables with dimension depending on input + real*8, dimension(size(field,1)/2+1,size(field,2),size(field,3),3), depend(field) :: xi + end function math_curlFFT - subroutine shape_compare(res,geomdim,defgrad,nodes,centroids,shape_mismatch) ! in :math:math.f90 + function math_divergenceFFT(geomdim,field) ! in :math:math.f90 ! input variables - integer, dimension(3), intent(in) :: res - real*8, dimension(3), intent(in) :: geomdim - real*8, dimension(res[0], res[1], res[2], 3,3),intent(in), depend(res[0],res[1],res[2]) :: defgrad - real*8, dimension(res[0]+1,res[1]+1,res[2]+1,3), intent(in), depend(res[0],res[1],res[2]) :: nodes - real*8, dimension(res[0], res[1], res[2], 3), intent(in), depend(res[0],res[1],res[2]) :: centroids - ! output variables - real*8, dimension(res[0], res[1], res[2]), intent(out), depend(res[0],res[1],res[2])) :: shape_mismatch - end subroutine shape_compare + real*8, dimension(3), intent(in) :: geomdim + real*8, dimension(:,:,:,:,:), intent(in), :: field + ! function definition + real*8, dimension(size(field,1),size(field,2),size(field,3),size(field,4)), depend(field) :: math_divergenceFFT + ! variables with dimension depending on input + real*8, dimension(size(field,1)/2+1,size(field,2),size(field,3),3), depend(field) :: xi + end function math_divergenceFFT + function math_divergenceFDM(geomdim,order,field) ! in :math:math.f90 + ! input variables + real*8 dimension(3), intent(in) :: geomdim + integer, intent(in) :: order + real*8, dimension(:,:,:,:,:), intent(in), :: field + ! function definition + real*8, dimension(size(field,1),size(field,2),size(field,3),size(field,4)), depend(field) :: math_divergenceFDM + end function math_divergenceFDM + + subroutine math_nearestNeighborSearch(spatialDim,Favg,geomdim,queryPoints,domainPoints,querySet,domainSet,indices) ! in :math:math.f90 + ! input variables + integer, intent(in) :: spatialDim + real*8, dimension(3,3), intent(in) :: Favg + real*8, dimension(3), intent(in) :: geomdim + integer, intent(in) :: queryPoints + integer, intent(in) :: domainPoints + real*8, dimension(spatialDim,queryPoints), intent(in), depend(spatialDim,queryPoints) :: querySet + real*8, dimension(spatialDim,domainPoints), intent(in), depend(spatialDim,domainPoints) :: domainSet + ! output variable + integer, dimension(queryPoints), intent(out), depend(queryPoints) :: indices + ! depending on input + real*8, dimension(spatialDim,(3**spatialDim)*domainPoints), depend(spatialDim,domainPoints) :: domainSetLarge + end subroutine math_nearestNeighborSearch + + function math_tensorAvg(field) ! in :math:math.f90 + ! input variables + real*8 dimension(:,:,:,:,:), intent(in), :: field + ! function definition + real*8 dimension(3,3), :: math_tensorAvg + end function math_tensorAvg + + function math_logstrainSpat(F) ! in :math:math.f90 + ! input variables + real*8, dimension(:,:,:,:,:), intent(in) :: F + ! output variables + real*8, dimension(size(F,1),size(F,2),size(F,3),3,3), depend(F) :: math_logstrainSpat + end function math_logstrainSpat + + function math_logstrainMat(F) ! in :math:math.f90 + ! input variables + real*8, dimension(:,:,:,:,:), intent(in) :: F + ! function + real*8, dimension(size(F,1),size(F,2),size(F,3),3,3), depend(F) :: math_logstrainMat + end function math_logstrainMat + + function math_cauchy(F,P) ! in :math:math.f90 + ! input variables + real*8, dimension(:,:,:,:,:), intent(in) :: F + real*8, dimension(:,:,:,:,:), intent(in) :: P + ! function + real*8, dimension(size(F,1),size(F,2),size(F,3),3,3), depend(F) :: math_cauchy + end function math_cauchy + + end module math + + module fesolving + subroutine FE_init + end subroutine FE_init + end module fesolving + + module mesh ! in :mesh:mesh.f90 + subroutine mesh_init(ip,element) + integer, parameter :: ip = 1 + integer, parameter :: element = 1 + end subroutine mesh_init + + function mesh_regrid(adaptive,resNewInput,minRes) ! in :mesh:mesh.f90 + logical, intent(in) :: adaptive + integer, dimension(3) :: mesh_regrid + integer, dimension(3), intent(in), optional :: resNewInput = -1 + integer, dimension(3), intent(in), optional :: minRes = -1 + end function mesh_regrid + subroutine mesh_regular_grid(res,geomdim,defgrad_av,centroids,nodes) ! in :math:math.f90 ! input variables integer, dimension(3), intent(in) :: res @@ -101,116 +173,28 @@ python module core ! in ! output variables real*8, dimension(res[0], res[1], res[2], 3), intent(out), depend(res[0],res[1],res[2]) :: coords end subroutine deformed_fft - - subroutine curl_fft(res,geomdim,vec_tens,field,curl) ! in :math:math.f90 - ! input variables - integer, dimension(3), intent(in) :: res - real*8, dimension(3), intent(in) :: geomdim - integer, intent(in) :: vec_tens - real*8, dimension(res[0], res[1], res[2], vec_tens,3), intent(in), depend(res[0],res[1],res[2],vec_tens) :: field - ! output variables - real*8, dimension(res[0], res[1], res[2], vec_tens,3), intent(out), depend(res[0],res[1],res[2],vec_tens) :: curl - ! variables with dimension depending on input - real*8, dimension(res[0]/2+1,res[1],res[2],3), depend(res[0],res[1],res[2]) :: xi - end subroutine curl_fft - - subroutine divergence_fft(res,geomdim,vec_tens,field,divergence) ! in :math:math.f90 - ! input variables - integer, dimension(3), intent(in) :: res - real*8, dimension(3), intent(in) :: geomdim - integer, intent(in) :: vec_tens - real*8, dimension(res[0], res[1], res[2], vec_tens,3), intent(in), depend(res[0],res[1],res[2],vec_tens) :: field - ! output variables - real*8, dimension(res[0], res[1], res[2], vec_tens), intent(out), depend(res[0],res[1],res[2],vec_tens) :: divergence - ! variables with dimension depending on input - real*8, dimension(res[0]/2+1,res[1],res[2],3), depend(res[0],res[1],res[2],3) :: xi - end subroutine divergence_fft - - subroutine divergence_fdm(res,geomdim,vec_tens,order,field,divergence) ! in :math:math.f90 - ! input variables - integer dimension(3), intent(in) :: res - real*8 dimension(3), intent(in) :: geomdim - integer intent(in) :: vec_tens - integer, intent(in) :: order - real*8 dimension(res[0], res[1], res[2], vec_tens,3), intent(in), depend(res[0],res[1],res[2],vec_tens) :: field - ! output variables - real*8 dimension(res[0], res[1], res[2], vec_tens), intent(out), depend(res[0],res[1],res[2],vec_tens) :: divergence - end subroutine divergence_fdm - - subroutine tensor_avg(res,tensor,avg) ! 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 - ! output variables - real*8 dimension(3,3), intent(out) :: avg - end subroutine tensor_avg - - subroutine logstrain_mat(res,defgrad,logstrain_field) ! 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]) :: 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 - end subroutine logstrain_mat - - subroutine logstrain_spat(res,defgrad,logstrain_field) ! 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]) :: 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 - end subroutine logstrain_spat - - subroutine calculate_cauchy(res,defgrad,p_stress,c_stress) ! 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]) :: defgrad - real*8, dimension(res[0],res[1],res[2],3,3), intent(in), depend(res[0],res[1],res[2]) :: p_stress - ! output variables - real*8, dimension(res[0],res[1],res[2],3,3), intent(out), depend(res[0],res[1],res[2]) :: c_stress - end subroutine calculate_cauchy - - 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 - ! 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(spatialDim,Favg,geomdim,queryPoints,domainPoints,querySet,domainSet,indices) ! in :math:math.f90 - ! input variables - integer, intent(in) :: spatialDim - real*8, dimension(3,3), intent(in) :: Favg - real*8, dimension(3), intent(in) :: geomdim - integer, intent(in) :: queryPoints - integer, intent(in) :: domainPoints - real*8, dimension(spatialDim,queryPoints), intent(in), depend(spatialDim,queryPoints) :: querySet - real*8, dimension(spatialDim,domainPoints), intent(in), depend(spatialDim,domainPoints) :: domainSet - ! output variable - integer, dimension(queryPoints), intent(out), depend(queryPoints) :: indices - ! depending on input - real*8, dimension(spatialDim,(3**spatialDim)*domainPoints), depend(spatialDim,domainPoints) :: domainSetLarge - end subroutine math_nearestNeighborSearch - end module math + subroutine volume_compare(res,geomdim,defgrad,nodes,volume_mismatch) ! in :math:math.f90 + ! input variables + integer, dimension(3), intent(in) :: res + real*8, dimension(3), intent(in) :: geomdim + real*8, dimension(res[0], res[1], res[2], 3,3),intent(in), depend(res[0],res[1],res[2]) :: defgrad + real*8, dimension(res[0]+1,res[1]+1,res[2]+1,3), intent(in), depend(res[0],res[1],res[2]) :: nodes + ! output variables + real*8, dimension(res[0], res[1], res[2]), intent(out), depend(res[0],res[1],res[2])) :: volume_mismatch + end subroutine volume_compare - module fesolving - subroutine FE_init - end subroutine FE_init - end module fesolving - - module mesh ! in :mesh:mesh.f90 - subroutine mesh_init(ip,element) - integer, parameter :: ip = 1 - integer, parameter :: element = 1 - end subroutine mesh_init - - function mesh_regrid(adaptive,resNewInput,minRes) ! in :mesh:mesh.f90 - logical, intent(in) :: adaptive - integer, dimension(3) :: mesh_regrid - integer, dimension(3), intent(in), optional :: resNewInput = -1 - integer, dimension(3), intent(in), optional :: minRes = -1 - end function mesh_regrid + subroutine shape_compare(res,geomdim,defgrad,nodes,centroids,shape_mismatch) ! in :math:math.f90 + ! input variables + integer, dimension(3), intent(in) :: res + real*8, dimension(3), intent(in) :: geomdim + real*8, dimension(res[0], res[1], res[2], 3,3),intent(in), depend(res[0],res[1],res[2]) :: defgrad + real*8, dimension(res[0]+1,res[1]+1,res[2]+1,3), intent(in), depend(res[0],res[1],res[2]) :: nodes + real*8, dimension(res[0], res[1], res[2], 3), intent(in), depend(res[0],res[1],res[2]) :: centroids + ! output variables + real*8, dimension(res[0], res[1], res[2]), intent(out), depend(res[0],res[1],res[2])) :: shape_mismatch + end subroutine shape_compare + end module mesh end interface end python module core diff --git a/code/math.f90 b/code/math.f90 index dbdccde6a..18f3dae0f 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -143,7 +143,45 @@ real(pReal), dimension(4,36), parameter, private :: & qsort, & math_range, & math_identity2nd, & - math_civita + math_civita, & + math_identity4th, & + math_vectorproduct, & + math_tensorproduct, & + math_mul3x3, & + math_mul6x6, & + math_mul33xx33, & + math_mul3333xx33, & + math_mul3333xx3333, & + math_mul33x33, & + math_mul66x66, & + math_mul99x99, & + math_mul33x3, & + math_mul33x3_complex, & + math_mul66x6 , & + math_qConj, & + math_qRot, & + math_transpose33, & + math_inv33, & + math_invert33, & + math_invSym3333, & + math_invert, & + math_symmetric33, & + math_symmetric66, & + math_skew33, & + math_deviatoric33, & + math_equivStrain33 +#ifdef Spectral + public :: math_curlFFT, & + math_divergenceFFT, & + math_divergenceFDM, & + math_tensorAvg, & + math_logstrainSpat, & + math_logstrainMat, & + math_cauchy, & + math_nearestNeighborSearch + private :: mesh_location, & + mesh_index +#endif private :: math_partition, & math_delta, & @@ -881,9 +919,9 @@ function math_invSym3333(A) end function math_invSym3333 -!************************************************************************** -! Gauss elimination to invert matrix of arbitrary dimension -!************************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief Gauss elimination to invert matrix of arbitrary dimension +!-------------------------------------------------------------------------------------------------- pure subroutine math_invert(dimen,A, InvA, AnzNegEW, error) ! Invertieren einer dimen x dimen - Matrix @@ -911,13 +949,12 @@ pure subroutine math_invert(dimen,A, InvA, AnzNegEW, error) end subroutine math_invert -! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +!-------------------------------------------------------------------------------------------------- +!> @brief Solves a linear EQS A * X = B with the GAUSS-Algorithm +! For numerical stabilization using a pivot search in rows and columns +!-------------------------------------------------------------------------------------------------- pure subroutine Gauss (dimen,A,B,LogAbsDetA,NegHDK,error) -! Solves a linear EQS A * X = B with the GAUSS-Algorithm -! For numerical stabilization using a pivot search in rows and columns -! ! input parameters ! A(dimen,dimen) = matrix A ! B(dimen,dimen) = right side B @@ -1155,42 +1192,10 @@ pure function math_equivStrain33(m) end function math_equivStrain33 -!******************************************************************** -subroutine math_equivStrain33_field(res,tensor,vm) -!******************************************************************** -!calculate von Mises equivalent of tensor field -! - implicit none - ! input variables - integer(pInt), intent(in), dimension(3) :: res - real(pReal), intent(in), dimension(res(1),res(2),res(3),3,3) :: tensor - ! output variables - real(pReal), intent(out), dimension(res(1),res(2),res(3)) :: vm - ! other variables - integer(pInt) :: i, j, k - real(pReal), dimension(3,3) :: deviator, delta = 0.0_pReal - real(pReal) :: J_2 - delta(1,1) = 1.0_pReal - delta(2,2) = 1.0_pReal - delta(3,3) = 1.0_pReal - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - deviator = tensor(i,j,k,1:3,1:3) - 1.0_pReal/3.0_pReal*tensor(i,j,k,1,1)*tensor(i,j,k,2,2)*tensor(i,j,k,3,3)*delta - J_2 = deviator(1,1)*deviator(2,2)& - + deviator(2,2)*deviator(3,3)& - + deviator(1,1)*deviator(3,3)& - - (deviator(1,2))**2.0_pReal& - - (deviator(2,3))**2.0_pReal& - - (deviator(1,3))**2.0_pReal - vm(i,j,k) = sqrt(3.0_pReal*J_2) - enddo; enddo; enddo - -end subroutine math_equivStrain33_field - - -!******************************************************************** -! determinant of a 33 matrix -!******************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief determinant of a 33 matrix +!-------------------------------------------------------------------------------------------------- pure function math_det33(m) implicit none @@ -1638,9 +1643,9 @@ pure function math_AxisAngleToR(axis,omega) end function math_AxisAngleToR -!**************************************************************** -! quaternion (w+ix+jy+kz) from axis and angle (in radians) -!**************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief quaternion (w+ix+jy+kz) from axis and angle (in radians) +!-------------------------------------------------------------------------------------------------- pure function math_AxisAngleToQuaternion(axis,omega) implicit none @@ -1750,9 +1755,9 @@ pure function math_QuaternionToAxisAngle(Q) end function math_QuaternionToAxisAngle -!******************************************************************** -! Rodrigues vector (x, y, z) from unit quaternion (w+ix+jy+kz) -!******************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief Rodrigues vector (x, y, z) from unit quaternion (w+ix+jy+kz) +!-------------------------------------------------------------------------------------------------- pure function math_QuaternionToRodrig(Q) use prec, only: DAMASK_NaN @@ -2048,9 +2053,9 @@ pure function math_symmetricEulers(sym,Euler) end function math_symmetricEulers -!******************************************************************** -! draw a random sample from Gauss variable -!******************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief draw a random sample from Gauss variable +!-------------------------------------------------------------------------------------------------- function math_sampleGaussVar(meanvalue, stddev, width) implicit none @@ -2299,8 +2304,8 @@ pure subroutine math_hi(M,HI1M,HI2M,HI3M) end subroutine math_hi -!******************************************************************************* -! GET_SEED returns a seed for the random number generator. +!-------------------------------------------------------------------------------------------------- +!> @brief GET_SEED returns a seed for the random number generator. ! ! The seed depends on the current time, and ought to be (slightly) ! different every millisecond. Once the seed is obtained, a random @@ -2315,7 +2320,7 @@ end subroutine math_hi ! ! Modified: 29 April 2005 ! Author: Franz Roters -! +!-------------------------------------------------------------------------------------------------- subroutine get_seed(seed) implicit none @@ -2424,7 +2429,7 @@ end subroutine halton ! ! Modified: 29 April 2005 ! Author: Franz Roters - +!-------------------------------------------------------------------------------------------------- subroutine halton_memory (action_halton, name_halton, ndim, value_halton) implicit none @@ -2591,7 +2596,7 @@ end subroutine halton_seed_set ! ! Modified: 29 April 2005 ! Author: Franz RotersA - +!-------------------------------------------------------------------------------------------------- subroutine i_to_halton (seed, base, ndim, r) use IO, only: IO_error @@ -2918,491 +2923,14 @@ pure function math_rotate_forward3333(tensor,rot_tensor) end function math_rotate_forward3333 -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -! Functions below are taken from the old postprocessingMath.f90 -! mostly they are used in combination with f2py to build fortran -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -! put the next two funtions into mesh? -function mesh_location(idx,resolution) - ! small helper functions for indexing - ! CAREFULL, index and location runs from 0 to N-1 (python style) - - integer(pInt), intent(in) :: idx - integer(pInt), intent(in), dimension(3) :: resolution - integer(pInt), dimension(3) :: mesh_location - mesh_location = (/modulo(idx/ resolution(3) / resolution(2),resolution(1)), & - modulo(idx/ resolution(3), resolution(2)), & - modulo(idx, resolution(3))/) - -end function mesh_location - - - function mesh_index(location,resolution) - ! small helper functions for indexing - ! CAREFULL, index and location runs from 0 to N-1 (python style) - integer(pInt), intent(in), dimension(3) :: resolution, location - integer(pInt) :: mesh_index - - mesh_index = modulo(location(3), resolution(3)) +& - (modulo(location(2), resolution(2)))*resolution(3) +& - (modulo(location(1), resolution(1)))*resolution(3)*resolution(2) - -end function mesh_index - - -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -subroutine volume_compare(res,geomdim,defgrad,nodes,volume_mismatch) -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -! Routine to calculate the mismatch between volume of reconstructed (compatible -! cube and determinant of defgrad at the FP - - use debug, only: debug_math, & - debug_level, & - debug_levelBasic - - implicit none - ! input variables - integer(pInt), intent(in), dimension(3) :: res - real(pReal), intent(in), dimension(3) :: geomdim - real(pReal), intent(in), dimension(res(1), res(2), res(3), 3,3) :: defgrad - real(pReal), intent(in), dimension(res(1)+1_pInt,res(2)+1_pInt,res(3)+1_pInt,3) :: nodes - ! output variables - real(pReal), intent(out), dimension(res(1), res(2), res(3)) :: volume_mismatch - ! other variables - real(pReal), dimension(8,3) :: coords - integer(pInt) i,j,k - real(pReal) vol_initial - - if (iand(debug_level(debug_math),debug_levelBasic) /= 0_pInt) then - print*, 'Calculating volume mismatch' - print '(a,3(e12.5))', ' Dimension: ', geomdim - print '(a,3(i5))', ' Resolution:', res - endif - - vol_initial = geomdim(1)*geomdim(2)*geomdim(3)/(real(res(1)*res(2)*res(3), pReal)) - do k = 1_pInt,res(3) - do j = 1_pInt,res(2) - do i = 1_pInt,res(1) - coords(1,1:3) = nodes(i, j, k ,1:3) - coords(2,1:3) = nodes(i+1_pInt,j, k ,1:3) - coords(3,1:3) = nodes(i+1_pInt,j+1_pInt,k ,1:3) - coords(4,1:3) = nodes(i, j+1_pInt,k ,1:3) - coords(5,1:3) = nodes(i, j, k+1_pInt,1:3) - coords(6,1:3) = nodes(i+1_pInt,j, k+1_pInt,1:3) - coords(7,1:3) = nodes(i+1_pInt,j+1_pInt,k+1_pInt,1:3) - coords(8,1:3) = nodes(i, j+1_pInt,k+1_pInt,1:3) - volume_mismatch(i,j,k) = abs(math_volTetrahedron(coords(7,1:3),coords(1,1:3),coords(8,1:3),coords(4,1:3))) & - + abs(math_volTetrahedron(coords(7,1:3),coords(1,1:3),coords(8,1:3),coords(5,1:3))) & - + abs(math_volTetrahedron(coords(7,1:3),coords(1,1:3),coords(3,1:3),coords(4,1:3))) & - + abs(math_volTetrahedron(coords(7,1:3),coords(1,1:3),coords(3,1:3),coords(2,1:3))) & - + abs(math_volTetrahedron(coords(7,1:3),coords(5,1:3),coords(2,1:3),coords(6,1:3))) & - + abs(math_volTetrahedron(coords(7,1:3),coords(5,1:3),coords(2,1:3),coords(1,1:3))) - volume_mismatch(i,j,k) = volume_mismatch(i,j,k)/math_det33(defgrad(i,j,k,1:3,1:3)) - enddo; enddo; enddo - volume_mismatch = volume_mismatch/vol_initial - -end subroutine volume_compare - - -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -subroutine shape_compare(res,geomdim,defgrad,nodes,centroids,shape_mismatch) -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -! Routine to calculate the mismatch between the vectors from the central point to -! the corners of reconstructed (combatible) volume element and the vectors calculated by deforming -! the initial volume element with the current deformation gradient - - use debug, only: debug_math, & - debug_level, & - debug_levelBasic - - implicit none - ! input variables - integer(pInt), intent(in), dimension(3) :: res - real(pReal), intent(in), dimension(3) :: geomdim - real(pReal), intent(in), dimension(res(1), res(2), res(3), 3,3) :: defgrad - real(pReal), intent(in), dimension(res(1)+1_pInt,res(2)+1_pInt,res(3)+1_pInt,3) :: nodes - real(pReal), intent(in), dimension(res(1), res(2), res(3), 3) :: centroids - ! output variables - real(pReal), intent(out), dimension(res(1), res(2), res(3)) :: shape_mismatch - ! other variables - real(pReal), dimension(8,3) :: coords_initial - integer(pInt) i,j,k - - if (iand(debug_level(debug_math),debug_levelBasic) /= 0_pInt) then - print*, 'Calculating shape mismatch' - print '(a,3(e12.5))', ' Dimension: ', geomdim - print '(a,3(i5))', ' Resolution:', res - endif - - coords_initial(1,1:3) = (/-geomdim(1)/2.0_pReal/real(res(1),pReal),& - -geomdim(2)/2.0_pReal/real(res(2),pReal),& - -geomdim(3)/2.0_pReal/real(res(3),pReal)/) - coords_initial(2,1:3) = (/+geomdim(1)/2.0_pReal/real(res(1),pReal),& - -geomdim(2)/2.0_pReal/real(res(2),pReal),& - -geomdim(3)/2.0_pReal/real(res(3),pReal)/) - coords_initial(3,1:3) = (/+geomdim(1)/2.0_pReal/real(res(1),pReal),& - +geomdim(2)/2.0_pReal/real(res(2),pReal),& - -geomdim(3)/2.0_pReal/real(res(3),pReal)/) - coords_initial(4,1:3) = (/-geomdim(1)/2.0_pReal/real(res(1),pReal),& - +geomdim(2)/2.0_pReal/real(res(2),pReal),& - -geomdim(3)/2.0_pReal/real(res(3),pReal)/) - coords_initial(5,1:3) = (/-geomdim(1)/2.0_pReal/real(res(1),pReal),& - -geomdim(2)/2.0_pReal/real(res(2),pReal),& - +geomdim(3)/2.0_pReal/real(res(3),pReal)/) - coords_initial(6,1:3) = (/+geomdim(1)/2.0_pReal/real(res(1),pReal),& - -geomdim(2)/2.0_pReal/real(res(2),pReal),& - +geomdim(3)/2.0_pReal/real(res(3),pReal)/) - coords_initial(7,1:3) = (/+geomdim(1)/2.0_pReal/real(res(1),pReal),& - +geomdim(2)/2.0_pReal/real(res(2),pReal),& - +geomdim(3)/2.0_pReal/real(res(3),pReal)/) - coords_initial(8,1:3) = (/-geomdim(1)/2.0_pReal/real(res(1),pReal),& - +geomdim(2)/2.0_pReal/real(res(2),pReal),& - +geomdim(3)/2.0_pReal/real(res(3),pReal)/) - do i=1_pInt,8_pInt - enddo - do k = 1_pInt,res(3) - do j = 1_pInt,res(2) - do i = 1_pInt,res(1) - shape_mismatch(i,j,k) = & - sqrt(sum((nodes(i, j, k, 1:3) - centroids(i,j,k,1:3)& - - matmul(defgrad(i,j,k,1:3,1:3), coords_initial(1,1:3)))**2.0_pReal))& - + sqrt(sum((nodes(i+1_pInt,j, k, 1:3) - centroids(i,j,k,1:3)& - - matmul(defgrad(i,j,k,1:3,1:3), coords_initial(2,1:3)))**2.0_pReal))& - + sqrt(sum((nodes(i+1_pInt,j+1_pInt,k, 1:3) - centroids(i,j,k,1:3)& - - matmul(defgrad(i,j,k,1:3,1:3), coords_initial(3,1:3)))**2.0_pReal))& - + sqrt(sum((nodes(i, j+1_pInt,k, 1:3) - centroids(i,j,k,1:3)& - - matmul(defgrad(i,j,k,1:3,1:3), coords_initial(4,1:3)))**2.0_pReal))& - + sqrt(sum((nodes(i, j, k+1_pInt,1:3) - centroids(i,j,k,1:3)& - - matmul(defgrad(i,j,k,1:3,1:3), coords_initial(5,1:3)))**2.0_pReal))& - + sqrt(sum((nodes(i+1_pInt,j, k+1_pInt,1:3) - centroids(i,j,k,1:3)& - - matmul(defgrad(i,j,k,1:3,1:3), coords_initial(6,1:3)))**2.0_pReal))& - + sqrt(sum((nodes(i+1_pInt,j+1_pInt,k+1_pInt,1:3) - centroids(i,j,k,1:3)& - - matmul(defgrad(i,j,k,1:3,1:3), coords_initial(7,1:3)))**2.0_pReal))& - + sqrt(sum((nodes(i, j+1_pInt,k+1_pInt,1:3) - centroids(i,j,k,1:3)& - - matmul(defgrad(i,j,k,1:3,1:3), coords_initial(8,1:3)))**2.0_pReal)) - enddo; enddo; enddo - -end subroutine shape_compare - - -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -subroutine mesh_regular_grid(res,geomdim,defgrad_av,centroids,nodes) -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -! Routine to build mesh of (distorted) cubes for given coordinates (= center of the cubes) -! - use debug, only: debug_math, & - debug_level, & - debug_levelBasic - - implicit none - ! input variables - integer(pInt), intent(in), dimension(3) :: res - real(pReal), intent(in), dimension(3) :: geomdim - real(pReal), intent(in), dimension(3,3) :: defgrad_av - real(pReal), intent(in), dimension(res(1), res(2), res(3), 3) :: centroids - ! output variables - real(pReal),intent(out), dimension(res(1)+1_pInt,res(2)+1_pInt,res(3)+1_pInt,3) :: nodes - ! variables with dimension depending on input - real(pReal), dimension(res(1)+2_pInt,res(2)+2_pInt,res(3)+2_pInt,3) :: wrappedCentroids - ! other variables - integer(pInt) :: i,j,k,n - integer(pInt), dimension(3), parameter :: diag = 1_pInt - integer(pInt), dimension(3) :: shift = 0_pInt, lookup = 0_pInt, me = 0_pInt - integer(pInt), dimension(3,8) :: neighbor = reshape((/ & - 0_pInt, 0_pInt, 0_pInt, & - 1_pInt, 0_pInt, 0_pInt, & - 1_pInt, 1_pInt, 0_pInt, & - 0_pInt, 1_pInt, 0_pInt, & - 0_pInt, 0_pInt, 1_pInt, & - 1_pInt, 0_pInt, 1_pInt, & - 1_pInt, 1_pInt, 1_pInt, & - 0_pInt, 1_pInt, 1_pInt & - /), & - (/3,8/)) - - if (iand(debug_level(debug_math),debug_levelBasic) /= 0_pInt) then - print*, 'Meshing cubes around centroids' - print '(a,3(e12.5))', ' Dimension: ', geomdim - print '(a,3(i5))', ' Resolution:', res - endif - - nodes = 0.0_pReal - wrappedCentroids = 0.0_pReal - wrappedCentroids(2_pInt:res(1)+1_pInt,2_pInt:res(2)+1_pInt,2_pInt:res(3)+1_pInt,1:3) = centroids - - do k = 0_pInt,res(3)+1_pInt - do j = 0_pInt,res(2)+1_pInt - do i = 0_pInt,res(1)+1_pInt - if (k==0_pInt .or. k==res(3)+1_pInt .or. & ! z skin - j==0_pInt .or. j==res(2)+1_pInt .or. & ! y skin - i==0_pInt .or. i==res(1)+1_pInt ) then ! x skin - me = (/i,j,k/) ! me on skin - shift = sign(abs(res+diag-2_pInt*me)/(res+diag),res+diag-2_pInt*me) - lookup = me-diag+shift*res - wrappedCentroids(i+1_pInt,j+1_pInt,k+1_pInt,1:3) = & - centroids(lookup(1)+1_pInt,lookup(2)+1_pInt,lookup(3)+1_pInt,1:3) - & - matmul(defgrad_av, shift*geomdim) - endif - enddo; enddo; enddo - do k = 0_pInt,res(3) - do j = 0_pInt,res(2) - do i = 0_pInt,res(1) - do n = 1_pInt,8_pInt - nodes(i+1_pInt,j+1_pInt,k+1_pInt,1:3) = & - nodes(i+1_pInt,j+1_pInt,k+1_pInt,1:3) + wrappedCentroids(i+1_pInt+neighbor(1_pInt,n), & - j+1_pInt+neighbor(2,n), & - k+1_pInt+neighbor(3,n),1:3) - enddo; enddo; enddo; enddo - nodes = nodes/8.0_pReal - -end subroutine mesh_regular_grid - - -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -subroutine deformed_linear(res,geomdim,defgrad_av,defgrad,coord) -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -! Routine to calculate coordinates in current configuration for given defgrad -! using linear interpolation (blurres out high frequency defomation) -! - implicit none - ! input variables - integer(pInt), intent(in), dimension(3) :: res - real(pReal), intent(in), dimension(3) :: geomdim - real(pReal), intent(in), dimension(3,3) :: defgrad_av - real(pReal), intent(in), dimension( res(1),res(2),res(3),3,3) :: defgrad - ! output variables - real(pReal), intent(out), dimension( res(1),res(2),res(3),3) :: coord - ! variables with dimension depending on input - real(pReal), dimension( 8,res(1),res(2),res(3),3) :: coord_avgOrder - ! other variables - real(pReal), dimension(3) :: myStep, fones = 1.0_pReal, parameter_coords, negative, positive, offset_coords - integer(pInt), dimension(3) :: rear, init, ones = 1_pInt, oppo, me - integer(pInt) i, j, k, s, o - integer(pInt), dimension(3,8) :: corner = reshape([ & - 0_pInt, 0_pInt, 0_pInt,& - 1_pInt, 0_pInt, 0_pInt,& - 1_pInt, 1_pInt, 0_pInt,& - 0_pInt, 1_pInt, 0_pInt,& - 1_pInt, 1_pInt, 1_pInt,& - 0_pInt, 1_pInt, 1_pInt,& - 0_pInt, 0_pInt, 1_pInt,& - 1_pInt, 0_pInt, 1_pInt & - ],[3,8]) - integer(pInt), dimension(3,8) :: step = reshape([& - 1_pInt, 1_pInt, 1_pInt,& - -1_pInt, 1_pInt, 1_pInt,& - -1_pInt,-1_pInt, 1_pInt,& - 1_pInt,-1_pInt, 1_pInt,& - -1_pInt,-1_pInt,-1_pInt,& - 1_pInt,-1_pInt,-1_pInt,& - 1_pInt, 1_pInt,-1_pInt,& - -1_pInt, 1_pInt,-1_pInt & - ], [3,8]) - integer(pInt), dimension(3,6) :: order = reshape([ & - 1_pInt, 2_pInt, 3_pInt,& - 1_pInt, 3_pInt, 2_pInt,& - 2_pInt, 1_pInt, 3_pInt,& - 2_pInt, 3_pInt, 1_pInt,& - 3_pInt, 1_pInt, 2_pInt,& - 3_pInt, 2_pInt, 1_pInt & - ], [3,6]) - - coord_avgOrder = 0.0_pReal - - do s = 0_pInt, 7_pInt ! corners (from 0 to 7) - init = corner(:,s+1_pInt)*(res-ones) +ones - oppo = corner(:,mod((s+4_pInt),8_pInt)+1_pInt)*(res-ones) +ones - - do o=1_pInt,6_pInt ! orders (from 1 to 6) - coord = 0_pReal - do k = init(order(3,o)), oppo(order(3,o)), step(order(3,o),s+1_pInt) - rear(order(2,o)) = init(order(2,o)) - do j = init(order(2,o)), oppo(order(2,o)), step(order(2,o),s+1_pInt) - rear(order(1,o)) = init(order(1,o)) - do i = init(order(1,o)), oppo(order(1,o)), step(order(1,o),s+1_pInt) - me(order(1:3,o)) = [i,j,k] - if ( all(me==init)) then - coord(me(1),me(2),me(3),1:3) = geomdim * (matmul(defgrad_av,real(corner(1:3,s+1),pReal)) + & - matmul(defgrad(me(1),me(2),me(3),1:3,1:3),0.5_pReal*real(step(1:3,s+1_pInt)/res,pReal))) - else - myStep = (me-rear)*geomdim/res - coord(me(1),me(2),me(3),1:3) = coord(rear(1),rear(2),rear(3),1:3) + & - 0.5_pReal*matmul(defgrad(me(1),me(2),me(3),1:3,1:3) + & - defgrad(rear(1),rear(2),rear(3),1:3,1:3),myStep) - endif - rear = me - enddo; enddo; enddo - coord_avgOrder(s+1_pInt,1:res(1),1:res(2),1:res(3),1:3) = & - coord_avgOrder(s+1_pInt,1:res(1),1:res(2),1:res(3),1:3) + coord/6.0_pReal - enddo - offset_coords = coord_avgOrder(s+1,1,1,1,1:3) - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - coord_avgOrder(s+1,i,j,k,1:3) = coord_avgOrder(s+1,i,j,k,1:3) - offset_coords - enddo; enddo; enddo - enddo - - do k = 0_pInt, res(3)-1_pInt - do j = 0_pInt, res(2)-1_pInt - do i = 0_pInt, res(1)-1_pInt - parameter_coords = (2.0_pReal*real([i,j,k]+1,pReal)-real(res,pReal))/(real(res,pReal)) - positive = fones + parameter_coords - negative = fones - parameter_coords - coord(i+1_pInt,j+1_pInt,k+1_pInt,1:3)& - =(coord_avgOrder(1,i+1_pInt,j+1_pInt,k+1_pInt,1:3) *negative(1)*negative(2)*negative(3)& - + coord_avgOrder(2,i+1_pInt,j+1_pInt,k+1_pInt,1:3) *positive(1)*negative(2)*negative(3)& - + coord_avgOrder(3,i+1_pInt,j+1_pInt,k+1_pInt,1:3) *positive(1)*positive(2)*negative(3)& - + coord_avgOrder(4,i+1_pInt,j+1_pInt,k+1_pInt,1:3) *negative(1)*positive(2)*negative(3)& - + coord_avgOrder(5,i+1_pInt,j+1_pInt,k+1_pInt,1:3) *positive(1)*positive(2)*positive(3)& - + coord_avgOrder(6,i+1_pInt,j+1_pInt,k+1_pInt,1:3) *negative(1)*positive(2)*positive(3)& - + coord_avgOrder(7,i+1_pInt,j+1_pInt,k+1_pInt,1:3) *negative(1)*negative(2)*positive(3)& - + coord_avgOrder(8,i+1_pInt,j+1_pInt,k+1_pInt,1:3) *positive(1)*negative(2)*positive(3))*0.125_pReal - enddo; enddo; enddo - - offset_coords = matmul(defgrad(1,1,1,1:3,1:3),geomdim/real(res, pReal)/2.0_pReal) - coord(1,1,1,1:3) - do k = 1_pInt, res(3) - do j = 1_pInt, res(2) - do i = 1_pInt, res(1) - coord(i,j,k,1:3) = coord(i,j,k,1:3)+ offset_coords - enddo; enddo; enddo - -end subroutine deformed_linear #ifdef Spectral -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -subroutine deformed_fft(res,geomdim,defgrad_av,scaling,defgrad,coords) -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -! Routine to calculate coordinates in current configuration for given defgrad -! using integration in Fourier space (more accurate than deformed(...)) -! - use IO, only: IO_error - use numerics, only: fftw_timelimit, fftw_planner_flag - use debug, only: debug_math, & - debug_level, & - debug_levelBasic - - implicit none - ! input variables - integer(pInt), intent(in), dimension(3) :: res - real(pReal), intent(in), dimension(3) :: geomdim - real(pReal), intent(in), dimension(3,3) :: defgrad_av - real(pReal), intent(in) :: scaling - real(pReal), intent(in), dimension(res(1), res(2),res(3),3,3) :: defgrad - ! output variables - real(pReal), intent(out), dimension(res(1), res(2),res(3),3) :: coords -! allocatable arrays for fftw c routines - type(C_PTR) :: fftw_forth, fftw_back - type(C_PTR) :: coords_fftw, defgrad_fftw - real(pReal), dimension(:,:,:,:,:), pointer :: defgrad_real - complex(pReal), dimension(:,:,:,:,:), pointer :: defgrad_fourier - real(pReal), dimension(:,:,:,:), pointer :: coords_real - complex(pReal), dimension(:,:,:,:), pointer :: coords_fourier - ! other variables - integer(pInt) :: i, j, k, m, res1_red - integer(pInt), dimension(3) :: k_s - real(pReal), dimension(3) :: step, offset_coords, integrator - - integrator = geomdim / 2.0_pReal / pi ! see notes where it is used - - if (iand(debug_level(debug_math),debug_levelBasic) /= 0_pInt) then - print*, 'Restore geometry using FFT-based integration' - print '(a,3(e12.5))', ' Dimension: ', geomdim - print '(a,3(i5))', ' Resolution:', res - endif - - res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c) - step = geomdim/real(res, pReal) - - if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt) - call fftw_set_timelimit(fftw_timelimit) - defgrad_fftw = fftw_alloc_complex(int(res1_red *res(2)*res(3)*9_pInt,C_SIZE_T)) !C_SIZE_T is of type integer(8) - call c_f_pointer(defgrad_fftw, defgrad_real, [res(1)+2_pInt,res(2),res(3),3_pInt,3_pInt]) - call c_f_pointer(defgrad_fftw, defgrad_fourier,[res1_red ,res(2),res(3),3_pInt,3_pInt]) - coords_fftw = fftw_alloc_complex(int(res1_red *res(2)*res(3)*3_pInt,C_SIZE_T)) !C_SIZE_T is of type integer(8) - call c_f_pointer(coords_fftw, coords_real, [res(1)+2_pInt,res(2),res(3),3_pInt]) - call c_f_pointer(coords_fftw, coords_fourier, [res1_red ,res(2),res(3),3_pInt]) - fftw_forth = fftw_plan_many_dft_r2c(3_pInt,(/res(3),res(2) ,res(1)/),9_pInt,& ! dimensions , length in each dimension in reversed order - defgrad_real,(/res(3),res(2) ,res(1)+2_pInt/),& ! input data , physical length in each dimension in reversed order - 1_pInt, res(3)*res(2)*(res(1)+2_pInt),& ! striding , product of physical lenght in the 3 dimensions - defgrad_fourier,(/res(3),res(2) ,res1_red/),& - 1_pInt, res(3)*res(2)* res1_red,fftw_planner_flag) - - fftw_back = fftw_plan_many_dft_c2r(3_pInt,(/res(3),res(2) ,res(1)/),3_pInt,& - coords_fourier,(/res(3),res(2) ,res1_red/),& - 1_pInt, res(3)*res(2)* res1_red,& - coords_real,(/res(3),res(2) ,res(1)+2_pInt/),& - 1_pInt, res(3)*res(2)*(res(1)+2_pInt),fftw_planner_flag) - - - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - defgrad_real(i,j,k,1:3,1:3) = defgrad(i,j,k,1:3,1:3) ! ensure that data is aligned properly (fftw_alloc) - enddo; enddo; enddo - - call fftw_execute_dft_r2c(fftw_forth, defgrad_real, defgrad_fourier) - - !remove highest frequency in each direction - if(res(1)>1_pInt) & - defgrad_fourier( res(1)/2_pInt+1_pInt,1:res(2) ,1:res(3) ,& - 1:3,1:3) = cmplx(0.0_pReal,0.0_pReal,pReal) - if(res(2)>1_pInt) & - defgrad_fourier(1:res1_red ,res(2)/2_pInt+1_pInt,1:res(3) ,& - 1:3,1:3) = cmplx(0.0_pReal,0.0_pReal,pReal) - if(res(3)>1_pInt) & - defgrad_fourier(1:res1_red ,1:res(2) ,res(3)/2_pInt+1_pInt,& - 1:3,1:3) = cmplx(0.0_pReal,0.0_pReal,pReal) - - coords_fourier = cmplx(0.0_pReal,0.0_pReal,pReal) - do k = 1_pInt, res(3) - k_s(3) = k-1_pInt - if(k > res(3)/2_pInt+1_pInt) k_s(3) = k_s(3)-res(3) - do j = 1_pInt, res(2) - k_s(2) = j-1_pInt - if(j > res(2)/2_pInt+1_pInt) k_s(2) = k_s(2)-res(2) - do i = 1_pInt, res1_red - k_s(1) = i-1_pInt - do m = 1_pInt,3_pInt - coords_fourier(i,j,k,m) = sum(defgrad_fourier(i,j,k,m,1:3)*cmplx(0.0_pReal,real(k_s,pReal)*integrator,pReal)) - enddo - if (k_s(3) /= 0_pInt .or. k_s(2) /= 0_pInt .or. k_s(1) /= 0_pInt) & - coords_fourier(i,j,k,1:3) = coords_fourier(i,j,k,1:3) / real(-sum(k_s*k_s),pReal) -! if(i/=1_pInt) coords_fourier(i,j,k,1:3) = coords_fourier(i,j,k,1:3)& ! substituting division by (on the fly calculated) xi * 2pi * img by multiplication with reversed img/real part -! - defgrad_fourier(i,j,k,1:3,1)*cmplx(0.0_pReal,integrator(1)/real(k_s(1),pReal),pReal) -! if(j/=1_pInt) coords_fourier(i,j,k,1:3) = coords_fourier(i,j,k,1:3)& -! - defgrad_fourier(i,j,k,1:3,2)*cmplx(0.0_pReal,integrator(2)/real(k_s(2),pReal),pReal) -! if(k/=1_pInt) coords_fourier(i,j,k,1:3) = coords_fourier(i,j,k,1:3)& -! - defgrad_fourier(i,j,k,1:3,3)*cmplx(0.0_pReal,integrator(3)/real(k_s(3),pReal),pReal) - enddo; enddo; enddo - - call fftw_execute_dft_c2r(fftw_back,coords_fourier,coords_real) - coords_real = coords_real/real(res(1)*res(2)*res(3),pReal) - - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - coords(i,j,k,1:3) = coords_real(i,j,k,1:3) ! ensure that data is aligned properly (fftw_alloc) - enddo; enddo; enddo - - offset_coords = matmul(defgrad(1,1,1,1:3,1:3),step/2.0_pReal) - scaling*coords(1,1,1,1:3) - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - coords(i,j,k,1:3) = scaling*coords(i,j,k,1:3) + offset_coords + matmul(defgrad_av,& - (/step(1)*real(i-1_pInt,pReal),& - step(2)*real(j-1_pInt,pReal),& - step(3)*real(k-1_pInt,pReal)/)) - - enddo; enddo; enddo - - call fftw_destroy_plan(fftw_forth) - call fftw_destroy_plan(fftw_back) - call fftw_free(defgrad_fftw) - call fftw_free(coords_fftw) - -end subroutine deformed_fft - - -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -subroutine curl_fft(res,geomdim,vec_tens,field,curl) -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -! calculates curl field using differentation in Fourier space +!-------------------------------------------------------------------------------------------------- +!> @brief calculates curl field using differentation in Fourier space ! use vec_tens to decide if tensor (3) or vector (1) +!-------------------------------------------------------------------------------------------------- +function math_curlFFT(geomdim,field) use IO, only: IO_error use numerics, only: fftw_timelimit, fftw_planner_flag @@ -3412,14 +2940,12 @@ subroutine curl_fft(res,geomdim,vec_tens,field,curl) implicit none ! input variables - integer(pInt), intent(in), dimension(3) :: res real(pReal), intent(in), dimension(3) :: geomdim - integer(pInt), intent(in) :: vec_tens - real(pReal), intent(in), dimension(res(1), res(2),res(3),vec_tens,3) :: field - ! output variables - real(pReal), intent(out), dimension(res(1), res(2),res(3),vec_tens,3) :: curl + real(pReal), intent(in), dimension(:,:,:,:,:) :: field + ! function + real(pReal), dimension(size(field,1),size(field,2),size(field,3),size(field,4),size(field,5)) :: math_curlFFT ! variables with dimension depending on input - real(pReal), dimension(res(1)/2_pInt+1_pInt,res(2),res(3),3) :: xi + real(pReal), dimension(size(field,1)/2_pInt+1_pInt,size(field,2),size(field,3),3) :: xi ! allocatable arrays for fftw c routines type(C_PTR) :: fftw_forth, fftw_back type(C_PTR) :: field_fftw, curl_fftw @@ -3429,8 +2955,12 @@ subroutine curl_fft(res,geomdim,vec_tens,field,curl) complex(pReal), dimension(:,:,:,:,:), pointer :: curl_fourier ! other variables integer(pInt) i, j, k, l, res1_red - integer(pInt), dimension(3) :: k_s + integer(pInt), dimension(3) :: k_s,res real(pReal) :: wgt + integer(pInt) :: vec_tens + + res = [size(field,1),size(field,2),size(field,3)] + vec_tens = size(field,4) if (iand(debug_level(debug_math),debug_levelBasic) /= 0_pInt) then print*, 'Calculating curl of vector/tensor field' @@ -3504,27 +3034,22 @@ subroutine curl_fft(res,geomdim,vec_tens,field,curl) call fftw_execute_dft_c2r(fftw_back, curl_fourier, curl_real) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - curl(i,j,k,1:vec_tens,1:3) = curl_real(i,j,k,1:vec_tens,1:3) ! ensure that data is aligned properly (fftw_alloc) + math_curlFFT(i,j,k,1:vec_tens,1:3) = curl_real(i,j,k,1:vec_tens,1:3) ! ensure that data is aligned properly (fftw_alloc) enddo; enddo; enddo - curl = curl * wgt - call fftw_destroy_plan(fftw_forth); call fftw_destroy_plan(fftw_back) - call c_f_pointer(C_NULL_PTR, field_real, [res(1)+2_pInt,res(2),res(3),vec_tens,3_pInt]) ! let all pointers point on NULL-Type - call c_f_pointer(C_NULL_PTR, field_fourier,[res1_red ,res(2),res(3),vec_tens,3_pInt]) - call c_f_pointer(C_NULL_PTR, curl_real, [res(1)+2_pInt,res(2),res(3),vec_tens,3_pInt]) - call c_f_pointer(C_NULL_PTR, curl_fourier, [res1_red ,res(2),res(3),vec_tens,3_pInt]) - if(.not. (c_associated(C_LOC(field_real(1,1,1,1,1))) .and. c_associated(C_LOC(field_fourier(1,1,1,1,1)))))& ! Check if pointers are deassociated and free memory - call fftw_free(field_fftw) ! This procedure ensures that optimization do not mix-up lines, because a - if(.not.(c_associated(C_LOC(curl_real(1,1,1,1,1))) .and. c_associated(C_LOC(curl_fourier(1,1,1,1,1)))))& ! simple fftw_free(field_fftw) could be done immediately after the last line where field_fftw appears, e.g: - call fftw_free(curl_fftw) ! call c_f_pointer(field_fftw, field_fourier, [res1_red ,res(2),res(3),vec_tens,3]) -end subroutine curl_fft + math_curlFFT = math_curlFFT * wgt + call fftw_destroy_plan(fftw_forth) + call fftw_destroy_plan(fftw_back) + call fftw_free(field_fftw) + call fftw_free(curl_fftw) +end function math_curlFFT -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -subroutine divergence_fft(res,geomdim,vec_tens,field,divergence) -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -! calculates divergence field using integration in Fourier space +!-------------------------------------------------------------------------------------------------- +!> @brief calculates divergence field using integration in Fourier space ! use vec_tens to decide if tensor (3) or vector (1) +!-------------------------------------------------------------------------------------------------- +function math_divergenceFFT(geomdim,field) use IO, only: IO_error use numerics, only: fftw_timelimit, fftw_planner_flag @@ -3534,14 +3059,13 @@ subroutine divergence_fft(res,geomdim,vec_tens,field,divergence) implicit none ! input variables - integer(pInt), intent(in), dimension(3) :: res + real(pReal), intent(in), dimension(3) :: geomdim - integer(pInt), intent(in) :: vec_tens - real(pReal), intent(in), dimension(res(1), res(2),res(3),vec_tens,3) :: field - ! output variables - real(pReal), intent(out), dimension(res(1), res(2),res(3),vec_tens) :: divergence + real(pReal), intent(in), dimension(:,:,:,:,:) :: field + ! function + real(pReal), dimension(size(field,1),size(field,2),size(field,3),size(field,4)) :: math_divergenceFFT ! variables with dimension depending on input - real(pReal), dimension(res(1)/2_pInt+1_pInt,res(2),res(3),3) :: xi + real(pReal), dimension(size(field,1)/2_pInt+1_pInt,size(field,2),size(field,3),3) :: xi ! allocatable arrays for fftw c routines type(C_PTR) :: fftw_forth, fftw_back type(C_PTR) :: field_fftw, divergence_fftw @@ -3552,7 +3076,11 @@ subroutine divergence_fft(res,geomdim,vec_tens,field,divergence) ! other variables integer(pInt) :: i, j, k, l, res1_red real(pReal) :: wgt - integer(pInt), dimension(3) :: k_s + integer(pInt), dimension(3) :: k_s, res + integer(pInt) :: vec_tens + + res = [size(field,1),size(field,2),size(field,3)] + vec_tens = size(field,4) if (iand(debug_level(debug_math),debug_levelBasic) /= 0_pInt) then print '(a)', 'Calculating divergence of tensor/vector field using FFT' @@ -3619,50 +3147,47 @@ if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt) call fftw_execute_dft_c2r(fftw_back, divergence_fourier, divergence_real) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - divergence(i,j,k,1:vec_tens) = divergence_real(i,j,k,1:vec_tens) ! ensure that data is aligned properly (fftw_alloc) + math_divergenceFFT(i,j,k,1:vec_tens) = divergence_real(i,j,k,1:vec_tens) ! ensure that data is aligned properly (fftw_alloc) enddo; enddo; enddo - divergence = divergence * wgt - call fftw_destroy_plan(fftw_forth); call fftw_destroy_plan(fftw_back) - call c_f_pointer(C_NULL_PTR, field_real, [res(1)+2_pInt,res(2),res(3),vec_tens,3_pInt]) ! let all pointers point on NULL-Type - call c_f_pointer(C_NULL_PTR, field_fourier, [res1_red ,res(2),res(3),vec_tens,3_pInt]) - call c_f_pointer(C_NULL_PTR, divergence_real, [res(1)+2_pInt,res(2),res(3),vec_tens]) - call c_f_pointer(C_NULL_PTR, divergence_fourier,[res1_red ,res(2),res(3),vec_tens]) - if(.not. (c_associated(C_LOC(field_real(1,1,1,1,1))) .and. c_associated(C_LOC(field_fourier(1,1,1,1,1)))))& ! Check if pointers are deassociated and free memory - call fftw_free(field_fftw) ! This procedure ensures that optimization do not mix-up lines, because a - if(.not.(c_associated(C_LOC(divergence_real(1,1,1,1))) .and. c_associated(C_LOC(divergence_fourier(1,1,1,1)))))& ! simple fftw_free(field_fftw) could be done immediately after the last line where field_fftw appears, e.g: - call fftw_free(divergence_fftw) ! call c_f_pointer(field_fftw, field_fourier, [res1_red ,res(2),res(3),vec_tens,3]) -end subroutine divergence_fft + math_divergenceFFT = math_divergenceFFT * wgt + call fftw_destroy_plan(fftw_forth) + call fftw_destroy_plan(fftw_back) + call fftw_free(field_fftw) + call fftw_free(divergence_fftw) +end function math_divergenceFFT -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -subroutine divergence_fdm(res,geomdim,vec_tens,order,field,divergence) -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -! calculates divergence field using FDM with variable accuracy +!-------------------------------------------------------------------------------------------------- +!> @brief calculates divergence field using FDM with variable accuracy ! use vec_tes to decide if tensor (3) or vector (1) +!-------------------------------------------------------------------------------------------------- +function math_divergenceFDM(geomdim,order,field) use debug, only: debug_math, & debug_level, & debug_levelBasic implicit none - integer(pInt), intent(in), dimension(3) :: res - integer(pInt), intent(in) :: vec_tens - integer(pInt), intent(inout) :: order + integer(pInt), intent(in) :: order real(pReal), intent(in), dimension(3) :: geomdim - real(pReal), intent(in), dimension(res(1),res(2),res(3),vec_tens,3) :: field - ! output variables - real(pReal), intent(out), dimension(res(1),res(2),res(3),vec_tens) :: divergence + real(pReal), intent(in), dimension(:,:,:,:,:) :: field + ! function + real(pReal), dimension(size(field,1),size(field,2),size(field,3),size(field,4)) :: math_divergenceFDM ! other variables integer(pInt), dimension(6,3) :: coordinates - integer(pInt) i, j, k, m, l - real(pReal), dimension(4,4), parameter :: FDcoefficient = reshape((/ & + integer(pInt) i, j, k, m, l, vec_tens + integer(pInt), dimension(3) :: res + real(pReal), dimension(4,4), parameter :: FDcoefficient = reshape([ & 1.0_pReal/2.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal,& !from http://en.wikipedia.org/wiki/Finite_difference_coefficients 2.0_pReal/3.0_pReal,-1.0_pReal/12.0_pReal, 0.0_pReal, 0.0_pReal,& 3.0_pReal/4.0_pReal,-3.0_pReal/20.0_pReal,1.0_pReal/ 60.0_pReal, 0.0_pReal,& - 4.0_pReal/5.0_pReal,-1.0_pReal/ 5.0_pReal,4.0_pReal/105.0_pReal,-1.0_pReal/280.0_pReal/),& - (/4,4/)) + 4.0_pReal/5.0_pReal,-1.0_pReal/ 5.0_pReal,4.0_pReal/105.0_pReal,-1.0_pReal/280.0_pReal],& + [4,4]) + + res = [size(field,1),size(field,2),size(field,3)] + vec_tens = size(field,4) if (iand(debug_level(debug_math),debug_levelBasic) /= 0_pInt) then print*, 'Calculating divergence of tensor/vector field using FDM' @@ -3670,10 +3195,9 @@ subroutine divergence_fdm(res,geomdim,vec_tens,order,field,divergence) print '(a,3(i5))', ' Resolution:', res endif - divergence = 0.0_pReal - order = order + 1_pInt + math_divergenceFDM = 0.0_pReal do k = 0_pInt, res(3)-1_pInt; do j = 0_pInt, res(2)-1_pInt; do i = 0_pInt, res(1)-1_pInt - do m = 1_pInt, order + do m = 1_pInt, order + 1_pInt coordinates(1,1:3) = mesh_location(mesh_index((/i+m,j,k/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/))& + (/1_pInt,1_pInt,1_pInt/) coordinates(2,1:3) = mesh_location(mesh_index((/i-m,j,k/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/))& @@ -3687,7 +3211,8 @@ subroutine divergence_fdm(res,geomdim,vec_tens,order,field,divergence) coordinates(6,1:3) = mesh_location(mesh_index((/i,j,k-m/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/))& + (/1_pInt,1_pInt,1_pInt/) do l = 1_pInt, vec_tens - divergence(i+1_pInt,j+1_pInt,k+1_pInt,l) = divergence(i+1_pInt,j+1_pInt,k+1_pInt,l) + FDcoefficient(m,order) * & + math_divergenceFDM(i+1_pInt,j+1_pInt,k+1_pInt,l) = math_divergenceFDM(i+1_pInt,j+1_pInt,k+1_pInt,l) & + + FDcoefficient(m,order) * & ((field(coordinates(1,1),coordinates(1,2),coordinates(1,3),l,1)- & field(coordinates(2,1),coordinates(2,2),coordinates(2,3),l,1))*real(res(1),pReal)/geomdim(1) +& (field(coordinates(3,1),coordinates(3,2),coordinates(3,3),l,2)- & @@ -3698,125 +3223,13 @@ subroutine divergence_fdm(res,geomdim,vec_tens,order,field,divergence) enddo enddo; enddo; enddo -end subroutine divergence_fdm -#endif +end function math_divergenceFDM - -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -subroutine tensor_avg(res,tensor,avg) -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -!calculate average of tensor field -! - implicit none - ! input variables - integer(pInt), intent(in), dimension(3) :: res - real(pReal), intent(in), dimension(res(1),res(2),res(3),3,3) ::tensor - ! output variables - real(pReal), intent(out), dimension(3,3) :: avg - ! other variables - real(pReal) wgt - integer(pInt) m,n - - wgt = 1.0_pReal/real(res(1)*res(2)*res(3), pReal) - - do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt - avg(m,n) = sum(tensor(1:res(1),1:res(2),1:res(3),m,n)) * wgt - enddo; enddo - -end subroutine tensor_avg - -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -subroutine logstrain_spat(res,defgrad,logstrain_field) -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -!calculate logarithmic strain in spatial configuration for given defgrad field -! - implicit none - ! input variables - integer(pInt), intent(in), dimension(3) :: res - real(pReal), intent(in), dimension(res(1),res(2),res(3),3,3) :: defgrad - ! output variables - real(pReal), intent(out), dimension(res(1),res(2),res(3),3,3) :: logstrain_field - ! other variables - real(pReal), dimension(3,3) :: temp33_Real, temp33_Real2 - real(pReal), dimension(3,3,3) :: eigenvectorbasis - real(pReal), dimension(3) :: eigenvalue - integer(pInt) :: i, j, k - logical :: errmatinv - - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - call math_pDecomposition(defgrad(i,j,k,1:3,1:3),temp33_Real2,temp33_Real,errmatinv) !store R in temp33_Real - temp33_Real2 = math_inv33(temp33_Real) - temp33_Real = math_mul33x33(defgrad(i,j,k,1:3,1:3),temp33_Real2) ! v = F o inv(R), store in temp33_Real2 - call math_spectral1(temp33_Real,eigenvalue(1), eigenvalue(2), eigenvalue(3),& - eigenvectorbasis(1,1:3,1:3),eigenvectorbasis(2,1:3,1:3),eigenvectorbasis(3,1:3,1:3)) - eigenvalue = log(sqrt(eigenvalue)) - logstrain_field(i,j,k,1:3,1:3) = eigenvalue(1)*eigenvectorbasis(1,1:3,1:3)+& - eigenvalue(2)*eigenvectorbasis(2,1:3,1:3)+& - eigenvalue(3)*eigenvectorbasis(3,1:3,1:3) - enddo; enddo; enddo - -end subroutine logstrain_spat - -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -subroutine logstrain_mat(res,defgrad,logstrain_field) -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -!calculate logarithmic strain in material configuration for given defgrad field -! - implicit none - ! input variables - integer(pInt), intent(in), dimension(3) :: res - real(pReal), intent(in), dimension(res(1),res(2),res(3),3,3) :: defgrad - ! output variables - real(pReal), intent(out), dimension(res(1),res(2),res(3),3,3) :: logstrain_field - ! other variables - real(pReal), dimension(3,3) :: temp33_Real, temp33_Real2 - real(pReal), dimension(3,3,3) :: eigenvectorbasis - real(pReal), dimension(3) :: eigenvalue - integer(pInt) :: i, j, k - logical :: errmatinv - - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - call math_pDecomposition(defgrad(i,j,k,1:3,1:3),temp33_Real,temp33_Real2,errmatinv) !store U in temp33_Real - call math_spectral1(temp33_Real,eigenvalue(1), eigenvalue(2), eigenvalue(3),& - eigenvectorbasis(1,1:3,1:3),eigenvectorbasis(2,1:3,1:3),eigenvectorbasis(3,1:3,1:3)) - eigenvalue = log(sqrt(eigenvalue)) - logstrain_field(i,j,k,1:3,1:3) = eigenvalue(1)*eigenvectorbasis(1,1:3,1:3)+& - eigenvalue(2)*eigenvectorbasis(2,1:3,1:3)+& - eigenvalue(3)*eigenvectorbasis(3,1:3,1:3) - enddo; enddo; enddo - -end subroutine logstrain_mat - -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -subroutine calculate_cauchy(res,defgrad,p_stress,c_stress) -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -!calculate cauchy stress for given PK1 stress and defgrad field -! - implicit none - ! input variables - integer(pInt), intent(in), dimension(3) :: res - real(pReal), intent(in), dimension(res(1),res(2),res(3),3,3) :: defgrad - real(pReal), intent(in), dimension(res(1),res(2),res(3),3,3) :: p_stress - ! output variables - real(pReal), intent(out), dimension(res(1),res(2),res(3),3,3) :: c_stress - ! other variables - real(pReal) :: jacobi - integer(pInt) :: i, j, k - - c_stress = 0.0_pReal - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - jacobi = math_det33(defgrad(i,j,k,1:3,1:3)) - c_stress(i,j,k,1:3,1:3) = matmul(p_stress(i,j,k,1:3,1:3),transpose(defgrad(i,j,k,1:3,1:3)))/jacobi - enddo; enddo; enddo - -end subroutine calculate_cauchy - -#ifdef Spectral -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +!-------------------------------------------------------------------------------------------------- +!> @brief Obtain the nearest neighbor in domain set for all points in querySet +!-------------------------------------------------------------------------------------------------- subroutine math_nearestNeighborSearch(spatialDim, Favg, geomdim, queryPoints, domainPoints, querySet, domainSet, indices) -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -! Obtain the nearest neighbor in domain set for all points in querySet -! + use kdtree2_module use IO, only: & IO_error @@ -3837,10 +3250,10 @@ subroutine math_nearestNeighborSearch(spatialDim, Favg, geomdim, queryPoints, do integer(pInt) :: i,j, l,m,n type(kdtree2), pointer :: tree type(kdtree2_result), dimension(1) :: Results - + if (size(querySet(:,1)) /= spatialDim) call IO_error(407_pInt,ext_msg='query set') if (size(domainSet(:,1)) /= spatialDim) call IO_error(407_pInt,ext_msg='domain set') - + i = 0_pInt if(spatialDim == 2_pInt) then @@ -3862,12 +3275,156 @@ subroutine math_nearestNeighborSearch(spatialDim, Favg, geomdim, queryPoints, do tree => kdtree2_create(domainSetLarge,sort=.true.,rearrange=.true.) do j = 1_pInt, queryPoints - call kdtree2_n_nearest(tp=tree, qv=querySet(1:spatialDim,j),nn=1_pInt, results = Results) + call kdtree2_n_nearest(tp=tree, qv=querySet(1:spatialDim,j),nn=1_pInt, results = Results) indices(j) = Results(1)%idx enddo indices = indices -1_pInt ! let them run from 0 to domainPoints -1 - + end subroutine math_nearestNeighborSearch + +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! Functions below are taken from the old postprocessingMath.f90 +! mostly they are used in combination with f2py to build fortran +! they now reside in mesh.f90 and should be removed from here soon +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +! put the next two funtions into mesh? +function mesh_location(idx,resolution) + ! small helper functions for indexing + ! CAREFULL, index and location runs from 0 to N-1 (python style) + + integer(pInt), intent(in) :: idx + integer(pInt), intent(in), dimension(3) :: resolution + integer(pInt), dimension(3) :: mesh_location + mesh_location = (/modulo(idx/ resolution(3) / resolution(2),resolution(1)), & + modulo(idx/ resolution(3), resolution(2)), & + modulo(idx, resolution(3))/) + +end function mesh_location + + + function mesh_index(location,resolution) + ! small helper functions for indexing + ! CAREFULL, index and location runs from 0 to N-1 (python style) + integer(pInt), intent(in), dimension(3) :: resolution, location + integer(pInt) :: mesh_index + + mesh_index = modulo(location(3), resolution(3)) +& + (modulo(location(2), resolution(2)))*resolution(3) +& + (modulo(location(1), resolution(1)))*resolution(3)*resolution(2) + +end function mesh_index #endif + +!-------------------------------------------------------------------------------------------------- +!> @brief calculate average of tensor field +!-------------------------------------------------------------------------------------------------- +function math_tensorAvg(field) + implicit none + ! input variables + real(pReal), intent(in), dimension(:,:,:,:,:) :: field + ! output variables + real(pReal), dimension(3,3) :: math_tensorAvg + ! other variables + real(pReal) :: wgt + integer(pInt) :: m,n + integer(pInt), dimension(3) :: res + + res = [size(field,1),size(field,2),size(field,3)] + + wgt = 1.0_pReal/real(res(1)*res(2)*res(3), pReal) + + do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt + math_tensorAvg(m,n) = sum(field(1:res(1),1:res(2),1:res(3),m,n)) * wgt + enddo; enddo + +end function math_tensorAvg + +!-------------------------------------------------------------------------------------------------- +!> @brief calculate logarithmic strain in spatial configuration for given F field +!-------------------------------------------------------------------------------------------------- +function math_logstrainSpat(F) + implicit none + ! input variables + real(pReal), intent(in), dimension(:,:,:,:,:) :: F + ! output variables + real(pReal) , dimension(size(F,1),size(F,2),size(F,3),3,3) :: math_logstrainSpat + ! other variables + integer(pInt), dimension(3) :: res + real(pReal), dimension(3,3) :: temp33_Real, temp33_Real2 + real(pReal), dimension(3,3,3) :: eigenvectorbasis + real(pReal), dimension(3) :: eigenvalue + integer(pInt) :: i, j, k + logical :: errmatinv + + res = [size(F,1),size(F,2),size(F,3)] + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + call math_pDecomposition(F(i,j,k,1:3,1:3),temp33_Real2,temp33_Real,errmatinv) !store R in temp33_Real + temp33_Real2 = math_inv33(temp33_Real) + temp33_Real = math_mul33x33(F(i,j,k,1:3,1:3),temp33_Real2) ! v = F o inv(R), store in temp33_Real2 + call math_spectral1(temp33_Real,eigenvalue(1), eigenvalue(2), eigenvalue(3),& + eigenvectorbasis(1,1:3,1:3),eigenvectorbasis(2,1:3,1:3),eigenvectorbasis(3,1:3,1:3)) + eigenvalue = log(sqrt(eigenvalue)) + math_logstrainSpat(i,j,k,1:3,1:3) = eigenvalue(1)*eigenvectorbasis(1,1:3,1:3)+& + eigenvalue(2)*eigenvectorbasis(2,1:3,1:3)+& + eigenvalue(3)*eigenvectorbasis(3,1:3,1:3) + enddo; enddo; enddo + +end function math_logstrainSpat + +!-------------------------------------------------------------------------------------------------- +!> @brief calculate logarithmic strain in material configuration for given F field +!-------------------------------------------------------------------------------------------------- +function math_logstrainMat(F) + implicit none + ! input variables + real(pReal), intent(in), dimension(:,:,:,:,:) :: F + ! output variables + real(pReal) , dimension(size(F,1),size(F,2),size(F,3),3,3) :: math_logstrainMat + ! other variables + integer(pInt), dimension(3) :: res + real(pReal), dimension(3,3) :: temp33_Real, temp33_Real2 + real(pReal), dimension(3,3,3) :: eigenvectorbasis + real(pReal), dimension(3) :: eigenvalue + integer(pInt) :: i, j, k + logical :: errmatinv + + res = [size(F,1),size(F,2),size(F,3)] + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + call math_pDecomposition(F(i,j,k,1:3,1:3),temp33_Real,temp33_Real2,errmatinv) !store U in temp33_Real + call math_spectral1(temp33_Real,eigenvalue(1), eigenvalue(2), eigenvalue(3),& + eigenvectorbasis(1,1:3,1:3),eigenvectorbasis(2,1:3,1:3),eigenvectorbasis(3,1:3,1:3)) + eigenvalue = log(sqrt(eigenvalue)) + math_logstrainMat(i,j,k,1:3,1:3) = eigenvalue(1)*eigenvectorbasis(1,1:3,1:3)+& + eigenvalue(2)*eigenvectorbasis(2,1:3,1:3)+& + eigenvalue(3)*eigenvectorbasis(3,1:3,1:3) + enddo; enddo; enddo + +end function math_logstrainMat + +!-------------------------------------------------------------------------------------------------- +!> @brief calculate cauchy stress for given PK1 stress and F field +!-------------------------------------------------------------------------------------------------- +function math_cauchy(F,P) + implicit none + ! input variables + real(pReal), intent(in), dimension(:,:,:,:,:) :: F + real(pReal), intent(in), dimension(:,:,:,:,:) :: P + ! output variables + real(pReal) , dimension(size(F,1),size(F,2),size(F,3),3,3) :: math_cauchy + ! other variables + integer(pInt), dimension(3) :: res + real(pReal) :: jacobi + integer(pInt) :: i, j, k + + res = [size(F,1),size(F,2),size(F,3)] + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + jacobi = math_det33(F(i,j,k,1:3,1:3)) + math_cauchy(i,j,k,1:3,1:3) = matmul(P(i,j,k,1:3,1:3),transpose(F(i,j,k,1:3,1:3)))/jacobi + enddo; enddo; enddo + +end function math_cauchy + end module math diff --git a/code/mesh.f90 b/code/mesh.f90 index e19e0ee8c..9eef62d99 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -28,6 +28,7 @@ !-------------------------------------------------------------------------------------------------- module mesh + use, intrinsic :: iso_c_binding use prec, only: pReal, pInt implicit none @@ -96,6 +97,7 @@ module mesh logical, private :: noPart !< for cases where the ABAQUS input file does not use part/assembly information #ifdef Spectral + include 'fftw3.f03' real(pReal), dimension(3), public :: geomdim, virt_dim ! physical dimension of volume element per direction integer(pInt), dimension(3), public :: res real(pReal), public :: wgt @@ -279,7 +281,12 @@ module mesh mesh_build_ipVolumes, & mesh_build_ipCoordinates #ifdef Spectral - public :: mesh_regrid + public :: mesh_regrid, & + mesh_regular_grid, & + deformed_linear, & + deformed_fft, & + volume_compare, & + shape_compare #endif private :: & @@ -858,6 +865,203 @@ function mesh_spectral_getHomogenization(fileUnit) call IO_error(error_ID = 842_pInt, ext_msg='mesh_spectral_getHomogenization') end function mesh_spectral_getHomogenization +!-------------------------------------------------------------------------------------------------- +!> @brief Count overall number of nodes and elements in mesh and stores them in +!! 'mesh_Nelems' and 'mesh_Nnodes' +!-------------------------------------------------------------------------------------------------- +subroutine mesh_spectral_count_nodesAndElements() + + implicit none + mesh_Nelems = res(1)*res(2)*res(3) + mesh_Nnodes = (1_pInt + res(1))*(1_pInt + res(2))*(1_pInt + res(3)) + +end subroutine mesh_spectral_count_nodesAndElements + + +!-------------------------------------------------------------------------------------------------- +!> @brief Count overall number of CP elements in mesh and stores them in 'mesh_NcpElems' +!-------------------------------------------------------------------------------------------------- +subroutine mesh_spectral_count_cpElements + + implicit none + + mesh_NcpElems = mesh_Nelems + +end subroutine mesh_spectral_count_cpElements + + +!-------------------------------------------------------------------------------------------------- +!> @brief Maps elements from FE ID to internal (consecutive) representation. +!! Allocates global array 'mesh_mapFEtoCPelem' +!-------------------------------------------------------------------------------------------------- +subroutine mesh_spectral_map_elements + + implicit none + integer(pInt) :: i + + allocate (mesh_mapFEtoCPelem(2_pInt,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt + + forall (i = 1_pInt:mesh_NcpElems) & + mesh_mapFEtoCPelem(1:2,i) = i + +end subroutine mesh_spectral_map_elements + + +!-------------------------------------------------------------------------------------------------- +!> @brief Maps node from FE ID to internal (consecutive) representation. +!! Allocates global array 'mesh_mapFEtoCPnode' +!-------------------------------------------------------------------------------------------------- +subroutine mesh_spectral_map_nodes + + implicit none + integer(pInt) :: i + + allocate (mesh_mapFEtoCPnode(2_pInt,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt + + forall (i = 1_pInt:mesh_Nnodes) & + mesh_mapFEtoCPnode(1:2,i) = i + +end subroutine mesh_spectral_map_nodes + + +!-------------------------------------------------------------------------------------------------- +!> @brief Gets maximum count of nodes, IPs, IP neighbors, and subNodes among cpElements. +!! Allocates global arrays 'mesh_maxNnodes', 'mesh_maxNips', mesh_maxNipNeighbors', +!! and mesh_maxNsubNodes +!-------------------------------------------------------------------------------------------------- +subroutine mesh_spectral_count_cpSizes + + implicit none + integer(pInt) :: t + + t = FE_mapElemtype('C3D8R') ! fake 3D hexahedral 8 node 1 IP element + + mesh_maxNnodes = FE_Nnodes(t) + mesh_maxNips = FE_Nips(t) + mesh_maxNipNeighbors = FE_NipNeighbors(t) + mesh_maxNsubNodes = FE_NsubNodes(t) + +end subroutine mesh_spectral_count_cpSizes + + +!-------------------------------------------------------------------------------------------------- +!> @brief Store x,y,z coordinates of all nodes in mesh. +!! Allocates global arrays 'mesh_node0' and 'mesh_node' +!-------------------------------------------------------------------------------------------------- +subroutine mesh_spectral_build_nodes() + + use IO, only: & + IO_error + + implicit none + integer(pInt) :: n + + allocate ( mesh_node0 (3,mesh_Nnodes) ); mesh_node0 = 0.0_pReal + allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0.0_pReal + + forall (n = 0_pInt:mesh_Nnodes-1_pInt) + mesh_node0(1,n+1_pInt) = & + geomdim(1) * real(mod(n,(res(1)+1_pInt) ),pReal) & + / real(res(1),pReal) + mesh_node0(2,n+1_pInt) = & + geomdim(2) * real(mod(n/(res(1)+1_pInt),(res(2)+1_pInt)),pReal) & + / real(res(2),pReal) + mesh_node0(3,n+1_pInt) = & + geomdim(3) * real(mod(n/(res(1)+1_pInt)/(res(2)+1_pInt),(res(3)+1_pInt)),pReal) & + / real(res(3),pReal) + end forall + + mesh_node = mesh_node0 !why? + +end subroutine mesh_spectral_build_nodes + + +!-------------------------------------------------------------------------------------------------- +!> @brief Store FEid, type, material, texture, and node list per element. +!! Allocates global array 'mesh_element' +!-------------------------------------------------------------------------------------------------- +subroutine mesh_spectral_build_elements(myUnit) + + use IO, only: & + IO_checkAndRewind, & + IO_lc, & + IO_stringValue, & + IO_stringPos, & + IO_error, & + IO_continuousIntValues, & + IO_intValue, & + IO_countContinuousIntValues + + implicit none + integer(pInt), intent(in) :: myUnit + + integer(pInt), parameter :: maxNchunks = 7_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos + integer(pInt) :: e, i, headerLength = 0_pInt, maxIntCount + integer(pInt), dimension(:), allocatable :: microstructures + integer(pInt), dimension(1,1) :: dummySet = 0_pInt + character(len=65536) :: line,keyword + character(len=64), dimension(1) :: dummyName = '' + + call IO_checkAndRewind(myUnit) + + read(myUnit,'(a65536)') line + myPos = IO_stringPos(line,2_pInt) + keyword = IO_lc(IO_StringValue(line,myPos,2_pInt)) + if (keyword(1:4) == 'head') then + headerLength = IO_intValue(line,myPos,1_pInt) + 1_pInt + else + call IO_error(error_ID=841_pInt, ext_msg='mesh_spectral_build_elements') + endif + + rewind(myUnit) + do i = 1_pInt, headerLength + read(myUnit,'(a65536)') line + enddo + + maxIntCount = 0_pInt + i = 1_pInt + + do while (i > 0_pInt) + i = IO_countContinuousIntValues(myUnit) + maxIntCount = max(maxIntCount, i) + enddo + + rewind (myUnit) + do i=1_pInt,headerLength ! skip header + read(myUnit,'(a65536)') line + enddo + + allocate (mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt + allocate (microstructures (1_pInt+maxIntCount)) ; microstructures = 2_pInt + + e = 0_pInt + do while (e < mesh_NcpElems .and. microstructures(1) > 0_pInt) ! fill expected number of elements, stop at end of data (or blank line!) + microstructures = IO_continuousIntValues(myUnit,maxIntCount,dummyName,dummySet,0_pInt) ! get affected elements + do i = 1_pInt,microstructures(1_pInt) + e = e+1_pInt ! valid element entry + mesh_element( 1,e) = e ! FE id + mesh_element( 2,e) = FE_mapElemtype('C3D8R') ! elem type + mesh_element( 3,e) = homog ! homogenization + mesh_element( 4,e) = microstructures(1_pInt+i) ! microstructure + mesh_element( 5,e) = e + (e-1_pInt)/res(1) + & + ((e-1_pInt)/(res(1)*res(2)))*(res(1)+1_pInt) ! base node + mesh_element( 6,e) = mesh_element(5,e) + 1_pInt + mesh_element( 7,e) = mesh_element(5,e) + res(1) + 2_pInt + mesh_element( 8,e) = mesh_element(5,e) + res(1) + 1_pInt + mesh_element( 9,e) = mesh_element(5,e) +(res(1) + 1_pInt) * (res(2) + 1_pInt) ! second floor base node + mesh_element(10,e) = mesh_element(9,e) + 1_pInt + mesh_element(11,e) = mesh_element(9,e) + res(1) + 2_pInt + mesh_element(12,e) = mesh_element(9,e) + res(1) + 1_pInt + mesh_maxValStateVar(1) = max(mesh_maxValStateVar(1),mesh_element(3,e)) !needed for statistics + mesh_maxValStateVar(2) = max(mesh_maxValStateVar(2),mesh_element(4,e)) + enddo + enddo + + deallocate(microstructures) + if (e /= mesh_NcpElems) call IO_error(880_pInt,e) + +end subroutine mesh_spectral_build_elements !-------------------------------------------------------------------------------------------------- @@ -880,7 +1084,6 @@ function mesh_regrid(adaptive,resNewInput,minRes) IO_error use math, only: & math_nearestNeighborSearch, & - deformed_FFT, & math_mul33x3 character(len=1024):: formatString, N_Digits logical, intent(in) :: adaptive ! if true, choose adaptive grid based on resNewInput, otherwise keep it constant @@ -1347,203 +1550,454 @@ function mesh_regrid(adaptive,resNewInput,minRes) end function mesh_regrid -!-------------------------------------------------------------------------------------------------- -!> @brief Count overall number of nodes and elements in mesh and stores them in -!! 'mesh_Nelems' and 'mesh_Nnodes' -!-------------------------------------------------------------------------------------------------- -subroutine mesh_spectral_count_nodesAndElements() +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +subroutine mesh_regular_grid(res,geomdim,defgrad_av,centroids,nodes) +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! Routine to build mesh of (distorted) cubes for given coordinates (= center of the cubes) +! + use debug, only: debug_math, & + debug_level, & + debug_levelBasic implicit none - mesh_Nelems = res(1)*res(2)*res(3) - mesh_Nnodes = (1_pInt + res(1))*(1_pInt + res(2))*(1_pInt + res(3)) + ! input variables + integer(pInt), intent(in), dimension(3) :: res + real(pReal), intent(in), dimension(3) :: geomdim + real(pReal), intent(in), dimension(3,3) :: defgrad_av + real(pReal), intent(in), dimension(res(1), res(2), res(3), 3) :: centroids + ! output variables + real(pReal),intent(out), dimension(res(1)+1_pInt,res(2)+1_pInt,res(3)+1_pInt,3) :: nodes + ! variables with dimension depending on input + real(pReal), dimension(res(1)+2_pInt,res(2)+2_pInt,res(3)+2_pInt,3) :: wrappedCentroids + ! other variables + integer(pInt) :: i,j,k,n + integer(pInt), dimension(3), parameter :: diag = 1_pInt + integer(pInt), dimension(3) :: shift = 0_pInt, lookup = 0_pInt, me = 0_pInt + integer(pInt), dimension(3,8) :: neighbor = reshape((/ & + 0_pInt, 0_pInt, 0_pInt, & + 1_pInt, 0_pInt, 0_pInt, & + 1_pInt, 1_pInt, 0_pInt, & + 0_pInt, 1_pInt, 0_pInt, & + 0_pInt, 0_pInt, 1_pInt, & + 1_pInt, 0_pInt, 1_pInt, & + 1_pInt, 1_pInt, 1_pInt, & + 0_pInt, 1_pInt, 1_pInt & + /), & + (/3,8/)) -end subroutine mesh_spectral_count_nodesAndElements + if (iand(debug_level(debug_math),debug_levelBasic) /= 0_pInt) then + print*, 'Meshing cubes around centroids' + print '(a,3(e12.5))', ' Dimension: ', geomdim + print '(a,3(i5))', ' Resolution:', res + endif + nodes = 0.0_pReal + wrappedCentroids = 0.0_pReal + wrappedCentroids(2_pInt:res(1)+1_pInt,2_pInt:res(2)+1_pInt,2_pInt:res(3)+1_pInt,1:3) = centroids -!-------------------------------------------------------------------------------------------------- -!> @brief Count overall number of CP elements in mesh and stores them in 'mesh_NcpElems' -!-------------------------------------------------------------------------------------------------- -subroutine mesh_spectral_count_cpElements + do k = 0_pInt,res(3)+1_pInt + do j = 0_pInt,res(2)+1_pInt + do i = 0_pInt,res(1)+1_pInt + if (k==0_pInt .or. k==res(3)+1_pInt .or. & ! z skin + j==0_pInt .or. j==res(2)+1_pInt .or. & ! y skin + i==0_pInt .or. i==res(1)+1_pInt ) then ! x skin + me = (/i,j,k/) ! me on skin + shift = sign(abs(res+diag-2_pInt*me)/(res+diag),res+diag-2_pInt*me) + lookup = me-diag+shift*res + wrappedCentroids(i+1_pInt,j+1_pInt,k+1_pInt,1:3) = & + centroids(lookup(1)+1_pInt,lookup(2)+1_pInt,lookup(3)+1_pInt,1:3) - & + matmul(defgrad_av, shift*geomdim) + endif + enddo; enddo; enddo + do k = 0_pInt,res(3) + do j = 0_pInt,res(2) + do i = 0_pInt,res(1) + do n = 1_pInt,8_pInt + nodes(i+1_pInt,j+1_pInt,k+1_pInt,1:3) = & + nodes(i+1_pInt,j+1_pInt,k+1_pInt,1:3) + wrappedCentroids(i+1_pInt+neighbor(1_pInt,n), & + j+1_pInt+neighbor(2,n), & + k+1_pInt+neighbor(3,n),1:3) + enddo; enddo; enddo; enddo + nodes = nodes/8.0_pReal - implicit none +end subroutine mesh_regular_grid - mesh_NcpElems = mesh_Nelems -end subroutine mesh_spectral_count_cpElements +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +subroutine deformed_linear(res,geomdim,defgrad_av,defgrad,coord) +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! Routine to calculate coordinates in current configuration for given defgrad +! using linear interpolation (blurres out high frequency defomation) +! + implicit none + ! input variables + integer(pInt), intent(in), dimension(3) :: res + real(pReal), intent(in), dimension(3) :: geomdim + real(pReal), intent(in), dimension(3,3) :: defgrad_av + real(pReal), intent(in), dimension( res(1),res(2),res(3),3,3) :: defgrad + ! output variables + real(pReal), intent(out), dimension( res(1),res(2),res(3),3) :: coord + ! variables with dimension depending on input + real(pReal), dimension( 8,res(1),res(2),res(3),3) :: coord_avgOrder + ! other variables + real(pReal), dimension(3) :: myStep, fones = 1.0_pReal, parameter_coords, negative, positive, offset_coords + integer(pInt), dimension(3) :: rear, init, ones = 1_pInt, oppo, me + integer(pInt) i, j, k, s, o + integer(pInt), dimension(3,8) :: corner = reshape([ & + 0_pInt, 0_pInt, 0_pInt,& + 1_pInt, 0_pInt, 0_pInt,& + 1_pInt, 1_pInt, 0_pInt,& + 0_pInt, 1_pInt, 0_pInt,& + 1_pInt, 1_pInt, 1_pInt,& + 0_pInt, 1_pInt, 1_pInt,& + 0_pInt, 0_pInt, 1_pInt,& + 1_pInt, 0_pInt, 1_pInt & + ],[3,8]) + integer(pInt), dimension(3,8) :: step = reshape([& + 1_pInt, 1_pInt, 1_pInt,& + -1_pInt, 1_pInt, 1_pInt,& + -1_pInt,-1_pInt, 1_pInt,& + 1_pInt,-1_pInt, 1_pInt,& + -1_pInt,-1_pInt,-1_pInt,& + 1_pInt,-1_pInt,-1_pInt,& + 1_pInt, 1_pInt,-1_pInt,& + -1_pInt, 1_pInt,-1_pInt & + ], [3,8]) + integer(pInt), dimension(3,6) :: order = reshape([ & + 1_pInt, 2_pInt, 3_pInt,& + 1_pInt, 3_pInt, 2_pInt,& + 2_pInt, 1_pInt, 3_pInt,& + 2_pInt, 3_pInt, 1_pInt,& + 3_pInt, 1_pInt, 2_pInt,& + 3_pInt, 2_pInt, 1_pInt & + ], [3,6]) + + coord_avgOrder = 0.0_pReal + + do s = 0_pInt, 7_pInt ! corners (from 0 to 7) + init = corner(:,s+1_pInt)*(res-ones) +ones + oppo = corner(:,mod((s+4_pInt),8_pInt)+1_pInt)*(res-ones) +ones + + do o=1_pInt,6_pInt ! orders (from 1 to 6) + coord = 0_pReal + do k = init(order(3,o)), oppo(order(3,o)), step(order(3,o),s+1_pInt) + rear(order(2,o)) = init(order(2,o)) + do j = init(order(2,o)), oppo(order(2,o)), step(order(2,o),s+1_pInt) + rear(order(1,o)) = init(order(1,o)) + do i = init(order(1,o)), oppo(order(1,o)), step(order(1,o),s+1_pInt) + me(order(1:3,o)) = [i,j,k] + if ( all(me==init)) then + coord(me(1),me(2),me(3),1:3) = geomdim * (matmul(defgrad_av,real(corner(1:3,s+1),pReal)) + & + matmul(defgrad(me(1),me(2),me(3),1:3,1:3),0.5_pReal*real(step(1:3,s+1_pInt)/res,pReal))) + else + myStep = (me-rear)*geomdim/res + coord(me(1),me(2),me(3),1:3) = coord(rear(1),rear(2),rear(3),1:3) + & + 0.5_pReal*matmul(defgrad(me(1),me(2),me(3),1:3,1:3) + & + defgrad(rear(1),rear(2),rear(3),1:3,1:3),myStep) + endif + rear = me + enddo; enddo; enddo + coord_avgOrder(s+1_pInt,1:res(1),1:res(2),1:res(3),1:3) = & + coord_avgOrder(s+1_pInt,1:res(1),1:res(2),1:res(3),1:3) + coord/6.0_pReal + enddo + offset_coords = coord_avgOrder(s+1,1,1,1,1:3) + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + coord_avgOrder(s+1,i,j,k,1:3) = coord_avgOrder(s+1,i,j,k,1:3) - offset_coords + enddo; enddo; enddo + enddo + + do k = 0_pInt, res(3)-1_pInt + do j = 0_pInt, res(2)-1_pInt + do i = 0_pInt, res(1)-1_pInt + parameter_coords = (2.0_pReal*real([i,j,k]+1,pReal)-real(res,pReal))/(real(res,pReal)) + positive = fones + parameter_coords + negative = fones - parameter_coords + coord(i+1_pInt,j+1_pInt,k+1_pInt,1:3)& + =(coord_avgOrder(1,i+1_pInt,j+1_pInt,k+1_pInt,1:3) *negative(1)*negative(2)*negative(3)& + + coord_avgOrder(2,i+1_pInt,j+1_pInt,k+1_pInt,1:3) *positive(1)*negative(2)*negative(3)& + + coord_avgOrder(3,i+1_pInt,j+1_pInt,k+1_pInt,1:3) *positive(1)*positive(2)*negative(3)& + + coord_avgOrder(4,i+1_pInt,j+1_pInt,k+1_pInt,1:3) *negative(1)*positive(2)*negative(3)& + + coord_avgOrder(5,i+1_pInt,j+1_pInt,k+1_pInt,1:3) *positive(1)*positive(2)*positive(3)& + + coord_avgOrder(6,i+1_pInt,j+1_pInt,k+1_pInt,1:3) *negative(1)*positive(2)*positive(3)& + + coord_avgOrder(7,i+1_pInt,j+1_pInt,k+1_pInt,1:3) *negative(1)*negative(2)*positive(3)& + + coord_avgOrder(8,i+1_pInt,j+1_pInt,k+1_pInt,1:3) *positive(1)*negative(2)*positive(3))*0.125_pReal + enddo; enddo; enddo + + offset_coords = matmul(defgrad(1,1,1,1:3,1:3),geomdim/real(res, pReal)/2.0_pReal) - coord(1,1,1,1:3) + do k = 1_pInt, res(3) + do j = 1_pInt, res(2) + do i = 1_pInt, res(1) + coord(i,j,k,1:3) = coord(i,j,k,1:3)+ offset_coords + enddo; enddo; enddo + +end subroutine deformed_linear -!-------------------------------------------------------------------------------------------------- -!> @brief Maps elements from FE ID to internal (consecutive) representation. -!! Allocates global array 'mesh_mapFEtoCPelem' -!-------------------------------------------------------------------------------------------------- -subroutine mesh_spectral_map_elements +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +subroutine deformed_fft(res,geomdim,defgrad_av,scaling,defgrad,coords) +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! Routine to calculate coordinates in current configuration for given defgrad +! using integration in Fourier space (more accurate than deformed(...)) +! + use IO, only: IO_error + use numerics, only: fftw_timelimit, fftw_planner_flag + use debug, only: debug_math, & + debug_level, & + debug_levelBasic + use math, only: PI implicit none - integer(pInt) :: i + ! input variables + integer(pInt), intent(in), dimension(3) :: res + real(pReal), intent(in), dimension(3) :: geomdim + real(pReal), intent(in), dimension(3,3) :: defgrad_av + real(pReal), intent(in) :: scaling + real(pReal), intent(in), dimension(res(1), res(2),res(3),3,3) :: defgrad + ! output variables + real(pReal), intent(out), dimension(res(1), res(2),res(3),3) :: coords +! allocatable arrays for fftw c routines + type(C_PTR) :: fftw_forth, fftw_back + type(C_PTR) :: coords_fftw, defgrad_fftw + real(pReal), dimension(:,:,:,:,:), pointer :: defgrad_real + complex(pReal), dimension(:,:,:,:,:), pointer :: defgrad_fourier + real(pReal), dimension(:,:,:,:), pointer :: coords_real + complex(pReal), dimension(:,:,:,:), pointer :: coords_fourier + ! other variables + integer(pInt) :: i, j, k, m, res1_red + integer(pInt), dimension(3) :: k_s + real(pReal), dimension(3) :: step, offset_coords, integrator - allocate (mesh_mapFEtoCPelem(2_pInt,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt + integrator = geomdim / 2.0_pReal / pi ! see notes where it is used - forall (i = 1_pInt:mesh_NcpElems) & - mesh_mapFEtoCPelem(1:2,i) = i - -end subroutine mesh_spectral_map_elements - - -!-------------------------------------------------------------------------------------------------- -!> @brief Maps node from FE ID to internal (consecutive) representation. -!! Allocates global array 'mesh_mapFEtoCPnode' -!-------------------------------------------------------------------------------------------------- -subroutine mesh_spectral_map_nodes - - implicit none - integer(pInt) :: i - - allocate (mesh_mapFEtoCPnode(2_pInt,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt - - forall (i = 1_pInt:mesh_Nnodes) & - mesh_mapFEtoCPnode(1:2,i) = i - -end subroutine mesh_spectral_map_nodes - - -!-------------------------------------------------------------------------------------------------- -!> @brief Gets maximum count of nodes, IPs, IP neighbors, and subNodes among cpElements. -!! Allocates global arrays 'mesh_maxNnodes', 'mesh_maxNips', mesh_maxNipNeighbors', -!! and mesh_maxNsubNodes -!-------------------------------------------------------------------------------------------------- -subroutine mesh_spectral_count_cpSizes - - implicit none - integer(pInt) :: t - - t = FE_mapElemtype('C3D8R') ! fake 3D hexahedral 8 node 1 IP element - - mesh_maxNnodes = FE_Nnodes(t) - mesh_maxNips = FE_Nips(t) - mesh_maxNipNeighbors = FE_NipNeighbors(t) - mesh_maxNsubNodes = FE_NsubNodes(t) - -end subroutine mesh_spectral_count_cpSizes - - -!-------------------------------------------------------------------------------------------------- -!> @brief Store x,y,z coordinates of all nodes in mesh. -!! Allocates global arrays 'mesh_node0' and 'mesh_node' -!-------------------------------------------------------------------------------------------------- -subroutine mesh_spectral_build_nodes() - - use IO, only: & - IO_error - - implicit none - integer(pInt) :: n - - allocate ( mesh_node0 (3,mesh_Nnodes) ); mesh_node0 = 0.0_pReal - allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0.0_pReal - - forall (n = 0_pInt:mesh_Nnodes-1_pInt) - mesh_node0(1,n+1_pInt) = & - geomdim(1) * real(mod(n,(res(1)+1_pInt) ),pReal) & - / real(res(1),pReal) - mesh_node0(2,n+1_pInt) = & - geomdim(2) * real(mod(n/(res(1)+1_pInt),(res(2)+1_pInt)),pReal) & - / real(res(2),pReal) - mesh_node0(3,n+1_pInt) = & - geomdim(3) * real(mod(n/(res(1)+1_pInt)/(res(2)+1_pInt),(res(3)+1_pInt)),pReal) & - / real(res(3),pReal) - end forall - - mesh_node = mesh_node0 !why? - -end subroutine mesh_spectral_build_nodes - - -!-------------------------------------------------------------------------------------------------- -!> @brief Store FEid, type, material, texture, and node list per element. -!! Allocates global array 'mesh_element' -!-------------------------------------------------------------------------------------------------- -subroutine mesh_spectral_build_elements(myUnit) - - use IO, only: & - IO_checkAndRewind, & - IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_error, & - IO_continuousIntValues, & - IO_intValue, & - IO_countContinuousIntValues - - implicit none - integer(pInt), intent(in) :: myUnit - - integer(pInt), parameter :: maxNchunks = 7_pInt - integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos - integer(pInt) :: e, i, headerLength = 0_pInt, maxIntCount - integer(pInt), dimension(:), allocatable :: microstructures - integer(pInt), dimension(1,1) :: dummySet = 0_pInt - character(len=65536) :: line,keyword - character(len=64), dimension(1) :: dummyName = '' - - call IO_checkAndRewind(myUnit) - - read(myUnit,'(a65536)') line - myPos = IO_stringPos(line,2_pInt) - keyword = IO_lc(IO_StringValue(line,myPos,2_pInt)) - if (keyword(1:4) == 'head') then - headerLength = IO_intValue(line,myPos,1_pInt) + 1_pInt - else - call IO_error(error_ID=841_pInt, ext_msg='mesh_spectral_build_elements') + if (iand(debug_level(debug_math),debug_levelBasic) /= 0_pInt) then + print*, 'Restore geometry using FFT-based integration' + print '(a,3(e12.5))', ' Dimension: ', geomdim + print '(a,3(i5))', ' Resolution:', res endif - rewind(myUnit) - do i = 1_pInt, headerLength - read(myUnit,'(a65536)') line - enddo + res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c) + step = geomdim/real(res, pReal) - maxIntCount = 0_pInt - i = 1_pInt + if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt) + call fftw_set_timelimit(fftw_timelimit) + defgrad_fftw = fftw_alloc_complex(int(res1_red *res(2)*res(3)*9_pInt,C_SIZE_T)) !C_SIZE_T is of type integer(8) + call c_f_pointer(defgrad_fftw, defgrad_real, [res(1)+2_pInt,res(2),res(3),3_pInt,3_pInt]) + call c_f_pointer(defgrad_fftw, defgrad_fourier,[res1_red ,res(2),res(3),3_pInt,3_pInt]) + coords_fftw = fftw_alloc_complex(int(res1_red *res(2)*res(3)*3_pInt,C_SIZE_T)) !C_SIZE_T is of type integer(8) + call c_f_pointer(coords_fftw, coords_real, [res(1)+2_pInt,res(2),res(3),3_pInt]) + call c_f_pointer(coords_fftw, coords_fourier, [res1_red ,res(2),res(3),3_pInt]) + fftw_forth = fftw_plan_many_dft_r2c(3_pInt,(/res(3),res(2) ,res(1)/),9_pInt,& ! dimensions , length in each dimension in reversed order + defgrad_real,(/res(3),res(2) ,res(1)+2_pInt/),& ! input data , physical length in each dimension in reversed order + 1_pInt, res(3)*res(2)*(res(1)+2_pInt),& ! striding , product of physical lenght in the 3 dimensions + defgrad_fourier,(/res(3),res(2) ,res1_red/),& + 1_pInt, res(3)*res(2)* res1_red,fftw_planner_flag) - do while (i > 0_pInt) - i = IO_countContinuousIntValues(myUnit) - maxIntCount = max(maxIntCount, i) - enddo - - rewind (myUnit) - do i=1_pInt,headerLength ! skip header - read(myUnit,'(a65536)') line - enddo - - allocate (mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt - allocate (microstructures (1_pInt+maxIntCount)) ; microstructures = 2_pInt + fftw_back = fftw_plan_many_dft_c2r(3_pInt,(/res(3),res(2) ,res(1)/),3_pInt,& + coords_fourier,(/res(3),res(2) ,res1_red/),& + 1_pInt, res(3)*res(2)* res1_red,& + coords_real,(/res(3),res(2) ,res(1)+2_pInt/),& + 1_pInt, res(3)*res(2)*(res(1)+2_pInt),fftw_planner_flag) - e = 0_pInt - do while (e < mesh_NcpElems .and. microstructures(1) > 0_pInt) ! fill expected number of elements, stop at end of data (or blank line!) - microstructures = IO_continuousIntValues(myUnit,maxIntCount,dummyName,dummySet,0_pInt) ! get affected elements - do i = 1_pInt,microstructures(1_pInt) - e = e+1_pInt ! valid element entry - mesh_element( 1,e) = e ! FE id - mesh_element( 2,e) = FE_mapElemtype('C3D8R') ! elem type - mesh_element( 3,e) = homog ! homogenization - mesh_element( 4,e) = microstructures(1_pInt+i) ! microstructure - mesh_element( 5,e) = e + (e-1_pInt)/res(1) + & - ((e-1_pInt)/(res(1)*res(2)))*(res(1)+1_pInt) ! base node - mesh_element( 6,e) = mesh_element(5,e) + 1_pInt - mesh_element( 7,e) = mesh_element(5,e) + res(1) + 2_pInt - mesh_element( 8,e) = mesh_element(5,e) + res(1) + 1_pInt - mesh_element( 9,e) = mesh_element(5,e) +(res(1) + 1_pInt) * (res(2) + 1_pInt) ! second floor base node - mesh_element(10,e) = mesh_element(9,e) + 1_pInt - mesh_element(11,e) = mesh_element(9,e) + res(1) + 2_pInt - mesh_element(12,e) = mesh_element(9,e) + res(1) + 1_pInt - mesh_maxValStateVar(1) = max(mesh_maxValStateVar(1),mesh_element(3,e)) !needed for statistics - mesh_maxValStateVar(2) = max(mesh_maxValStateVar(2),mesh_element(4,e)) + + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + defgrad_real(i,j,k,1:3,1:3) = defgrad(i,j,k,1:3,1:3) ! ensure that data is aligned properly (fftw_alloc) + enddo; enddo; enddo + + call fftw_execute_dft_r2c(fftw_forth, defgrad_real, defgrad_fourier) + + !remove highest frequency in each direction + if(res(1)>1_pInt) & + defgrad_fourier( res(1)/2_pInt+1_pInt,1:res(2) ,1:res(3) ,& + 1:3,1:3) = cmplx(0.0_pReal,0.0_pReal,pReal) + if(res(2)>1_pInt) & + defgrad_fourier(1:res1_red ,res(2)/2_pInt+1_pInt,1:res(3) ,& + 1:3,1:3) = cmplx(0.0_pReal,0.0_pReal,pReal) + if(res(3)>1_pInt) & + defgrad_fourier(1:res1_red ,1:res(2) ,res(3)/2_pInt+1_pInt,& + 1:3,1:3) = cmplx(0.0_pReal,0.0_pReal,pReal) + + coords_fourier = cmplx(0.0_pReal,0.0_pReal,pReal) + do k = 1_pInt, res(3) + k_s(3) = k-1_pInt + if(k > res(3)/2_pInt+1_pInt) k_s(3) = k_s(3)-res(3) + do j = 1_pInt, res(2) + k_s(2) = j-1_pInt + if(j > res(2)/2_pInt+1_pInt) k_s(2) = k_s(2)-res(2) + do i = 1_pInt, res1_red + k_s(1) = i-1_pInt + do m = 1_pInt,3_pInt + coords_fourier(i,j,k,m) = sum(defgrad_fourier(i,j,k,m,1:3)*cmplx(0.0_pReal,real(k_s,pReal)*integrator,pReal)) + enddo + if (k_s(3) /= 0_pInt .or. k_s(2) /= 0_pInt .or. k_s(1) /= 0_pInt) & + coords_fourier(i,j,k,1:3) = coords_fourier(i,j,k,1:3) / real(-sum(k_s*k_s),pReal) +! if(i/=1_pInt) coords_fourier(i,j,k,1:3) = coords_fourier(i,j,k,1:3)& ! substituting division by (on the fly calculated) xi * 2pi * img by multiplication with reversed img/real part +! - defgrad_fourier(i,j,k,1:3,1)*cmplx(0.0_pReal,integrator(1)/real(k_s(1),pReal),pReal) +! if(j/=1_pInt) coords_fourier(i,j,k,1:3) = coords_fourier(i,j,k,1:3)& +! - defgrad_fourier(i,j,k,1:3,2)*cmplx(0.0_pReal,integrator(2)/real(k_s(2),pReal),pReal) +! if(k/=1_pInt) coords_fourier(i,j,k,1:3) = coords_fourier(i,j,k,1:3)& +! - defgrad_fourier(i,j,k,1:3,3)*cmplx(0.0_pReal,integrator(3)/real(k_s(3),pReal),pReal) + enddo; enddo; enddo + + call fftw_execute_dft_c2r(fftw_back,coords_fourier,coords_real) + coords_real = coords_real/real(res(1)*res(2)*res(3),pReal) + + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + coords(i,j,k,1:3) = coords_real(i,j,k,1:3) ! ensure that data is aligned properly (fftw_alloc) + enddo; enddo; enddo + + offset_coords = matmul(defgrad(1,1,1,1:3,1:3),step/2.0_pReal) - scaling*coords(1,1,1,1:3) + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + coords(i,j,k,1:3) = scaling*coords(i,j,k,1:3) + offset_coords + matmul(defgrad_av,& + (/step(1)*real(i-1_pInt,pReal),& + step(2)*real(j-1_pInt,pReal),& + step(3)*real(k-1_pInt,pReal)/)) + + enddo; enddo; enddo + + call fftw_destroy_plan(fftw_forth) + call fftw_destroy_plan(fftw_back) + call fftw_free(defgrad_fftw) + call fftw_free(coords_fftw) + +end subroutine deformed_fft + + +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +subroutine volume_compare(res,geomdim,defgrad,nodes,volume_mismatch) +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! Routine to calculate the mismatch between volume of reconstructed (compatible +! cube and determinant of defgrad at the FP + + use debug, only: debug_math, & + debug_level, & + debug_levelBasic + use math, only: PI, & + math_det33, & + math_volTetrahedron + + implicit none + ! input variables + integer(pInt), intent(in), dimension(3) :: res + real(pReal), intent(in), dimension(3) :: geomdim + real(pReal), intent(in), dimension(res(1), res(2), res(3), 3,3) :: defgrad + real(pReal), intent(in), dimension(res(1)+1_pInt,res(2)+1_pInt,res(3)+1_pInt,3) :: nodes + ! output variables + real(pReal), intent(out), dimension(res(1), res(2), res(3)) :: volume_mismatch + ! other variables + real(pReal), dimension(8,3) :: coords + integer(pInt) i,j,k + real(pReal) vol_initial + + if (iand(debug_level(debug_math),debug_levelBasic) /= 0_pInt) then + print*, 'Calculating volume mismatch' + print '(a,3(e12.5))', ' Dimension: ', geomdim + print '(a,3(i5))', ' Resolution:', res + endif + + vol_initial = geomdim(1)*geomdim(2)*geomdim(3)/(real(res(1)*res(2)*res(3), pReal)) + do k = 1_pInt,res(3) + do j = 1_pInt,res(2) + do i = 1_pInt,res(1) + coords(1,1:3) = nodes(i, j, k ,1:3) + coords(2,1:3) = nodes(i+1_pInt,j, k ,1:3) + coords(3,1:3) = nodes(i+1_pInt,j+1_pInt,k ,1:3) + coords(4,1:3) = nodes(i, j+1_pInt,k ,1:3) + coords(5,1:3) = nodes(i, j, k+1_pInt,1:3) + coords(6,1:3) = nodes(i+1_pInt,j, k+1_pInt,1:3) + coords(7,1:3) = nodes(i+1_pInt,j+1_pInt,k+1_pInt,1:3) + coords(8,1:3) = nodes(i, j+1_pInt,k+1_pInt,1:3) + volume_mismatch(i,j,k) = abs(math_volTetrahedron(coords(7,1:3),coords(1,1:3),coords(8,1:3),coords(4,1:3))) & + + abs(math_volTetrahedron(coords(7,1:3),coords(1,1:3),coords(8,1:3),coords(5,1:3))) & + + abs(math_volTetrahedron(coords(7,1:3),coords(1,1:3),coords(3,1:3),coords(4,1:3))) & + + abs(math_volTetrahedron(coords(7,1:3),coords(1,1:3),coords(3,1:3),coords(2,1:3))) & + + abs(math_volTetrahedron(coords(7,1:3),coords(5,1:3),coords(2,1:3),coords(6,1:3))) & + + abs(math_volTetrahedron(coords(7,1:3),coords(5,1:3),coords(2,1:3),coords(1,1:3))) + volume_mismatch(i,j,k) = volume_mismatch(i,j,k)/math_det33(defgrad(i,j,k,1:3,1:3)) + enddo; enddo; enddo + volume_mismatch = volume_mismatch/vol_initial + +end subroutine volume_compare + + +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +subroutine shape_compare(res,geomdim,defgrad,nodes,centroids,shape_mismatch) +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! Routine to calculate the mismatch between the vectors from the central point to +! the corners of reconstructed (combatible) volume element and the vectors calculated by deforming +! the initial volume element with the current deformation gradient + + use debug, only: debug_math, & + debug_level, & + debug_levelBasic + + implicit none + ! input variables + integer(pInt), intent(in), dimension(3) :: res + real(pReal), intent(in), dimension(3) :: geomdim + real(pReal), intent(in), dimension(res(1), res(2), res(3), 3,3) :: defgrad + real(pReal), intent(in), dimension(res(1)+1_pInt,res(2)+1_pInt,res(3)+1_pInt,3) :: nodes + real(pReal), intent(in), dimension(res(1), res(2), res(3), 3) :: centroids + ! output variables + real(pReal), intent(out), dimension(res(1), res(2), res(3)) :: shape_mismatch + ! other variables + real(pReal), dimension(8,3) :: coords_initial + integer(pInt) i,j,k + + if (iand(debug_level(debug_math),debug_levelBasic) /= 0_pInt) then + print*, 'Calculating shape mismatch' + print '(a,3(e12.5))', ' Dimension: ', geomdim + print '(a,3(i5))', ' Resolution:', res + endif + + coords_initial(1,1:3) = (/-geomdim(1)/2.0_pReal/real(res(1),pReal),& + -geomdim(2)/2.0_pReal/real(res(2),pReal),& + -geomdim(3)/2.0_pReal/real(res(3),pReal)/) + coords_initial(2,1:3) = (/+geomdim(1)/2.0_pReal/real(res(1),pReal),& + -geomdim(2)/2.0_pReal/real(res(2),pReal),& + -geomdim(3)/2.0_pReal/real(res(3),pReal)/) + coords_initial(3,1:3) = (/+geomdim(1)/2.0_pReal/real(res(1),pReal),& + +geomdim(2)/2.0_pReal/real(res(2),pReal),& + -geomdim(3)/2.0_pReal/real(res(3),pReal)/) + coords_initial(4,1:3) = (/-geomdim(1)/2.0_pReal/real(res(1),pReal),& + +geomdim(2)/2.0_pReal/real(res(2),pReal),& + -geomdim(3)/2.0_pReal/real(res(3),pReal)/) + coords_initial(5,1:3) = (/-geomdim(1)/2.0_pReal/real(res(1),pReal),& + -geomdim(2)/2.0_pReal/real(res(2),pReal),& + +geomdim(3)/2.0_pReal/real(res(3),pReal)/) + coords_initial(6,1:3) = (/+geomdim(1)/2.0_pReal/real(res(1),pReal),& + -geomdim(2)/2.0_pReal/real(res(2),pReal),& + +geomdim(3)/2.0_pReal/real(res(3),pReal)/) + coords_initial(7,1:3) = (/+geomdim(1)/2.0_pReal/real(res(1),pReal),& + +geomdim(2)/2.0_pReal/real(res(2),pReal),& + +geomdim(3)/2.0_pReal/real(res(3),pReal)/) + coords_initial(8,1:3) = (/-geomdim(1)/2.0_pReal/real(res(1),pReal),& + +geomdim(2)/2.0_pReal/real(res(2),pReal),& + +geomdim(3)/2.0_pReal/real(res(3),pReal)/) + do i=1_pInt,8_pInt enddo - enddo + do k = 1_pInt,res(3) + do j = 1_pInt,res(2) + do i = 1_pInt,res(1) + shape_mismatch(i,j,k) = & + sqrt(sum((nodes(i, j, k, 1:3) - centroids(i,j,k,1:3)& + - matmul(defgrad(i,j,k,1:3,1:3), coords_initial(1,1:3)))**2.0_pReal))& + + sqrt(sum((nodes(i+1_pInt,j, k, 1:3) - centroids(i,j,k,1:3)& + - matmul(defgrad(i,j,k,1:3,1:3), coords_initial(2,1:3)))**2.0_pReal))& + + sqrt(sum((nodes(i+1_pInt,j+1_pInt,k, 1:3) - centroids(i,j,k,1:3)& + - matmul(defgrad(i,j,k,1:3,1:3), coords_initial(3,1:3)))**2.0_pReal))& + + sqrt(sum((nodes(i, j+1_pInt,k, 1:3) - centroids(i,j,k,1:3)& + - matmul(defgrad(i,j,k,1:3,1:3), coords_initial(4,1:3)))**2.0_pReal))& + + sqrt(sum((nodes(i, j, k+1_pInt,1:3) - centroids(i,j,k,1:3)& + - matmul(defgrad(i,j,k,1:3,1:3), coords_initial(5,1:3)))**2.0_pReal))& + + sqrt(sum((nodes(i+1_pInt,j, k+1_pInt,1:3) - centroids(i,j,k,1:3)& + - matmul(defgrad(i,j,k,1:3,1:3), coords_initial(6,1:3)))**2.0_pReal))& + + sqrt(sum((nodes(i+1_pInt,j+1_pInt,k+1_pInt,1:3) - centroids(i,j,k,1:3)& + - matmul(defgrad(i,j,k,1:3,1:3), coords_initial(7,1:3)))**2.0_pReal))& + + sqrt(sum((nodes(i, j+1_pInt,k+1_pInt,1:3) - centroids(i,j,k,1:3)& + - matmul(defgrad(i,j,k,1:3,1:3), coords_initial(8,1:3)))**2.0_pReal)) + enddo; enddo; enddo - deallocate(microstructures) - if (e /= mesh_NcpElems) call IO_error(880_pInt,e) - -end subroutine mesh_spectral_build_elements +end subroutine shape_compare #endif diff --git a/lib/damask/__init__.py b/lib/damask/__init__.py index 563e9225e..35858f910 100644 --- a/lib/damask/__init__.py +++ b/lib/damask/__init__.py @@ -26,22 +26,22 @@ try: core.numerics.init = core.numerics.numerics_init core.debug.init = core.debug.debug_init core.math.init = core.math.math_init - #core.math.volumeMismatch = core.math.math_volumeMismatch - #core.math.shapeMismatch = core.math.math_shapeMismatch - #core.math.deformedCoordsLin = core.math.math_deformedCoordsLin - #core.math.deformedCoordsFFT = core.math.math_deformedCoordsFFT - #core.math.curlFFT = core.math.math_curlFFT - #core.math.divergenceFFT = core.math.math_divergenceFFT - #core.math.divergenceFDM = core.math.math_divergenceFDM - #core.math.tensorAvg = core.math.math_tensorAvg - #core.math.logStrainSpat = core.math.math_logStrainSpat - #core.math.logStrainMat = core.math.math_logStrainMat - #core.math.cauchyStress = core.math.math_cauchyStress + core.math.curlFFT = core.math.math_curlFFT + core.math.divergenceFFT = core.math.math_divergenceFFT + core.math.divergenceFDM = core.math.math_divergenceFDM #core.math.periodicNearestNeighbor = core.math.math_periodicNearestNeighbor + core.math.tensorAvg = core.math.math_tensorAvg + core.math.logstrainSpat = core.math.math_logstrainSpat + core.math.logstrainMat = core.math.math_logstrainMat + core.math.cauchy = core.math.math_cauchy core.FEsolving.init = core.FEsolving.FE_init core.mesh.init = core.mesh.mesh_init core.mesh.regrid = core.mesh.mesh_regrid + #core.mesh.volumeMismatch = core.mesh.mesh_volumeMismatch + #core.mesh.shapeMismatch = core.mesh.mesh_shapeMismatch #core.mesh.nodesAroundCentroids = core.mesh.mesh_spectral_nodesAroundCentroids + #core.mesh.deformedCoordsLin = core.mesh.mesh_deformedCoordsLin + #core.mesh.deformedCoordsFFT = core.mesh.mesh_deformedCoordsFFT except ImportError: sys.stderr.write('\nWARNING: Core module (Fortran code) not available, try to run setup_processing.py\nError Message when importing core.so: \n\n') core = None # from http://www.python.org/dev/peps/pep-0008/ diff --git a/processing/post/3Dvisualize.py b/processing/post/3Dvisualize.py index c14481620..d728dbbb1 100755 --- a/processing/post/3Dvisualize.py +++ b/processing/post/3Dvisualize.py @@ -430,8 +430,7 @@ for filename in args: values = numpy.array(sorted([map(transliterateToFloat,line.split()[:maxcol]) for line in content[headrow+1:]], - key=lambda x:(x[locol+0],x[locol+1],x[locol+2])), # sort with z as fastest and x as slowest index - 'd') + key=lambda x:(x[locol+0],x[locol+1],x[locol+2])),'d') # sort with z as fastest and x as slowest index N = len(values) grid = [{},{},{}] @@ -457,20 +456,20 @@ for filename in args: if options.undeformed: defgrad_av = numpy.eye(3) else: - defgrad_av = damask.core.math.tensor_avg(res,numpy.reshape(values[:,column['tensor'][options.defgrad]: - column['tensor'][options.defgrad]+9], - (res[0],res[1],res[2],3,3))) + defgrad_av = damask.core.math.tensorAvg(numpy.reshape(values[:,column['tensor'][options.defgrad]: + column['tensor'][options.defgrad]+9], + (res[0],res[1],res[2],3,3))) if options.linearreconstruction: - centroids = damask.core.math.deformed_linear(res,dim,defgrad_av, + centroids = damask.core.mesh.deformed_linear(res,dim,defgrad_av, numpy.reshape(values[:,column['tensor'][options.defgrad]: column['tensor'][options.defgrad]+9], (res[0],res[1],res[2],3,3))) else: - centroids = damask.core.math.deformed_fft(res,dim,defgrad_av,options.scaling, + centroids = damask.core.mesh.deformed_fft(res,dim,defgrad_av,options.scaling, numpy.reshape(values[:,column['tensor'][options.defgrad]: column['tensor'][options.defgrad]+9], (res[0],res[1],res[2],3,3))) - ms = damask.core.math.mesh_regular_grid(res,dim,defgrad_av,centroids) + ms = damask.core.mesh.mesh_regular_grid(res,dim,defgrad_av,centroids) fields = {\ 'tensor': {},\ 'vector': {},\ diff --git a/processing/post/addCompatibilityMismatch.py b/processing/post/addCompatibilityMismatch.py index 4e8e3d6fd..4ed54ffed 100755 --- a/processing/post/addCompatibilityMismatch.py +++ b/processing/post/addCompatibilityMismatch.py @@ -144,11 +144,11 @@ for file in files: idx += 1 defgrad[x,y,z] = numpy.array(map(float,table.data[column:column+9]),'d').reshape(3,3) - defgrad_av = damask.core.math.tensor_avg(res,defgrad) - centroids = damask.core.math.deformed_fft(res,geomdim,defgrad_av,1.0,defgrad) - nodes = damask.core.math.mesh_regular_grid(res,geomdim,defgrad_av,centroids) - if not options.noShape: shape_mismatch = damask.core.math.shape_compare( res,geomdim,defgrad,nodes,centroids) - if not options.noVolume: volume_mismatch = damask.core.math.volume_compare(res,geomdim,defgrad,nodes) + defgrad_av = damask.core.math.tensorAvg(defgrad) + centroids = damask.core.mesh.deformed_fft(res,geomdim,defgrad_av,1.0,defgrad) + nodes = damask.core.mesh.mesh_regular_grid(res,geomdim,defgrad_av,centroids) + if not options.noShape: shape_mismatch = damask.core.mesh.shape_compare( res,geomdim,defgrad,nodes,centroids) + if not options.noVolume: volume_mismatch = damask.core.mesh.volume_compare(res,geomdim,defgrad,nodes) # ------------------------------------------ process data --------------------------------------- @@ -169,4 +169,4 @@ for file in files: file['input'].close() # close input ASCII table if file['name'] != 'STDIN': file['output'].close # close output ASCII table - os.rename(file['name']+'_tmp',file['name']) # overwrite old one with tmp new \ No newline at end of file + os.rename(file['name']+'_tmp',file['name']) # overwrite old one with tmp new diff --git a/processing/post/addCurl.py b/processing/post/addCurl.py index eff147cc4..d379f534a 100755 --- a/processing/post/addCurl.py +++ b/processing/post/addCurl.py @@ -170,7 +170,7 @@ for file in files: for datatype,labels in active.items(): # loop over vector,tensor for label in labels: # loop over all requested curls - curl[datatype][label] = damask.core.math.curl_fft(resolution,dimension,datainfo[datatype]['len']//3,values[datatype][label]) + curl[datatype][label] = damask.core.math.curlFFT(dimension,values[datatype][label]) # ------------------------------------------ process data --------------------------------------- diff --git a/processing/post/addDeformedConfiguration.py b/processing/post/addDeformedConfiguration.py index c5fc75edd..c2b0a0ae2 100755 --- a/processing/post/addDeformedConfiguration.py +++ b/processing/post/addDeformedConfiguration.py @@ -133,8 +133,8 @@ for file in files: defgrad[x,y,z] = numpy.array(map(float,table.data[column:column+9]),'d').reshape(3,3) # ------------------------------------------ process value field ---------------------------- - defgrad_av = damask.core.math.tensor_avg(res,defgrad) - centroids = damask.core.math.deformed_fft(res,geomdim,defgrad_av,1.0,defgrad) + defgrad_av = damask.core.math.tensorAvg(defgrad) + centroids = damask.core.mesh.deformed_fft(res,geomdim,defgrad_av,1.0,defgrad) # ------------------------------------------ process data --------------------------------------- diff --git a/processing/post/addDivergence.py b/processing/post/addDivergence.py index 2d985a18a..9b9cf11fd 100755 --- a/processing/post/addDivergence.py +++ b/processing/post/addDivergence.py @@ -187,9 +187,9 @@ for file in files: for label in labels: # loop over all requested divergencies for accuracy in options.accuracy: if accuracy == 'FFT': - divergence[datatype][label][accuracy] = damask.core.math.divergence_fft(resolution,dimension,datainfo[datatype]['len']//3,values[datatype][label]) + divergence[datatype][label][accuracy] = damask.core.math.divergenceFFT(dimension,values[datatype][label]) else: - divergence[datatype][label][accuracy] = damask.core.math.divergence_fdm(resolution,dimension,datainfo[datatype]['len']//3,eval(accuracy)//2-1,values[datatype][label]) + divergence[datatype][label][accuracy] = damask.core.math.divergenceFDM(dimension,eval(accuracy)//2-1,values[datatype][label]) # ------------------------------------------ process data --------------------------------------- table.data_rewind() diff --git a/processing/post/spectral_buildElements.py b/processing/post/spectral_buildElements.py index 0764c28c7..b4959a637 100755 --- a/processing/post/spectral_buildElements.py +++ b/processing/post/spectral_buildElements.py @@ -136,9 +136,9 @@ for file in files: # ------------------------------------------ process value field ---------------------------- - defgrad_av = damask.core.math.tensor_avg(res,defgrad) - centroids = damask.core.math.deformed_fft(res,geomdim,defgrad_av,1.0,defgrad) - nodes = damask.core.math.mesh_regular_grid(res,geomdim,defgrad_av,centroids) + defgrad_av = damask.core.math.tensorAvg(defgrad) + centroids = damask.core.mesh.deformed_fft(res,geomdim,defgrad_av,1.0,defgrad) + nodes = damask.core.mesh.mesh_regular_grid(res,geomdim,defgrad_av,centroids) # ------------------------------------------ process data ---------------------------------------