discretization_grid.f90

Go to the documentation of this file.

# 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/grid/discretization_grid.f90"
# 1 "<built-in>"
# 1 "<command-line>"
# 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/grid/discretization_grid.f90"
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
module discretization_grid
 
# 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/discretization_grid.f90" 2
  use petscsys
 
  use prec
  use system_routines
  use damask_interface
  use io
  use debug
  use numerics
  use results
  use discretization
  use geometry_plastic_nonlocal
  use fesolving
 
  implicit none
  private
 
  integer,     dimension(3), public, protected :: &
    grid
  integer,                   public, protected :: &
    grid3, &                                                                                        !< (local) grid in 3rd direction
    grid3offset
  real(preal), dimension(3), public, protected :: &
    geomsize
  real(preal),               public, protected :: &
    size3, &                                                                                        !< (local) size in 3rd direction
    size3offset
 
  public :: &
    discretization_grid_init
 
contains
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine discretization_grid_init
 
  include 'fftw3-mpi.f03'
  real(preal), dimension(3) :: &
    mysize, &                                                                                       !< domain size of this process
    origin
  integer,     dimension(3) :: &
    mygrid
 
  integer,     dimension(:),   allocatable :: &
    microstructureat, &
    homogenizationat
 
  integer :: j
  integer(C_INTPTR_T) :: &
    devnull, z, z_offset
 
  write(6,'(/,a)') ' <<<+-  discretization_grid init  -+>>>'; flush(6)
 
  call readgeom(grid,geomsize,origin,microstructureat,homogenizationat)
 
!--------------------------------------------------------------------------------------------------
! grid solver specific quantities
  if(worldsize>grid(3)) call io_error(894, ext_msg='number of processes exceeds grid(3)')
 
  call fftw_mpi_init
  devnull = fftw_mpi_local_size_3d(int(grid(3),c_intptr_t), &
                                   int(grid(2),c_intptr_t), &
                                   int(grid(1),c_intptr_t)/2+1, &
                                   petsc_comm_world, &
                                   z, &                                                             ! domain grid size along z
                                   z_offset)                                                        ! domain grid offset along z
  grid3       = int(z)
  grid3offset = int(z_offset)
  size3       = geomsize(3)*real(grid3,preal)      /real(grid(3),preal)
  size3offset = geomsize(3)*real(grid3offset,preal)/real(grid(3),preal)
  mygrid = [grid(1:2),grid3]
  mysize = [geomsize(1:2),size3]
 
!--------------------------------------------------------------------------------------------------
! general discretization
  microstructureat = microstructureat(product(grid(1:2))*grid3offset+1: &
                                      product(grid(1:2))*(grid3offset+grid3))                       ! reallocate/shrink in case of MPI
  homogenizationat = homogenizationat(product(grid(1:2))*grid3offset+1: &
                                      product(grid(1:2))*(grid3offset+grid3))                       ! reallocate/shrink in case of MPI
 
  call discretization_init(homogenizationat,microstructureat, &
                           ipcoordinates0(mygrid,mysize,grid3offset), &
                           nodes0(mygrid,mysize,grid3offset),&
                           merge((grid(1)+1) * (grid(2)+1) * (grid3+1),&                            ! write bottom layer
                                 (grid(1)+1) * (grid(2)+1) *  grid3,&                               ! do not write bottom layer (is top of rank-1)
                                 worldrank<1))
 
  fesolving_execelem = [1,product(mygrid)]                                                          ! parallel loop bounds set to comprise all elements
  fesolving_execip   = [1,1]                                                                        ! parallel loop bounds set to comprise the only IP
 
!--------------------------------------------------------------------------------------------------
! store geometry information for post processing
  call results_openjobfile
  call results_closegroup(results_addgroup('geometry'))
  call results_addattribute('grid',  grid,    'geometry')
  call results_addattribute('size',  geomsize,'geometry')
  call results_addattribute('origin',origin,  'geometry')
  call results_closejobfile
 
!--------------------------------------------------------------------------------------------------
! geometry information required by the nonlocal CP model
  call geometry_plastic_nonlocal_setipvolume(reshape([(product(mysize/real(mygrid,preal)),j=1,product(mygrid))], &
                                                     [1,product(mygrid)]))
  call geometry_plastic_nonlocal_setiparea        (cellsurfacearea(mysize,mygrid))
  call geometry_plastic_nonlocal_setipareanormal  (cellsurfacenormal(product(mygrid)))
  call geometry_plastic_nonlocal_setipneighborhood(ipneighborhood(mygrid))
 
