spectral_utilities.f90

Go to the documentation of this file.

# 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/grid/spectral_utilities.f90"
# 1 "<built-in>"
# 1 "<command-line>"
# 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/grid/spectral_utilities.f90"
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
module spectral_utilities
  use, intrinsic :: iso_c_binding
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h" 1
!
!
!  Part of the base include file for Fortran use of 1.
!  Note: This file should contain only define statements and
!  not the declaration of variables.
 
! No spaces for #defines as some compilers (PGI) also adds
! those additional spaces during preprocessing - bad for fixed format
!
 
 
 
# 1 "/opt/petsc-3.10.3/Intel-18.4-IntelMPI-2018/include/petscconf.h" 1
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
# 13 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h" 2
 
 
 
 
# 1 "/opt/petsc-3.10.3/include/petscversion.h" 1
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
# 17 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscviewer.h" 1
!
!  Include file for Fortran use of the PetscViewer package in 1
!
 
 
 
 
 
 
 
 
 
 
 
 
 
# 31 "/opt/petsc-3.10.3/include/petsc/finclude/petscviewer.h"
 
# 18 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscerror.h" 1
 
!
!  Include file for Fortran error codes
!    These are also in include/petscerror.h
!
 
 
 
# 24 "/opt/petsc-3.10.3/include/petsc/finclude/petscerror.h"
 
# 38 "/opt/petsc-3.10.3/include/petsc/finclude/petscerror.h"
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
# 19 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petsclog.h" 1
!
!  No includes needed for logging
# 20 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscbag.h" 1
!
!
!  Include file for Fortran use of the Bag package in 1
!
 
 
 
 
 
!
!  End of Fortran include file for the IS package in 1
 
# 21 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h" 2
 
!
! The real*8,complex*16 notatiton is used so that the
! 1 double/complex variables are not affected by
! compiler options like -r4,-r8, sometimes invoked
! by the user. NAG compiler does not like integer*4,real*8
 
# 41 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h"
 
 
 
 
 
 
 
 
 
# 63 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h"
 
 
 
 
 
 
!
 
 
 
 
 
!
# 85 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h"
!
 
 
 
 
 
 
!
 
 
 
 
!
 
!
 
 
!
# 128 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h"
 
# 150 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h"
!
!     Macro for templating between real and complex
!
# 174 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h"
 
 
 
 
 
 
 
 
 
!
!    Allows the matrix Fortran Kernels to work with single precision
!    matrix data structures
!
 
!
!     PetscLogDouble variables are used to contain double precision numbers
!     that are not used in the numerical computations, but rather in logging,
!     timing etc.
!
 
 
!
!     Macros for error checking
!
 
 
 
 
 
 
 
 
# 215 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h"
 
 
 
 
 
 
 
 
 
 
 
 
# 9 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/grid/spectral_utilities.f90" 2
  use petscsys
 
  use prec
  use damask_interface
  use math
  use rotations
  use io
  use discretization_grid
  use numerics
  use debug
  use config
  use discretization
  use homogenization
 
  implicit none
  private
 
  include 'fftw3-mpi.f03'
 
!--------------------------------------------------------------------------------------------------
! field labels information
  enum, bind(c); enumerator :: &
    field_undefined_id, &
    field_mech_id, &
    field_thermal_id, &
    field_damage_id
  end enum
 
!--------------------------------------------------------------------------------------------------
! grid related information information
  real(preal), protected,  public                :: wgt
  integer,     protected,  public                :: grid1red
  real(preal), protected,  public,  dimension(3) :: scaledgeomsize
 
!--------------------------------------------------------------------------------------------------
! variables storing information for spectral method and FFTW
 
  real   (c_double),        public,  dimension(:,:,:,:,:),     pointer     :: tensorfield_real
  complex(C_DOUBLE_COMPLEX),public,  dimension(:,:,:,:,:),     pointer     :: tensorfield_fourier
  real(c_double),           public,  dimension(:,:,:,:),       pointer     :: vectorfield_real
  complex(C_DOUBLE_COMPLEX),public,  dimension(:,:,:,:),       pointer     :: vectorfield_fourier
  real(c_double),           public,  dimension(:,:,:),         pointer     :: scalarfield_real
  complex(C_DOUBLE_COMPLEX),public,  dimension(:,:,:),         pointer     :: scalarfield_fourier
  complex(pReal),           private, dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat
  complex(pReal),           private, dimension(:,:,:,:),       allocatable :: xi1st
  complex(pReal),           private, dimension(:,:,:,:),       allocatable :: xi2nd
  real(preal),              private, dimension(3,3,3,3)                    :: c_ref
 
 
!--------------------------------------------------------------------------------------------------
! plans for FFTW
  type(c_ptr),   private :: &
    plantensorforth, &                                                                              !< FFTW MPI plan P(x) to P(k)
    plantensorback, &                                                                               !< FFTW MPI plan F(k) to F(x)
    planvectorforth, &                                                                              !< FFTW MPI plan v(x) to v(k)
    planvectorback, &                                                                               !< FFTW MPI plan v(k) to v(x)
    planscalarforth, &                                                                              !< FFTW MPI plan s(x) to s(k)
    planscalarback
 
!--------------------------------------------------------------------------------------------------
! variables controlling debugging
  logical, private :: &
    debuggeneral, &                                                                                 !< general debugging of spectral solver
    debugrotation, &                                                                                !< also printing out results in lab frame
    debugpetsc
 
!--------------------------------------------------------------------------------------------------
! derived types
  type, public :: tsolutionstate                                                                    
    integer :: &
       iterationsneeded  = 0
    logical :: &
       converged         = .true., &
       stagconverged     = .true., &
       termill           = .false.
  end type tsolutionstate
 
  type, public :: tboundarycondition                                                                
    real(preal), dimension(3,3) :: values      = 0.0_preal, &
                                   maskfloat   = 0.0_preal
    logical,     dimension(3,3) :: masklogical = .false.
    character(len=pStringLen)   :: mytype      = 'None'
  end type tboundarycondition
 
  type, public :: tloadcase
    type(rotation)               :: rot
    type(tboundarycondition) ::     stress, &                                                       !< stress BC
                                    deformation
    real(preal) ::                  time                   = 0.0_preal                              
    integer ::                      incs                   = 0, &                                   !< number of increments
                                    outputfrequency        = 1, &                                   
                                    restartfrequency       = huge(0), &                             
                                    logscale               = 0                                      
    logical ::                      followformertrajectory = .true.                                 
    integer(kind(FIELD_UNDEFINED_ID)), allocatable :: id(:)
  end type tloadcase
 
  type, public :: tsolutionparams                                                                   
    real(preal), dimension(3,3) :: stress_mask, stress_bc
    type(rotation)              :: rotation_bc
    real(preal) :: timeinc
    real(preal) :: timeincold
  end type tsolutionparams
 
  type, private :: tnumerics
    real(preal) :: &
      fftw_timelimit
    integer :: &
      divergence_correction
    logical :: &
      memory_efficient
    character(len=pStringLen) :: &
      spectral_derivative, &                                                                        !< approximation used for derivatives in Fourier space
      fftw_plan_mode, &                                                                             !< FFTW plan mode, see www.fftw.org
      petsc_options
  end type tnumerics
 
  type(tnumerics) :: num                                                                            ! numerics parameters. Better name?
 
  enum, bind(c); enumerator :: &
    derivative_continuous_id, &
    derivative_central_diff_id, &
    derivative_fwbw_diff_id
  end enum
 
  integer(kind(DERIVATIVE_CONTINUOUS_ID)) :: &
    spectral_derivative_id
 
  public :: &
    utilities_init, &
    utilities_updategamma, &
    utilities_ffttensorforward, &
    utilities_ffttensorbackward, &
    utilities_fftvectorforward, &
    utilities_fftvectorbackward, &
    utilities_fftscalarforward, &
    utilities_fftscalarbackward, &
    utilities_fouriergammaconvolution, &
    utilities_fouriergreenconvolution, &
    utilities_divergencerms, &
    utilities_curlrms, &
    utilities_fourierscalargradient, &
    utilities_fouriervectordivergence, &
    utilities_fouriervectorgradient, &
    utilities_fouriertensordivergence, &
    utilities_maskedcompliance, &
    utilities_constitutiveresponse, &
    utilities_calculaterate, &
    utilities_forwardfield, &
    utilities_updatecoords, &
    utilities_savereferencestiffness, &
    field_undefined_id, &
    field_mech_id, &
    field_thermal_id, &
    field_damage_id
 
