grid_mech_FEM.f90

Go to the documentation of this file.

# 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/grid/grid_mech_FEM.f90"
# 1 "<built-in>"
# 1 "<command-line>"
# 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/grid/grid_mech_FEM.f90"
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
module grid_mech_fem
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscsnes.h" 1
!
!  Include file for Fortran use of the SNES package in 1
!
 
 
 
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscksp.h" 1
!
!
!  Include file for Fortran use of the KSP package in 1
!
 
 
 
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscpc.h" 1
!
!
!  Include file for Fortran use of the PC (preconditioner) package in 1
!
 
 
 
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscmat.h" 1
!
!
!  Include file for Fortran use of the Mat package in 1
!
 
 
 
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscvec.h" 1
!
!
!  Include file for Fortran use of the Vec package in 1
!
 
 
 
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscao.h" 1
!
!
!  Include file for Fortran use of the AO (application ordering) package in 1
!
 
 
 
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscis.h" 1
!
!
!  Include file for Fortran use of the IS (index set) package in 1
!
 
 
 
 
# 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 "/opt/petsc-3.10.3/include/petsc/finclude/petscis.h" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscviewer.h" 1
!
!  Include file for Fortran use of the type(tPetscViewer) package in 1
!
# 10 "/opt/petsc-3.10.3/include/petsc/finclude/petscis.h" 2
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
# 9 "/opt/petsc-3.10.3/include/petsc/finclude/petscao.h" 2
 
 
 
 
 
 
 
 
 
# 9 "/opt/petsc-3.10.3/include/petsc/finclude/petscvec.h" 2
 
 
 
 
 
 
 
 
 
 
 
 
 
# 30 "/opt/petsc-3.10.3/include/petsc/finclude/petscvec.h"
 
 
 
# 9 "/opt/petsc-3.10.3/include/petsc/finclude/petscmat.h" 2
 
 
 
 
 
# 49 "/opt/petsc-3.10.3/include/petsc/finclude/petscmat.h"
 
 
 
 
 
 
 
 
 
 
 
 
 
# 71 "/opt/petsc-3.10.3/include/petsc/finclude/petscmat.h"
!
!  Matrix types
!
# 115 "/opt/petsc-3.10.3/include/petsc/finclude/petscmat.h"
 
# 140 "/opt/petsc-3.10.3/include/petsc/finclude/petscmat.h"
 
!
! MatMFFDType values
!
 
 
 
!
! MatSolverTypes
!
# 168 "/opt/petsc-3.10.3/include/petsc/finclude/petscmat.h"
 
!
! GPU Storage Formats for CUSPARSE
!
 
 
 
!
! sparsity reducing ordering for STRUMPACK
!
 
 
# 9 "/opt/petsc-3.10.3/include/petsc/finclude/petscpc.h" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscdm.h" 1
 
!
!  Include file for Fortran use of the DM package in 1
!
 
 
 
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscis.h" 1
!
!
!  Include file for Fortran use of the type(tIS) (index set) package in 1
!
# 9 "/opt/petsc-3.10.3/include/petsc/finclude/petscdm.h" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscvec.h" 1
!
!
!  Include file for Fortran use of the type(tVec) package in 1
!
# 10 "/opt/petsc-3.10.3/include/petsc/finclude/petscdm.h" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscmat.h" 1
!
!
!  Include file for Fortran use of the type(tMat) package in 1
!
# 11 "/opt/petsc-3.10.3/include/petsc/finclude/petscdm.h" 2
 
 
 
 
 
 
 
 
 
 
 
# 32 "/opt/petsc-3.10.3/include/petsc/finclude/petscdm.h"
 
 
# 47 "/opt/petsc-3.10.3/include/petsc/finclude/petscdm.h"
 
# 10 "/opt/petsc-3.10.3/include/petsc/finclude/petscpc.h" 2
 
 
 
# 26 "/opt/petsc-3.10.3/include/petsc/finclude/petscpc.h"
!
! GAMG types
!
 
 
!
!  Various preconditioners
!
# 74 "/opt/petsc-3.10.3/include/petsc/finclude/petscpc.h"
 
# 9 "/opt/petsc-3.10.3/include/petsc/finclude/petscksp.h" 2
 
 
 
 
# 21 "/opt/petsc-3.10.3/include/petsc/finclude/petscksp.h"
!
!  Various Krylov subspace methods
!
# 57 "/opt/petsc-3.10.3/include/petsc/finclude/petscksp.h"
!
!  Various Initial guesses for Krylov subspace methods
!
# 8 "/opt/petsc-3.10.3/include/petsc/finclude/petscsnes.h" 2
 
 
 
 
 
