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

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
This commit is contained in:
Martin Diehl 2012-08-27 08:04:47 +00:00
parent 55a5112f36
commit 96ba5ecae4
16 changed files with 1072 additions and 1085 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
code/DAMASK_spectral.f90
View File

@ -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

1
code/DAMASK_spectral_Driver.f90
View File

@ -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, &

14
code/DAMASK_spectral_SolverAL.f90
View File

@ -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

12
code/DAMASK_spectral_SolverBasic.f90
View File

@ -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

13
code/DAMASK_spectral_SolverBasicPETSC.f90
View File

@ -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

5
code/DAMASK_spectral_Utilities.f90
View File

@ -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

234
code/damask.core.pyf
View File

@ -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

997
code/math.f90
View File

File diff suppressed because it is too large Load Diff

810
code/mesh.f90
View File

@ -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

22
lib/damask/__init__.py
View File

@ -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/

15
processing/post/3Dvisualize.py
View File

@ -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': {},\

12
processing/post/addCompatibilityMismatch.py
View File

@ -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
os.rename(file['name']+'_tmp',file['name']) # overwrite old one with tmp new

2
processing/post/addCurl.py
View File

@ -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 ---------------------------------------

4
processing/post/addDeformedConfiguration.py
View File

@ -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 ---------------------------------------

4
processing/post/addDivergence.py
View File

@ -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()

6
processing/post/spectral_buildElements.py
View File

@ -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 ---------------------------------------