contains
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine utilities_init
 
  integer(kind=selected_int_kind(5)) :: ierr
  integer        :: i, j, k, &
    fftw_planner_flag
  integer, dimension(3) :: k_s
  type(c_ptr) :: &
    tensorfield, &                                                                                  !< field containing data for FFTW in real and fourier space (in place)
    vectorfield, &                                                                                  !< field containing data for FFTW in real space when debugging FFTW (no in place)
    scalarfield
  integer(C_INTPTR_T), dimension(3) :: gridfftw
  integer(C_INTPTR_T) :: alloc_local, local_k, local_k_offset
  integer(C_INTPTR_T), parameter :: &
    scalarsize = 1_c_intptr_t, &
    vecsize    = 3_c_intptr_t, &
    tensorsize = 9_c_intptr_t
 
  write(6,'(/,a)') ' <<<+-  spectral_utilities init  -+>>>'
 
  write(6,'(/,a)') ' Diehl, Diploma Thesis TU München, 2010'
  write(6,'(a)')   ' https://doi.org/10.13140/2.1.3234.3840'
 
  write(6,'(/,a)') ' Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013'
  write(6,'(a)')   ' https://doi.org/10.1016/j.ijplas.2012.09.012'
 
  write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
  write(6,'(a)')   ' https://doi.org/10.1016/j.ijplas.2014.02.006'
 
  write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, 2019'
  write(6,'(a)')   ' https://doi.org/10.1007/978-981-10-6855-3_80'
 
!--------------------------------------------------------------------------------------------------
! set debugging parameters
  debuggeneral    = iand(debug_level(debug_spectral),debug_levelbasic)         /= 0
  debugrotation   = iand(debug_level(debug_spectral),debug_spectralrotation)   /= 0
  debugpetsc      = iand(debug_level(debug_spectral),debug_spectralpetsc)      /= 0
 
  if(debugpetsc) write(6,'(3(/,a),/)') &
                 ' Initializing PETSc with debug options: ', &
                 trim(petscdebug), &
                 ' add more using the PETSc_Options keyword in numerics.config '; flush(6)
 
  call petscoptionsclear(petsc_null_options,ierr)
  if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  if(debugpetsc) call petscoptionsinsertstring(petsc_null_options,trim(petscdebug),ierr)
  if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call petscoptionsinsertstring(petsc_null_options,trim(petsc_options),ierr)
  if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
 
  grid1red = grid(1)/2 + 1
  wgt = 1.0/real(product(grid),preal)
 
  write(6,'(/,a,3(i12  ))')  ' grid     a b c: ', grid
  write(6,'(a,3(es12.5))')   ' size     x y z: ', geomsize
 
  num%memory_efficient      = config_numerics%getInt   ('memory_efficient',       defaultval=1) > 0
  num%FFTW_timelimit        = config_numerics%getFloat ('fftw_timelimit',         defaultval=-1.0_preal)
  num%divergence_correction = config_numerics%getInt   ('divergence_correction',  defaultval=2)
  num%spectral_derivative   = config_numerics%getString('spectral_derivative',    defaultval='continuous')
  num%FFTW_plan_mode        = config_numerics%getString('fftw_plan_mode',         defaultval='FFTW_MEASURE')
 
  if (num%divergence_correction < 0 .or. num%divergence_correction > 2) &
    call io_error(301,ext_msg='divergence_correction')
 
  select case (num%spectral_derivative)
    case ('continuous')
      spectral_derivative_id = derivative_continuous_id
    case ('central_difference')
      spectral_derivative_id = derivative_central_diff_id
    case ('fwbw_difference')
      spectral_derivative_id = derivative_fwbw_diff_id
    case default
      call io_error(892,ext_msg=trim(num%spectral_derivative))
  end select
 
!--------------------------------------------------------------------------------------------------
! scale dimension to calculate either uncorrected, dimension-independent, or dimension- and
! resolution-independent divergence
  if (num%divergence_correction == 1) then
    do j = 1, 3
     if (j /= minloc(geomsize,1) .and. j /= maxloc(geomsize,1)) &
       scaledgeomsize = geomsize/geomsize(j)
    enddo
  elseif (num%divergence_correction == 2) then
    do j = 1, 3
     if (      j /= int(minloc(geomsize/real(grid,preal),1)) &
         .and. j /= int(maxloc(geomsize/real(grid,preal),1))) &
       scaledgeomsize = geomsize/geomsize(j)*real(grid(j),preal)
    enddo
  else
    scaledgeomsize = geomsize
  endif
 
 
  select case(io_lc(num%FFTW_plan_mode))                                                            ! setting parameters for the plan creation of FFTW. Basically a translation from fftw3.f
    case('fftw_estimate')                                                                           ! ordered from slow execution (but fast plan creation) to fast execution
      fftw_planner_flag = fftw_estimate
    case('fftw_measure')
      fftw_planner_flag = fftw_measure
    case('fftw_patient')
      fftw_planner_flag = fftw_patient
    case('fftw_exhaustive')
      fftw_planner_flag = fftw_exhaustive
    case default
      call io_warning(warning_id=47,ext_msg=trim(io_lc(num%FFTW_plan_mode)))
      fftw_planner_flag = fftw_measure
  end select
 
!--------------------------------------------------------------------------------------------------
! general initialization of FFTW (see manual on fftw.org for more details)
  if (preal /= c_double .or. kind(1) /= c_int) call io_error(0,ext_msg='Fortran to C')              ! check for correct precision in C
  call fftw_set_timelimit(num%FFTW_timelimit)                                                       ! set timelimit for plan creation
 
  if (debuggeneral) write(6,'(/,a)') ' FFTW initialized'; flush(6)
 
