mesh_cube_map.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!! 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_cube_map_oct_m
  use accel_oct_m
  use debug_oct_m
  use global_oct_m
  use index_oct_m
  use mesh_oct_m
  use messages_oct_m
  use profiling_oct_m
  use types_oct_m
 
  implicit none
 
  private
  public ::                               &
    mesh_cube_map_t,                      &
    mesh_cube_map_init,                   &
    mesh_cube_map_end
 
  type mesh_cube_map_t
    ! Components are public by default
    integer              :: nmap
    integer, allocatable :: map(:, :)
    type(accel_mem_t)    :: map_buffer
    logical              :: is_trivial = .false. 
  end type mesh_cube_map_t
 
  integer, public, parameter :: MCM_POINT = 4, mcm_count = 5
 
contains
 
  ! ---------------------------------------------------------
 
  subroutine mesh_cube_map_init(this, mesh, np)
    type(mesh_cube_map_t),         intent(out) :: this
    class(mesh_t),                 intent(in)  :: mesh
    integer,                       intent(in)  :: np
 
    integer :: i1(1:3), i2(1:3)
    integer :: step
    integer :: ip
    integer :: im, n1, n2, n3, expected_y, expected_z
 
    push_sub(mesh_cube_map_init)
 
    if (mesh%idx%dim <= 3) then
      do step = 1, 2
        if (step == 2) then
          safe_allocate(this%map(1:5, 1:this%nmap))
        end if
 
        this%nmap = 0
        i2 = 0
        do ip = 1, np
 
          i1 = 0
          call mesh_local_index_to_coords(mesh, ip, i1)
 
          if (any(i1(2:3) /= i2(2:3)) .or. i1(1) /= i2(1) + 1 .or. ip == 1) then
            this%nmap = this%nmap + 1
            if (step == 2) then
              this%map(1:3, this%nmap) = i1(1:3)
              this%map(mesh%idx%dim + 1:3, this%nmap) = 0
              this%map(mcm_point, this%nmap) = ip
              this%map(mcm_count, this%nmap) = 1
            end if
          else
            if (step == 2) this%map(mcm_count, this%nmap) = this%map(mcm_count, this%nmap) + 1
          end if
          i2 = i1
        end do
      end do
 
      ! Check whether the mesh ordering matches an X-fastest identity over the
      ! full inner dimensions (mesh%idx%ll). When true, the scatter degenerates
      ! to a flat memcopy and callers can short-circuit through a cheap path.
      this%is_trivial = .false.
      if (mesh%idx%dim == 3 .and. np == product(mesh%idx%ll(1:mesh%idx%dim)) &
        .and. this%nmap > 0) then
        n1 = mesh%idx%ll(1)
        n2 = mesh%idx%ll(2)
        n3 = mesh%idx%ll(3)
        if (this%nmap == n2 * n3) then
          ! Use segment 1 as the reference origin and verify every other
          ! segment follows the X-fastest, Y-medium, Z-slowest pattern.
          this%is_trivial = .true.
          do im = 1, this%nmap
            expected_y = this%map(2, 1) + mod(im - 1, n2)
            expected_z = this%map(3, 1) + (im - 1) / n2
            if (this%map(mcm_count, im) /= n1 .or. &
              this%map(mcm_point, im) /= (im - 1) * n1 + 1 .or. &
              this%map(1, im) /= this%map(1, 1) .or. &
              this%map(2, im) /= expected_y .or. &
              this%map(3, im) /= expected_z) then
              this%is_trivial = .false.
              exit
            end if
          end do
        end if
      end if
 
      if (accel_is_enabled()) then
        call accel_create_buffer(this%map_buffer, accel_mem_read_only, type_integer, this%nmap*5)
        call accel_write_buffer(this%map_buffer, 5, this%nmap, this%map)
      end if
 
    end if
 
    pop_sub(mesh_cube_map_init)
  end subroutine mesh_cube_map_init
 
  ! ---------------------------------------------------------
 
  subroutine mesh_cube_map_end(this)
    type(mesh_cube_map_t), intent(inout) :: this
 
    push_sub(mesh_cube_map_end)
 
    if (allocated(this%map)) then
 
      safe_deallocate_a(this%map)
 
      if (accel_is_enabled()) call accel_free_buffer(this%map_buffer)
 
    end if
 
    pop_sub(mesh_cube_map_end)
  end subroutine mesh_cube_map_end
 
  ! ---------------------------------------------------------
 
end module mesh_cube_map_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: