mesh.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2011 M. Marques, A. Castro, A. Rubio, G. Bertsch, M. Oliveira
!! Copyright (C) 2021 S. Ohlmann
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
 
module mesh_oct_m
  use basis_set_abst_oct_m
  use box_oct_m
  use comm_oct_m
  use coordinate_system_oct_m
  use debug_oct_m
  use global_oct_m
  use iihash_oct_m
  use index_oct_m
  use io_oct_m
  use io_binary_oct_m
  use messages_oct_m
  use mpi_oct_m
  use mpi_lib_oct_m
  use namespace_oct_m
  use par_vec_oct_m
  use partition_oct_m
  use parser_oct_m
  use profiling_oct_m
  use space_oct_m
  use symmetries_oct_m
  use symm_op_oct_m
  use unit_oct_m
  use unit_system_oct_m
 
  implicit none
 
  private
  public ::                         &
    mesh_t,                         &
    mesh_plane_t,                   &
    mesh_line_t,                    &
    mesh_check_dump_compatibility,  &
    mesh_end,                       &
    mesh_double_box,                &
    mesh_r,                         &
    mesh_gcutoff,                   &
    mesh_write_info,                &
    mesh_nearest_point,             &
    mesh_discretize_values_to_mesh, &
    mesh_periodic_point,            &
    mesh_periodic_point_global,     &
    mesh_global_memory,             &
    mesh_local_memory,              &
    mesh_x_global,                  &
    mesh_write_fingerprint,         &
    mesh_read_fingerprint,          &
    mesh_check_symmetries,          &
    mesh_global_index_to_coords,    &
    mesh_global_index_from_coords,  &
    mesh_local_index_to_coords,     &
    mesh_local_index_from_coords,   &
    mesh_local2global,              &
    mesh_global2local,              &
    mesh_minmaxloc
 
  type, extends(basis_set_abst_t) :: mesh_t
    ! Components are public by default
    class(box_t),               pointer :: box
    class(coordinate_system_t), pointer :: coord_system
    type(index_t)                       :: idx
    logical :: use_curvilinear
 
    real(real64), allocatable :: spacing(:)
 
    ! When running serially, the local number of points is
    ! equal to the global number of points.
    ! Otherwise, the next two are different on each node.
    integer      :: np
    integer      :: np_part
    integer(int64)  :: np_global
    integer(int64)  :: np_part_global
    logical         :: parallel_in_domains
    type(mpi_grp_t) :: mpi_grp
    type(par_vec_t) :: pv
    type(partition_t) :: partition
 
    real(real64),   allocatable :: x(:,:)
    real(real64),   allocatable :: x_t(:,:)
    real(real64),   allocatable :: chi(:,:)
    real(real64)                :: volume_element
    real(real64),   allocatable :: vol_pp(:)
    real(real64),   allocatable :: jacobian_inverse(:,:,:)
 
    logical            :: masked_periodic_boundaries
    character(len=256) :: periodic_boundary_mask
  contains
    procedure :: end => mesh_end                
    procedure :: init => mesh_init              
    procedure :: write_info => mesh_write_info  
    procedure :: dmesh_allreduce_0, zmesh_allreduce_0, imesh_allreduce_0, lmesh_allreduce_0
    procedure :: dmesh_allreduce_1, zmesh_allreduce_1, imesh_allreduce_1, lmesh_allreduce_1
    procedure :: dmesh_allreduce_2, zmesh_allreduce_2, imesh_allreduce_2, lmesh_allreduce_2
    procedure :: dmesh_allreduce_3, zmesh_allreduce_3, imesh_allreduce_3, lmesh_allreduce_3
    procedure :: dmesh_allreduce_4, zmesh_allreduce_4, imesh_allreduce_4, lmesh_allreduce_4
    procedure :: dmesh_allreduce_5, zmesh_allreduce_5, imesh_allreduce_5, lmesh_allreduce_5
    procedure :: dmesh_allreduce_6, zmesh_allreduce_6, imesh_allreduce_6, lmesh_allreduce_6
    generic :: allreduce => dmesh_allreduce_0, zmesh_allreduce_0, imesh_allreduce_0, lmesh_allreduce_0
    generic :: allreduce => dmesh_allreduce_1, zmesh_allreduce_1, imesh_allreduce_1, lmesh_allreduce_1
    generic :: allreduce => dmesh_allreduce_2, zmesh_allreduce_2, imesh_allreduce_2, lmesh_allreduce_2
    generic :: allreduce => dmesh_allreduce_3, zmesh_allreduce_3, imesh_allreduce_3, lmesh_allreduce_3
    generic :: allreduce => dmesh_allreduce_4, zmesh_allreduce_4, imesh_allreduce_4, lmesh_allreduce_4
    generic :: allreduce => dmesh_allreduce_5, zmesh_allreduce_5, imesh_allreduce_5, lmesh_allreduce_5
    generic :: allreduce => dmesh_allreduce_6, zmesh_allreduce_6, imesh_allreduce_6, lmesh_allreduce_6
    procedure :: red_min => mesh_red_min      
  end type mesh_t
 
  type mesh_plane_t
    ! Components are public by default
    real(real64) :: n(3)
    real(real64) :: u(3), v(3)
    real(real64) :: origin(3)
    real(real64) :: spacing
    integer :: nu, mu, nv, mv
  end type mesh_plane_t
 
  type mesh_line_t
    ! Components are public by default
    real(real64) :: n(2)
    real(real64) :: u(2)
    real(real64) :: origin(2)
    real(real64) :: spacing
    integer :: nu, mu
  end type mesh_line_t
 