!--------------------------------------------------------------------------------------------------
! MPI allocation
  gridfftw = int(grid,c_intptr_t)
  alloc_local = fftw_mpi_local_size_3d(gridfftw(3), gridfftw(2), gridfftw(1)/2 +1, &
                                       petsc_comm_world, local_k, local_k_offset)
  allocate (xi1st(3,grid1red,grid(2),grid3),source = cmplx(0.0_preal,0.0_preal,preal))             ! frequencies for first derivatives, only half the size for first dimension
  allocate (xi2nd(3,grid1red,grid(2),grid3),source = cmplx(0.0_preal,0.0_preal,preal))             ! frequencies for second derivatives, only half the size for first dimension
 
  tensorfield = fftw_alloc_complex(tensorsize*alloc_local)
  call c_f_pointer(tensorfield, tensorfield_real,    [3_c_intptr_t,3_c_intptr_t, &
                   2_c_intptr_t*(gridfftw(1)/2_c_intptr_t + 1_c_intptr_t),gridfftw(2),local_k])     ! place a pointer for a real tensor representation
  call c_f_pointer(tensorfield, tensorfield_fourier, [3_c_intptr_t,3_c_intptr_t, &
                   gridfftw(1)/2_c_intptr_t + 1_c_intptr_t ,              gridfftw(2),local_k])     ! place a pointer for a fourier tensor representation
 
  vectorfield = fftw_alloc_complex(vecsize*alloc_local)
  call c_f_pointer(vectorfield, vectorfield_real,   [3_c_intptr_t,&
                   2_c_intptr_t*(gridfftw(1)/2_c_intptr_t + 1_c_intptr_t),gridfftw(2),local_k])     ! place a pointer for a real vector representation
  call c_f_pointer(vectorfield, vectorfield_fourier,[3_c_intptr_t,&
                   gridfftw(1)/2_c_intptr_t + 1_c_intptr_t,               gridfftw(2),local_k])     ! place a pointer for a fourier vector representation
 
  scalarfield = fftw_alloc_complex(scalarsize*alloc_local)                                          ! allocate data for real representation (no in place transform)
  call c_f_pointer(scalarfield,    scalarfield_real, &
                   [2_c_intptr_t*(gridfftw(1)/2_c_intptr_t + 1),gridfftw(2),local_k])               ! place a pointer for a real scalar representation
  call c_f_pointer(scalarfield, scalarfield_fourier, &
                    [             gridfftw(1)/2_c_intptr_t + 1 ,gridfftw(2),local_k])               ! place a pointer for a fourier scarlar representation
 
!--------------------------------------------------------------------------------------------------
! tensor MPI fftw plans
  plantensorforth = fftw_mpi_plan_many_dft_r2c(3, [gridfftw(3),gridfftw(2),gridfftw(1)], &          ! dimension, logical length in each dimension in reversed order
                                        tensorsize, fftw_mpi_default_block, fftw_mpi_default_block, &! no. of transforms, default iblock and oblock
                                                       tensorfield_real, tensorfield_fourier, &     ! input data, output data
                                                               petsc_comm_world, fftw_planner_flag) ! use all processors, planer precision
  if (.not. c_associated(plantensorforth)) call io_error(810, ext_msg='planTensorForth')
  plantensorback  = fftw_mpi_plan_many_dft_c2r(3, [gridfftw(3),gridfftw(2),gridfftw(1)], &          ! dimension, logical length in each dimension in reversed order
                                       tensorsize,  fftw_mpi_default_block, fftw_mpi_default_block, &! no. of transforms, default iblock and oblock
                                                        tensorfield_fourier,tensorfield_real, &     ! input data, output data
                                                               petsc_comm_world, fftw_planner_flag) ! all processors, planer precision
  if (.not. c_associated(plantensorback)) call io_error(810, ext_msg='planTensorBack')
 
!--------------------------------------------------------------------------------------------------
! vector MPI fftw plans
  planvectorforth = fftw_mpi_plan_many_dft_r2c(3, [gridfftw(3),gridfftw(2),gridfftw(1)], &          ! dimension, logical length in each dimension in reversed order
                                           vecsize, fftw_mpi_default_block, fftw_mpi_default_block,&! no. of transforms, default iblock and oblock
                                                       vectorfield_real, vectorfield_fourier, &     ! input data, output data
                                                               petsc_comm_world, fftw_planner_flag) ! use all processors, planer precision
  if (.not. c_associated(planvectorforth)) call io_error(810, ext_msg='planVectorForth')
  planvectorback  = fftw_mpi_plan_many_dft_c2r(3, [gridfftw(3),gridfftw(2),gridfftw(1)],  &         ! dimension, logical length in each dimension in reversed order
                                         vecsize, fftw_mpi_default_block, fftw_mpi_default_block, & ! no. of transforms, default iblock and oblock
                                                      vectorfield_fourier,vectorfield_real, &       ! input data, output data
                                                               petsc_comm_world, fftw_planner_flag) ! all processors, planer precision
  if (.not. c_associated(planvectorback)) call io_error(810, ext_msg='planVectorBack')
 
!--------------------------------------------------------------------------------------------------
! scalar MPI fftw plans
  planscalarforth = fftw_mpi_plan_many_dft_r2c(3, [gridfftw(3),gridfftw(2),gridfftw(1)], &          ! dimension, logical length in each dimension in reversed order
                                       scalarsize, fftw_mpi_default_block, fftw_mpi_default_block, &! no. of transforms, default iblock and oblock
                                                       scalarfield_real, scalarfield_fourier, &     ! input data, output data
                                                               petsc_comm_world, fftw_planner_flag) ! use all processors, planer precision
  if (.not. c_associated(planscalarforth)) call io_error(810, ext_msg='planScalarForth')
  planscalarback  = fftw_mpi_plan_many_dft_c2r(3, [gridfftw(3),gridfftw(2),gridfftw(1)], &          ! dimension, logical length in each dimension in reversed order, no. of transforms
                                       scalarsize, fftw_mpi_default_block, fftw_mpi_default_block, &! no. of transforms, default iblock and oblock
                                                       scalarfield_fourier,scalarfield_real, &      ! input data, output data
                                                               petsc_comm_world, fftw_planner_flag) ! use all processors, planer precision
  if (.not. c_associated(planscalarback)) call io_error(810, ext_msg='planScalarBack')
 
!--------------------------------------------------------------------------------------------------
! calculation of discrete angular frequencies, ordered as in FFTW (wrap around)
  do k = grid3offset+1, grid3offset+grid3
    k_s(3) = k - 1
    if(k > grid(3)/2 + 1) k_s(3) = k_s(3) - grid(3)                                                 ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
      do j = 1, grid(2)
        k_s(2) = j - 1
        if(j > grid(2)/2 + 1) k_s(2) = k_s(2) - grid(2)                                             ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
          do i = 1, grid1red
            k_s(1) = i - 1                                                                          ! symmetry, junst running from 0,1,...,N/2,N/2+1
            xi2nd(1:3,i,j,k-grid3offset) = utilities_getfreqderivative(k_s)
            where(mod(grid,2)==0 .and. [i,j,k] == grid/2+1 .and. &
                  spectral_derivative_id == derivative_continuous_id)                               ! for even grids, set the Nyquist Freq component to 0.0
              xi1st(1:3,i,j,k-grid3offset) = cmplx(0.0_preal,0.0_preal,preal)
            elsewhere
              xi1st(1:3,i,j,k-grid3offset) = xi2nd(1:3,i,j,k-grid3offset)
            endwhere
  enddo; enddo; enddo
 
  if(num%memory_efficient) then                                                                     ! allocate just single fourth order tensor
    allocate (gamma_hat(3,3,3,3,1,1,1), source = cmplx(0.0_preal,0.0_preal,preal))
  else                                                                                              ! precalculation of gamma_hat field
    allocate (gamma_hat(3,3,3,3,grid1red,grid(2),grid3), source = cmplx(0.0_preal,0.0_preal,preal))
  endif
 