# 30 "/opt/petsc-3.10.3/include/petsc/finclude/petscsnes.h"
 
!
!  character*(80)
!
# 47 "/opt/petsc-3.10.3/include/petsc/finclude/petscsnes.h"
 
!
! character*(80)
!
 
 
 
 
 
 
 
!
! integer(kind=selected_int_kind(5))
!
 
 
 
 
 
 
!
!  character*(80)
!
# 78 "/opt/petsc-3.10.3/include/petsc/finclude/petscsnes.h"
 
 
 
# 9 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/grid/grid_mech_FEM.f90" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscdmda.h" 1
 
!
!  Include file for Fortran use of the 'da' (distributed array) package in 1
!
 
 
 
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscis.h" 1
!
!
!  Include file for Fortran use of the type(tIS) (index set) package in 1
!
# 9 "/opt/petsc-3.10.3/include/petsc/finclude/petscdmda.h" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscvec.h" 1
!
!
!  Include file for Fortran use of the type(tVec) package in 1
!
# 10 "/opt/petsc-3.10.3/include/petsc/finclude/petscdmda.h" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscmat.h" 1
!
!
!  Include file for Fortran use of the type(tMat) package in 1
!
# 11 "/opt/petsc-3.10.3/include/petsc/finclude/petscdmda.h" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscdm.h" 1
 
!
!  Include file for Fortran use of the type(tDM) package in 1
!
# 12 "/opt/petsc-3.10.3/include/petsc/finclude/petscdmda.h" 2
 
 
 
 
 
!
!   DMDA_LOCAL_INFO_SIZE is one large than the size incase the DA is larger than an integer (on 64 bit systems).
!   non-int fields are not accessiable from fortran.
!
# 39 "/opt/petsc-3.10.3/include/petsc/finclude/petscdmda.h"
 
 
 
 
 
 
 
 
 
 
 
 
 
# 10 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/grid/grid_mech_FEM.f90" 2
  use petscdmda
  use petscsnes
 
  use prec
  use damask_interface
  use hdf5_utilities
  use math
  use spectral_utilities
  use fesolving
  use numerics
  use homogenization
  use discretization
  use discretization_grid
  use debug
 
  implicit none
  private
 
!--------------------------------------------------------------------------------------------------
! derived types
  type(tsolutionparams), private :: params
 
!--------------------------------------------------------------------------------------------------
! 1 data
  type(tdm),   private :: mech_grid
  type(tsnes), private :: mech_snes
  type(tvec),  private :: solution_current, solution_lastinc, solution_rate
 
!--------------------------------------------------------------------------------------------------
! common pointwise data
  real(preal), private, dimension(:,:,:,:,:), allocatable ::  f, p_current, f_lastinc
  real(preal), private :: detj
  real(preal), private, dimension(3)   :: delta
  real(preal), private, dimension(3,8) :: bmat
  real(preal), private, dimension(8,8) :: hgmat
  integer(kind=selected_int_kind(5)),    private :: xstart,ystart,zstart,xend,yend,zend
 
!--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc.
  real(preal), private, dimension(3,3) :: &
    f_aimdot = 0.0_preal, &                                                                         
    f_aim = math_i3, &                                                                              
    f_aim_lastiter = math_i3, &
    f_aim_lastinc  = math_i3, &                                                                     
    p_av = 0.0_preal                                                                                
 
  character(len=pStringLen), private :: incinfo
 
  real(preal), private, dimension(3,3,3,3) :: &
    c_volavg = 0.0_preal, &                                                                         
    c_volavglastinc = 0.0_preal, &                                                                  
    s = 0.0_preal                                                                                   
 
  real(preal), private :: &
    err_bc
 
  integer, private :: &
    totaliter = 0
 
  public :: &
    grid_mech_fem_init, &
    grid_mech_fem_solution, &
    grid_mech_fem_forward, &
    grid_mech_fem_updatecoords, &
    grid_mech_fem_restartwrite
 
contains
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_fem_init
 
  real(preal) :: hgcoeff = 0.0e-2_preal
  integer(kind=selected_int_kind(5)), dimension(0:worldsize-1) :: localk
  real(preal), dimension(3,3) :: &
    temp33_real = 0.0_preal
  real(preal), dimension(4,8) :: &
    hgcomp = reshape([ 1.0_preal, 1.0_preal, 1.0_preal,-1.0_preal, &
                       1.0_preal,-1.0_preal,-1.0_preal, 1.0_preal, &
                      -1.0_preal, 1.0_preal,-1.0_preal, 1.0_preal, &
                      -1.0_preal,-1.0_preal, 1.0_preal,-1.0_preal, &
                      -1.0_preal,-1.0_preal, 1.0_preal, 1.0_preal, &
                      -1.0_preal, 1.0_preal,-1.0_preal,-1.0_preal, &
                       1.0_preal,-1.0_preal,-1.0_preal,-1.0_preal, &
                       1.0_preal, 1.0_preal, 1.0_preal, 1.0_preal], [4,8])
  integer(kind=selected_int_kind(5)) :: ierr
  integer(HID_T) :: filehandle, grouphandle
  character(len=pStringLen) :: filename
  real(preal), dimension(3,3,3,3) :: devnull
  real(kind=selected_real_kind(10)), pointer, dimension(:,:,:,:) :: &
  u_current,u_lastinc
 
  write(6,'(/,a)') ' <<<+-  grid_mech_FEM init  -+>>>'; flush(6)
 