!--------------------------------------------------------------------------------------------------
! sanity checks for debugging
  if (debug_e < 1 .or. debug_e > product(mygrid)) call io_error(602,ext_msg='element')              ! selected element does not exist
  if (debug_i /= 1)                               call io_error(602,ext_msg='IP')                   ! selected IP does not exist
 
end subroutine discretization_grid_init
 
 
!--------------------------------------------------------------------------------------------------
! supposed to be called only once!
!--------------------------------------------------------------------------------------------------
subroutine readgeom(grid,geomSize,origin,microstructure,homogenization)
 
  integer,     dimension(3), intent(out) :: &
    grid                                                                                            ! grid (for all processes!)
  real(pReal), dimension(3), intent(out) :: &
    geomSize, &                                                                                     ! size (for all processes!)
    origin                                                                                          ! origin (for all processes!)
  integer,     dimension(:), intent(out), allocatable :: &
    microstructure, &
    homogenization
 
  character(len=:),      allocatable :: rawData
  character(len=65536)               :: line
  integer, allocatable, dimension(:) :: chunkPos
  integer :: &
    h =- 1, &
    headerlength = -1, &                                                                            
    filelength, &                                                                                   
    fileunit, &
    startpos, endpos, &
    mystat, &
    l, &                                                                                            
    c, &                                                                                            
    o, &                                                                                            
    e, &                                                                                            
    i, j
 
  grid = -1
  geomsize = -1.0_preal
 
!--------------------------------------------------------------------------------------------------
! read raw data as stream
  inquire(file = trim(geometryfile), size=filelength)
  open(newunit=fileunit, file=trim(geometryfile), access='stream',&
       status='old', position='rewind', action='read',iostat=mystat)
  if(mystat /= 0) call io_error(100,ext_msg=trim(geometryfile))
  allocate(character(len=fileLength)::rawData)
  read(fileunit) rawdata
  close(fileunit)
 
!--------------------------------------------------------------------------------------------------
! get header length
  endpos = index(rawdata,io_eol)
  if(endpos <= index(rawdata,'head')) then                                                          ! ToDo: Should be 'header'
    startpos = len(rawdata)
    call io_error(error_id=841, ext_msg='readGeom')
  else
    chunkpos = io_stringpos(rawdata(1:endpos))
    if (chunkpos(1) < 2) call io_error(error_id=841, ext_msg='readGeom')
    headerlength = io_intvalue(rawdata(1:endpos),chunkpos,1)
    startpos = endpos + 1
  endif
 
!--------------------------------------------------------------------------------------------------
! read and interprete header
  origin = 0.0_preal
  l = 0
  do while (l < headerlength .and. startpos < len(rawdata))
    endpos = startpos + index(rawdata(startpos:),io_eol) - 1
    if (endpos < startpos) endpos = len(rawdata)                                                    ! end of file without new line
    line = rawdata(startpos:endpos)
    startpos = endpos + 1
    l = l + 1
 
    chunkpos = io_stringpos(trim(line))
    if (chunkpos(1) < 2) cycle                                                                      ! need at least one keyword value pair
 
    select case (io_lc(io_stringvalue(trim(line),chunkpos,1)) )
      case ('grid')
        if (chunkpos(1) > 6) then
          do j = 2,6,2
            select case (io_lc(io_stringvalue(line,chunkpos,j)))
              case('a')
                grid(1) = io_intvalue(line,chunkpos,j+1)
              case('b')
                grid(2) = io_intvalue(line,chunkpos,j+1)
              case('c')
                grid(3) = io_intvalue(line,chunkpos,j+1)
            end select
          enddo
        endif
 
      case ('size')
        if (chunkpos(1) > 6) then
          do j = 2,6,2
            select case (io_lc(io_stringvalue(line,chunkpos,j)))
              case('x')
                geomsize(1) = io_floatvalue(line,chunkpos,j+1)
              case('y')
                geomsize(2) = io_floatvalue(line,chunkpos,j+1)
              case('z')
                geomsize(3) = io_floatvalue(line,chunkpos,j+1)
            end select
          enddo
        endif
 
      case ('origin')
        if (chunkpos(1) > 6) then
          do j = 2,6,2
            select case (io_lc(io_stringvalue(line,chunkpos,j)))
              case('x')
                origin(1) = io_floatvalue(line,chunkpos,j+1)
              case('y')
                origin(2) = io_floatvalue(line,chunkpos,j+1)
              case('z')
                origin(3) = io_floatvalue(line,chunkpos,j+1)
            end select
          enddo
        endif
 
      case ('homogenization')
        if (chunkpos(1) > 1) h = io_intvalue(line,chunkpos,2)
 
    end select
 
  enddo
 