end subroutine utilities_init
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
subroutine utilities_updategamma(C)
 
  real(preal), intent(in), dimension(3,3,3,3) :: c
  complex(pReal),              dimension(3,3) :: temp33_complex, xidyad_cmplx
  real(preal),                 dimension(6,6) :: a, a_inv
  integer :: &
    i, j, k, &
    l, m, n, o
  logical :: err
 
  c_ref = c
 
  if(.not. num%memory_efficient) then
    gamma_hat =  cmplx(0.0_preal,0.0_preal,preal)                                                   ! for the singular point and any non invertible A
    do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid1red
      if (any([i,j,k] /= 1)) then                                                                   ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
        forall(l = 1:3, m = 1:3) &
          xidyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-grid3offset))*xi1st(m,i,j,k-grid3offset)
        forall(l = 1:3, m = 1:3) &
          temp33_complex(l,m) = sum(cmplx(c_ref(l,1:3,m,1:3),0.0_preal)*xidyad_cmplx)
        a(1:3,1:3) = real(temp33_complex);  a(4:6,4:6) =   real(temp33_complex)
        a(1:3,4:6) = aimag(temp33_complex); a(4:6,1:3) = -aimag(temp33_complex)
        if (abs(math_det33(a(1:3,1:3))) > 1e-16) then
          call math_invert(a_inv, err, a)
          temp33_complex = cmplx(a_inv(1:3,1:3),a_inv(1:3,4:6),preal)
          forall(l=1:3, m=1:3, n=1:3, o=1:3) &
            gamma_hat(l,m,n,o,i,j,k-grid3offset) = temp33_complex(l,n)* &
                                                   conjg(-xi1st(o,i,j,k-grid3offset))*xi1st(m,i,j,k-grid3offset)
        endif
      endif
    enddo; enddo; enddo
  endif
 
end subroutine utilities_updategamma
 
 
!--------------------------------------------------------------------------------------------------
! to 0.0
!--------------------------------------------------------------------------------------------------
subroutine utilities_ffttensorforward
 
  tensorfield_real(1:3,1:3,grid(1)+1:grid1red*2,:,:) = 0.0_preal
  call fftw_mpi_execute_dft_r2c(plantensorforth,tensorfield_real,tensorfield_fourier)
 
end subroutine utilities_ffttensorforward
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine utilities_ffttensorbackward
 
  call fftw_mpi_execute_dft_c2r(plantensorback,tensorfield_fourier,tensorfield_real)
  tensorfield_real = tensorfield_real * wgt                                                         ! normalize the result by number of elements
 
end subroutine utilities_ffttensorbackward
 
!--------------------------------------------------------------------------------------------------
! to 0.0
!--------------------------------------------------------------------------------------------------
subroutine utilities_fftscalarforward
 
  scalarfield_real(grid(1)+1:grid1red*2,:,:) = 0.0_preal
  call fftw_mpi_execute_dft_r2c(planscalarforth,scalarfield_real,scalarfield_fourier)
 
end subroutine utilities_fftscalarforward
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine utilities_fftscalarbackward
 
  call fftw_mpi_execute_dft_c2r(planscalarback,scalarfield_fourier,scalarfield_real)
  scalarfield_real = scalarfield_real * wgt                                                         ! normalize the result by number of elements
 
end subroutine utilities_fftscalarbackward
 
 
!--------------------------------------------------------------------------------------------------
! to 0.0
!--------------------------------------------------------------------------------------------------
subroutine utilities_fftvectorforward
 
  vectorfield_real(1:3,grid(1)+1:grid1red*2,:,:) = 0.0_preal
  call fftw_mpi_execute_dft_r2c(planvectorforth,vectorfield_real,vectorfield_fourier)
 
end subroutine utilities_fftvectorforward
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine utilities_fftvectorbackward
 
  call fftw_mpi_execute_dft_c2r(planvectorback,vectorfield_fourier,vectorfield_real)
  vectorfield_real = vectorfield_real * wgt                                                         ! normalize the result by number of elements
 
end subroutine utilities_fftvectorbackward
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine utilities_fouriergammaconvolution(fieldAim)
 
  real(preal), intent(in), dimension(3,3) :: fieldaim
  complex(pReal),          dimension(3,3) :: temp33_complex, xidyad_cmplx
  real(preal),             dimension(6,6) :: a, a_inv
 
  integer :: &
    i, j, k, &
    l, m, n, o
  logical :: err
 
 
  write(6,'(/,a)') ' ... doing gamma convolution ...............................................'
  flush(6)
 
!--------------------------------------------------------------------------------------------------
! do the actual spectral method calculation (mechanical equilibrium)
  memoryefficient: if(num%memory_efficient) then
    do k = 1, grid3; do j = 1, grid(2); do i = 1, grid1red
      if (any([i,j,k+grid3offset] /= 1)) then                                                       ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
        forall(l = 1:3, m = 1:3) &
          xidyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k))*xi1st(m,i,j,k)
        forall(l = 1:3, m = 1:3) &
          temp33_complex(l,m) = sum(cmplx(c_ref(l,1:3,m,1:3),0.0_preal)*xidyad_cmplx)
        a(1:3,1:3) =  real(temp33_complex); a(4:6,4:6) =   real(temp33_complex)
        a(1:3,4:6) = aimag(temp33_complex); a(4:6,1:3) = -aimag(temp33_complex)
        if (abs(math_det33(a(1:3,1:3))) > 1e-16) then
          call math_invert(a_inv, err, a)
          temp33_complex = cmplx(a_inv(1:3,1:3),a_inv(1:3,4:6),preal)
          forall(l=1:3, m=1:3, n=1:3, o=1:3) &
            gamma_hat(l,m,n,o,1,1,1) =  temp33_complex(l,n)*conjg(-xi1st(o,i,j,k))*xi1st(m,i,j,k)
        else
          gamma_hat(1:3,1:3,1:3,1:3,1,1,1) = cmplx(0.0_preal,0.0_preal,preal)
        endif
        forall(l = 1:3, m = 1:3) &
          temp33_complex(l,m) = sum(gamma_hat(l,m,1:3,1:3,1,1,1)*tensorfield_fourier(1:3,1:3,i,j,k))
        tensorfield_fourier(1:3,1:3,i,j,k) = temp33_complex
      endif
    enddo; enddo; enddo
  else memoryefficient
    do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid1red
      forall(l = 1:3, m = 1:3) &
        temp33_complex(l,m) = sum(gamma_hat(l,m,1:3,1:3,i,j,k) * tensorfield_fourier(1:3,1:3,i,j,k))
      tensorfield_fourier(1:3,1:3,i,j,k) = temp33_complex
    enddo; enddo; enddo
  endif memoryefficient
 
  if (grid3offset == 0) tensorfield_fourier(1:3,1:3,1,1,1) = cmplx(fieldaim/wgt,0.0_preal,preal)
 
end subroutine utilities_fouriergammaconvolution
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine utilities_fouriergreenconvolution(D_ref, mobility_ref, deltaT)
 
  real(preal), dimension(3,3), intent(in) :: d_ref
  real(preal),                 intent(in) :: mobility_ref, deltat
  complex(pReal)                          :: greenop_hat
  integer                                 :: i, j, k
 
