scale_meshutil_3d.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Mesh / utility for 3D mesh
!!
!! @par Description
!!          A module to provide subroutines for generating 3D mesh 
!!
!! @author Yuta Kawai, Team SCALE
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_meshutil_3d
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
 
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !  
  public :: meshutil3d_gencubedomain
  public :: meshutil3d_genconnectivity
  public :: meshutil3d_buildinteriormap
  public :: meshutil3d_genpatchboundarymap
  public :: meshutil3d_buildglobalmap
 
contains
 
!OCL SERIAL
  subroutine meshutil3d_gencubedomain( pos_v, EToV,        &
    Ke_x, xmin, xmax, Ke_y, ymin, ymax, Ke_z, zmin, zmax,  &
    Fz )
 
    implicit none
 
    integer, intent(in) :: ke_x
    real(rp), intent(in) :: xmin, xmax
    integer, intent(in) :: ke_y
    real(rp), intent(in) :: ymin, ymax
    integer, intent(in) :: ke_z
    real(rp), intent(in) :: zmin, zmax
    real(rp), intent(in), optional :: fz(ke_z+1)
    
    real(rp), intent(out) :: pos_v((ke_x+1)*(ke_y+1)*(ke_z+1),3)
    integer, intent(out) :: etov(ke_x*ke_y*ke_z,8)
 
    integer :: i, j, k
    integer :: n
    integer :: nvx, nvy, nvz
    !-----------------------------------------------------------------------------
 
    nvx = ke_x + 1
    nvy = ke_y + 1
    nvz = ke_z + 1
 
    !$omp parallel private(i,j,k,n)
 
    !$omp do collapse(2) 
    do k=1, nvz
    do j=1, nvy
    do i=1, nvx
      n = i + ( j - 1 ) * nvx + ( k - 1 ) *nvx*nvy
      pos_v(n,1) = ( xmax - xmin ) * dble(i - 1)/dble(nvx - 1) + xmin
      pos_v(n,2) = ( ymax - ymin ) * dble(j - 1)/dble(nvy - 1) + ymin
      if ( present(fz) ) then
        pos_v(n,3) = fz(k)
      else
        pos_v(n,3) = (zmax - zmin)*dble(k - 1)/dble(nvz - 1) + zmin
      end if
    end do
    end do
    end do
 
    !$omp do collapse(2)
    do k=1, ke_z
    do j=1, ke_y
    do i=1, ke_x
      n = i + (j-1)*ke_x + (k-1)*ke_x*ke_y
      etov(n,1:4) = (k-1)*nvx*nvy + (j-1)*nvx + &
            & (/ i, i+1,              &
            &    i + nvx, i+1 + nvx /)
      etov(n,5:8) = k    *nvx*nvy + (j-1)*nvx + &
            & (/ i, i+1,              &
            &    i + nvx, i+1 + nvx /)              
    end do
    end do
    end do
 
    !$omp end parallel
    
    !---
  !!$
  !!$    write(*,*) "-- vx, vy --"
  !!$    do j=1, NvY
  !!$    do i=1, NvX
  !!$       n = i + (j-1)*NvY
  !!$       write(*,*) i, j, " :", mesh%vx(n), mesh%vy(n)
  !!$    end do
  !!$    end do
  !!$
  !!$    write(*,*) "-- EToV --"
  !!$    do j=1, Ke_y
  !!$    do i=1, Ke_x
  !!$       n = i + (j-1)*Ke_y
  !!$       write(*,*) i, j, " :", mesh%EToV(n,:)
  !!$    end do
  !!$    end do
 
    return
  end subroutine meshutil3d_gencubedomain
 