!--------------------------------------------------------------------------------------------------
! set default and user defined options for 1
  call petscoptionsinsertstring(petsc_null_options,'-mech_snes_type newtonls -mech_ksp_type fgmres &
                                &-mech_ksp_max_it 25 -mech_pc_type ml -mech_mg_levels_ksp_type chebyshev',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
 
!--------------------------------------------------------------------------------------------------
! allocate global fields
  allocate(f(3,3,grid(1),grid(2),grid3),source = 0.0_preal)
  allocate(p_current(3,3,grid(1),grid(2),grid3),source = 0.0_preal)
  allocate(f_lastinc(3,3,grid(1),grid(2),grid3),source = 0.0_preal)
 
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of 1
  call snescreate(petsc_comm_world,mech_snes,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call snessetoptionsprefix(mech_snes,'mech_',ierr);if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  localk            = 0
  localk(worldrank) = grid3
  call mpi_allreduce(mpi_in_place,localk,worldsize,mpi_integer,mpi_sum,petsc_comm_world,ierr)
  call dmdacreate3d(petsc_comm_world, &
         dm_boundary_periodic, dm_boundary_periodic, dm_boundary_periodic, &
         dmda_stencil_box, &
         grid(1),grid(2),grid(3), &
         1, 1, worldsize, &
         3, 1, &
         [grid(1)],[grid(2)],localk, &
         mech_grid,ierr)
  if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call dmdasetuniformcoordinates(mech_grid,0.0_preal,geomsize(1),0.0_preal,geomsize(2),0.0_preal,geomsize(3),ierr)
  if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call snessetdm(mech_snes,mech_grid,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call dmsetfromoptions(mech_grid,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call dmsetup(mech_grid,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call dmcreateglobalvector(mech_grid,solution_current,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call dmcreateglobalvector(mech_grid,solution_lastinc,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call dmcreateglobalvector(mech_grid,solution_rate   ,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call dmsnessetfunctionlocal(mech_grid,formresidual,petsc_null_snes,ierr)
  if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call dmsnessetjacobianlocal(mech_grid,formjacobian,petsc_null_snes,ierr)
  if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call snessetconvergencetest(mech_snes,converged,petsc_null_snes,petsc_null_function,ierr)         ! specify custom convergence check function "_converged"
  if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call snessetmaxlinearsolvefailures(mech_snes, huge(1), ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif                       ! ignore linear solve failures
  call snessetfromoptions(mech_snes,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif                                            ! pull it all together with additional cli arguments
 
!--------------------------------------------------------------------------------------------------
! init fields
  call vecset(solution_current,0.0_preal,ierr);if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call vecset(solution_lastinc,0.0_preal,ierr);if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call vecset(solution_rate   ,0.0_preal,ierr);if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call dmdavecgetarrayf90(mech_grid,solution_current,u_current,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call dmdavecgetarrayf90(mech_grid,solution_lastinc,u_lastinc,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
 
  call dmdagetcorners(mech_grid,xstart,ystart,zstart,xend,yend,zend,ierr)                           ! local grid extent
  if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  xend = xstart+xend-1
  yend = ystart+yend-1
  zend = zstart+zend-1
  delta = geomsize/real(grid,preal)                                                                 ! grid spacing
  detj = product(delta)                                                                             ! cell volume
 
  bmat = reshape(real([-1.0_preal/delta(1),-1.0_preal/delta(2),-1.0_preal/delta(3), &
                        1.0_preal/delta(1),-1.0_preal/delta(2),-1.0_preal/delta(3), &
                       -1.0_preal/delta(1), 1.0_preal/delta(2),-1.0_preal/delta(3), &
                        1.0_preal/delta(1), 1.0_preal/delta(2),-1.0_preal/delta(3), &
                       -1.0_preal/delta(1),-1.0_preal/delta(2), 1.0_preal/delta(3), &
                        1.0_preal/delta(1),-1.0_preal/delta(2), 1.0_preal/delta(3), &
                       -1.0_preal/delta(1), 1.0_preal/delta(2), 1.0_preal/delta(3), &
                        1.0_preal/delta(1), 1.0_preal/delta(2), 1.0_preal/delta(3)],preal), [3,8])/4.0_preal ! shape function derivative matrix
 
  hgmat = matmul(transpose(hgcomp),hgcomp) &
        * hgcoeff*(delta(1)*delta(2) + delta(2)*delta(3) + delta(3)*delta(1))/16.0_preal            ! hourglass stabilization matrix
 
!--------------------------------------------------------------------------------------------------
! init fields
  restartread: if (interface_restartinc > 0) then
    write(6,'(/,a,i0,a)') ' reading restart data of increment ', interface_restartinc, ' from file'
 
    write(filename,'(a,a,i0,a)') trim(getsolverjobname()),'_',worldrank,'.hdf5'
    filehandle  = hdf5_openfile(filename)
    grouphandle = hdf5_opengroup(filehandle,'solver')
 
    call hdf5_read(grouphandle,f_aim,        'F_aim')
    call hdf5_read(grouphandle,f_aim_lastinc,'F_aim_lastInc')
    call hdf5_read(grouphandle,f_aimdot,     'F_aimDot')
    call hdf5_read(grouphandle,f,            'F')
    call hdf5_read(grouphandle,f_lastinc,    'F_lastInc')
    call hdf5_read(grouphandle,u_current,    'u')
    call hdf5_read(grouphandle,u_lastinc,    'u_lastInc')
 
  elseif (interface_restartinc == 0) then restartread
    f_lastinc = spread(spread(spread(math_i3,3,grid(1)),4,grid(2)),5,grid3)                         ! initialize to identity
    f         = spread(spread(spread(math_i3,3,grid(1)),4,grid(2)),5,grid3)
  endif restartread
  materialpoint_f0 = reshape(f_lastinc, [3,3,1,product(grid(1:2))*grid3])                           ! set starting condition for materialpoint_stressAndItsTangent
  call utilities_updatecoords(f)
  call utilities_constitutiveresponse(p_current,temp33_real,c_volavg,devnull, &                     ! stress field, stress avg, global average of stiffness and (min+max)/2
                                      f, &                                                          ! target F
                                      0.0_preal)                                                    ! time increment
  call dmdavecrestorearrayf90(mech_grid,solution_current,u_current,ierr)
  if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call dmdavecrestorearrayf90(mech_grid,solution_lastinc,u_lastinc,ierr)
  if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
 
  restartread2: if (interface_restartinc > 0) then
    write(6,'(/,a,i0,a)') ' reading more restart data of increment ', interface_restartinc, ' from file'
    call hdf5_read(grouphandle,c_volavg,       'C_volAvg')
    call hdf5_read(grouphandle,c_volavglastinc,'C_volAvgLastInc')
 
    call hdf5_closegroup(grouphandle)
    call hdf5_closefile(filehandle)
 
  endif restartread2
 
end subroutine grid_mech_fem_init
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
function grid_mech_fem_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution)
 
!--------------------------------------------------------------------------------------------------
! input data for solution
  character(len=*),            intent(in) :: &
    incinfoin
  real(preal),                 intent(in) :: &
    timeinc, &                                                                                      !< time increment of current solution
    timeinc_old
  type(tboundarycondition),    intent(in) :: &
    stress_bc
  type(rotation),              intent(in) :: &
    rotation_bc
  type(tsolutionstate)                    :: &
    solution
!--------------------------------------------------------------------------------------------------
! 1 Data
  integer(kind=selected_int_kind(5)) :: ierr
  integer(kind=selected_int_kind(5)) :: reason
 
  incinfo = incinfoin
 
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
  s = utilities_maskedcompliance(rotation_bc,stress_bc%maskLogical,c_volavg)
!--------------------------------------------------------------------------------------------------
! set module wide available data
  params%stress_mask = stress_bc%maskFloat
  params%stress_BC   = stress_bc%values
  params%rotation_BC = rotation_bc
  params%timeinc     = timeinc
  params%timeincOld  = timeinc_old
 
!--------------------------------------------------------------------------------------------------
! solve BVP
  call snessolve(mech_snes,petsc_null_vec,solution_current,ierr);if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
 
!--------------------------------------------------------------------------------------------------
! check convergence
  call snesgetconvergedreason(mech_snes,reason,ierr);if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
 
  solution%converged = reason > 0
  solution%iterationsNeeded = totaliter
  solution%termIll = terminallyill
  terminallyill = .false.
 
end function grid_mech_fem_solution
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_fem_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,&
                                 deformation_BC,stress_BC,rotation_BC)
 
  logical,                     intent(in) :: &
    cutback, &
    guess
  real(preal),                 intent(in) :: &
    timeinc_old, &
    timeinc, &
    loadcasetime
  type(tboundarycondition),    intent(in) :: &
    stress_bc, &
    deformation_bc
  type(rotation),              intent(in) :: &
    rotation_bc
  integer(kind=selected_int_kind(5)) :: ierr
  real(kind=selected_real_kind(10)), pointer, dimension(:,:,:,:) :: &
    u_current,u_lastinc
 
  call dmdavecgetarrayf90(mech_grid,solution_current,u_current,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call dmdavecgetarrayf90(mech_grid,solution_lastinc,u_lastinc,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
 
  if (cutback) then
    c_volavg = c_volavglastinc
  else
    c_volavglastinc    = c_volavg
 
    f_aimdot = merge(stress_bc%maskFloat*(f_aim-f_aim_lastinc)/timeinc_old, 0.0_preal, guess)
    f_aim_lastinc = f_aim
 
    !--------------------------------------------------------------------------------------------------
    ! calculate rate for aim
    if     (deformation_bc%myType=='l') then                                                        ! calculate F_aimDot from given L and current F
      f_aimdot = &
      f_aimdot + deformation_bc%maskFloat * matmul(deformation_bc%values, f_aim_lastinc)
    elseif(deformation_bc%myType=='fdot') then                                                      ! F_aimDot is prescribed
      f_aimdot = &
      f_aimdot + deformation_bc%maskFloat * deformation_bc%values
    elseif (deformation_bc%myType=='f') then                                                        ! aim at end of load case is prescribed
      f_aimdot = &
      f_aimdot + deformation_bc%maskFloat * (deformation_bc%values - f_aim_lastinc)/loadcasetime
    endif
 
    if (guess) then
      call vecwaxpy(solution_rate,-1.0_preal,solution_lastinc,solution_current,ierr)
      if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
      call vecscale(solution_rate,1.0_preal/timeinc_old,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
    else
      call vecset(solution_rate,0.0_preal,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
    endif
    call veccopy(solution_current,solution_lastinc,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
 
    f_lastinc = f
 
    materialpoint_f0 = reshape(f, [3,3,1,product(grid(1:2))*grid3])
  endif
 
!--------------------------------------------------------------------------------------------------
! update average and local deformation gradients
  f_aim = f_aim_lastinc + f_aimdot * timeinc
  call vecaxpy(solution_current,timeinc,solution_rate,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
 
  call dmdavecrestorearrayf90(mech_grid,solution_current,u_current,ierr);if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call dmdavecrestorearrayf90(mech_grid,solution_lastinc,u_lastinc,ierr);if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
 
end subroutine grid_mech_fem_forward
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_fem_updatecoords
 
  call utilities_updatecoords(f)
 
end subroutine grid_mech_fem_updatecoords
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_fem_restartwrite
 
  integer(kind=selected_int_kind(5)) :: ierr
  integer(HID_T) :: filehandle, grouphandle
  real(kind=selected_real_kind(10)), dimension(:,:,:,:), pointer :: u_current,u_lastinc
  character(len=pStringLen) :: filename
 
  call dmdavecgetarrayf90(mech_grid,solution_current,u_current,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call dmdavecgetarrayf90(mech_grid,solution_lastinc,u_lastinc,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
 
  write(6,'(a)') ' writing solver data required for restart to file'; flush(6)
 
  write(filename,'(a,a,i0,a)') trim(getsolverjobname()),'_',worldrank,'.hdf5'
  filehandle  = hdf5_openfile(filename,'w')
  grouphandle = hdf5_addgroup(filehandle,'solver')
 
  call hdf5_write(grouphandle,f_aim,        'F_aim')
  call hdf5_write(grouphandle,f_aim_lastinc,'F_aim_lastInc')
  call hdf5_write(grouphandle,f_aimdot,     'F_aimDot')
  call hdf5_write(grouphandle,f,            'F')
  call hdf5_write(grouphandle,f_lastinc,    'F_lastInc')
  call hdf5_write(grouphandle,u_current,    'u')
  call hdf5_write(grouphandle,u_lastinc,    'u_lastInc')
 
  call hdf5_write(grouphandle,c_volavg,       'C_volAvg')
  call hdf5_write(grouphandle,c_volavglastinc,'C_volAvgLastInc')
 
  call hdf5_closegroup(grouphandle)
  call hdf5_closefile(filehandle)
 
  call dmdavecrestorearrayf90(mech_grid,solution_current,u_current,ierr);if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call dmdavecrestorearrayf90(mech_grid,solution_lastinc,u_lastinc,ierr);if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
 
end subroutine grid_mech_fem_restartwrite
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,ierr)
 
  type(tsnes) :: snes_local
   integer(kind=selected_int_kind(5)),  intent(in) :: PETScIter
   real(kind=selected_real_kind(10)), intent(in) :: &
     devnull1, &
     devnull2, &
     fnorm
  integer(kind=selected_int_kind(5)) :: reason
  integer(kind=selected_int_kind(10)) :: dummy
  integer(kind=selected_int_kind(5)) :: ierr
  real(pReal) :: &
    err_div, &
    divTol, &
    BCTol
 
  err_div = fnorm*sqrt(wgt)*geomsize(1)/scaledgeomsize(1)/detj
  divtol = max(maxval(abs(p_av))*err_div_tolrel   ,err_div_tolabs)
  bctol  = max(maxval(abs(p_av))*err_stress_tolrel,err_stress_tolabs)
 
  if ((totaliter >= itmin .and. &
                            all([ err_div/divtol, &
                                  err_bc /bctol       ] < 1.0_preal)) &
              .or.    terminallyill) then
    reason = 1
  elseif (totaliter >= itmax) then
    reason = -1
  else
    reason = 0
  endif
 
!--------------------------------------------------------------------------------------------------
! report
  write(6,'(1/,a)') ' ... reporting .............................................................'
  write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', &
          err_div/divtol,  ' (',err_div,' / m, tol = ',divtol,')'
  write(6,'(a,f12.2,a,es8.2,a,es9.2,a)')    ' error stress BC  = ', &
          err_bc/bctol,    ' (',err_bc, ' Pa,  tol = ',bctol,')'
  write(6,'(/,a)') ' ==========================================================================='
  flush(6)
 
end subroutine converged
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine formresidual(da_local,x_local, &
                        f_local,dummy,ierr)
 
  type(tdm)                   :: da_local
  type(tvec)                  :: x_local, f_local
  real(kind=selected_real_kind(10)), pointer,dimension(:,:,:,:) :: x_scal, f_scal
  real(kind=selected_real_kind(10)), dimension(8,3) :: x_elem,  f_elem
  integer(kind=selected_int_kind(5))             :: i, ii, j, jj, k, kk, ctr, ele
  real(pReal), dimension(3,3) :: &
    deltaF_aim
  integer(kind=selected_int_kind(5)) :: &
    PETScIter, &
    nfuncs
  integer(kind=selected_int_kind(10)) :: dummy
  integer(kind=selected_int_kind(5)) :: ierr
  real(pReal), dimension(3,3,3,3) :: devNull
 
 
  call snesgetnumberfunctionevals(mech_snes,nfuncs,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call snesgetiterationnumber(mech_snes,petsciter,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
 
  if (nfuncs == 0 .and. petsciter == 0) totaliter = -1                                              ! new increment
 
!--------------------------------------------------------------------------------------------------
! begin of new iteration
  newiteration: if (totaliter <= petsciter) then
    totaliter = totaliter + 1
    write(6,'(1x,a,3(a,i0))') trim(incinfo), ' @ Iteration ', itmin, '≤',totaliter+1, '≤', itmax
    if (iand(debug_level(debug_spectral),debug_spectralrotation) /= 0) &
      write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
              ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(f_aim,active=.true.))
    write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
              ' deformation gradient aim       =', transpose(f_aim)
    flush(6)
  endif newiteration
 
!--------------------------------------------------------------------------------------------------
! get deformation gradient
  call dmdavecgetarrayf90(da_local,x_local,x_scal,ierr);if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  do k = zstart, zend; do j = ystart, yend; do i = xstart, xend
    ctr = 0
    do kk = 0, 1; do jj = 0, 1; do ii = 0, 1
      ctr = ctr + 1
      x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk)
    enddo; enddo; enddo
    ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1
    f(1:3,1:3,ii,jj,kk) = params%rotation_BC%rotate(f_aim,active=.true.) + transpose(matmul(bmat,x_elem))
  enddo; enddo; enddo
  call dmdavecrestorearrayf90(da_local,x_local,x_scal,ierr);if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
 
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
  call utilities_constitutiveresponse(p_current,&
                                      p_av,c_volavg,devnull, &
                                      f,params%timeinc,params%rotation_BC)
  call mpi_allreduce(mpi_in_place,terminallyill,1,mpi_logical,mpi_lor,petsc_comm_world,ierr)
 
!--------------------------------------------------------------------------------------------------
! stress BC handling
  f_aim_lastiter = f_aim
  deltaf_aim = math_mul3333xx33(s, p_av - params%stress_BC)
  f_aim = f_aim - deltaf_aim
  err_bc = maxval(abs(params%stress_mask * (p_av - params%stress_BC)))                              ! mask = 0.0 when no stress bc
 
!--------------------------------------------------------------------------------------------------
! constructing residual
  call vecset(f_local,0.0_preal,ierr);if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call dmdavecgetarrayf90(da_local,f_local,f_scal,ierr);if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call dmdavecgetarrayf90(da_local,x_local,x_scal,ierr);if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  ele = 0
  do k = zstart, zend; do j = ystart, yend; do i = xstart, xend
    ctr = 0
    do kk = 0, 1; do jj = 0, 1; do ii = 0, 1
      ctr = ctr + 1
      x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk)
    enddo; enddo; enddo
    ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1
    ele = ele + 1
    f_elem = matmul(transpose(bmat),transpose(p_current(1:3,1:3,ii,jj,kk)))*detj + &
             matmul(hgmat,x_elem)*(materialpoint_dpdf(1,1,1,1,1,ele) + &
                                   materialpoint_dpdf(2,2,2,2,1,ele) + &
                                   materialpoint_dpdf(3,3,3,3,1,ele))/3.0_preal
    ctr = 0
    do kk = 0, 1; do jj = 0, 1; do ii = 0, 1
      ctr = ctr + 1
      f_scal(0:2,i+ii,j+jj,k+kk) = f_scal(0:2,i+ii,j+jj,k+kk) + f_elem(ctr,1:3)
    enddo; enddo; enddo
  enddo; enddo; enddo
  call dmdavecrestorearrayf90(da_local,x_local,x_scal,ierr);if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call dmdavecrestorearrayf90(da_local,f_local,f_scal,ierr);if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
 
!--------------------------------------------------------------------------------------------------
! applying boundary conditions
  call dmdavecgetarrayf90(da_local,f_local,f_scal,ierr);if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  if (zstart == 0) then
    f_scal(0:2,xstart,ystart,zstart) = 0.0
    f_scal(0:2,xend+1,ystart,zstart) = 0.0
    f_scal(0:2,xstart,yend+1,zstart) = 0.0
    f_scal(0:2,xend+1,yend+1,zstart) = 0.0
  endif
  if (zend + 1 == grid(3)) then
    f_scal(0:2,xstart,ystart,zend+1) = 0.0
    f_scal(0:2,xend+1,ystart,zend+1) = 0.0
    f_scal(0:2,xstart,yend+1,zend+1) = 0.0
    f_scal(0:2,xend+1,yend+1,zend+1) = 0.0
  endif
  call dmdavecrestorearrayf90(da_local,f_local,f_scal,ierr);if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
 
end subroutine formresidual
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine formjacobian(da_local,x_local,Jac_pre,Jac,dummy,ierr)
 
  type(tdm)                                   :: da_local
  type(tvec)                                  :: x_local, coordinates
  type(tmat)                                  :: Jac_pre, Jac
  integer(kind=selected_int_kind(5)),dimension(4,24)           :: row, col
  real(kind=selected_real_kind(10)),pointer,dimension(:,:,:,:) :: x_scal
  real(kind=selected_real_kind(10)),dimension(24,24)         :: k_ele
  real(kind=selected_real_kind(10)),dimension(9,24)          :: bmatfull
  integer(kind=selected_int_kind(5))                             :: i, ii, j, jj, k, kk, ctr, ele
  integer(kind=selected_int_kind(5)),dimension(3),parameter      :: rows = [0, 1, 2]
  real(kind=selected_real_kind(10))                          :: diag
  integer(kind=selected_int_kind(10))                          :: dummy
  type(tmatnullspace)                         :: matnull
  integer(kind=selected_int_kind(5))                       :: ierr
 
  bmatfull = 0.0
  bmatfull(1:3,1 :8 ) = bmat
  bmatfull(4:6,9 :16) = bmat
  bmatfull(7:9,17:24) = bmat
  call matsetoption(jac,mat_keep_nonzero_pattern,petsc_true,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call matsetoption(jac,mat_new_nonzero_allocation_err,petsc_false,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call matzeroentries(jac,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  ele = 0
  do k = zstart, zend; do j = ystart, yend; do i = xstart, xend
    ctr = 0
    do kk = 0, 1; do jj = 0, 1; do ii = 0, 1
      ctr = ctr + 1
      col(3,ctr   ) = i+ii
      col(2,ctr   ) = j+jj
      col(1,ctr   ) = k+kk
      col(4,ctr   ) = 0
      col(3,ctr+8 ) = i+ii
      col(2,ctr+8 ) = j+jj
      col(1,ctr+8 ) = k+kk
      col(4,ctr+8 ) = 1
      col(3,ctr+16) = i+ii
      col(2,ctr+16) = j+jj
      col(1,ctr+16) = k+kk
      col(4,ctr+16) = 2
    enddo; enddo; enddo
    row = col
    ele = ele + 1
    k_ele = 0.0
    k_ele(1 :8 ,1 :8 ) = hgmat*(materialpoint_dpdf(1,1,1,1,1,ele) + &
                                materialpoint_dpdf(2,2,2,2,1,ele) + &
                                materialpoint_dpdf(3,3,3,3,1,ele))/3.0_preal
    k_ele(9 :16,9 :16) = hgmat*(materialpoint_dpdf(1,1,1,1,1,ele) + &
                                materialpoint_dpdf(2,2,2,2,1,ele) + &
                                materialpoint_dpdf(3,3,3,3,1,ele))/3.0_preal
    k_ele(17:24,17:24) = hgmat*(materialpoint_dpdf(1,1,1,1,1,ele) + &
                                materialpoint_dpdf(2,2,2,2,1,ele) + &
                                materialpoint_dpdf(3,3,3,3,1,ele))/3.0_preal
    k_ele = k_ele + &
            matmul(transpose(bmatfull), &
                   matmul(reshape(reshape(materialpoint_dpdf(1:3,1:3,1:3,1:3,1,ele), &
                                          shape=[3,3,3,3], order=[2,1,4,3]),shape=[9,9]),bmatfull))*detj
    call matsetvaluesstencil(jac,24,row,24,col,k_ele,add_values,ierr)
    if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  enddo; enddo; enddo
  call matassemblybegin(jac,mat_final_assembly,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call matassemblyend(jac,mat_final_assembly,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call matassemblybegin(jac_pre,mat_final_assembly,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call matassemblyend(jac_pre,mat_final_assembly,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
 
!--------------------------------------------------------------------------------------------------
! applying boundary conditions
  diag = (c_volavg(1,1,1,1)/delta(1)**2.0_preal + &
          c_volavg(2,2,2,2)/delta(2)**2.0_preal + &
          c_volavg(3,3,3,3)/delta(3)**2.0_preal)*detj
  call matzerorowscolumns(jac,size(rows),rows,diag,petsc_null_vec,petsc_null_vec,ierr)
  if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call dmgetglobalvector(da_local,coordinates,ierr);if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call dmdavecgetarrayf90(da_local,coordinates,x_scal,ierr);if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  ele = 0
  do k = zstart, zend; do j = ystart, yend; do i = xstart, xend
    ele = ele + 1
    x_scal(0:2,i,j,k) = discretization_ipcoords(1:3,ele)
  enddo; enddo; enddo
  call dmdavecrestorearrayf90(da_local,coordinates,x_scal,ierr);if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif                       ! initialize to undeformed coordinates (ToDo: use ip coordinates)
  call matnullspacecreaterigidbody(coordinates,matnull,ierr);if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif                          ! get rigid body deformation modes
  call dmrestoreglobalvector(da_local,coordinates,ierr);if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call matsetnullspace(jac,matnull,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call matsetnearnullspace(jac,matnull,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call matnullspacedestroy(matnull,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
 
end subroutine formjacobian
 
end module grid_mech_fem