!--------------------------------------------------------------------------------------------------
! do the actual spectral method calculation
  do k = 1, grid3; do j = 1, grid(2) ;do i = 1, grid1red
    greenop_hat =  cmplx(1.0_preal,0.0_preal,preal)/ &
                   (cmplx(mobility_ref,0.0_preal,preal) + cmplx(deltat,0.0_preal)*&
                    sum(conjg(xi1st(1:3,i,j,k))* matmul(cmplx(d_ref,0.0_preal),xi1st(1:3,i,j,k))))
    scalarfield_fourier(i,j,k) = scalarfield_fourier(i,j,k)*greenop_hat
  enddo; enddo; enddo
 
end subroutine utilities_fouriergreenconvolution
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
real(pReal) function utilities_divergencerms()
 
  integer :: i, j, k, ierr
  complex(pReal), dimension(3)   :: rescaledgeom
 
  write(6,'(/,a)') ' ... calculating divergence ................................................'
  flush(6)
 
  rescaledgeom = cmplx(geomsize/scaledgeomsize,0.0_preal)
 
!--------------------------------------------------------------------------------------------------
! calculating RMS divergence criterion in Fourier space
  utilities_divergencerms = 0.0_preal
  do k = 1, grid3; do j = 1, grid(2)
    do i = 2, grid1red -1                                                                           ! Has somewhere a conj. complex counterpart. Therefore count it twice.
      utilities_divergencerms = utilities_divergencerms &
            + 2.0_preal*(sum(real(matmul(tensorfield_fourier(1:3,1:3,i,j,k),&                      ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2. do not take square root and square again
                                          conjg(-xi1st(1:3,i,j,k))*rescaledgeom))**2.0_preal)&      ! --> sum squared L_2 norm of vector
                        +sum(aimag(matmul(tensorfield_fourier(1:3,1:3,i,j,k),&
                                          conjg(-xi1st(1:3,i,j,k))*rescaledgeom))**2.0_preal))
    enddo
    utilities_divergencerms = utilities_divergencerms &                                             ! these two layers (DC and Nyquist) do not have a conjugate complex counterpart (if grid(1) /= 1)
               + sum( real(matmul(tensorfield_fourier(1:3,1:3,1       ,j,k), &
                                  conjg(-xi1st(1:3,1,j,k))*rescaledgeom))**2.0_preal) &
               + sum(aimag(matmul(tensorfield_fourier(1:3,1:3,1       ,j,k), &
                                  conjg(-xi1st(1:3,1,j,k))*rescaledgeom))**2.0_preal) &
               + sum( real(matmul(tensorfield_fourier(1:3,1:3,grid1red,j,k), &
                                  conjg(-xi1st(1:3,grid1red,j,k))*rescaledgeom))**2.0_preal) &
               + sum(aimag(matmul(tensorfield_fourier(1:3,1:3,grid1red,j,k), &
                                  conjg(-xi1st(1:3,grid1red,j,k))*rescaledgeom))**2.0_preal)
  enddo; enddo
  if(grid(1) == 1) utilities_divergencerms = utilities_divergencerms * 0.5_preal                    ! counted twice in case of grid(1) == 1
  call mpi_allreduce(mpi_in_place,utilities_divergencerms,1,mpi_double,mpi_sum,petsc_comm_world,ierr)
  if(ierr /=0) call io_error(894, ext_msg='utilities_divergenceRMS')
  utilities_divergencerms = sqrt(utilities_divergencerms) * wgt                                     ! RMS in real space calculated with Parsevals theorem from Fourier space
 
 
end function utilities_divergencerms
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
real(pReal) function utilities_curlrms()
 
  integer  ::  i, j, k, l, ierr
  complex(pReal), dimension(3,3) :: curl_fourier
  complex(pReal), dimension(3)   :: rescaledgeom
 
  write(6,'(/,a)') ' ... calculating curl ......................................................'
  flush(6)
 
  rescaledgeom = cmplx(geomsize/scaledgeomsize,0.0_preal)
 
!--------------------------------------------------------------------------------------------------
! calculating max curl criterion in Fourier space
  utilities_curlrms = 0.0_preal
 
  do k = 1, grid3; do j = 1, grid(2);
    do i = 2, grid1red - 1
      do l = 1, 3
        curl_fourier(l,1) = (+tensorfield_fourier(l,3,i,j,k)*xi1st(2,i,j,k)*rescaledgeom(2) &
                             -tensorfield_fourier(l,2,i,j,k)*xi1st(3,i,j,k)*rescaledgeom(3))
        curl_fourier(l,2) = (+tensorfield_fourier(l,1,i,j,k)*xi1st(3,i,j,k)*rescaledgeom(3) &
                             -tensorfield_fourier(l,3,i,j,k)*xi1st(1,i,j,k)*rescaledgeom(1))
        curl_fourier(l,3) = (+tensorfield_fourier(l,2,i,j,k)*xi1st(1,i,j,k)*rescaledgeom(1) &
                             -tensorfield_fourier(l,1,i,j,k)*xi1st(2,i,j,k)*rescaledgeom(2))
      enddo
      utilities_curlrms = utilities_curlrms &
                        +2.0_preal*sum(real(curl_fourier)**2.0_preal+aimag(curl_fourier)**2.0_preal)! Has somewhere a conj. complex counterpart. Therefore count it twice.
    enddo
    do l = 1, 3
       curl_fourier = (+tensorfield_fourier(l,3,1,j,k)*xi1st(2,1,j,k)*rescaledgeom(2) &
                       -tensorfield_fourier(l,2,1,j,k)*xi1st(3,1,j,k)*rescaledgeom(3))
       curl_fourier = (+tensorfield_fourier(l,1,1,j,k)*xi1st(3,1,j,k)*rescaledgeom(3) &
                       -tensorfield_fourier(l,3,1,j,k)*xi1st(1,1,j,k)*rescaledgeom(1))
       curl_fourier = (+tensorfield_fourier(l,2,1,j,k)*xi1st(1,1,j,k)*rescaledgeom(1) &
                       -tensorfield_fourier(l,1,1,j,k)*xi1st(2,1,j,k)*rescaledgeom(2))
    enddo
    utilities_curlrms = utilities_curlrms &
                      + sum(real(curl_fourier)**2.0_preal + aimag(curl_fourier)**2.0_preal)         ! this layer (DC) does not have a conjugate complex counterpart (if grid(1) /= 1)
    do l = 1, 3
      curl_fourier = (+tensorfield_fourier(l,3,grid1red,j,k)*xi1st(2,grid1red,j,k)*rescaledgeom(2) &
                      -tensorfield_fourier(l,2,grid1red,j,k)*xi1st(3,grid1red,j,k)*rescaledgeom(3))
      curl_fourier = (+tensorfield_fourier(l,1,grid1red,j,k)*xi1st(3,grid1red,j,k)*rescaledgeom(3) &
                      -tensorfield_fourier(l,3,grid1red,j,k)*xi1st(1,grid1red,j,k)*rescaledgeom(1))
      curl_fourier = (+tensorfield_fourier(l,2,grid1red,j,k)*xi1st(1,grid1red,j,k)*rescaledgeom(1) &
                      -tensorfield_fourier(l,1,grid1red,j,k)*xi1st(2,grid1red,j,k)*rescaledgeom(2))
    enddo
    utilities_curlrms = utilities_curlrms &
                      + sum(real(curl_fourier)**2.0_preal + aimag(curl_fourier)**2.0_preal)         ! this layer (Nyquist) does not have a conjugate complex counterpart (if grid(1) /= 1)
  enddo; enddo
 
  call mpi_allreduce(mpi_in_place,utilities_curlrms,1,mpi_double,mpi_sum,petsc_comm_world,ierr)
  if(ierr /=0) call io_error(894, ext_msg='utilities_curlRMS')
  utilities_curlrms = sqrt(utilities_curlrms) * wgt
  if(grid(1) == 1) utilities_curlrms = utilities_curlrms * 0.5_preal                                ! counted twice in case of grid(1) == 1
 
end function utilities_curlrms
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
function utilities_maskedcompliance(rot_BC,mask_stress,C)
 
  real(preal),                dimension(3,3,3,3) :: utilities_maskedcompliance
  real(preal),    intent(in), dimension(3,3,3,3) :: c
  type(rotation), intent(in)                     :: rot_bc
  logical,        intent(in), dimension(3,3)     :: mask_stress
 
  integer :: i, j
  logical, dimension(9)   :: mask_stressvector
  logical, dimension(9,9) :: mask
  real(preal), dimension(9,9) :: temp99_real
  integer :: size_reduced = 0
  real(preal),              dimension(:,:), allocatable ::  &
    s_reduced, &                                                                                    !< reduced compliance matrix (depending on number of stress BC)
    c_reduced, &                                                                                    !< reduced stiffness (depending on number of stress BC)
    stimesc
  logical :: errmatinv
  character(len=pStringLen):: formatstring
 
  mask_stressvector = reshape(transpose(mask_stress), [9])
  size_reduced = count(mask_stressvector)
  if(size_reduced > 0) then
    temp99_real = math_3333to99(rot_bc%rotate(c))
 
    if(debuggeneral) then
      write(6,'(/,a)') ' ... updating masked compliance ............................................'
      write(6,'(/,a,/,9(9(2x,f12.7,1x)/))',advance='no') ' Stiffness C (load) / GPa =',&
                                                   transpose(temp99_real)*1.0e-9_preal
      flush(6)
    endif
 
    do i = 1,9; do j = 1,9
      mask(i,j) = mask_stressvector(i) .and. mask_stressvector(j)
    enddo; enddo
    c_reduced = reshape(pack(temp99_real,mask),[size_reduced,size_reduced])
 
    allocate(s_reduced,mold = c_reduced)
    call math_invert(s_reduced, errmatinv, c_reduced)                                               ! invert reduced stiffness
    if (any(ieee_is_nan(s_reduced))) errmatinv = .true.
    if (errmatinv) call io_error(error_id=400,ext_msg='utilities_maskedCompliance')
 
!--------------------------------------------------------------------------------------------------
! check if inversion was successful
    stimesc = matmul(c_reduced,s_reduced)
    errmatinv = errmatinv .or. any(dneq(stimesc,math_identity2nd(size_reduced),1.0e-12_preal))
    if (debuggeneral .or. errmatinv) then
      write(formatstring, '(i2)') size_reduced
      formatstring = '(/,a,/,'//trim(formatstring)//'('//trim(formatstring)//'(2x,es9.2,1x)/))'
      write(6,trim(formatstring),advance='no') ' C * S (load) ', &
                                                             transpose(matmul(c_reduced,s_reduced))
      write(6,trim(formatstring),advance='no') ' S (load) ', transpose(s_reduced)
      if(errmatinv) call io_error(error_id=400,ext_msg='utilities_maskedCompliance')
    endif
    temp99_real = reshape(unpack(reshape(s_reduced,[size_reduced**2]),reshape(mask,[81]),0.0_preal),[9,9])
  else
    temp99_real = 0.0_preal
  endif
 
  utilities_maskedcompliance = math_99to3333(temp99_real)
 
  if(debuggeneral) then
    write(6,'(/,a,/,9(9(2x,f10.5,1x)/),/)',advance='no') &
      ' Masked Compliance (load) * GPa =', transpose(temp99_real)*1.0e9_preal
    flush(6)
  endif
 
end function utilities_maskedcompliance
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine utilities_fourierscalargradient()
 
  integer :: i, j, k
 
  do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid1red
    vectorfield_fourier(1:3,i,j,k) = scalarfield_fourier(i,j,k)*xi1st(1:3,i,j,k)                    ! ToDo: no -conjg?
  enddo; enddo; enddo
 
end subroutine utilities_fourierscalargradient
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine utilities_fouriervectordivergence()
 
  integer :: i, j, k
 
  do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid1red
    scalarfield_fourier(i,j,k) = sum(vectorfield_fourier(1:3,i,j,k)*conjg(-xi1st(1:3,i,j,k)))
  enddo; enddo; enddo
 
end subroutine utilities_fouriervectordivergence
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine utilities_fouriervectorgradient()
 
  integer :: i, j, k, m, n
 
  do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid1red
    do m = 1, 3; do n = 1, 3
      tensorfield_fourier(m,n,i,j,k) = vectorfield_fourier(m,i,j,k)*xi1st(n,i,j,k)
    enddo; enddo
  enddo; enddo; enddo
 
end subroutine utilities_fouriervectorgradient
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine utilities_fouriertensordivergence()
 
  integer :: i, j, k
 
  do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid1red
    vectorfield_fourier(:,i,j,k) = matmul(tensorfield_fourier(:,:,i,j,k),conjg(-xi1st(:,i,j,k)))
  enddo; enddo; enddo
 