!OCL SERIAL
  subroutine meshutil3d_genconnectivity( EToE, EToF, &
    EToV, Ne, Nfaces )
    
    use scale_quicksort, only: quicksort_exec_with_idx
    implicit none
 
    integer, intent(in) :: ne
    integer, intent(in) :: nfaces
    integer, intent(out) :: etoe(ne, nfaces)
    integer, intent(out) :: etof(ne, nfaces)
    integer, intent(in) :: etov(ne,8)
 
    integer :: nodes(ne*nfaces,4)
    integer(kind=8) :: face_ids(ne*nfaces)
    integer :: ke
    integer :: f
    integer :: n
    integer :: nnodes
    integer :: nnodes_row
    integer :: matchl(2,3), matchr(2,3)
 
    real(rp) :: etoe_1d(ne*nfaces)
    real(rp) :: etof_1d(ne*nfaces)
 
    integer :: spnodetonode(3,ne*nfaces)
    integer :: spnodetonoderowtmp(3,ne*nfaces)
    integer :: sort_indx(ne*nfaces)
    integer(kind=8) :: sorted_faceid(ne*nfaces)
 
    real(rp) :: vtmp(4)
    
    !-----------------------------------------------------------------------------
 
    nnodes = maxval( etov )
    nnodes_row = size(nodes,1)
    
    !---------
    do ke=1, ne
       nodes(ke     ,:) = etov(ke,(/ 1, 2, 6, 5 /))
       nodes(ke+ne  ,:) = etov(ke,(/ 2, 4, 8, 6 /))
       nodes(ke+2*ne,:) = etov(ke,(/ 4, 3, 7, 8 /))
       nodes(ke+3*ne,:) = etov(ke,(/ 3, 1, 5, 7 /))
       nodes(ke+4*ne,:) = etov(ke,(/ 1, 3, 4, 2 /))       
       nodes(ke+5*ne,:) = etov(ke,(/ 5, 6, 8, 7 /))       
    end do
   
    ! Sort
    do n=1, nnodes_row
      vtmp(:) = nodes(n,:)
      call bubblesort( vtmp )
      nodes(n,:) = vtmp(:)
      ! write(*,*) n, ":", nodes(n,:)
    end do
    nodes = nodes - 1
    !---------
 
    do ke=1, ne
       etoe(ke,:) = ke
       etof(ke,:) = (/ 1, 2, 3, 4, 5, 6 /)
    end do
 
    face_ids(:) = nodes(:,1) * nnodes**3 + nodes(:,2) * nnodes**2 &
                + nodes(:,3) * nnodes    + nodes(:,4) + 1
 
    do f=1, nfaces
    do ke=1, ne
      n = ke + (f-1)*ne
      spnodetonode(:,n) = (/ n, etoe(ke,f), etof(ke,f) /)
      sorted_faceid(n)  = face_ids(n)
      sort_indx(n)      = n
       ! write(*,*) "face_id, n, EToE, EToF:", spNodeToNode(:,n)
    end do
    end do
    
    !- sort row
    call quicksort_exec_with_idx( ne*nfaces, sorted_faceid, sort_indx )
    spnodetonoderowtmp(:,:) = spnodetonode(:,:)
 
    do n=1, nnodes_row
      spnodetonode(:,n) = spnodetonoderowtmp(:,sort_indx(n))
      !  write(*,'(a,4i10)') "(sorted) face_id, n, EToE, EToF:", spNodeToNode(:,n)
    end do
 
    etoe_1d(:) = -1
    etof_1d(:) = -1
    do n=1, nnodes_row-1
       if ( sorted_faceid(n) - sorted_faceid(n+1) == 0 ) then
          matchl(:,:) = transpose( spnodetonode(:,(/ n, n+1 /)) )
          matchr(:,:) = transpose( spnodetonode(:,(/ n+1, n /)) )
 
          etoe_1d(matchl(:,1)) = matchr(:,2)
          etof_1d(matchl(:,1)) = matchr(:,3)
       end if
    end do
 
    do f=1, nfaces
    do ke=1, ne
       n = ke + (f-1)*ne
       if ( etoe_1d(n) /= -1 ) then
          etoe(ke,f) = etoe_1d(n)
          etof(ke,f) = etof_1d(n)
       end if
    end do
    end do
 
    !-------------------
 
!    write(*,*) "EToE------"
!    do k=1, mesh%Ne
!      write(*,*) "k=", k, ":", mesh%EToE(k,:)
!    end do
!    write(*,*) "EToF------"
!    do k=1, mesh%Ne
!       write(*,*) "k=", k, ":", mesh%EToF(k,:)
!    end do
 
    return
  end subroutine meshutil3d_genconnectivity
 