!--------------------------------------------------------------------------------------------------
! sanity checks
  if(h < 1) &
    call io_error(error_id = 842, ext_msg='homogenization (readGeom)')
  if(any(grid < 1)) &
    call io_error(error_id = 842, ext_msg='grid (readGeom)')
  if(any(geomsize < 0.0_preal)) &
    call io_error(error_id = 842, ext_msg='size (readGeom)')
 
  allocate(microstructure(product(grid)), source = -1)                                              ! too large in case of MPI (shrink later, not very elegant)
  allocate(homogenization(product(grid)), source = h)                                               ! too large in case of MPI (shrink later, not very elegant)
 
!--------------------------------------------------------------------------------------------------
! read and interpret content
  e = 1
  do while (startpos < len(rawdata))
    endpos = startpos + index(rawdata(startpos:),io_eol) - 1
    if (endpos < startpos) endpos = len(rawdata)                                                    ! end of file without new line
    line = rawdata(startpos:endpos)
    startpos = endpos + 1
    l = l + 1
    chunkpos = io_stringpos(trim(line))
 
    nocompression: if (chunkpos(1) /= 3) then
      c = chunkpos(1)
      microstructure(e:e+c-1) =  [(io_intvalue(line,chunkpos,i+1), i=0, c-1)]
    else nocompression
      compression: if (io_lc(io_stringvalue(line,chunkpos,2))  == 'of') then
        c = io_intvalue(line,chunkpos,1)
        microstructure(e:e+c-1) = [(io_intvalue(line,chunkpos,3),i = 1,io_intvalue(line,chunkpos,1))]
      else         if (io_lc(io_stringvalue(line,chunkpos,2))  == 'to') then compression
        c = abs(io_intvalue(line,chunkpos,3) - io_intvalue(line,chunkpos,1)) + 1
        o = merge(+1, -1, io_intvalue(line,chunkpos,3) > io_intvalue(line,chunkpos,1))
        microstructure(e:e+c-1) = [(i, i = io_intvalue(line,chunkpos,1),io_intvalue(line,chunkpos,3),o)]
      else compression
        c = chunkpos(1)
        microstructure(e:e+c-1) = [(io_intvalue(line,chunkpos,i+1), i=0, c-1)]
      endif compression
    endif nocompression
 
    e = e+c
  end do
 
  if (e-1 /= product(grid)) call io_error(error_id = 843, el=e)
 
end subroutine readgeom
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
function ipcoordinates0(grid,geomSize,grid3Offset)
 
  integer,     dimension(3), intent(in) :: grid                                                     ! grid (for this process!)
  real(preal), dimension(3), intent(in) :: geomsize                                                 ! size (for this process!)
  integer,                   intent(in) :: grid3offset                                              ! grid(3) offset
 
  real(preal), dimension(3,product(grid))  :: ipcoordinates0
 
  integer :: &
    a,b,c, &
    i
 
  i = 0
  do c = 1, grid(3); do b = 1, grid(2); do a = 1, grid(1)
    i = i + 1
    ipcoordinates0(1:3,i) = geomsize/real(grid,preal) * (real([a,b,grid3offset+c],preal) -0.5_preal)
  enddo; enddo; enddo
 