end subroutine utilities_fouriertensordivergence
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine utilities_constitutiveresponse(P,P_av,C_volAvg,C_minmaxAvg,&
                                          F,timeinc,rotation_BC)
 
  real(preal),    intent(out), dimension(3,3,3,3)                   :: c_volavg, c_minmaxavg
  real(preal),    intent(out), dimension(3,3)                       :: p_av
  real(preal),    intent(out), dimension(3,3,grid(1),grid(2),grid3) :: p
  real(preal),    intent(in),  dimension(3,3,grid(1),grid(2),grid3) :: f
  real(preal),    intent(in)                                        :: timeinc
  type(rotation), intent(in),  optional                             :: rotation_bc
 
 
  integer :: &
    i,ierr
  real(preal), dimension(3,3,3,3) :: dpdf_max,      dpdf_min
  real(preal)                     :: dpdf_norm_max, dpdf_norm_min
  real(preal), dimension(2) :: valueandrank
 
  write(6,'(/,a)') ' ... evaluating constitutive response ......................................'
  flush(6)
 
  materialpoint_f  = reshape(f,[3,3,1,product(grid(1:2))*grid3])                                    ! set materialpoint target F to estimated field
 
  call materialpoint_stressanditstangent(.true.,timeinc)                                            ! calculate P field
 
  p = reshape(materialpoint_p, [3,3,grid(1),grid(2),grid3])
  p_av = sum(sum(sum(p,dim=5),dim=4),dim=3) * wgt                                                   ! average of P
  call mpi_allreduce(mpi_in_place,p_av,9,mpi_double,mpi_sum,petsc_comm_world,ierr)
  if (debugrotation) &
  write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',&
                                                      transpose(p_av)*1.e-6_preal
  if(present(rotation_bc)) &
    p_av = rotation_bc%rotate(p_av)
  write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress       / MPa =',&
                                                      transpose(p_av)*1.e-6_preal
  flush(6)
 
  dpdf_max = 0.0_preal
  dpdf_norm_max = 0.0_preal
  dpdf_min = huge(1.0_preal)
  dpdf_norm_min = huge(1.0_preal)
  do i = 1, product(grid(1:2))*grid3
    if (dpdf_norm_max < sum(materialpoint_dpdf(1:3,1:3,1:3,1:3,1,i)**2.0_preal)) then
      dpdf_max = materialpoint_dpdf(1:3,1:3,1:3,1:3,1,i)
      dpdf_norm_max = sum(materialpoint_dpdf(1:3,1:3,1:3,1:3,1,i)**2.0_preal)
    endif
    if (dpdf_norm_min > sum(materialpoint_dpdf(1:3,1:3,1:3,1:3,1,i)**2.0_preal)) then
      dpdf_min = materialpoint_dpdf(1:3,1:3,1:3,1:3,1,i)
      dpdf_norm_min = sum(materialpoint_dpdf(1:3,1:3,1:3,1:3,1,i)**2.0_preal)
    endif
  end do
 
  valueandrank = [dpdf_norm_max,real(worldrank,preal)]
  call mpi_allreduce(mpi_in_place,valueandrank,1, mpi_2double_precision, mpi_maxloc, petsc_comm_world, ierr)
  if (ierr /= 0) call io_error(894, ext_msg='MPI_Allreduce max')
  call mpi_bcast(dpdf_max,81,mpi_double,int(valueandrank(2)),petsc_comm_world, ierr)
  if (ierr /= 0) call io_error(894, ext_msg='MPI_Bcast max')
 
  valueandrank = [dpdf_norm_min,real(worldrank,preal)]
  call mpi_allreduce(mpi_in_place,valueandrank,1, mpi_2double_precision, mpi_minloc, petsc_comm_world, ierr)
  if (ierr /= 0) call io_error(894, ext_msg='MPI_Allreduce min')
  call mpi_bcast(dpdf_min,81,mpi_double,int(valueandrank(2)),petsc_comm_world, ierr)
  if (ierr /= 0) call io_error(894, ext_msg='MPI_Bcast min')
 
  c_minmaxavg = 0.5_preal*(dpdf_max + dpdf_min)
 
  c_volavg = sum(sum(materialpoint_dpdf,dim=6),dim=5)
  call mpi_allreduce(mpi_in_place,c_volavg,81,mpi_double,mpi_sum,petsc_comm_world,ierr)
  c_volavg = c_volavg * wgt
 
 