!OCL SERIAL
  subroutine bubblesort( array )
    implicit none
    real(rp), intent(inout) :: array(:)
    
    integer :: i, j, n
    real(rp) :: t
    !--------------------
 
    n = size(array)
    do i=1, n-1
      do j=i+1,n
        t = array(i)      
        if(t > array(j)) then
          array(i) = array(j)
          array(j) = t
        end if
      end do    
    end do
 
    return
  end subroutine bubblesort
 
!OCL SERIAL
  subroutine meshutil3d_buildinteriormap( VMapM, VMapP, MapM, MapP, &
    pos_en, pos_ev, EtoE, EtoF, EtoV, Fmask_h, Fmask_v,             &
    Ne, Nv, Np, Nfp_h, Nfp_v, NfpTot, Nfaces_h, Nfaces_v, Nfaces)
 
    implicit none
 
    integer, intent(in) :: ne
    integer, intent(in) :: nv
    integer, intent(in) :: np
    integer, intent(in) :: nfp_h
    integer, intent(in) :: nfp_v
    integer, intent(in) :: nfptot
    integer, intent(in) :: nfaces_h
    integer, intent(in) :: nfaces_v
    integer, intent(in) :: nfaces
 
    integer, intent(out) :: vmapm(nfptot,ne)
    integer, intent(out) :: vmapp(nfptot,ne)
    integer, intent(out) :: mapm(nfptot,ne)
    integer, intent(out) :: mapp(nfptot,ne)
 
    real(rp), intent(in) :: pos_en(np,ne,3)
    real(rp), intent(in) :: pos_ev(nv,3)
    integer, intent(in) :: etoe(ne,nfaces)
    integer, intent(in) :: etof(ne,nfaces)
    integer, intent(in) :: etov(ne,nv)
    integer, intent(in) :: fmask_h(nfp_h,nfaces_h)
    integer, intent(in) :: fmask_v(nfp_v,nfaces_v)
 
    integer :: ke, ke1, ke2
    integer :: f, f1, f2
    integer :: p
    integer :: n
    integer :: i
    integer :: idp, idm
    integer :: v1, v2
    integer :: nodeids(np,ne)
    real(rp) :: r_h(nfp_h,nfp_h,3,2)
    real(rp) :: r_v(nfp_v,nfp_v,3,2)
    real(rp) :: dist_h(nfp_h,nfp_h)
    real(rp) :: dist_v(nfp_v,nfp_v)
    real(rp) :: x(np*ne), y(np*ne), z(np*ne)
 
    integer :: vmapm_h(nfp_h,nfaces_h,ne)
    integer :: vmapp_h(nfp_h,nfaces_h,ne)
    integer :: vmapm_v(nfp_v,nfaces_v,ne)
    integer :: vmapp_v(nfp_v,nfaces_v,ne)
    integer :: mapm_h(nfp_h,nfaces_h,ne)
    integer :: mapp_h(nfp_h,nfaces_h,ne)
    integer :: mapm_v(nfp_v,nfaces_v,ne)
    integer :: mapp_v(nfp_v,nfaces_v,ne)
 
    integer :: mindist_indx(1)
    !-----------------------------------------------------------------------------
 
    !$omp parallel private(ke,f,p,n)
 
    !$omp do
    do ke=1, ne
    do p=1, np
      n = p + (ke-1)*np
      nodeids(p,ke) = n
      x(n) = pos_en(p,ke,1)
      y(n) = pos_en(p,ke,2)
      z(n) = pos_en(p,ke,3)
    end do
    end do
    !$omp end do
 
    !$omp do
    do ke=1, ne
      do f=1, nfaces_h
      do p=1, nfp_h
        n = p + (f-1)*nfp_h + (ke-1)*nfptot
        mapm_h(p,f,ke) = n
        mapp_h(p,f,ke) = n
        vmapm_h(p,f,ke) = nodeids(fmask_h(p,f),ke)
      end do
      end do
      do f=1, nfaces_v
      do p=1, nfp_v
        n = p + nfaces_h*nfp_h + (f-1)*nfp_v + (ke-1)*nfptot
        mapm_v(p,f,ke) = n
        mapp_v(p,f,ke) = n
        vmapm_v(p,f,ke) = nodeids(fmask_v(p,f),ke)
      end do
      end do    
    end do
    !$omp end do
 
    !$omp workshare
    vmapp_h(:,:,:) = -1
    vmapp_v(:,:,:) = -1
    !$omp end workshare
    !$omp end parallel
 
    !$omp parallel private( &
    !$omp ke1, f1, ke2, f2, v1, v2,                         &
    !$omp r_h, r_v, dist_h, dist_v, mindist_indx, idP, idM  )
 
    !$omp do
    do ke1=1, ne
    do f1=1, nfaces_h
      ke2 = etoe(ke1,f1); f2 = etof(ke1,f1)
 
      v1 = etov(ke1,f1); v2 = etov(ke1,1+mod(f1,nfaces_h))
 
      r_h(:,:,1,1) = spread( x(vmapm_h(:,f1,ke1)), 2, nfp_h )
      r_h(:,:,1,2) = spread( x(vmapm_h(:,f2,ke2)), 1, nfp_h )
      r_h(:,:,2,1) = spread( y(vmapm_h(:,f1,ke1)), 2, nfp_h )
      r_h(:,:,2,2) = spread( y(vmapm_h(:,f2,ke2)), 1, nfp_h )
      r_h(:,:,3,1) = spread( z(vmapm_h(:,f1,ke1)), 2, nfp_h )
      r_h(:,:,3,2) = spread( z(vmapm_h(:,f2,ke2)), 1, nfp_h )
      
      dist_h(:,:) =   (r_h(:,:,1,1) - r_h(:,:,1,2))**2  &
                    + (r_h(:,:,2,1) - r_h(:,:,2,2))**2  &
                    + (r_h(:,:,3,1) - r_h(:,:,3,2))**2
      do idm=1, nfp_h
        mindist_indx(:) = minloc(dist_h(idm,:))
        idp = mindist_indx(1)
        vmapp_h(idm,f1,ke1) = vmapm_h(idp,f2,ke2)
        mapp_h(idm,f1,ke1) = idp + (f2-1)*nfp_h + (ke2-1)*nfptot
      end do
    end do
    end do
    !omp end do
 
    !$omp do
    do ke1=1, ne
    do f1=1, nfaces_v
      ke2 = etoe(ke1,nfaces_h+f1); f2 = etof(ke1,nfaces_h+f1) - nfaces_h
      
      v1 = etov(ke1,1); v2 = etov(ke1,nfaces_h+1)
 
      r_v(:,:,1,1) = spread( x(vmapm_v(:,f1,ke1)), 2, nfp_v )
      r_v(:,:,1,2) = spread( x(vmapm_v(:,f2,ke2)), 1, nfp_v )
      r_v(:,:,2,1) = spread( y(vmapm_v(:,f1,ke1)), 2, nfp_v )
      r_v(:,:,2,2) = spread( y(vmapm_v(:,f2,ke2)), 1, nfp_v )
      r_v(:,:,3,1) = spread( z(vmapm_v(:,f1,ke1)), 2, nfp_v )
      r_v(:,:,3,2) = spread( z(vmapm_v(:,f2,ke2)), 1, nfp_v )
      
      dist_v(:,:) =   (r_v(:,:,1,1) - r_v(:,:,1,2))**2  &
                    + (r_v(:,:,2,1) - r_v(:,:,2,2))**2  &
                    + (r_v(:,:,3,1) - r_v(:,:,3,2))**2
      do idm=1, nfp_v
        mindist_indx(:) = minloc(dist_v(idm,:))
        idp = mindist_indx(1)
        vmapp_v(idm,f1,ke1) = vmapm_v(idp,f2,ke2)
        mapp_v(idm,f1,ke1) = idp + nfaces_h*nfp_h + (f2-1)*nfp_v + (ke2-1)*nfptot
      end do
    end do
    end do
    !omp end do
    !$omp end parallel
 
    !$omp parallel do private(ke,f,n,i)
    do ke=1, ne
      do f=1, nfaces_h
      do n=1, nfp_h
        i = n + (f-1)*nfp_h
        vmapm(i,ke) = vmapm_h(n,f,ke)
        mapm(i,ke) = mapm_h(n,f,ke)
        vmapp(i,ke) = vmapp_h(n,f,ke)
        mapp(i,ke) = mapp_h(n,f,ke)
      end do
      end do
      do f=1, nfaces_v
      do n=1, nfp_v
        i = n + nfaces_h*nfp_h + (f-1)*nfp_v
        vmapm(i,ke) = vmapm_v(n,f,ke)
        mapm(i,ke) = mapm_v(n,f,ke)
        vmapp(i,ke) = vmapp_v(n,f,ke)
        mapp(i,ke) = mapp_v(n,f,ke)
      end do
      end do
    end do
 
    !-----
  !    mapB_counter = 0
  !    do k=1,mesh%Ne
  !    do f=1,elem%Nfaces
  !    do p=1,elem%Nfp
  !      n = p + (f-1)*elem%Nfp + (k-1)*elem%Nfp*elem%Nfaces
  !      if (mesh%VMapM(p,f,k) == mesh%VMapP(p,f,k)) then
  !        mapB_counter = mapB_counter + 1
  !        mapB_tmp(mapB_counter) = n
  !
  !        vmapB_tmp(mapB_counter) = mesh%VMapM(p,f,k)
  !        mesh%VMapP(p,f,k) = elem%Np*mesh%NeE + mapB_counter
  !      end if
  !    end do
  !    end do
  !    end do
  !
  !    allocate( mesh%mapB(mapB_counter) )
  !    allocate( mesh%vmapB(mapB_counter) )
  !    mesh%mapB(:) = mapB_tmp(1:mapB_counter)
  !    mesh%vmapB(:) = vmapB_tmp(1:mapB_counter)
 
    !-------
  !    write(*,*) "Build MapInfo: "
  !    do k=mesh%NeS,mesh%NeE
  !       write(*,*) "k=", k, "---"
  !       write(*,*) " VMapM:", mesh%VMapM(:,:,k)
  !       write(*,*) " VMapP:", mesh%VMapP(:,:,k)
  !    end do
 
  !    write(*,*) "mapB:", mesh%mapB(:)
  !    write(*,*) "vmapB:", mesh%vmapB(:)
 
    return
  end subroutine meshutil3d_buildinteriormap
 
!OCL SERIAL
  subroutine meshutil3d_genpatchboundarymap(  VMapB, MapB, VMapP,                  &
    pos_en, xmin, xmax, ymin, ymax, zmin, zmax,                                    &
    Fmask_h, Fmask_v, Ne, Nv, Np, Nfp_h, Nfp_v, NfpTot, Nfaces_h, Nfaces_v, Nfaces )
 
    implicit none
 
    integer, intent(in) :: ne
    integer, intent(in) :: nv
    integer, intent(in) :: np
    integer, intent(in) :: nfp_h
    integer, intent(in) :: nfp_v
    integer, intent(in) :: nfptot
    integer, intent(in) :: nfaces_h
    integer, intent(in) :: nfaces_v
    integer, intent(in) :: nfaces
 
    integer, intent(inout), allocatable :: vmapb(:)
    integer, intent(inout), allocatable :: mapb(:)
    integer, intent(inout) :: vmapp(nfptot,ne)
 
    real(rp), intent(in) :: pos_en(np,ne,3)
    real(rp), intent(in) :: xmin, xmax
    real(rp), intent(in) :: ymin, ymax
    real(rp), intent(in) :: zmin, zmax
    integer, intent(in) :: fmask_h(nfp_h,nfaces_h)
    integer, intent(in) :: fmask_v(nfp_v,nfaces_v)
 
    integer :: ke, n
    integer :: b
    integer :: f
    integer :: i, j
    real(rp) :: x, y, z
 
    real(rp), parameter :: node_tol = 1.0e-12_rp
 
    integer :: elemids_h(nfp_h*ne,nfaces_h)
    real(rp) :: ordinfo_h(nfp_h*ne,nfaces_h)
    integer :: faceids_h(nfp_h*ne,nfaces_h)
    integer :: counterb_h(nfaces_h)
 
    integer :: elemids_v(nfp_v*ne,nfaces_v)
    real(rp) :: ordinfo_v(nfp_v*ne,nfaces_v)
    integer :: faceids_v(nfp_v*ne,nfaces_v)
    integer :: counterb_v(nfaces_v)
 
    integer :: mapb_counter
    real(rp) :: rdomx, rdomy, rdomz
    !-----------------------------------------------------------------------------
 
    counterb_h(:) = 0
    counterb_v(:) = 0
 
    rdomx = 1.0_rp/(xmax - xmin)
    rdomy = 1.0_rp/(ymax - ymin)
    rdomz = 1.0_rp/abs(zmax - zmin)
    
    do ke=1, ne
      do f=1, nfaces_h
        x = sum(pos_en(fmask_h(:,f),ke,1)) / dble(nfp_h)
        y = sum(pos_en(fmask_h(:,f),ke,2)) / dble(nfp_h)
 
        call eval_domain_boundary( &
          elemids_h, ordinfo_h, faceids_h, counterb_h, & ! (inout)
          1, y, ymin, x, ke, f, rdomy                  ) ! (in)
        call eval_domain_boundary( &
          elemids_h, ordinfo_h, faceids_h, counterb_h, & ! (inout)
          2, x, xmax, y, ke, f, rdomx                  ) ! (in)    
        call eval_domain_boundary( &
          elemids_h, ordinfo_h, faceids_h, counterb_h, & ! (inout)
          3, y, ymax, x, ke, f, rdomy                  ) ! (in)
        call eval_domain_boundary( &
          elemids_h, ordinfo_h, faceids_h, counterb_h, & ! (inout)
          4, x, xmin, y, ke, f, rdomx                  ) ! (in)
      end do
      do f=1, nfaces_v
        x = sum(pos_en(fmask_v(:,f),ke,1)) / dble(nfp_v)
        z = sum(pos_en(fmask_v(:,f),ke,3)) / dble(nfp_v)
 
        call eval_domain_boundary( &
          elemids_v, ordinfo_v, faceids_v, counterb_v, & ! (inout)
          1, z, zmin, x, ke, f, rdomz                  ) ! (in)
        call eval_domain_boundary( &
          elemids_v, ordinfo_v, faceids_v, counterb_v, & ! (inout)
          2, z, zmax, x, ke, f, rdomz                  ) ! (in)  
      end do    
    end do
 
 
    allocate( mapb(sum(counterb_h(:)*nfp_h)+sum(counterb_v(:)*nfp_v)) )
    allocate( vmapb(size(mapb)) )
 
    mapb_counter = 1
    do b = 1, nfaces_h
      ! write(*,*) "LocalMesh boundary ID:", b
      ! write(*,*) counterB_h(b)
      ! write(*,*) ordInfo_h(1:counterB_h(1),b)
      ! write(*,*) elemIds_h(1:counterB_h(1),b)
      ! write(*,*) faceIds_h(1:counterB_h(1),b)
 
      do i=1, counterb_h(b)
        ke = elemids_h(i,b); f = faceids_h(i,b)
        do j=1, nfp_h
          n = j + (f-1)*nfp_h
          vmapp(n,ke) = np*ne + mapb_counter
          vmapb(mapb_counter) = fmask_h(j,f) + (ke-1)*np
          mapb_counter = mapb_counter + 1
        end do
      end do
    end do
 
    do b = 1, nfaces_v
      do i=1, counterb_v(b)
        ke = elemids_v(i,b); f = faceids_v(i,b)
        do j=1, nfp_v
          n = j + nfp_h*nfaces_h + nfp_v*(f-1)
          vmapp(n,ke) = np*ne + mapb_counter
          vmapb(mapb_counter) = fmask_v(j,f) + (ke-1)*np
          mapb_counter = mapb_counter + 1
        end do
      end do
    end do
 
    ! write(*,*) "VMapP:-----"
    ! do b=1, 4
    !   do i=1, counterB(b)
    !     k = elemIds(i,b); f = faceIDs(i,b)
    !     write(*,*) "bid=", b, ":", mesh%VmapP(:,f,k)
    !    end do
    ! end do
    ! write(*,*) "-----"
    ! write(*,*) "VMapB:", mesh%VmapB(:)
    ! write(*,*) "NeA=", mesh%NeA
 
    return
 
  contains
!OCL SERIAL
    subroutine eval_domain_boundary( &
        elemIDs, ordInfo, faceIDs, counterB,              &
        domb_id, r, rbc, ord_info, k_, f_, normalized_fac )
      implicit none
 
      integer, intent(inout) :: elemids(:,:)
      real(rp), intent(inout) :: ordinfo(:,:)
      integer, intent(inout) :: counterb(:)
      integer, intent(inout) :: faceids(:,:)
      integer, intent(in) :: domb_id
      real(rp), intent(in) :: r
      real(rp), intent(in) :: rbc
      real(rp), intent(in) :: ord_info
      integer, intent(in) :: k_, f_
      real(rp), intent(in) :: normalized_fac
      !-------------------------------------------------------------
 
      if ( abs(r - rbc)*normalized_fac < node_tol ) then
        counterb(domb_id) = counterb(domb_id) + 1
        ordinfo(counterb(domb_id),domb_id) = ord_info
        elemids(counterb(domb_id),domb_id) = k_
        faceids(counterb(domb_id),domb_id) = f_
      end if
 
      return
    end subroutine eval_domain_boundary
  end subroutine meshutil3d_genpatchboundarymap
 
!OCL SERIAL
  subroutine meshutil3d_buildglobalmap( &
    panelID_table, pi_table, pj_table, pk_table,    &
    tileID_map, tileFaceID_map, tilePanelID_map,    &
    Ntile, NtileFace, NtileVertex,                  &
    isPeriodicX, isPeriodicY, isPeriodicZ,          &
    Ne_x, Ne_y, Ne_z )
    
    ! use scale_prc, only: PRC_isMaster
    implicit none
 
    integer, intent(in) :: ntile
    integer, intent(in) :: ntileface
    integer, intent(in) :: ntilevertex
    integer, intent(out) :: panelid_table(ntile)
    integer, intent(out) :: pi_table(ntile)
    integer, intent(out) :: pj_table(ntile)
    integer, intent(out) :: pk_table(ntile)
    integer, intent(out) :: tileid_map(ntileface,ntile)
    integer, intent(out) :: tilefaceid_map(ntileface,ntile)
    integer, intent(out) :: tilepanelid_map(ntileface,ntile)
    logical, intent(in) :: isperiodicx
    logical, intent(in) :: isperiodicy
    logical, intent(in) :: isperiodicz
    integer, intent(in) :: ne_x
    integer, intent(in) :: ne_y
    integer, intent(in) :: ne_z
 
    integer :: ntileperpanel
    integer :: nv_x, nv_y, nv_z
    integer, allocatable :: nodesid_3d(:,:,:)
    integer, allocatable :: etov(:,:)
    integer, allocatable :: etoe(:,:)
    integer, allocatable :: etof(:,:)
    integer :: i, j, k, f
    integer :: tileid, tileid_r
    integer :: counter
 
    !-----------------------------------------------------------------------------
 
    ntileperpanel = ntile/1
    
    nv_x = ne_x + 1
    nv_y = ne_y + 1
    nv_z = ne_z + 1
 
    allocate( nodesid_3d(nv_x,nv_y,nv_z) )
    allocate( etov(ntile,ntilevertex), etoe(ntile,ntileface), etof(ntile,ntileface) )
 
    counter = 0
    do k = 1, nv_z
    do j = 1, nv_y
    do i = 1, nv_x
      counter = counter + 1
      nodesid_3d(i,j,k) = counter
    end do
    end do
    end do
 
 
    !----
 
    tileid = 0
    do k = 1, ne_z
    do j = 1, ne_y
    do i = 1, ne_x
      tileid = tileid + 1
      panelid_table(tileid) = 1
      pi_table(tileid) = i; pj_table(tileid) = j; pk_table(tileid) = k
      etov(tileid,:) = (/ nodesid_3d(i,j  ,k  ), nodesid_3d(i+1,j  ,k  ),   &
                          nodesid_3d(i,j+1,k  ), nodesid_3d(i+1,j+1,k  ),   &
                          nodesid_3d(i,j  ,k+1), nodesid_3d(i+1,j  ,k+1),   &
                          nodesid_3d(i,j+1,k+1), nodesid_3d(i+1,j+1,k+1) /)
    end do
    end do
    end do
 
    call meshutil3d_genconnectivity( etoe, etof, &
      etov, ntile, ntileface )
    tileid_map(:,:) = transpose(etoe)
    tilefaceid_map(:,:) = transpose(etof)
 
    do tileid=1, ntile
    do f=1, ntileface
      tileid_r = tileid_map(f,tileid)
      tilepanelid_map(f,tileid) = panelid_table(tileid_r)
    end do
    end do
 
    if (isperiodicx) then
      do tileid=1, ntile
        if (pi_table(tileid) == 1 .and. tilefaceid_map(4,tileid) == 4) then
          tileid_map(4,tileid) = ne_x + (pj_table(tileid) - 1)*ne_x + (pk_table(tileid) - 1)*ne_x*ne_y
          tilefaceid_map(4,tileid) = 2
        end if
        if (pi_table(tileid) == ne_x .and. tilefaceid_map(2,tileid) == 2) then
          tileid_map(2,tileid) = 1 + (pj_table(tileid) - 1)*ne_x + (pk_table(tileid) - 1)*ne_x*ne_y
          tilefaceid_map(2,tileid) = 4
        end if
      end do
    end if
 
    if (isperiodicy) then
      do tileid=1, ntile
        if (pj_table(tileid) == 1 .and. tilefaceid_map(1,tileid) == 1) then
          tileid_map(1,tileid) = pi_table(tileid) + (ne_y - 1)*ne_x + (pk_table(tileid) - 1)*ne_x*ne_y
          tilefaceid_map(1,tileid) = 3
        end if
        if (pj_table(tileid) == ne_y .and. tilefaceid_map(3,tileid) == 3) then
          tileid_map(3,tileid) = pi_table(tileid) + (pk_table(tileid) - 1)*ne_x*ne_y
          tilefaceid_map(3,tileid) = 1
        end if
      end do
    end if
 
    if (isperiodicz) then
      do tileid=1, ntile
        if (pk_table(tileid) == 1 .and. tilefaceid_map(5,tileid) == 5) then
          tileid_map(5,tileid) = pi_table(tileid) + (pj_table(tileid) - 1)*ne_x  + (ne_z - 1)*ne_x*ne_y
          tilefaceid_map(5,tileid) = 6
        end if
        if (pk_table(tileid) == ne_z .and. tilefaceid_map(6,tileid) == 6) then
          tileid_map(6,tileid) = pi_table(tileid) + (pj_table(tileid) - 1)*ne_x
          tilefaceid_map(6,tileid) = 5
        end if
      end do
    end if
 
    return
    
    !--
    ! if (PRC_isMaster) then
    !   write(*,*) "TotTile", Ntile
    !   do tileID = 1, Ntile
    !     write(*,*) "tileID:", tileID, ", EtoV:", EtoV(tileID,:)
    !   end do
    !   write(*,*) "-----------"
    !   do tileID = 1, Ntile
    !     write(*,*) "tileID:", tileID, ", EtoE:", EtoE(tileID,:)
    !   end do
    !   write(*,*) "-----------"
    !   do tileID = 1, Ntile
    !     write(*,*) "tileID:", tileID, ", EtoF:", EtoF(tileID,:)
    !   end do
    ! end if
    
  end subroutine meshutil3d_buildglobalmap
 
end module scale_meshutil_3d