contains
 
  subroutine mesh_init(this)
    class(mesh_t), intent(inout) :: this
 
    push_sub(mesh_init)
 
    call this%set_time_dependent(.false.)
 
    pop_sub(mesh_init)
  end subroutine mesh_init
 
! ---------------------------------------------------------
  subroutine mesh_double_box(space, mesh, alpha, db)
    class(space_t),    intent(in)  :: space
    type(mesh_t),      intent(in)  :: mesh
    real(real64),      intent(in)  :: alpha
    integer,           intent(out) :: db(:)
 
    integer :: idir
 
    push_sub(mesh_double_box)
 
    db = 1
 
    ! double mesh with 2n points
    do idir = 1, space%periodic_dim
      db(idir) = mesh%idx%ll(idir)
    end do
    do idir = space%periodic_dim + 1, space%dim
      db(idir) = nint(alpha * (mesh%idx%ll(idir) - 1)) + 1
    end do
 
    pop_sub(mesh_double_box)
  end subroutine mesh_double_box
 
 
  ! ---------------------------------------------------------
  subroutine mesh_write_info(this, iunit, namespace)
    class(mesh_t),               intent(in) :: this
    integer,           optional, intent(in) :: iunit
    type(namespace_t), optional, intent(in) :: namespace
 
    integer :: ii
    real(real64) :: cutoff
 
    push_sub(mesh_write_info)
 
    write(message(1),'(3a)') '  Spacing [', trim(units_abbrev(units_out%length)), '] = ('
    do ii = 1, this%box%dim
      if (ii > 1) write(message(1), '(2a)') trim(message(1)), ','
      write(message(1), '(a,f6.3)') trim(message(1)), units_from_atomic(units_out%length, this%spacing(ii))
    end do
    write(message(1), '(5a,f12.5)') trim(message(1)), ') ', &
      '   volume/point [', trim(units_abbrev(units_out%length**this%box%dim)), '] = ',      &
      units_from_atomic(units_out%length**this%box%dim, this%vol_pp(1))
 
    write(message(2),'(a, i10)') '  # inner mesh = ', this%np_global
    write(message(3),'(a, i10)') '  # total mesh = ', this%np_part_global
 
    cutoff = mesh_gcutoff(this)**2 / m_two
    write(message(4),'(3a,f12.6,a,f12.6)') '  Grid Cutoff [', trim(units_abbrev(units_out%energy)),'] = ', &
      units_from_atomic(units_out%energy, cutoff), '    Grid Cutoff [Ry] = ', cutoff * m_two
    call messages_info(4, iunit=iunit, namespace=namespace)
 
    pop_sub(mesh_write_info)
  end subroutine mesh_write_info
 
 
  pure subroutine mesh_r(mesh, ip, rr, origin, coords)
    class(mesh_t), intent(in)  :: mesh
    integer,       intent(in)  :: ip
    real(real64),  intent(out) :: rr
    real(real64),  intent(in),  optional :: origin(:)
    real(real64),  intent(out), optional :: coords(:)
 
    real(real64) :: xx(1:mesh%box%dim)
 
    xx = mesh%x(:, ip)
    if (present(origin)) xx = xx - origin
    rr = norm2(xx)
 
    if (present(coords)) then
      coords = xx
    end if
 
  end subroutine mesh_r
 
  subroutine mesh_r_global(mesh, ipg, rr, origin, coords)
    class(mesh_t), intent(in)  :: mesh
    integer(int64),intent(in)  :: ipg
    real(real64),  intent(out) :: rr
    real(real64),  intent(in),  optional :: origin(:)
    real(real64),  intent(out), optional :: coords(:)
 
    real(real64) :: xx(1:mesh%box%dim)
 
    xx = mesh_x_global(mesh, ipg)
    if (present(origin)) xx = xx - origin
    rr = norm2(xx)
 
    if (present(coords)) then
      coords = xx
    end if
 
  end subroutine mesh_r_global
 
 
  integer function mesh_nearest_point(mesh, pos, dmin, rankmin) result(ind)
    class(mesh_t),intent(in)  :: mesh
    real(real64), intent(in)  :: pos(:)
    real(real64), intent(out) :: dmin
    integer,      intent(out) :: rankmin
    !                                          This returns 0 in the absence of domain decomposition
 
    real(real64) :: dd
    integer :: imin, ip
 
    push_sub(mesh_nearest_point)
 
    ! find the point of the grid that is closer to the atom
    dmin = m_zero
    do ip = 1, mesh%np
      dd = sum((pos - mesh%x(:, ip))**2)
      if ((dd < dmin) .or. (ip == 1)) then
        imin = ip
        dmin = dd
      end if
    end do
 
    call mesh_minmaxloc(mesh, dmin, rankmin, mpi_minloc)
    call mesh%mpi_grp%bcast(imin, 1, mpi_integer, rankmin)
 
    ind = imin
    pop_sub(mesh_nearest_point)
  end function mesh_nearest_point
 
  subroutine mesh_discretize_values_to_mesh(mesh, values)
    class(mesh_t),   intent(in  )  :: mesh
    real(real64),    intent(inout) :: values(:, :)
    !                                                   Out: Values discretized to mesh points
    integer        :: i, ip, ndim
    integer        :: process
    real(real64)   :: dummy
 
    push_sub(mesh_discretize_values_to_mesh)
 
    if (mesh%parallel_in_domains) then
      ndim = size(values, 1)
      do i = 1, size(values, 2)
        ip = mesh_nearest_point(mesh, values(:, i), dummy, process)
        ! ip is defined for `process`, only
        if (mesh%mpi_grp%rank == process) values(:, i) = mesh%x(:, ip)
        call mesh%mpi_grp%bcast(values(:, i), ndim, mpi_double_precision, process)
      enddo
    else
      do i = 1, size(values, 2)
        ip = mesh_nearest_point(mesh, values(:, i), dummy, process)
        values(:, i) = mesh%x(:, ip)
      enddo
    endif
 
    pop_sub(mesh_discretize_values_to_mesh)
 
  end subroutine mesh_discretize_values_to_mesh
 
  ! --------------------------------------------------------------
  ! --------------------------------------------------------------
  real(real64) function mesh_gcutoff(mesh) result(gmax)
    class(mesh_t), intent(in) :: mesh
 
    push_sub(mesh_gcutoff)
    gmax = m_pi / (maxval(mesh%spacing))
 
    pop_sub(mesh_gcutoff)
  end function mesh_gcutoff
 
  ! --------------------------------------------------------------
  subroutine mesh_write_fingerprint(mesh, dir, filename, mpi_grp, namespace, ierr)
    type(mesh_t),     intent(in)  :: mesh
    character(len=*), intent(in)  :: dir
    character(len=*), intent(in)  :: filename
    type(mpi_grp_t),  intent(in)  :: mpi_grp
    type(namespace_t),intent(in)  :: namespace
    integer,          intent(out) :: ierr
 
    integer :: iunit, ii
 
    push_sub(mesh_write_fingerprint)
 
    ierr = 0
 
    iunit = io_open(trim(dir)//"/"//trim(filename), namespace, action='write', &
      die=.false., grp=mpi_grp)
    if (iunit == -1) then
      message(1) = "Unable to open file '"//trim(dir)//"/"//trim(filename)//"'."
      call messages_warning(1, namespace=namespace)
      ierr = ierr + 1
    else
      if (mpi_grp%is_root()) then
        write(iunit, '(a20,i21)')  'np_part_global=     ', mesh%np_part_global
        write(iunit, '(a20,i21)')  'np_global=          ', mesh%np_global
        write(iunit, '(a20,i21)')  'algorithm=          ', 1
        write(iunit, '(a20,i21)')  'checksum=           ', mesh%idx%checksum
        write(iunit, '(a20,i21)')  'bits=               ', mesh%idx%bits
        write(iunit, '(a20,i21)')  'dim=                ', mesh%idx%dim
        write(iunit, '(a20,i21)')  'type=               ', mesh%idx%type
        do ii = 1, mesh%idx%dim
          write(iunit, '(a7,i2,a11,i21)') 'offset(',ii,')=         ', mesh%idx%offset(ii)
        end do
        do ii = 1, mesh%idx%dim
          write(iunit, '(a7,i2,a11,i21)') 'nn(',ii,')=             ', mesh%idx%nr(2, ii) - mesh%idx%nr(1, ii) + 1
        end do
      end if
      call io_close(iunit, grp=mpi_grp)
    end if
 
    pop_sub(mesh_write_fingerprint)
  end subroutine mesh_write_fingerprint
 
 
  ! -----------------------------------------------------------------------
  subroutine mesh_read_fingerprint(mesh, dir, filename, mpi_grp, namespace, &
    read_np_part, read_np, bits, type, offset, nn, ierr)
    type(mesh_t),     intent(in)  :: mesh
    character(len=*), intent(in)  :: dir
    character(len=*), intent(in)  :: filename
    type(mpi_grp_t),  intent(in)  :: mpi_grp
    type(namespace_t),intent(in)  :: namespace
    integer(int64),      intent(out) :: read_np_part
    integer(int64),      intent(out) :: read_np
    integer,          intent(out) :: bits
    integer,          intent(out) :: type
    integer,          intent(out) :: offset(1:mesh%idx%dim)
    integer,          intent(out) :: nn(1:mesh%idx%dim)
    integer,          intent(out) :: ierr
 
    character(len=20)  :: str
    character(len=100) :: lines(7)
    integer :: iunit, algorithm, dim, err, ii
    integer(int64) :: checksum
 
    push_sub(mesh_read_fingerprint)
 
    ierr = 0
 
    read_np_part = 0_int64
    read_np = 0_int64
 
    iunit = io_open(trim(dir)//"/"//trim(filename), namespace, action='read', &
      status='old', die=.false., grp=mpi_grp)
    if (iunit == -1) then
      ierr = ierr + 1
      message(1) = "Unable to open file '"//trim(dir)//"/"//trim(filename)//"'."
      call messages_warning(1, namespace=namespace)
    else
      call iopar_read(mpi_grp, iunit, lines, 7, err)
      if (err /= 0) then
        ierr = ierr + 4
      else
        read(lines(1), '(a20,i21)')  str, read_np_part
        read(lines(2), '(a20,i21)')  str, read_np
        read(lines(3), '(a20,i21)')  str, algorithm
        read(lines(4), '(a20,i21)')  str, checksum
        read(lines(5), '(a20,i21)')  str, bits
        read(lines(6), '(a20,i21)')  str, dim
        read(lines(7), '(a20,i21)')  str, type
        ! only allow restarting simulations with the same dimensions
        if (dim /= mesh%idx%dim) then
          ierr = ierr + 8
        else
          ! read offset, has dim lines
          call iopar_read(mpi_grp, iunit, lines, dim, err)
          if (err /= 0) then
            ierr = ierr + 4
          else
            do ii = 1, dim
              read(lines(ii), '(a20,i21)')  str, offset(ii)
            end do
          end if
 
          ! read nn, has dim lines
          call iopar_read(mpi_grp, iunit, lines, dim, err)
          if (err /= 0) then
            ierr = ierr + 4
          else
            do ii = 1, dim
              read(lines(ii), '(a20,i21)')  str, nn(ii)
            end do
          end if
        end if
 
        assert(read_np_part >= read_np)
 
        if (read_np_part == mesh%np_part_global &
          .and. read_np == mesh%np_global &
          .and. algorithm == 1 &
          .and. checksum == mesh%idx%checksum) then
          read_np_part = 0
          read_np = 0
        end if
      end if
 
      call io_close(iunit, grp=mpi_grp)
    end if
 
    pop_sub(mesh_read_fingerprint)
  end subroutine mesh_read_fingerprint
 
  ! ---------------------------------------------------------
  subroutine mesh_check_dump_compatibility(mesh, dir, filename, namespace, mpi_grp, grid_changed, grid_reordered, map, ierr)
    type(mesh_t),             intent(in)  :: mesh
    character(len=*),         intent(in)  :: dir
    character(len=*),         intent(in)  :: filename
    type(namespace_t),        intent(in)  :: namespace
    type(mpi_grp_t),          intent(in)  :: mpi_grp
    logical,                  intent(out) :: grid_changed
    logical,                  intent(out) :: grid_reordered
    integer(int64), allocatable, intent(out) :: map(:)
    integer,                  intent(out) :: ierr
 
    integer(int64) :: ipg, ipg_new, read_np_part, read_np
    integer :: err, idir
    integer :: bits, type, offset(mesh%idx%dim), point(mesh%idx%dim), nn(mesh%idx%dim)
    integer(int64), allocatable :: read_indices(:)
    type(index_t) :: idx_old
 
    push_sub(mesh_check_dump_compatibility)
 
    ierr = 0
 
    grid_changed = .false.
    grid_reordered = .false.
 
    ! Read the mesh fingerprint
    call mesh_read_fingerprint(mesh, dir, filename, mpi_grp, namespace, read_np_part, read_np, &
      bits, type, offset, nn, err)
    if (err /= 0) then
      ierr = 1
      message(1) = "Unable to read mesh fingerprint from '"//trim(dir)//"/"//trim(filename)//"'."
      call messages_warning(1, namespace=namespace)
 
    else if (read_np > 0) then
      if (.not. associated(mesh%box)) then
        ! We can only check the compatibility of two meshes that have different fingerprints if we also
        ! have the simulation box. In the case we do not, we will assume that the fingerprint is enough.
        ierr = ierr + 2
      else
        grid_changed = .true.
 
        ! perhaps only the order of the points changed, this can only
        ! happen if the number of points is the same and no points maps
        ! to zero (this is checked below)
        grid_reordered = (read_np == mesh%np_global)
 
        ! the grid is different, so we read the coordinates.
        safe_allocate(read_indices(1:read_np_part))
        call io_binary_read(trim(io_workpath(dir, namespace))//"/indices.obf", read_np_part, &
          read_indices, err)
        if (err /= 0) then
          ierr = ierr + 4
          message(1) = "Unable to read index map from '"//trim(dir)//"'."
          call messages_warning(1, namespace=namespace)
        else
          ! dummy index object
          call index_init(idx_old, mesh%idx%dim)
          idx_old%type = type
          idx_old%bits = bits
          idx_old%nr(1, :) = -offset
          idx_old%nr(2, :) = -offset + nn - 1
          idx_old%offset = offset
          idx_old%stride(1) = 1
          do idir = 2, mesh%idx%dim
            idx_old%stride(idir) = idx_old%stride(idir-1) * nn(idir-1)
          end do
          ! generate the map
          safe_allocate(map(1:read_np))
          do ipg = 1, read_np
            ! get nd-index from old 1d index
            call index_spatial_to_point(idx_old, mesh%idx%dim, read_indices(ipg), point)
            ! get new global index
            ipg_new = mesh_global_index_from_coords(mesh, point)
            map(ipg) = ipg_new
            ! ignore boundary points
            if (map(ipg) > mesh%np_global) map(ipg) = 0
            ! if the map is zero for one point, it is not a simple reordering
            if (map(ipg) == 0) grid_reordered = .false.
          end do
          call index_end(idx_old)
        end if
 
        safe_deallocate_a(read_indices)
      end if
    end if
 
    pop_sub(mesh_check_dump_compatibility)
  end subroutine mesh_check_dump_compatibility
 
 
  ! --------------------------------------------------------------
  recursive subroutine mesh_end(this)
    class(mesh_t), intent(inout)   :: this
 
    push_sub(mesh_end)
 
