scale_mesh_cubedom3d.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Mesh / Cubic 3D domain
!!
!! @par Description
!!      Manage mesh data of cubic 3D domain for element-based methods
!!
!! @author Yuta Kawai, Team SCALE
!<
#include "scaleFElib.h"
module scale_mesh_cubedom3d
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_io
  use scale_precision
  use scale_prc
 
  use scale_mesh_base3d, only: &
    meshbase3d, meshbase3d_init, meshbase3d_final, &
    meshbase3d_setgeometricinfo
 
  use scale_mesh_base2d, only: &
    meshbase2d, meshbase2d_init, meshbase2d_final, &
    meshbase2d_setgeometricinfo
 
  use scale_localmesh_3d, only: &
    localmesh3d
  use scale_element_base, only: elementbase3d
  use scale_element_hexahedral, only: hexahedralelement
 
  use scale_mesh_rectdom2d, only: &
    meshrectdom2d, meshrectdom2d_setuplocaldom
  
  use scale_element_quadrilateral, only: quadrilateralelement
 
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  ! 
  !> Derived type to manage a cubic 3D computational domain
  type, extends(meshbase3d), public :: meshcubedom3d
    integer :: negx   !< Number of elements in X direction for the global mesh
    integer :: negy   !< Number of elements in Y direction for the global mesh
    integer :: negz   !< Number of elements in Z direction for the global mesh
 
    integer :: nprcx  !< Number of processes in X direction for domain decomposition
    integer :: nprcy  !< Number of processes in Y direction for domain decomposition
    integer :: nprcz  !< Number of processes in Z direction for domain decomposition
    
    real(rp), public :: xmin_gl !< Minimum X coordinate of the global domain
    real(rp), public :: xmax_gl !< Maximum X coordinate of the global domain
    real(rp), public :: ymin_gl !< Minimum Y coordinate of the global domain
    real(rp), public :: ymax_gl !< Maximum Y coordinate of the global domain
    real(rp), public :: zmin_gl !< Minimum Z coordinate of the global domain
    real(rp), public :: zmax_gl !< Maximum Z coordinate of the global domain
 
    real(rp), allocatable :: fz(:)
 
    integer, allocatable :: rcdomijk2lcmeshid(:,:,:)
 
    logical :: isperiodicx  !< Flag whether the domain is periodic in X direction
    logical :: isperiodicy  !< Flag whether the domain is periodic in Y direction
    logical :: isperiodicz  !< Flag whether the domain is periodic in Z direction
 
    type(meshrectdom2d) :: mesh2d             !< 2D mesh for the horizontal plane (X-Y plane)
    type(quadrilateralelement) :: refelem2d   !< Reference element for the 2D mesh
  contains
    procedure :: init => meshcubedom3d_init
    procedure :: final => meshcubedom3d_final
    procedure :: generate => meshcubedom3d_generate
    procedure :: getmesh2d => meshcubedom3d_getmesh2d
    procedure :: set_geometric_with_vcoord => meshcubedom3d_set_geometric_with_vcoord
  end type meshcubedom3d
 
  public :: meshcubedom3d_coord_conv
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
 
contains
!> Initialize an object to manage a cubic 3D domain
!OCL SERIAL
  subroutine meshcubedom3d_init( this,                          &
    NeGX, NeGY, NeGZ,                                           &
    dom_xmin, dom_xmax, dom_ymin, dom_ymax, dom_zmin, dom_zmax, &
    isPeriodicX, isPeriodicY, isPeriodicZ,                      &
    refElem, NLocalMeshPerPrc, NprcX, NprcY,                    &
    nproc, myrank,                                              &
    FZ                                                          )
    
    implicit none
 
    class(meshcubedom3d), intent(inout) :: this
    integer, intent(in) :: NeGX          !< Number of elements in X direction for the global mesh
    integer, intent(in) :: NeGY          !< Number of elements in Y direction for the global mesh
    integer, intent(in) :: NeGZ          !< Number of elements in Z direction for the global mesh
    real(RP), intent(in) :: dom_xmin     !< Minimum X coordinate of the global domain
    real(RP), intent(in) :: dom_xmax     !< Maximum X coordinate of the global domain
    real(RP), intent(in) :: dom_ymin     !< Minimum Y coordinate of the global domain
    real(RP), intent(in) :: dom_ymax     !< Maximum Y coordinate of the global domain
    real(RP), intent(in) :: dom_zmin     !< Minimum Z coordinate of the global domain
    real(RP), intent(in) :: dom_zmax     !< Maximum Z coordinate of the global domain
    logical, intent(in) :: isPeriodicX   !< Flag whether the domain is periodic in X direction
    logical, intent(in) :: isPeriodicY   !< Flag whether the domain is periodic in Y direction
    logical, intent(in) :: isPeriodicZ   !< Flag whether the domain is periodic in Z direction
    type(hexahedralelement), intent(in), target :: refElem  !< Reference element for the 3D mesh
    integer, intent(in) :: NLocalMeshPerPrc  !< Number of local meshes managed by each process
    integer, intent(in) :: NprcX         !< Number of processes in X direction for domain decomposition
    integer, intent(in) :: NprcY         !< Number of processes in Y direction for domain decomposition
    integer, intent(in), optional :: nproc       !< Total number of processes (if not provided, it will be determined from the parallel environment)
    integer, intent(in), optional :: myrank      !< Rank of the current process (if not provided, it will be determined from the parallel environment)
    real(RP), intent(in), optional :: FZ(NeGZ+1) !< Optional array of Z coordinates for the mesh
 
    integer :: k
    real(RP) :: dz
    !-----------------------------------------------------------------------------
    
    this%NeGX = negx
    this%NeGY = negy
    this%NeGZ = negz
 
    this%xmin_gl       = dom_xmin
    this%xmax_gl       = dom_xmax
    this%ymin_gl       = dom_ymin
    this%ymax_gl       = dom_ymax
    this%zmin_gl       = dom_zmin
    this%zmax_gl       = dom_zmax
    this%dom_vol = (this%xmax_gl - this%xmin_gl) * (this%ymax_gl - this%ymin_gl) * (this%zmax_gl - this%zmin_gl)  
 
    this%isPeriodicX = isperiodicx
    this%isPeriodicY = isperiodicy
    this%isPeriodicZ = isperiodicz
 
    this%NprcX = nprcx
    this%NprcY = nprcy
    this%NprcZ = 1
 
    !- Fz
    allocate( this%FZ(this%NeGZ+1) )
    if ( present(fz) ) then
      this%FZ(:) = fz(:)
    else
      this%FZ(1          ) = dom_zmin
      this%FZ(this%NeGZ+1) = dom_zmax
      dz = (dom_zmax - dom_zmin) / dble(this%NeGZ)
      do k=2, this%NeGZ
        this%FZ(k) = this%FZ(k-1) + dz
      end do
    end if
 
    !--
    call meshbase3d_init( this, refelem, nlocalmeshperprc, 6, &
                          nproc, myrank                       )
 
    !--- 2D mesh
 
    call this%refElem2D%Init( this%refElem3D%PolyOrder_h, refelem%IsLumpedMatrix() )
 
    call this%mesh2D%Init( this%NeGX, this%NeGY, dom_xmin, dom_xmax, dom_ymin, dom_ymax, &
      isperiodicx, isperiodicy, this%refElem2D, nlocalmeshperprc, &
      nprcx, nprcy, nproc, myrank )
 
    return
  end subroutine meshcubedom3d_init
 
!> Finalize an object to manage a cubic 3D domain
!OCL SERIAL
  subroutine meshcubedom3d_final( this )
    implicit none
    class(meshcubedom3d), intent(inout) :: this
    !-----------------------------------------------------------------------------
  
    if (this%isGenerated) then
      if ( allocated(this%rcdomIJK2LCMeshID) ) then
        !$acc exit data delete( this%rcdomIJK2LCMeshID )
        deallocate( this%rcdomIJK2LCMeshID )
      end if
    else
      if ( allocated( this%FZ ) ) deallocate( this%FZ )
    end if
 
    call this%mesh2D%Final()
    call this%refElem2D%Final()
 
    call meshbase3d_final(this)
 
    return
  end subroutine meshcubedom3d_final
  
!OCL SERIAL
  subroutine meshcubedom3d_getmesh2d( this, ptr_mesh2D )
    implicit none
    class(meshcubedom3d), intent(in), target :: this
    class(meshbase2d), pointer, intent(out) :: ptr_mesh2D
    !-------------------------------------------------------
 
    ptr_mesh2d => this%mesh2D
    return
  end subroutine meshcubedom3d_getmesh2d
 
!> Generate a cubic 3D computational domain based on the parameters set in the initialization
!OCL SERIAL
  subroutine meshcubedom3d_generate( this )
    implicit none
 
    class(meshcubedom3d), intent(inout), target :: this
            
    integer :: n
    type(localmesh3d), pointer :: mesh
 
    integer :: tileID_table(this%LOCAL_MESH_NUM, this%PRC_NUM)
    integer :: panelID_table(this%LOCAL_MESH_NUM*this%PRC_NUM)
    integer :: pi_table(this%LOCAL_MESH_NUM*this%PRC_NUM)
    integer :: pj_table(this%LOCAL_MESH_NUM*this%PRC_NUM)
    integer :: pk_table(this%LOCAL_MESH_NUM*this%PRC_NUM)
 
    integer :: TILE_NUM_PER_PANEL
    integer :: tileID
    !-----------------------------------------------------------------------------
 
    tile_num_per_panel = this%LOCAL_MESH_NUM_global / 1
    
    
    !--- Construct the connectivity of patches  (only master node)
 
    call messhcubedom3d_assigndomid( this,    & ! (in)
      tileid_table, panelid_table,            & ! (out)
      pi_table, pj_table, pk_table )            ! (out)
    
    !--- Setup local meshes managed by my process
 
    do n=1, this%LOCAL_MESH_NUM
      mesh => this%lcmesh_list(n)
      tileid = tileid_table(n, mesh%PRC_myrank+1)
 
      call meshcubedom3d_setuplocaldom( mesh, &
        tileid,  panelid_table(tileid),                                                           &
        pi_table(tileid), pj_table(tileid), pk_table(tileid), this%NprcX, this%NprcY, this%NprcZ, &
        this%xmin_gl, this%xmax_gl, this%ymin_gl, this%ymax_gl, this%zmin_gl, this%zmax_gl,       &
        this%NeGX/this%NprcX, this%NeGY/this%NprcY, this%NeGZ/this%NprcZ, this%FZ(:)              )
 
      call meshrectdom2d_setuplocaldom( this%mesh2D%lcmesh_list(n), &
        tileid,  panelid_table(tileid),                             &
        pi_table(tileid), pj_table(tileid), this%NprcX, this%NprcY, &
        this%xmin_gl, this%xmax_gl, this%ymin_gl, this%ymax_gl,     &
        this%NeGX/this%NprcX, this%NeGY/this%NprcY )
 
      call mesh%SetLocalMesh2D( this%mesh2D%lcmesh_list(n) )
      !---
      ! write(*,*) "** my_rank=", mesh%PRC_myrank
      ! write(*,*) " tileID:", mesh%tileID
      ! write(*,*) " pnlID:", mesh%panelID, "-- i,j (within a panel)=", pi_table(tileID), pj_table(tileID)
      ! write(*,*) " local mesh:", n, "( total", this%LOCAL_MESH_NUM, ")"
      ! write(*,*) " panel_connect:", this%tilePanelID_globalMap(:,mesh%tileID)
      ! write(*,*) " tile_connect:", this%tileID_globalMap(:,mesh%tileID)
      ! write(*,*) " face_connect:", this%tileFaceID_globalMap(:,mesh%tileID)
      ! write(*,*) " domain size"
      ! write(*,*) "   NeX, NeY:", mesh%NeX, mesh%NeY
      ! write(*,*) "   [X], [Y]:",  mesh%xmin, mesh%xmax, ":", mesh%ymin, mesh%ymax
    end do
 
    ! To set rcdomIJP2LCMeshID, call AssignDomID for 2D mesh
    call this%mesh2D%AssignDomID( & 
      tileid_table, panelid_table,   & ! (out)
      pi_table, pj_table             ) ! (out)
 
    this%isGenerated = .true.
    this%mesh2D%isGenerated = .true.
    
    return
  end subroutine meshcubedom3d_generate
 
  !> Set geometric information of the local mesh with the vertical coordinate transformation
!OCL SERIAL
  subroutine meshcubedom3d_set_geometric_with_vcoord(this, lcdomID, GsqrtV_lc, zlev_lc, G13_lc, G23_lc)
    implicit none
    class(meshcubedom3d), intent(inout), target :: this
    integer, intent(in) :: lcdomID
    real(RP), intent(in) :: GsqrtV_lc(this%refElem3D%Np,this%lcmesh_list(lcdomID)%NeA)
    real(RP), intent(in) :: zlev_lc(this%refElem3D%Np,this%lcmesh_list(lcdomID)%NeA)
    real(RP), intent(in) :: G13_lc(this%refElem3D%Np,this%lcmesh_list(lcdomID)%NeA)
    real(RP), intent(in) :: G23_lc(this%refElem3D%Np,this%lcmesh_list(lcdomID)%NeA)
 
    integer :: ke, p
    class(localmesh3d), pointer :: lcmesh
    integer :: Np
    !-------------------------------------------------------
 
    lcmesh => this%lcmesh_list(lcdomid)
    np = lcmesh%refElem%Np
 
    !$omp parallel private(ke)
    !$omp do
    !$acc parallel present(lcmesh%zlev, lcmesh%Gsqrt, lcmesh%GI3, zlev_lc, GsqrtV_lc, G13_lc, G23_lc)
    !$acc loop gang
    do ke=lcmesh%NeS, lcmesh%NeE
      !$acc loop vector
      do p=1, np
        lcmesh%zlev(p,ke) = zlev_lc(p,ke)
      end do
    end do
    !$omp do
    !$acc loop gang
    do ke=lcmesh%NeS, lcmesh%NeA
      !$acc loop vector
      do p=1, np
        lcmesh%Gsqrt(p,ke) = gsqrtv_lc(p,ke)
        lcmesh%GI3(p,ke,1) = g13_lc(p,ke)
        lcmesh%GI3(p,ke,2) = g23_lc(p,ke)
      end do
    end do
    !$acc end parallel
    !$omp end parallel
    !$acc update host(lcmesh%zlev, lcmesh%Gsqrt, lcmesh%GI3)
 
    return
  end subroutine meshcubedom3d_set_geometric_with_vcoord
 
  !- private ------------------------------------------------------
 
!OCL SERIAL
  subroutine meshcubedom3d_setuplocaldom( lcmesh, &
    tileID, panelID,                                            &
    i, j, k, NprcX, NprcY, NprcZ,                               &
    dom_xmin, dom_xmax, dom_ymin, dom_ymax, dom_zmin, dom_zmax, &
    NeX, NeY, NeZ,                                              &
    FZ                                                          )
 
    use scale_prof
    use scale_meshutil_3d, only: &
      meshutil3d_genconnectivity,   &
      meshutil3d_gencubedomain,     &
      meshutil3d_buildinteriormap,  &
      meshutil3d_genpatchboundarymap
 
    use scale_localmesh_base, only: bctype_interior
    implicit none
      
    type(localmesh3d), intent(inout) :: lcmesh
    integer, intent(in) :: tileID
    integer, intent(in) :: panelID
    integer, intent(in) :: i, j, k
    integer, intent(in) :: NprcX, NprcY, NprcZ
    real(RP), intent(in) :: dom_xmin, dom_xmax
    real(RP), intent(in) :: dom_ymin, dom_ymax
    real(RP), intent(in) :: dom_Zmin, dom_zmax
    integer, intent(in) :: NeX, NeY, NeZ
    real(RP), intent(in) :: FZ(NeZ*NprcZ+1)
    
    class(elementbase3d), pointer :: elem
    real(RP) :: delx, dely
    real(RP) :: FZ_lc(NeZ+1)
 
    integer :: ii, jj, kk
    integer :: ke
    !-----------------------------------------------------------------------------
 
    elem => lcmesh%refElem3D
    
    lcmesh%tileID = tileid
    lcmesh%panelID = panelid
    !$acc update device(lcmesh%tileID, lcmesh%panelID)
    
    !--
    lcmesh%Ne   = nex * ney * nez
    lcmesh%Nv  = (nex + 1)*(ney + 1)*(nez + 1)
    lcmesh%NeS = 1
    lcmesh%NeE = lcmesh%Ne
    lcmesh%NeA = lcmesh%Ne + 2*(nex + ney)*nez + 2*nex*ney
    !$acc update device(lcmesh%Ne, lcmesh%Nv, lcmesh%NeS, lcmesh%NeE, lcmesh%NeA)
 
    lcmesh%NeX = nex
    lcmesh%NeY = ney
    lcmesh%NeZ = nez
    !$acc update device(lcmesh%NeX, lcmesh%NeY, lcmesh%NeZ)
 
    lcmesh%Ne2D  = nex * ney
    lcmesh%Ne2DA = nex * ney + 2*(nex + ney)
    !$acc update device(lcmesh%Ne2D, lcmesh%Ne2DA)
 
    !--
    delx = (dom_xmax - dom_xmin)/dble(nprcx)
    dely = (dom_ymax - dom_ymin)/dble(nprcy)
    fz_lc(:) = fz((k-1)*nez+1:k*nez+1)
    lcmesh%xmin = dom_xmin + (i-1)*delx
    lcmesh%xmax = dom_xmin +  i   *delx
    lcmesh%ymin = dom_ymin + (j-1)*dely
    lcmesh%ymax = dom_ymin +  j   *dely
    lcmesh%zmin = fz_lc(1)
    lcmesh%zmax = fz_lc(nez+1)
    !$acc update device(lcmesh%xmin, lcmesh%xmax, lcmesh%ymin, lcmesh%ymax, lcmesh%zmin, lcmesh%zmax)
    
    !-
    allocate( lcmesh%pos_ev(lcmesh%Nv,3) )
    allocate( lcmesh%EToV(lcmesh%Ne,elem%Nv) )
    allocate( lcmesh%EToE(lcmesh%Ne,elem%Nfaces) )
    allocate( lcmesh%EToF(lcmesh%Ne,elem%Nfaces) )
    allocate( lcmesh%BCType(lcmesh%refElem%Nfaces,lcmesh%Ne) )
    allocate( lcmesh%VMapM(elem%NfpTot, lcmesh%Ne) )
    allocate( lcmesh%VMapP(elem%NfpTot, lcmesh%Ne) )
    allocate( lcmesh%MapM(elem%NfpTot, lcmesh%Ne) )
    allocate( lcmesh%MapP(elem%NfpTot, lcmesh%Ne) )
    !$acc enter data create( lcmesh%pos_ev, lcmesh%EToV, lcmesh%EToE, lcmesh%EToF, lcmesh%BCType, &
    !$acc   lcmesh%VMapM, lcmesh%VMapP, lcmesh%MapM, lcmesh%MapP )
    
    allocate( lcmesh%EMap3Dto2D(lcmesh%Ne) )
    !$acc enter data create( lcmesh%EMap3Dto2D )
 
    lcmesh%BCType(:,:) = bctype_interior
    !$acc update device(lcmesh%BCType)
    
    !----
 
    call meshutil3d_gencubedomain( lcmesh%pos_ev, lcmesh%EToV,   & ! (out)
        lcmesh%NeX, lcmesh%xmin, lcmesh%xmax,                    & ! (in)
        lcmesh%NeY, lcmesh%ymin, lcmesh%ymax,                    & ! (in) 
        lcmesh%NeZ, lcmesh%zmin, lcmesh%zmax, fz=fz_lc           ) ! (in) 
    !$acc update device(lcmesh%pos_ev, lcmesh%EToV)
    
    !---
    call meshbase3d_setgeometricinfo( lcmesh, meshcubedom3d_coord_conv, meshcubedom3d_calc_normal )
 
    !---
    call meshutil3d_genconnectivity( lcmesh%EToE, lcmesh%EToF, & ! (out)
        lcmesh%EToV, lcmesh%Ne, elem%Nfaces )                    ! (in)
    !$acc update device(lcmesh%EToE, lcmesh%EToF)
    
    !---
    call meshutil3d_buildinteriormap( lcmesh%VmapM, lcmesh%VMapP, lcmesh%MapM, lcmesh%MapP,           & ! (out)
      lcmesh%pos_en, lcmesh%pos_ev, lcmesh%EToE, lcmesh%EtoF, lcmesh%EtoV,                            & ! (in)
      elem%Fmask_h, elem%Fmask_v, lcmesh%Ne, lcmesh%Nv, elem%Np, elem%Nfp_h, elem%Nfp_v, elem%NfpTot, & ! (in)
      elem%Nfaces_h, elem%Nfaces_v, elem%Nfaces )                                                       ! (in)
 
    call meshutil3d_genpatchboundarymap( lcmesh%VMapB, lcmesh%MapB, lcmesh%VMapP,                       & !(out)
      lcmesh%pos_en, lcmesh%xmin, lcmesh%xmax, lcmesh%ymin, lcmesh%ymax, lcmesh%zmin, lcmesh%zmax,      & ! (in)
      elem%Fmask_h, elem%Fmask_v, lcmesh%Ne, lcmesh%Nv, elem%Np, elem%Nfp_h, elem%Nfp_v, elem%NfpTot,   & ! (in)
      elem%Nfaces_h, elem%Nfaces_v, elem%Nfaces )                                                         ! (in)
    !$acc update device(lcmesh%VMapM, lcmesh%VMapP, lcmesh%MapM, lcmesh%MapP)
    !$acc enter data copyin(lcmesh%VMapB, lcmesh%MapB)
    
    !---
    !$omp parallel do collapse(2) private(ii,ke)
    do kk=1, lcmesh%NeZ
    do jj=1, lcmesh%NeY
    do ii=1, lcmesh%NeX
      ke = ii + (jj-1) * lcmesh%NeX + (kk-1) * lcmesh%NeX * lcmesh%NeY
      lcmesh%EMap3Dto2D(ke) = ii + (jj-1) * lcmesh%NeX
    end do
    end do
    end do
    !$acc update device(lcmesh%EMap3Dto2D)
 
    return
  end subroutine meshcubedom3d_setuplocaldom
 
!OCL SERIAL
  subroutine messhcubedom3d_assigndomid( this, &
    tileID_table, panelID_table,               &
    pi_table, pj_table, pk_table )
  
    use scale_meshutil_3d, only: &       
      meshutil3d_buildglobalmap    
    implicit none
 
    type(meshcubedom3d), target, intent(inout) :: this    
    integer, intent(out) :: tileID_table(this%LOCAL_MESH_NUM, this%PRC_NUM)
    integer, intent(out) :: panelID_table(this%LOCAL_MESH_NUM*this%PRC_NUM)
    integer, intent(out) :: pi_table(this%LOCAL_MESH_NUM*this%PRC_NUM)
    integer, intent(out) :: pj_table(this%LOCAL_MESH_NUM*this%PRC_NUM)
    integer, intent(out) :: pk_table(this%LOCAL_MESH_NUM*this%PRC_NUM)
    
    integer :: n
    integer :: p
    integer :: tileID
    integer :: is_lc, js_lc, ks_lc
    integer :: ilc_count, jlc_count, klc_count
    integer :: ilc, jlc, klc
 
    type(localmesh3d), pointer :: lcmesh
    !-----------------------------------------------------------------------------
    
    call meshutil3d_buildglobalmap( &
      panelid_table, pi_table, pj_table, pk_table,                                          & ! (out)
      this%tileID_globalMap, this%tileFaceID_globalMap, this%tilePanelID_globalMap,         & ! (out)
      this%LOCAL_MESH_NUM_global, 6, 8,                                                     & ! (in)
      this%isPeriodicX, this%isPeriodicY, this%isPeriodicZ,                                 & ! (in)
      this%NprcX, this%NprcY, this%NprcZ )                                                    ! (in)
    !$acc update device(this%tileID_globalMap, this%tileFaceID_globalMap, this%tilePanelID_globalMap)
 
    !----
 
    do p=1, this%PRC_NUM
    do n=1, this%LOCAL_MESH_NUM
      tileid = n + (p-1)*this%LOCAL_MESH_NUM
      lcmesh => this%lcmesh_list(n)
      !-
      tileid_table(n,p)                   = tileid
      this%tileID_global2localMap(tileid) = n
      this%PRCRank_globalMap(tileid)      = p - 1
 
      !-
      if ( this%PRCRank_globalMap(tileid) == lcmesh%PRC_myrank ) then
        if (n==1) then
          is_lc = pi_table(tileid); ilc_count = 1
          js_lc = pj_table(tileid); jlc_count = 1
          ks_lc = pk_table(tileid); klc_count = 1
        end if
        if(is_lc < pi_table(tileid)) ilc_count = ilc_count + 1
        if(js_lc < pj_table(tileid)) jlc_count = jlc_count + 1
        if(ks_lc < pk_table(tileid)) klc_count = klc_count + 1
      end if 
    end do
    end do
    !$acc update device(this%tileID_global2localMap, this%PRCRank_globalMap)
 
    allocate( this%rcdomIJK2LCMeshID(ilc_count,jlc_count,klc_count) )
    do klc=1, klc_count
    do jlc=1, jlc_count
    do ilc=1, ilc_count
      this%rcdomIJK2LCMeshID(ilc,jlc,klc) = ilc + (jlc - 1)*ilc_count + (klc - 1)*ilc_count*jlc_count
    end do
    end do
    end do
    !$acc enter data copyin(this%rcdomIJK2LCMeshID)
 
    return
  end subroutine messhcubedom3d_assigndomid
  
!OCL SERIAL
  subroutine meshcubedom3d_coord_conv( x, y, z, xX, xY, xZ, yX, yY, yZ, zX, zY, zZ, &
    vx, vy, vz, elem )
 
    implicit none
 
    type(elementbase3d), intent(in) :: elem
    real(rp), intent(out) :: x(elem%np), y(elem%np), z(elem%np)
    real(rp), intent(out) :: xx(elem%np), xy(elem%np), xz(elem%np)
    real(rp), intent(out) :: yx(elem%np), yy(elem%np), yz(elem%np)
    real(rp), intent(out) :: zx(elem%np), zy(elem%np), zz(elem%np)
    real(rp), intent(in) :: vx(elem%nv), vy(elem%nv), vz(elem%nv)
 
    !-------------------------------------------------
    
    x(:) = vx(1) + 0.5_rp*(elem%x1(:) + 1.0_rp)*(vx(2) - vx(1))
    y(:) = vy(1) + 0.5_rp*(elem%x2(:) + 1.0_rp)*(vy(3) - vy(1))
    z(:) = vz(1) + 0.5_rp*(elem%x3(:) + 1.0_rp)*(vz(5) - vz(1))
 
    xx(:) = 0.5_rp*(vx(2) - vx(1)) !matmul(refElem%Dx1,mesh%x1(:,n))
    xy(:) = 0.0_rp                 !matmul(refElem%Dx2,mesh%x1(:,n))
    xz(:) = 0.0_rp                 !matmul(refElem%Dx3,mesh%x1(:,n))
    yx(:) = 0.0_rp                 !matmul(refElem%Dx1,mesh%x2(:,n))
    yy(:) = 0.5_rp*(vy(3) - vy(1)) !matmul(refElem%Dx2,mesh%x2(:,n))
    yz(:) = 0.0_rp                 !matmul(refElem%Dx3,mesh%x2(:,n))
    zx(:) = 0.0_rp                 !matmul(refElem%Dx1,mesh%x3(:,n))
    zy(:) = 0.0_rp                 !matmul(refElem%Dx2,mesh%x3(:,n))
    zz(:) = 0.5_rp*(vz(5) - vz(1)) !matmul(refElem%Dx3,mesh%x3(:,n))
 
    return
  end subroutine meshcubedom3d_coord_conv  
 
!OCL SERIAL
  subroutine meshcubedom3d_calc_normal( normal_fn, &
    Escale_f, fid_h, fid_v, elem )
 
    implicit none
 
    type(elementbase3d), intent(in) :: elem
    real(rp), intent(out) :: normal_fn(elem%nfptot,3)
    integer, intent(in) :: fid_h(elem%nfp_h,elem%nfaces_h)
    integer, intent(in) :: fid_v(elem%nfp_v,elem%nfaces_v)        
    real(rp), intent(in) :: escale_f(elem%nfptot,3,3)
 
    integer :: d
    !-------------------------------------------------
 
    do d=1, 3
      normal_fn(fid_h(:,1),d) = - escale_f(fid_h(:,1),2,d)
      normal_fn(fid_h(:,2),d) = + escale_f(fid_h(:,2),1,d)
      normal_fn(fid_h(:,3),d) = + escale_f(fid_h(:,3),2,d)
      normal_fn(fid_h(:,4),d) = - escale_f(fid_h(:,4),1,d)
 
      normal_fn(fid_v(:,1),d) = - escale_f(fid_v(:,1),3,d)
      normal_fn(fid_v(:,2),d) = + escale_f(fid_v(:,2),3,d)    
    end do
 
    return
  end subroutine meshcubedom3d_calc_normal 
 
end module scale_mesh_cubedom3d