end function ipcoordinates0
 
 
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
pure function nodes0(grid,geomSize,grid3Offset)
 
  integer,     dimension(3), intent(in) :: grid                                                     ! grid (for this process!)
  real(preal), dimension(3), intent(in) :: geomsize                                                 ! size (for this process!)
  integer,                   intent(in) :: grid3offset                                              ! grid(3) offset
 
  real(preal), dimension(3,product(grid+1)) :: nodes0
 
  integer :: &
    a,b,c, &
    n
 
  n = 0
  do c = 0, grid3; do b = 0, grid(2); do a = 0, grid(1)
    n = n + 1
    nodes0(1:3,n) = geomsize/real(grid,preal) * real([a,b,grid3offset+c],preal)
  enddo; enddo; enddo
 
end function nodes0
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
pure function cellsurfacearea(geomSize,grid)
 
  real(preal), dimension(3), intent(in) :: geomsize                                                 ! size (for this process!)
  integer,     dimension(3), intent(in) :: grid                                                     ! grid (for this process!)
 
  real(preal), dimension(6,1,product(grid)) :: cellsurfacearea
 
  cellsurfacearea(1:2,1,:) = geomsize(2)/real(grid(2)) * geomsize(3)/real(grid(3))
  cellsurfacearea(3:4,1,:) = geomsize(3)/real(grid(3)) * geomsize(1)/real(grid(1))
  cellsurfacearea(5:6,1,:) = geomsize(1)/real(grid(1)) * geomsize(2)/real(grid(2))
 
end function cellsurfacearea
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
pure function cellsurfacenormal(nElems)
 
  integer, intent(in) :: nelems
 
  real(preal), dimension(3,6,1,nElems) :: cellsurfacenormal
 
  cellsurfacenormal(1:3,1,1,:) = spread([+1.0_preal, 0.0_preal, 0.0_preal],2,nelems)
  cellsurfacenormal(1:3,2,1,:) = spread([-1.0_preal, 0.0_preal, 0.0_preal],2,nelems)
  cellsurfacenormal(1:3,3,1,:) = spread([ 0.0_preal,+1.0_preal, 0.0_preal],2,nelems)
  cellsurfacenormal(1:3,4,1,:) = spread([ 0.0_preal,-1.0_preal, 0.0_preal],2,nelems)
  cellsurfacenormal(1:3,5,1,:) = spread([ 0.0_preal, 0.0_preal,+1.0_preal],2,nelems)
  cellsurfacenormal(1:3,6,1,:) = spread([ 0.0_preal, 0.0_preal,-1.0_preal],2,nelems)
 
end function cellsurfacenormal
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
pure function ipneighborhood(grid)
 
  integer, dimension(3), intent(in) :: grid                                                         ! grid (for this process!)
 
  integer, dimension(3,6,1,product(grid)) :: ipneighborhood
 
  integer :: &
   x,y,z, &
   e
 
  e = 0
  do z = 0,grid(3)-1; do y = 0,grid(2)-1; do x = 0,grid(1)-1
    e = e + 1
    ! element ID
    ipneighborhood(1,1,1,e) = z * grid(1) * grid(2) &
                            + y * grid(1) &
                            + modulo(x+1,grid(1)) &
                            + 1
    ipneighborhood(1,2,1,e) = z * grid(1) * grid(2) &
                            + y * grid(1) &
                            + modulo(x-1,grid(1)) &
                            + 1
    ipneighborhood(1,3,1,e) = z * grid(1) * grid(2) &
                            + modulo(y+1,grid(2)) * grid(1) &
                            + x &
                            + 1
    ipneighborhood(1,4,1,e) = z * grid(1) * grid(2) &
                            + modulo(y-1,grid(2)) * grid(1) &
                            + x &
                            + 1
    ipneighborhood(1,5,1,e) = modulo(z+1,grid(3)) * grid(1) * grid(2) &
                            + y * grid(1) &
                            + x &
                            + 1
    ipneighborhood(1,6,1,e) = modulo(z-1,grid(3)) * grid(1) * grid(2) &
                            + y * grid(1) &
                            + x &
                            + 1
    ! IP ID
    ipneighborhood(2,:,1,e) = 1
 
    ! face ID
    ipneighborhood(3,1,1,e) = 2
    ipneighborhood(3,2,1,e) = 1
    ipneighborhood(3,3,1,e) = 4
    ipneighborhood(3,4,1,e) = 3
    ipneighborhood(3,5,1,e) = 6
    ipneighborhood(3,6,1,e) = 5
 
  enddo; enddo; enddo
 
end function ipneighborhood
 
 
end module discretization_grid