end subroutine utilities_constitutiveresponse
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
pure function utilities_calculaterate(heterogeneous,field0,field,dt,avRate)
 
  real(preal), intent(in), dimension(3,3) :: &
    avrate
  real(preal), intent(in) :: &
    dt
  logical, intent(in) :: &
    heterogeneous
  real(preal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: &
    field0, &                                                                                       !< data of previous step
    field
  real(preal),             dimension(3,3,grid(1),grid(2),grid3) :: &
    utilities_calculaterate
 
  if (heterogeneous) then
    utilities_calculaterate = (field-field0) / dt
  else
    utilities_calculaterate = spread(spread(spread(avrate,3,grid(1)),4,grid(2)),5,grid3)
  endif
 
end function utilities_calculaterate
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
function utilities_forwardfield(timeinc,field_lastInc,rate,aim)
 
  real(preal), intent(in) :: &
    timeinc
  real(preal), intent(in),           dimension(3,3,grid(1),grid(2),grid3) :: &
    field_lastinc, &                                                                                !< initial field
    rate
  real(preal), intent(in), optional, dimension(3,3) :: &
    aim
  real(preal),                       dimension(3,3,grid(1),grid(2),grid3) :: &
    utilities_forwardfield
  real(preal),                       dimension(3,3)                       :: fielddiff
  integer(kind=selected_int_kind(5)) :: ierr
 
  utilities_forwardfield = field_lastinc + rate*timeinc
  if (present(aim)) then                                                                            
    fielddiff = sum(sum(sum(utilities_forwardfield,dim=5),dim=4),dim=3)*wgt
    call mpi_allreduce(mpi_in_place,fielddiff,9,mpi_double,mpi_sum,petsc_comm_world,ierr)
    fielddiff = fielddiff - aim
    utilities_forwardfield = utilities_forwardfield - &
                     spread(spread(spread(fielddiff,3,grid(1)),4,grid(2)),5,grid3)
  endif
 
end function utilities_forwardfield
 
 
!--------------------------------------------------------------------------------------------------
! standard approach
!--------------------------------------------------------------------------------------------------
pure function utilities_getfreqderivative(k_s)
 
  integer, intent(in),  dimension(3) :: k_s
  complex(pReal),       dimension(3) :: utilities_getfreqderivative
 
  select case (spectral_derivative_id)
    case (derivative_continuous_id)
      utilities_getfreqderivative = cmplx(0.0_preal, 2.0_preal*pi*real(k_s,preal)/geomsize,preal)
 
    case (derivative_central_diff_id)
      utilities_getfreqderivative = cmplx(0.0_preal, sin(2.0_preal*pi*real(k_s,preal)/real(grid,preal)), preal)/ &
                                    cmplx(2.0_preal*geomsize/real(grid,preal), 0.0_preal, preal)
 
    case (derivative_fwbw_diff_id)
      utilities_getfreqderivative(1) = &
                               cmplx(cos(2.0_preal*pi*real(k_s(1),preal)/real(grid(1),preal)) - 1.0_preal, &
                                     sin(2.0_preal*pi*real(k_s(1),preal)/real(grid(1),preal)), preal)* &
                               cmplx(cos(2.0_preal*pi*real(k_s(2),preal)/real(grid(2),preal)) + 1.0_preal, &
                                     sin(2.0_preal*pi*real(k_s(2),preal)/real(grid(2),preal)), preal)* &
                               cmplx(cos(2.0_preal*pi*real(k_s(3),preal)/real(grid(3),preal)) + 1.0_preal, &
                                     sin(2.0_preal*pi*real(k_s(3),preal)/real(grid(3),preal)), preal)/ &
                               cmplx(4.0_preal*geomsize(1)/real(grid(1),preal), 0.0_preal, preal)
      utilities_getfreqderivative(2) = &
                               cmplx(cos(2.0_preal*pi*real(k_s(1),preal)/real(grid(1),preal)) + 1.0_preal, &
                                     sin(2.0_preal*pi*real(k_s(1),preal)/real(grid(1),preal)), preal)* &
                               cmplx(cos(2.0_preal*pi*real(k_s(2),preal)/real(grid(2),preal)) - 1.0_preal, &
                                     sin(2.0_preal*pi*real(k_s(2),preal)/real(grid(2),preal)), preal)* &
                               cmplx(cos(2.0_preal*pi*real(k_s(3),preal)/real(grid(3),preal)) + 1.0_preal, &
                                     sin(2.0_preal*pi*real(k_s(3),preal)/real(grid(3),preal)), preal)/ &
                               cmplx(4.0_preal*geomsize(2)/real(grid(2),preal), 0.0_preal, preal)
      utilities_getfreqderivative(3) = &
                               cmplx(cos(2.0_preal*pi*real(k_s(1),preal)/real(grid(1),preal)) + 1.0_preal, &
                                     sin(2.0_preal*pi*real(k_s(1),preal)/real(grid(1),preal)), preal)* &
                               cmplx(cos(2.0_preal*pi*real(k_s(2),preal)/real(grid(2),preal)) + 1.0_preal, &
                                     sin(2.0_preal*pi*real(k_s(2),preal)/real(grid(2),preal)), preal)* &
                               cmplx(cos(2.0_preal*pi*real(k_s(3),preal)/real(grid(3),preal)) - 1.0_preal, &
                                     sin(2.0_preal*pi*real(k_s(3),preal)/real(grid(3),preal)), preal)/ &
                               cmplx(4.0_preal*geomsize(3)/real(grid(3),preal), 0.0_preal, preal)
  end select
 
end function utilities_getfreqderivative
 
 
!--------------------------------------------------------------------------------------------------
! using integration in Fourier space. Similar as in mesh.f90, but using data already defined for
! convolution
!--------------------------------------------------------------------------------------------------
subroutine utilities_updatecoords(F)
 
  real(preal),   dimension(3,3,grid(1),grid(2),grid3), intent(in) :: f
  real(preal),   dimension(3,  grid(1),grid(2),grid3)             :: ipcoords
  real(preal),   dimension(3,  grid(1),grid(2),grid3+2)           :: ipfluct_padded                  ! Fluctuations of cell center displacement (padded along z for MPI)
  real(preal),   dimension(3,  grid(1)+1,grid(2)+1,grid3+1)       :: nodecoords
  integer :: &
    i,j,k,n, &
    rank_t, &
    rank_b, &
    c, r, &
    ierr
  integer, dimension(MPI_STATUS_SIZE) :: &
    s
  real(preal),   dimension(3)   :: step
  real(preal),   dimension(3,3) :: favg
  integer,       dimension(3) :: me
  integer, dimension(3,8) :: &
    neighbor = reshape([ &
                        0, 0, 0, &
                        1, 0, 0, &
                        1, 1, 0, &
                        0, 1, 0, &
                        0, 0, 1, &
                        1, 0, 1, &
                        1, 1, 1, &
                        0, 1, 1  ], [3,8])
 
  step = geomsize/real(grid, preal)
 !--------------------------------------------------------------------------------------------------
 ! integration in Fourier space to get fluctuations of cell center discplacements
  tensorfield_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = f
  call utilities_ffttensorforward()
 
  do k = 1, grid3; do j = 1, grid(2); do i = 1, grid1red
    if(any([i,j,k+grid3offset] /= 1)) then
      vectorfield_fourier(1:3,i,j,k) = matmul(tensorfield_fourier(1:3,1:3,i,j,k),xi2nd(1:3,i,j,k)) &
                                     / sum(conjg(-xi2nd(1:3,i,j,k))*xi2nd(1:3,i,j,k)) * cmplx(wgt,0.0,preal)
    else
      vectorfield_fourier(1:3,i,j,k) = cmplx(0.0,0.0,preal)
    endif
  enddo; enddo; enddo
 
  call fftw_mpi_execute_dft_c2r(planvectorback,vectorfield_fourier,vectorfield_real)
 
 !--------------------------------------------------------------------------------------------------
 ! average F
  if (grid3offset == 0) favg = real(tensorfield_fourier(1:3,1:3,1,1,1),preal)*wgt
  call mpi_bcast(favg,9,mpi_double,0,petsc_comm_world,ierr)
  if(ierr /=0) call io_error(894, ext_msg='update_IPcoords/MPI_Bcast')
 
 !--------------------------------------------------------------------------------------------------
 ! pad cell center fluctuations along z-direction (needed when running MPI simulation)
  ipfluct_padded(1:3,1:grid(1),1:grid(2),2:grid3+1) = vectorfield_real(1:3,1:grid(1),1:grid(2),1:grid3)
  c = product(shape(ipfluct_padded(:,:,:,1)))                                                        
  rank_t = modulo(worldrank+1,worldsize)
  rank_b = modulo(worldrank-1,worldsize)
 
  ! send bottom layer to process below
  call mpi_isend(ipfluct_padded(:,:,:,2),      c,mpi_double,rank_b,0,petsc_comm_world,r,ierr)
  if(ierr /=0) call io_error(894, ext_msg='update_IPcoords/MPI_Isend')
  call mpi_irecv(ipfluct_padded(:,:,:,grid3+2),c,mpi_double,rank_t,0,petsc_comm_world,r,ierr)
  if(ierr /=0) call io_error(894, ext_msg='update_IPcoords/MPI_Irecv')
  call mpi_wait(r,s,ierr)
 
  ! send top layer to process above
  if(ierr /=0) call io_error(894, ext_msg='update_IPcoords/MPI_Wait')
  call mpi_isend(ipfluct_padded(:,:,:,grid3+1),c,mpi_double,rank_t,0,petsc_comm_world,r,ierr)
  if(ierr /=0) call io_error(894, ext_msg='update_IPcoords/MPI_Isend')
  call mpi_irecv(ipfluct_padded(:,:,:,1),      c,mpi_double,rank_b,0,petsc_comm_world,r,ierr)
  if(ierr /=0) call io_error(894, ext_msg='update_IPcoords/MPI_Irecv')
  call mpi_wait(r,s,ierr)
 
 !--------------------------------------------------------------------------------------------------
 ! calculate nodal displacements
  nodecoords = 0.0_preal
  do k = 0,grid3; do j = 0,grid(2); do i = 0,grid(1)
    nodecoords(1:3,i+1,j+1,k+1) = matmul(favg,step*(real([i,j,k+grid3offset],preal)))
    averagefluct: do n = 1,8
      me = [i+neighbor(1,n),j+neighbor(2,n),k+neighbor(3,n)]
      nodecoords(1:3,i+1,j+1,k+1) = nodecoords(1:3,i+1,j+1,k+1) &
                                  + ipfluct_padded(1:3,modulo(me(1)-1,grid(1))+1,modulo(me(2)-1,grid(2))+1,me(3)+1)*0.125_preal
    enddo averagefluct
  enddo; enddo; enddo
 
 !--------------------------------------------------------------------------------------------------
 ! calculate cell center displacements
  do k = 1,grid3; do j = 1,grid(2); do i = 1,grid(1)
    ipcoords(1:3,i,j,k) = vectorfield_real(1:3,i,j,k) &
                        + matmul(favg,step*(real([i,j,k+grid3offset],preal)-0.5_preal))
  enddo; enddo; enddo
 
  call discretization_setnodecoords(reshape(nodecoords,[3,(grid(1)+1)*(grid(2)+1)*(grid3+1)]))
  call discretization_setipcoords  (reshape(ipcoords,  [3,grid(1)*grid(2)*grid3]))
 
end subroutine utilities_updatecoords
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
subroutine utilities_savereferencestiffness
 
  integer :: &
    fileunit
 
  if (worldrank == 0) then
    write(6,'(a)') ' writing reference stiffness data required for restart to file'; flush(6)
    fileunit = io_open_binary(trim(getsolverjobname())//'.C_ref','w')
    write(fileunit) c_ref
    close(fileunit)
  endif
 
end subroutine utilities_savereferencestiffness
 
end module spectral_utilities