#ifdef HAVE_MPI
    call lmpi_destroy_shared_memory_window(this%idx%window_grid_to_spatial)
    call lmpi_destroy_shared_memory_window(this%idx%window_spatial_to_grid)
    nullify(this%idx%grid_to_spatial_global)
    nullify(this%idx%spatial_to_grid_global)
#else
    safe_deallocate_p(this%idx%grid_to_spatial_global)
    safe_deallocate_p(this%idx%spatial_to_grid_global)
#endif
 
    safe_deallocate_a(this%x)
    safe_deallocate_a(this%chi)
    safe_deallocate_a(this%vol_pp)
    safe_deallocate_a(this%jacobian_inverse)
 
    if (this%parallel_in_domains) then
      call par_vec_end(this%pv)
      call partition_end(this%partition)
    end if
 
    call index_end(this%idx)
    safe_deallocate_a(this%spacing)
 
    pop_sub(mesh_end)
  end subroutine mesh_end
 
 
  ! ---------------------------------------------------------
  integer(int64) function mesh_periodic_point(mesh, space, ip) result(ipg)
    class(mesh_t), intent(in)    :: mesh
    class(space_t),intent(in)    :: space
    integer,       intent(in)    :: ip
 
    integer :: ix(space%dim), nr(2, space%dim), idim
    real(real64) :: xx(space%dim), rr, ufn_re, ufn_im
 
    ! no push_sub, called too frequently
 
    call mesh_local_index_to_coords(mesh, ip, ix)
    nr(1, :) = mesh%idx%nr(1, :) + mesh%idx%enlarge(:)
    nr(2, :) = mesh%idx%nr(2, :) - mesh%idx%enlarge(:)
 
    do idim = 1, space%periodic_dim
      do while (ix(idim) < nr(1, idim))
        ix(idim) = ix(idim) + mesh%idx%ll(idim)
      end do
      do while (ix(idim) > nr(2, idim))
        ix(idim) = ix(idim) - mesh%idx%ll(idim)
      end do
    end do
 
    ipg = mesh_global_index_from_coords(mesh, ix)
    assert(ipg > 0)
 
    if (mesh%masked_periodic_boundaries) then
      call mesh_r(mesh, ip, rr, coords = xx)
      call parse_expression(ufn_re, ufn_im, space%dim, xx, rr, m_zero, mesh%periodic_boundary_mask)
      if (int(ufn_re) == 0) ipg = mesh_local2global(mesh, ip) ! Nothing will be done
    end if
 
  end function mesh_periodic_point
 
  integer(int64) function mesh_periodic_point_global(mesh, space, ipg_in) result(ipg)
    class(mesh_t), intent(in)    :: mesh
    class(space_t),intent(in)    :: space
    integer(int64),intent(in)    :: ipg_in
 
    integer :: ix(space%dim), nr(2, space%dim), idim
    real(real64) :: xx(space%dim), rr, ufn_re, ufn_im
 
    ! no push_sub, called too frequently
 
    call mesh_global_index_to_coords(mesh, ipg_in, ix)
    nr(1, :) = mesh%idx%nr(1, :) + mesh%idx%enlarge(:)
    nr(2, :) = mesh%idx%nr(2, :) - mesh%idx%enlarge(:)
 
    do idim = 1, space%periodic_dim
      if (ix(idim) < nr(1, idim)) ix(idim) = ix(idim) + mesh%idx%ll(idim)
      if (ix(idim) > nr(2, idim)) ix(idim) = ix(idim) - mesh%idx%ll(idim)
    end do
 
    ipg = mesh_global_index_from_coords(mesh, ix)
    assert(ipg > 0)
 
    if (mesh%masked_periodic_boundaries) then
      call mesh_r_global(mesh, ipg_in, rr, coords = xx)
      call parse_expression(ufn_re, ufn_im, space%dim, xx, rr, m_zero, mesh%periodic_boundary_mask)
      if (int(ufn_re) == 0) ipg = ipg_in ! Nothing will be done
    end if
 
  end function mesh_periodic_point_global
 
 
 
  ! ---------------------------------------------------------
  real(real64) pure function mesh_global_memory(mesh) result(memory)
    use iso_c_binding, only: c_sizeof, c_long_long
    class(mesh_t), intent(in) :: mesh
 
    ! 2 global index arrays
    memory = c_sizeof(c_long_long) * real(mesh%np_part_global, real64)  * 2
 
  end function mesh_global_memory
 
 
  ! ---------------------------------------------------------
  real(real64) pure function mesh_local_memory(mesh) result(memory)
    use iso_c_binding, only: c_sizeof, c_long_long
    class(mesh_t), intent(in) :: mesh
 
    memory = m_zero
 
    ! x
    memory = memory + real64 * real(mesh%np_part, real64)  * mesh%idx%dim
    ! local index arrays
    memory = memory + c_sizeof(c_long_long) * real(mesh%np_part, real64)  * 2
  end function mesh_local_memory
 
 
  function mesh_x_global(mesh, ipg) result(xx)
    class(mesh_t),    intent(in) :: mesh
    integer(int64),   intent(in) :: ipg
    real(real64)                 :: xx(1:mesh%box%dim)
 
    real(real64) :: chi(1:mesh%box%dim)
    integer :: ix(1:mesh%box%dim)
 
    ! no push_sub because function is called too frequently
 
    call mesh_global_index_to_coords(mesh, ipg, ix)
    chi = ix * mesh%spacing
    xx = mesh%coord_system%to_cartesian(chi)
 
  end function mesh_x_global
 
 
  ! ---------------------------------------------------------
  subroutine mesh_check_symmetries(mesh, symm, periodic_dim)
    class(mesh_t),       intent(in) :: mesh
    type(symmetries_t),  intent(in) :: symm
    integer,             intent(in) :: periodic_dim
 
    integer :: iop, ip, idim, nops, ix(1:3)
    real(real64) :: destpoint(1:3), srcpoint(1:3), lsize(1:3), offset(1:3)
    real(real64), parameter :: tol_spacing = 1e-12_real64
 
    !If all the axis have the same spacing and the same length
    !the grid is by obviously symmetric
    !Indeed, reduced coordinates are proportional to the point index
    !and the reduced rotation are integer matrices
    !The result of the product is also proportional to an integer
    !and therefore belong to the grid.
    if (mesh%idx%ll(1) == mesh%idx%ll(2)  .and. &
      mesh%idx%ll(2)   == mesh%idx%ll(3)  .and. &
      abs(mesh%spacing(1) - mesh%spacing(2)) < tol_spacing .and. &
      abs(mesh%spacing(2) - mesh%spacing(3)) < tol_spacing) return
 
    push_sub(mesh_check_symmetries)
 
    message(1) = "Checking if the real-space grid is symmetric"
    call messages_info(1)
 
    lsize(1:3) = real(mesh%idx%ll(1:3), real64)
    offset(1:3) = real(mesh%idx%nr(1, 1:3) + mesh%idx%enlarge(1:3), real64)
 
    nops = symmetries_number(symm)
 
    do ip = 1, mesh%np
      !We use floating point coordinates to check if the symmetric point
      !belong to the grid.
      !If yes, it should have integer reduced coordinates
      call mesh_local_index_to_coords(mesh, ip, ix)
      destpoint(1:3) = real(ix(1:3), real64) - offset(1:3)
      ! offset moves corner of cell to origin, in integer mesh coordinates
      assert(all(destpoint >= 0))
      assert(all(destpoint < lsize))
 
      ! move to center of cell in real coordinates
      destpoint = destpoint - real(int(lsize)/2, real64)
 
      !convert to proper reduced coordinates
      do idim = 1, 3
        destpoint(idim) = destpoint(idim)/lsize(idim)
      end do
 
      ! iterate over all points that go to this point by a symmetry operation
      do iop = 1, nops
        srcpoint = symm_op_apply_red(symm%ops(iop), destpoint)
 
        !We now come back to what should be an integer, if the symmetric point beloings to the grid
        do idim = 1, 3
          srcpoint(idim) = srcpoint(idim)*lsize(idim)
        end do
 
        ! move back to reference to origin at corner of cell
        srcpoint = srcpoint + real(int(lsize)/2, real64)
 
        ! apply periodic boundary conditions in periodic directions
        do idim = 1, periodic_dim
          if (nint(srcpoint(idim)) < 0 .or. nint(srcpoint(idim)) >= lsize(idim)) then
            srcpoint(idim) = modulo(srcpoint(idim)+m_half*symprec, lsize(idim))
          end if
        end do
        assert(all(srcpoint >= -symprec))
        assert(all(srcpoint < lsize))
 
        srcpoint(1:3) = srcpoint(1:3) + offset(1:3)
 
        if (any(srcpoint-anint(srcpoint)> symprec*m_two)) then
          message(1) = "The real-space grid breaks at least one of the symmetries of the system."
          message(2) = "Change your spacing or use SymmetrizeDensity=no."
          call messages_fatal(2)
        end if
      end do
    end do
 
    pop_sub(mesh_check_symmetries)
  end subroutine
 
  integer(int64) function mesh_global_index_from_coords(mesh, ix) result(index)
    class(mesh_t), intent(in)    :: mesh
    integer,       intent(in)    :: ix(:)
 
    index = index_from_coords(mesh%idx, ix)
  end function mesh_global_index_from_coords
 
  subroutine mesh_global_index_to_coords(mesh, ipg, ix)
    type(mesh_t),  intent(in)    :: mesh
    integer(int64),   intent(in)    :: ipg
    integer,       intent(out)   :: ix(:)
 
    call index_to_coords(mesh%idx, ipg, ix)
  end subroutine mesh_global_index_to_coords
 
  integer function mesh_local_index_from_coords(mesh, ix) result(ip)
    type(mesh_t),  intent(in)    :: mesh
    integer,       intent(in)    :: ix(:)
 
    integer(int64) :: ipg
 
    ipg = index_from_coords(mesh%idx, ix)
    ip = mesh_global2local(mesh, ipg)
  end function mesh_local_index_from_coords
 
  subroutine mesh_local_index_to_coords(mesh, ip, ix)
    type(mesh_t),  intent(in)    :: mesh
    integer,       intent(in)    :: ip
    integer,       intent(out)   :: ix(:)
 
    integer(int64) :: ipg
 
    ipg = mesh_local2global(mesh, ip)
    call index_to_coords(mesh%idx, ipg, ix)
  end subroutine mesh_local_index_to_coords
 
  integer(int64) function mesh_local2global(mesh, ip) result(ipg)
    type(mesh_t),  intent(in)    :: mesh
    integer,       intent(in)    :: ip
 
    ipg = par_vec_local2global(mesh%pv, ip)
  end function mesh_local2global
 
  integer function mesh_global2local(mesh, ipg) result(ip)
    type(mesh_t),  intent(in)    :: mesh
    integer(int64),   intent(in)    :: ipg
 
    ip = par_vec_global2local(mesh%pv, ipg)
  end function mesh_global2local
 
  !-----------------------------------------------------------------------------
  subroutine mesh_minmaxloc(mesh, min_or_max, rank_min_or_max, op)
    type(mesh_t),  intent(in)    :: mesh
    real(real64),  intent(inout) :: min_or_max
    integer,       intent(out)   :: rank_min_or_max
    type(mpi_op),  intent(in)    :: op
 
    real(real64) :: loc_in(2), loc_out(2)
 
    push_sub(mesh_minmaxloc)
 
    assert(op == mpi_minloc .or. op == mpi_maxloc)
 
    rank_min_or_max = 0
    if (mesh%parallel_in_domains) then
      loc_in(1) = min_or_max
      loc_in(2) = mesh%mpi_grp%rank
      call mesh%mpi_grp%allreduce(loc_in, loc_out, 1, mpi_2double_precision, op)
      min_or_max = loc_out(1)
      rank_min_or_max = nint(loc_out(2))
    end if
 
    pop_sub(mesh_minmaxloc)
  end subroutine mesh_minmaxloc
 
  function mesh_red_min(mesh, space)
    class(mesh_t),  intent(in)    :: mesh
    type(space_t),  intent(in)    :: space
    real(real64) :: mesh_red_min(space%dim)
 
    integer :: idir
 
    push_sub(mesh_red_min)
 
    mesh_red_min = m_zero
    do idir = 1, space%periodic_dim
      mesh_red_min(idir) = real(mesh%idx%nr(1, idir) + mesh%idx%enlarge(idir), real64) &
        / real(mesh%idx%ll(idir), real64)
    end do
 
    pop_sub(mesh_red_min)
  end function mesh_red_min
 
 
 
#include "undef.F90"
#include "real.F90"
#include "mesh_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "mesh_inc.F90"
 
#include "undef.F90"
#include "integer.F90"
#include "mesh_inc.F90"
 
#include "undef.F90"
#include "integer8.F90"
#include "mesh_inc.F90"
 
end module mesh_oct_m
 
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: