regridding.F90

Go to the documentation of this file.

!! Copyright (C) 2022, 2024 F. Bonafé, 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 regridding_oct_m
  use affine_coordinates_oct_m
  use coordinate_system_oct_m
  use debug_oct_m
  use global_oct_m
  use grid_oct_m
  use iihash_oct_m
  use lalg_adv_oct_m
  use math_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use namespace_oct_m
  use parser_oct_m
  use partition_oct_m
  use partition_transfer_oct_m
  use profiling_oct_m
  use space_oct_m
  use stencil_oct_m
  use utils_oct_m
 
  implicit none
  public :: regridding_t
 
  type :: regridding_t
    private
    class(mesh_t), pointer :: mesh_in, mesh_out
    type(partition_transfer_t) :: partition_transfer
    integer :: nsend, nrecv, dim
    integer,     allocatable :: order_in(:), order_out(:)
    integer(int64), allocatable :: order_in_global(:), order_out_global(:)
    logical, allocatable :: overlap_map(:)
    logical :: do_restriction, do_prolongation
    type(stencil_t) :: transfer_stencil
    real(real64), allocatable :: weights(:)
    real(real64), allocatable :: weights_generic(:, :)
    integer, allocatable :: stencil_points(:, :)
    integer, allocatable :: eta(:)
    integer :: interpolation_level
    integer :: scale_norms
    logical :: generic_interpolation
  contains
    procedure :: dregridding_do_transfer_1, zregridding_do_transfer_1
    procedure :: dregridding_do_transfer_2, zregridding_do_transfer_2
    generic :: do_transfer => dregridding_do_transfer_1, zregridding_do_transfer_1
    generic :: do_transfer => dregridding_do_transfer_2, zregridding_do_transfer_2
    final :: regridding_finalize
  end type regridding_t
 
  interface regridding_t
    procedure regridding_init
  end interface regridding_t
 
  integer, parameter :: &
    NEAREST_NEIGHBOR = 0, &
    linear = 1, &
    scale_none = 0, &
    scale_norm2 = 1
 
contains
 
  ! ---------------------------------------------------------
  function regridding_init(mesh_out, mesh_in, space_in, namespace) result(this)
    class(mesh_t), target, intent(in) :: mesh_out
    class(mesh_t), target, intent(in) :: mesh_in
    class(space_t),        intent(in) :: space_in
    type(namespace_t),     intent(in) :: namespace
    class(regridding_t), pointer :: this
 
    integer :: ip_in, ip_out, index(1:space_in%dim), idim, ii, is, size_array
    integer(int64), allocatable :: global_indices_in(:), order_in(:), order_out(:)
    integer(int64), allocatable :: order_in_original(:), order_out_original(:)
    integer(int64), allocatable :: global_indices_in_tmp(:)
    integer(int64), allocatable :: global_indices_in_original(:)
    integer, allocatable :: partition_in(:)
    integer(int64) :: ipg_coarse, ipg_in
    real(real64) :: spacing_ratio(1:space_in%dim)
    real(real64) :: x_fine(1:space_in%dim), x_coarse(1:space_in%dim), x_primitive(1:space_in%dim)
    real(real64) :: x_out(1:space_in%dim), x_in(1:space_in%dim), rmin, r_current
    real(real64), allocatable :: M(:, :), p(:)
    real(real64) :: scaling(1:space_in%dim, 1:space_in%dim), rotation(1:space_in%dim, 1:space_in%dim)
    integer :: i_coarse, ip_fine, ip_mindist, cube_size
    integer :: index_stencil(1:space_in%dim), shift(1:space_in%dim)
    logical :: same_eta, on_coarse_grid, found, on_in_grid
    type(lihash_t) :: global_indices
    class(mesh_t), pointer :: mesh_coarse, mesh_fine
    class(coordinate_system_t), pointer :: coarse_coord, fine_coord
 
    push_sub_with_profile(regridding_init)
 
    safe_allocate(this)
 
    nullify(coarse_coord)
    nullify(fine_coord)
 
    !%Variable RegriddingInterpolationLevel
    !%Type integer
    !%Default regridding_linear
    !%Section Mesh
    !%Description
    !% Choose the interpolation order for the regridding. 0 is equivalent to nearest neighbor
    !% interpolation. The default is linear interpolation.
    !%Option regridding_nearest_neighbor 0
    !% Use the nearest neighbor for the regridding. This is faster than the linear interpolation.
    !%Option regridding_linear 1
    !% Use a trilinear interpolation. This is implemented similar to restriction and prolongation
    !% operations in multigrid methods. This ensures that both directions are adjoint to each other.
    !%End
    call parse_variable(namespace, "RegriddingInterpolationLevel", linear, this%interpolation_level)
 
    !%Variable RegriddingRescale
    !%Type integer
    !%Default scale_norm2
    !%Section Mesh
    !%Description
    !% Rescale the regridded quantities. Not using a rescaling can lead to bad results if the
    !% ratio of the grid spacings is large. The default is to rescale by the 2-norm of the quantity
    !% except for generic interpolations (i.e. between curvilinear or non-commensurate grids).
    !%Option scale_none 0
    !% Do not rescale the regridded quantities.
    !%Option scale_norm2 1
    !% Scale the regridded quantities by the 2-norm of the quantity on the overlapping
    !% region of the grid.
    !%End
    call parse_variable(namespace, "RegriddingRescale", scale_norm2, this%scale_norms)
 
    ! check some conditions which are not yet supported
    if (.not. (mesh_out%coord_system%dim == space_in%dim)) then
      message(1) = "For regridding, both grids need to have the same space dimension."
      call messages_fatal(1, namespace=namespace)
    end if
    ! use generic interpolation if the grids are not commensurate or curvilinear
    this%generic_interpolation = mesh_out%coord_system%local_basis .or. &
      mesh_in%coord_system%local_basis .or. &
      .not. (same_type_as(mesh_out%coord_system, mesh_in%coord_system))
 
    this%dim = space_in%dim
    spacing_ratio(:) = mesh_out%spacing(1:this%dim)/mesh_in%spacing(1:this%dim)
 
    this%do_restriction = all(spacing_ratio > m_one)
    this%do_prolongation = all(spacing_ratio < m_one)
    ! invert spacing ratio for prolongations
    if (this%do_prolongation) spacing_ratio = m_one/spacing_ratio
    ! get the integer ratio of the spacings
    safe_allocate(this%eta(1:this%dim))
    do idim = 1, this%dim
      this%eta(idim) = nint(spacing_ratio(idim))
    end do
    if (any(abs((spacing_ratio - this%eta)/spacing_ratio) > 10.0_real64*m_epsilon)) then
      this%generic_interpolation = .true.
    end if
    same_eta = .true.
    do idim = 2, this%dim
      same_eta = same_eta .and. this%eta(idim) == this%eta(idim-1)
    end do
    if (.not. same_eta) then
      this%generic_interpolation = .true.
    end if
    ! use always generic interpolation for nearest neighbor interpolation
    if (this%interpolation_level == nearest_neighbor) then
      this%generic_interpolation = .true.
    end if
 
    if (this%generic_interpolation) then
      ! Generic interpolation samples the input at the points of the output and does not
      ! integrate, so a quantity from a finer input mesh silently loses the values between
      ! those points. Warn in that case; anisotropic ratios reach this path by construction.
      if (any(mesh_out%spacing(1:this%dim) > mesh_in%spacing(1:this%dim))) then
        write(message(1), '(a)') "Regridding falls back to generic interpolation between meshes"
        write(message(2), '(a)') "whose spacings differ non-uniformly, and the input mesh is finer"
        write(message(3), '(a)') "in at least one direction. The transfer interpolates (samples) the input at"
        write(message(4), '(a)') "the output points instead of averaging over output cells."
        call messages_warning(4, namespace=namespace)
      end if
      this%do_prolongation = .true.
      this%do_restriction = .false.
      this%scale_norms = scale_none
    end if
 
    this%mesh_in => mesh_in
    this%mesh_out => mesh_out
    if (.not. this%do_prolongation) then
      mesh_fine => this%mesh_in
      mesh_coarse => this%mesh_out
    else
      mesh_coarse => this%mesh_in
      mesh_fine => this%mesh_out
    end if
 
    ! collect all locally available points in mesh_coarse that are also in mesh_fine
    call lihash_init(global_indices)
    size_array = mesh_fine%np
    safe_allocate(global_indices_in_tmp(size_array))
    safe_allocate(global_indices_in_original(size_array))
 
    if (.not. this%generic_interpolation) then
      ii = 0
      do ip_fine = 1, mesh_fine%np
        call mesh_local_index_to_coords(mesh_fine, ip_fine, index)
        on_coarse_grid = .true.
        if (this%do_restriction .or. this%do_prolongation) then
          do idim = 1, this%dim
            on_coarse_grid = on_coarse_grid .and. mod(index(idim), this%eta(idim)) == 0
          end do
        end if
        if (on_coarse_grid) then
          ! translate between fine and coarse grid
          index(:) = index(:) / this%eta(:)
          ipg_coarse = mesh_global_index_from_coords(mesh_coarse, index)
          call insert_global_point(mesh_coarse, ipg_coarse, ii)
        else if (this%do_restriction .or. this%do_prolongation) then
          ! now get all coarse points that surround the fine point
          ! start always from the same corner of the cube
          ! need a cube of size of the interpolation level, so shift to the corner of that cube
          index(:) = floor(real(index(:)) / this%eta(:)) - (this%interpolation_level-1)/2
          cube_size = 1+this%interpolation_level
          do is = 1, cube_size**this%dim
            ! get index of coarse point in the surrounding parallelepiped
            ! using the corresponding index as number in base cube_size
            call convert_to_base(is-1, cube_size, shift)
            index_stencil(:) = index(:) + shift(:)
            ipg_coarse = mesh_global_index_from_coords(mesh_coarse, index_stencil)
            call insert_global_point(mesh_coarse, ipg_coarse, ii)
          end do
        end if
      end do
      call lihash_end(global_indices)
    else
      ii = 0
      do ip_out = 1, mesh_out%np
        call mesh_local_index_to_coords(mesh_out, ip_out, index)
        x_out = mesh_out%coord_system%to_cartesian(real(index, real64)*mesh_out%spacing)
        x_in = mesh_in%coord_system%from_cartesian(x_out)/mesh_in%spacing
        on_in_grid = all(abs((x_in - nint(x_in))) < 10.0_real64*m_epsilon)
        if (on_in_grid) then
          ! translate between out and in grid
          ipg_in = mesh_global_index_from_coords(mesh_in, nint(x_in))
          call insert_global_point(mesh_in, ipg_in, ii)
        else
          ! now get all points in the input mesh that surround the point in the output mesh
          ! start always from the same corner of the cube
          ! need a cube of size of the interpolation level, so shift to the corner of that cube
          index(:) = floor(x_in) - (max(this%interpolation_level, 1)-1)/2
          cube_size = 1+max(this%interpolation_level, 1)
          do is = 1, cube_size**this%dim
            ! get index of coarse point in the surrounding parallelepiped
            ! using the corresponding index as number in base cube_size
            call convert_to_base(is-1, cube_size, shift)
            index_stencil(:) = index(:) + shift(:)
            ipg_in = mesh_global_index_from_coords(mesh_in, index_stencil)
            call insert_global_point(mesh_in, ipg_in, ii)
          end do
        end if
      end do
      call lihash_end(global_indices)
    end if
 
    safe_allocate(global_indices_in(ii))
    safe_allocate(partition_in(ii))
    global_indices_in(1:ii) = global_indices_in_tmp(1:ii)
    safe_deallocate_a(global_indices_in_tmp)
 
    if (mesh_coarse%parallel_in_domains) then
      ! determine where the points of the coarse mesh are stored
      call partition_get_partition_number(mesh_coarse%partition, ii, global_indices_in, partition_in)
    else
      partition_in = 1
    end if
 
    ! Init partition transfer
    ! need inverse direction for prolongations
    call partition_transfer_init(this%partition_transfer, ii, global_indices_in, &
      mesh_in%mpi_grp, mesh_out%mpi_grp, partition_in, &
      this%nsend, this%nrecv, order_in, order_out, inverse=this%do_prolongation)
 
    ! we need to transfer the global indices of the boundary point and the corresponding
    ! inner point for periodic meshes
    if (space_in%is_periodic()) then
      global_indices_in(1:ii) = global_indices_in_original(1:ii)
      call partition_transfer_init(this%partition_transfer, ii, global_indices_in, &
        mesh_in%mpi_grp, mesh_out%mpi_grp, partition_in, &
        this%nsend, this%nrecv, order_in_original, order_out_original, inverse=this%do_prolongation)
    else
      order_in_original = order_in
      order_out_original = order_out
    end if
 
    ! convert from global to local indices
    safe_allocate(this%order_in(1:this%nsend))
    safe_allocate(this%order_in_global(1:this%nsend))
    safe_allocate(this%order_out(1:this%nrecv))
    safe_allocate(this%order_out_global(1:this%nrecv))
 
    ! get the mapping for mesh_in in the order given by the global indices of mesh_out
    do ip_in = 1, this%nsend
      if (this%do_restriction) then
        ! in this case, we have an index on mesh_out
        this%order_in_global(ip_in) = order_in_original(ip_in)
        this%order_in(ip_in) = 0
      else if (this%do_prolongation) then
        ! in this case, we have an index on mesh_in
        this%order_in_global(ip_in) = order_in(ip_in)
        this%order_in(ip_in) = mesh_global2local(mesh_in, this%order_in_global(ip_in))
      else
        ! convert back from the global index of mesh_out to a local index of mesh_in
        call mesh_global_index_to_coords(mesh_out, order_in_original(ip_in), index)
        this%order_in_global(ip_in) = mesh_global_index_from_coords(mesh_in, index)
        this%order_in(ip_in) = mesh_global2local(mesh_in, this%order_in_global(ip_in))
      end if
      if (.not. this%do_restriction) then
        if (this%order_in(ip_in) == 0) then
          write(message(1),'(a,i10,a,i10)') "Error in regridding part 1: mesh point ", &
            this%order_in(ip_in), " is not stored in partition ", mesh_in%pv%partno
          call messages_fatal(1, namespace=namespace)
        end if
      end if
    end do
 
    ! for mapping back to the global grid after the transfer, convert to local indices of mesh_out
    do ip_out = 1, this%nrecv
      if (.not. this%do_prolongation) then
        this%order_out_global(ip_out) = order_out(ip_out)
        this%order_out(ip_out) = mesh_global2local(mesh_out, order_out(ip_out))
      else
        ! store the global index of mesh_in
        this%order_out_global(ip_out) = order_out_original(ip_out)
        this%order_out(ip_out) = 0
      end if
      if (.not. this%do_prolongation) then
        if (this%order_out(ip_out) == 0) then
          write(message(1),'(a,i10,a,i10)') "Error in regridding part 2: mesh point ", &
            this%order_out(ip_out), " is not stored in partition ", mesh_out%pv%partno
          call messages_fatal(1, namespace=namespace)
        end if
      end if
    end do
 
    ! remove rotations for affine coordinates because they can lead to singularities in the
    ! computation of the weights; the weights are independent of rotations
    ! remove the rotation for the coarse mesh and use the same rotation matrix for the fine mesh
    select type(coord_system => mesh_coarse%coord_system)
    class is (affine_coordinates_t)
      scaling = lalg_remove_rotation(this%dim, coord_system%basis%vectors)
      coarse_coord => affine_coordinates_t(namespace, this%dim, scaling)
      ! get rotation matrix from polar decomposition A = U P -> U = A P^{-1}
      call lalg_inverse(this%dim, scaling, "dir")
      rotation = matmul(coord_system%basis%vectors, scaling)
    class default
      coarse_coord => mesh_coarse%coord_system
      rotation = diagonal_matrix(this%dim, m_one)
    end select
    select type(coord_system => mesh_fine%coord_system)
    class is (affine_coordinates_t)
      ! here we compute the new basis as U^T A
      fine_coord => affine_coordinates_t(namespace, this%dim, &
        matmul(transpose(rotation), coord_system%basis%vectors))
    class default
      fine_coord => mesh_fine%coord_system
    end select
    ! create transfer stencil
    if (.not. this%generic_interpolation) then
      if (this%do_restriction .or. this%do_prolongation) then
        ! determine weights by requiring a set of polynomials to be exactly reproduced
        ! on the corners of the polyhedron generated by neighbouring grid points
        call get_transfer_stencil(this%transfer_stencil, this%dim, this%eta(1), this%interpolation_level)
        safe_allocate(this%weights(1:this%transfer_stencil%size))
        ! get M matrix
        cube_size = 1+this%interpolation_level
        safe_allocate(m(1:cube_size**this%dim, 1:cube_size**this%dim))
        safe_allocate(p(1:cube_size**this%dim))
        do ii = 1, cube_size**this%dim
          ! get index of coarse point in the surrounding parallelepiped
          ! using the corresponding index as number in base cube_size
          call convert_to_base(ii-1, cube_size, shift)
          x_primitive = (shift - (this%interpolation_level-1)/2)*mesh_coarse%spacing
          x_coarse = coarse_coord%to_cartesian(x_primitive)
          ! fill one row of M with the polynomials
          m(ii, :) = evaluate_polynomials(x_coarse)
        end do
        call lalg_inverse(cube_size**this%dim, m, "dir")
        do ii = 1, this%transfer_stencil%size
          ! Cartesian coordinates of point in fine mesh
          shift = -floor(real(this%transfer_stencil%points(:, ii), real64)/real(this%eta, real64))
          x_fine = real(this%transfer_stencil%points(:, ii), real64) + shift*this%eta
          x_fine = fine_coord%to_cartesian(x_fine * mesh_fine%spacing)
          ! get polynomials
          p = evaluate_polynomials(x_fine)
          ! determine coefficients for interpolation
          p = matmul(p, m)
          ! take the coefficient for the correct coarse point
          call convert_from_base(shift+(this%interpolation_level-1)/2, cube_size, i_coarse)
          i_coarse = i_coarse + 1
          this%weights(ii) = p(i_coarse)
        end do
        ! for the restriction, we need to take into account the ratio of volume elements
        ! like this, restriction is adjoint to prolongation
        if (this%do_restriction) then
          this%weights(:) = this%weights(:) * mesh_fine%volume_element/mesh_coarse%volume_element
        end if
        safe_deallocate_a(m)
        safe_deallocate_a(p)
      end if
    else
      ! compute interpolation weights for each point
      ! we always need at least a cube of size 2
      cube_size = 1+max(this%interpolation_level, 1)
      safe_allocate(this%weights_generic(1:cube_size**this%dim, 1:mesh_out%np))
      this%weights_generic = m_zero
      safe_allocate(this%stencil_points(1:cube_size**this%dim, 1:mesh_out%np))
      this%stencil_points = -1
 
      call lihash_init(global_indices)
      ! create hash map for the global indices of the input mesh
      do ip_out = 1, this%nrecv
        ! it is called ip_out because it is the output of the partition_transfer
        call lihash_insert(global_indices, this%order_out_global(ip_out), ip_out)
      end do
 
      safe_allocate(m(1:cube_size**this%dim, 1:cube_size**this%dim))
      safe_allocate(p(1:cube_size**this%dim))
      do ip_out = 1, mesh_out%np
        call mesh_local_index_to_coords(mesh_out, ip_out, index)
        x_out = fine_coord%to_cartesian(real(index, real64)*mesh_out%spacing)
        x_in = coarse_coord%from_cartesian(x_out)/mesh_in%spacing
        on_in_grid = all(abs((x_in - nint(x_in))) < 10.0_real64*m_epsilon)
        if (on_in_grid) then
          ! the points coincide, so use only one point in the stencil with a weight of 1
          ipg_in = mesh_global_index_from_coords(mesh_in, nint(x_in))
          this%stencil_points(1, ip_out) = lihash_lookup(global_indices, ipg_in, found)
          this%weights_generic(1, ip_out) = m_one
        else
          ! now get all in points that surround the out point
          ! start always from the same corner of the cube
          ! need a cube of size of the interpolation level, so shift to the corner of that cube
          index(:) = floor(x_in) - (max(this%interpolation_level, 1)-1)/2
          rmin = m_huge
          ip_mindist = -1
          do is = 1, cube_size**this%dim
            ! get index of coarse point in the surrounding parallelepiped
            ! using the corresponding index as number in base cube_size
            call convert_to_base(is-1, cube_size, shift)
            index_stencil(:) = index(:) + shift(:)
            ipg_in = mesh_global_index_from_coords(mesh_in, index_stencil)
            ! the lookup will return -1 for points that are not locally available in the input mesh
            ! those will be ignored during the regridding
            ! the points we need from other processes are already in this hash table
            this%stencil_points(is, ip_out) = lihash_lookup(global_indices, ipg_in, found)
            x_in = coarse_coord%to_cartesian(index_stencil*mesh_in%spacing)
            ! compute M
            ! fill one row of M with the polynomials
            m(is, :) = evaluate_polynomials(x_in)
            ! determine closest point for nearest neighbor interpolation
            r_current = norm2(x_in - x_out)
            ! avoid rounding errors by requesting the distance to be smaller
            ! by a certain amount
            if (r_current < rmin - 10.*m_epsilon .and. found) then
              rmin = r_current
              ip_mindist = this%stencil_points(is, ip_out)
            end if
          end do
          if (this%interpolation_level == nearest_neighbor) then
            ! only use one point for th nearest-neighbor interpolation
            this%weights_generic(1, ip_out) = m_one
            this%stencil_points(:, ip_out) = -1
            this%stencil_points(1, ip_out) = ip_mindist
          else
            call lalg_inverse(cube_size**this%dim, m, "dir")
            ! get polynomials at point to be interpolated
            p = evaluate_polynomials(x_out)
            ! determine coefficients for interpolation
            this%weights_generic(:, ip_out) = matmul(p, m)
          end if
        end if
      end do
      safe_deallocate_a(m)
      safe_deallocate_a(p)
      call lihash_end(global_indices)
    end if
 
    select case (this%scale_norms)
    case(scale_none)
      ! do nothing
    case(scale_norm2)
      ! create overlap map
      safe_allocate(this%overlap_map(1:this%mesh_in%np))
      this%overlap_map = this%mesh_out%box%contains_points(this%mesh_in%np, this%mesh_in%x_t)
    end select
 
    safe_deallocate_a(partition_in)
    safe_deallocate_a(global_indices_in)
    safe_deallocate_a(order_in)
    safe_deallocate_a(order_out)
 
    ! deallocate coordinate systems correctly
    select type(coord_system => mesh_coarse%coord_system)
    class is (affine_coordinates_t)
      safe_deallocate_p(coarse_coord)
    class default
      nullify(coarse_coord)
    end select
    select type(coord_system => mesh_fine%coord_system)
    class is (affine_coordinates_t)
      safe_deallocate_p(fine_coord)
    class default
      nullify(fine_coord)
    end select
 
    pop_sub_with_profile(regridding_init)
  contains
    subroutine insert_global_point(mesh, ipg, ii)
      class(mesh_t),  intent(in)    :: mesh
      integer(int64), intent(in)    :: ipg
      integer,        intent(inout) :: ii
 
      integer :: ip_tmp
      integer(int64) :: ipg_tmp
 
      ! take boundary into account for periodic meshes
      if ((ipg > 0 .and. ipg <= mesh%np_global) .or. &
        (ipg > 0 .and. space_in%is_periodic())) then
        ip_tmp = lihash_lookup(global_indices, ipg, found)
        if (.not. found) then
          if (ipg > mesh%np_global) then
            ipg_tmp = mesh_periodic_point_global(mesh, space_in, ipg)
            ! early exit if periodic copy is not inside the mesh
            ! this can happen for partly periodic systems
            if (ipg_tmp > mesh%np_global) return
          end if
          ii = ii + 1
          ! enlarge array if necessary
          if (ii >= size_array) then
            size_array = size_array * 2
            call make_array_larger(global_indices_in_tmp, size_array)
            call make_array_larger(global_indices_in_original, size_array)
          end if
          if (ipg <= mesh%np_global) then
            global_indices_in_tmp(ii) = ipg
            global_indices_in_original(ii) = ipg
          else if (ipg > mesh%np_global) then
            ! save the index of the corresponding inner point
            global_indices_in_tmp(ii) = ipg_tmp
            global_indices_in_original(ii) = ipg
          end if
          call lihash_insert(global_indices, ipg, ii)
        end if
      end if
    end subroutine insert_global_point
 
    function evaluate_polynomials(x)
      real(real64), intent(in) :: x(1:this%dim)
      real(real64) :: evaluate_polynomials(1:(1+this%interpolation_level)**this%dim)
      integer :: cube_size, i, j, shift(1:this%dim)
 
      evaluate_polynomials = m_one
      cube_size = 1+this%interpolation_level
      do i = 1, cube_size**this%dim
        ! get index of coarse point in the surrounding parallelepiped
        ! using the corresponding index as number in base cube_size
        call convert_to_base(i-1, cube_size, shift)
        ! evaluate polynomials with the corresponding powers
        do j = 1, this%dim
          evaluate_polynomials(i) = evaluate_polynomials(i) * x(j)**shift(j)
        end do
      end do
    end function evaluate_polynomials
 
    subroutine get_transfer_stencil(this, dim, eta, order)
      type(stencil_t), intent(inout) :: this
      integer,         intent(in)    :: dim
      integer,         intent(in)    :: eta
      integer,         intent(in)    :: order
 
      integer :: i, offset, cube_size, shift(1:dim)
 
      push_sub(get_transfer_stencil)
 
      cube_size = eta*(order + 1) - 1
      call stencil_allocate(this, dim, cube_size**dim)
 
      offset = eta * (order/2 + 1) - 1
      do i = 1, cube_size**dim
        call convert_to_base(i-1, cube_size, shift)
        this%points(:, i) = shift(:) - offset
      end do
 
      call stencil_init_center(this)
 
      pop_sub(get_transfer_stencil)
    end subroutine get_transfer_stencil
  end function regridding_init
 
  subroutine regridding_finalize(this)
    type(regridding_t), intent(inout) :: this
 
    push_sub(regridding_finalize)
 
    call partition_transfer_end(this%partition_transfer)
    safe_deallocate_a(this%eta)
    safe_deallocate_a(this%order_in)
    safe_deallocate_a(this%order_in_global)
    safe_deallocate_a(this%order_out)
    safe_deallocate_a(this%order_out_global)
    safe_deallocate_a(this%weights)
    safe_deallocate_a(this%weights_generic)
    safe_deallocate_a(this%stencil_points)
    select case (this%scale_norms)
    case(scale_none)
      ! do nothing
    case(scale_norm2)
      safe_deallocate_a(this%overlap_map)
    end select
 
    pop_sub(regridding_finalize)
  end subroutine regridding_finalize
 
#include "real.F90"
#include "regridding_inc.F90"
#include "undef.F90"
 
#include "complex.F90"
#include "regridding_inc.F90"
#include "undef.F90"
 
end module regridding_oct_m