scale_meshutil_2d.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Mesh / utility for 2D mesh
!!
!! @par Description
!!          A module useful for generating 2D mesh 
!!
!! @author Yuta Kawai, Team SCALE
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_meshutil_2d
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
 
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !  
  public :: meshutil2d_genrectdomain
  public :: meshutil2d_genconnectivity
  public :: meshutil2d_buildinteriormap
  public :: meshutil2d_genpatchboundarymap
  public :: meshutil2d_genpatchboundarymap_wide
  public :: meshutil2d_buildglobalmap
 
contains
!OCL SERIAL
  subroutine meshutil2d_genrectdomain( pos_v, EToV, &
    Ke_x, xmin, xmax, Ke_y, ymin, ymax )
 
    implicit none
 
    integer, intent(in) :: ke_x
    real(rp), intent(in) :: xmin, xmax
    integer, intent(in) :: ke_y
    real(rp), intent(in) :: ymin, ymax
 
    real(rp), intent(out) :: pos_v((ke_x+1)*(ke_y+1),2)
    integer, intent(out) :: etov(ke_x*ke_y,4)
 
    integer :: i, j
    integer :: n
    integer :: nvx, nvy
    !-----------------------------------------------------------------------------
 
    nvx = ke_x + 1
    nvy = ke_y + 1
 
    do j=1, nvy
    do i=1, nvx
        n = i + (j-1)*nvx
        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
    end do
    end do
 
    do j=1, ke_y
    do i=1, ke_x
        n = i + (j-1)*ke_x
        etov(n,1:4) = (j-1)*nvx + &
            & (/ i, i+1,              &
            &    i + nvx, i+1 + nvx /)
    end do
    end do
 
      !---
  !!$
  !!$    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 meshutil2d_genrectdomain
 
!OCL SERIAL
  subroutine meshutil2d_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,4)
 
    integer :: nodes(ne*nfaces,2)
    integer(kind=8) :: face_ids(ne*nfaces)
    integer :: ke
    integer :: f
    integer :: n
    integer :: n1, n2
    integer :: nnodes
    integer :: tmp
    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)
    !-----------------------------------------------------------------------------
 
    nnodes = maxval( etov )
    nnodes_row = size(nodes,1)
 
    !---------
    do ke=1, ne
       nodes(ke     ,:) = etov(ke,(/ 1, 2 /))
       nodes(ke+ne  ,:) = etov(ke,(/ 2, 4 /))
       nodes(ke+2*ne,:) = etov(ke,(/ 4, 3 /))
       nodes(ke+3*ne,:) = etov(ke,(/ 3, 1 /))
    end do
   
    ! Sort
    do n=1, nnodes_row
       if (nodes(n,1) > nodes(n,2) ) then
          tmp = nodes(n,1);
          nodes(n,1) = nodes(n,2)
          nodes(n,2) = tmp
       end if
       ! write(*,*) n, ":", nodes(n,:)
    end do
    nodes = nodes - 1
    !---------
 
    do ke=1, ne
       etoe(ke,:) = ke
       etof(ke,:) = (/ 1, 2, 3, 4 /)
    end do
 
    face_ids(:) = nodes(:,1)*nnodes + nodes(:,2) + 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 meshutil2d_genconnectivity
 
!OCL SERIAL
  subroutine meshutil2d_buildinteriormap( VMapM, VMapP, MapM, MapP, &
    pos_en, pos_ev, EtoE, EtoF, EtoV, Fmask, Ne, Np, Nfp, Nfaces, Nv)
 
    implicit none
 
    integer, intent(in) :: ne
    integer, intent(in) :: np
    integer, intent(in) :: nfp
    integer, intent(in) :: nfaces
    integer, intent(in) :: nv
 
    integer, intent(out) :: vmapm(nfp,nfaces,ne)
    integer, intent(out) :: vmapp(nfp,nfaces,ne)
    integer, intent(out) :: mapm(nfp,nfaces,ne)
    integer, intent(out) :: mapp(nfp,nfaces,ne)
 
    real(rp), intent(in) :: pos_en(np,ne,2)
    real(rp), intent(in) :: pos_ev(nv,2)
    integer, intent(in) :: etoe(ne, nfaces)
    integer, intent(in) :: etof(ne, nfaces)
    integer, intent(in) :: etov(ne, 4)
    integer, intent(in) :: fmask(nfp,4)
 
 
    integer :: ke, ke1, ke2
    integer :: f, f1, f2
    integer :: p
    integer :: n
    integer :: idp, idm
    integer :: v1, v2
    integer :: nodeids(np,ne)
    real(rp) :: x1(nfp,nfp), x2(nfp,nfp)
    real(rp) :: y1(nfp,nfp), y2(nfp,nfp)
    real(rp) :: dist(nfp,nfp)
    real(rp) :: x(np*ne), y(np*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)
    end do
    end do
    !$omp end do
    !$omp do collapse(2)
    do ke=1, ne
    do f=1, nfaces
    do p=1, nfp
      n = p + (f-1)*nfp + (ke-1)*nfp*nfaces
      mapm(p,f,ke) = n
      mapp(p,f,ke) = n
      vmapm(p,f,ke) = nodeids(fmask(p,f),ke)
    end do
    end do
    end do
 
    !$omp workshare
    vmapp(:,:,:) = -1
    !$omp end workshare
    !$omp end parallel
 
    !$omp parallel private( &
    !$omp ke1, f1, ke2, f2, v1, v2,                     &
    !$omp x1, x2, y1, y2, dist, mindist_indx, idP, idM  )
    
    !$omp do
    do ke1=1, ne
    do f1=1, nfaces
      ke2 = etoe(ke1,f1); f2 = etof(ke1,f1)
 
      v1 = etov(ke1,f1); v2 = etov(ke1,1+mod(f1,nfaces))
 
      x1(:,:) = spread( x(vmapm(:,f1,ke1)), 2, nfp )
      x2(:,:) = spread( x(vmapm(:,f2,ke2)), 1, nfp )
      y1(:,:) = spread( y(vmapm(:,f1,ke1)), 2, nfp )
      y2(:,:) = spread( y(vmapm(:,f2,ke2)), 1, nfp )
 
      dist(:,:) = (x1(:,:) - x2(:,:))**2 &
                + (y1(:,:) - y2(:,:))**2
      do idm=1, nfp
        mindist_indx(:) = minloc(dist(idm,:))
        idp = mindist_indx(1)
        vmapp(idm,f1,ke1) = vmapm(idp,f2,ke2)
        mapp(idm,f1,ke1) = idp + (f2-1)*nfp + (ke2-1)*nfp*nfaces
      end do
    end do
    end do
    !$omp end do
    !$omp end parallel
 
    !-----
  !    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(:)
 
  end subroutine meshutil2d_buildinteriormap
 
!OCL SERIAL
  subroutine meshutil2d_genpatchboundarymap(  VMapB, MapB, VMapP, &
    pos_en, xmin, xmax, ymin, ymax, Fmask, Ne, Np, Nfp, Nfaces, Nv)
 
    implicit none
 
    integer, intent(in) :: ne
    integer, intent(in) :: np
    integer, intent(in) :: nfp
    integer, intent(in) :: nfaces
    integer, intent(in) :: nv
 
    integer, intent(inout), allocatable :: vmapb(:)
    integer, intent(inout), allocatable :: mapb(:)
    integer, intent(inout) :: vmapp(nfp,nfaces,ne)
 
    real(rp), intent(in) :: pos_en(np,ne,2)
    real(rp), intent(in) :: xmin, xmax
    real(rp), intent(in) :: ymin, ymax
    integer, intent(in) :: fmask(nfp,4)
 
 
    integer :: ke
    integer :: b
    integer :: f
    integer :: i, j
    real(rp) :: x, y
 
    real(rp), parameter :: node_tol = 1.0e-12_rp
 
    real(rp) :: ordinfo(nfp*ne,4)
    integer :: elemids(nfp*ne,4)
    integer :: faceids(nfp*ne,4)
    integer :: counterb(4)
    integer :: mapb_counter
    real(rp) :: rdomx, rdomy
    !-----------------------------------------------------------------------------
 
 
    counterb(:) = 0
    rdomx = 1.0_rp/(xmax - xmin)
    rdomy = 1.0_rp/(ymax - ymin)
 
    do ke=1, ne
    do f=1, nfaces
      x = sum(pos_en(fmask(:,f),ke,1)/dble(nfp))
      y = sum(pos_en(fmask(:,f),ke,2)/dble(nfp))
 
      call eval_domain_boundary(1, y, ymin, x, ke, f, rdomy)
      call eval_domain_boundary(2, x, xmax, y, ke, f, rdomx)
      call eval_domain_boundary(3, y, ymax, x, ke, f, rdomy)
      call eval_domain_boundary(4, x, xmin, y, ke, f, rdomx)
    end do
    end do
 
    allocate( mapb(sum(counterb*nfp))    )
    allocate( vmapb(size(mapb)) )
 
    mapb_counter = 1
    do b = 1, 4
    !  write(*,*) "LocalMesh boundary ID:", b
    !  write(*,*) counterB(b)
    !  write(*,*) ordInfo(1:counterB(b),b)
    !  write(*,*) elemIds(1:counterB(b),b)
    !  write(*,*) faceIds(1:counterB(b),b)
 
      do i=1, counterb(b)
        ke = elemids(i,b); f = faceids(i,b)
        do j=1, nfp
          vmapp(j,f,ke) = np*ne + mapb_counter
          vmapb(mapb_counter) = fmask(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
     subroutine eval_domain_boundary(domb_id, r, rbc, ord_info, ke_, f_, normalized_fac)
        integer, intent(in) :: domb_id
        real(rp), intent(in) :: r
        real(rp), intent(in) :: rbc
        real(rp), intent(in) :: ord_info
        integer, intent(in) :: ke_, 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) = ke_
          faceids(counterb(domb_id),domb_id) = f_
        end if
 
        return
     end subroutine eval_domain_boundary
     
  end subroutine meshutil2d_genpatchboundarymap
 
 
!OCL SERIAL
  subroutine meshutil2d_genpatchboundarymap_wide( VMapB2, &
    VMapB, HaloSize,                                      &
    NeX, NeY, Nfp )
    use scale_prc, only: prc_abort
    implicit none
    integer, intent(in) :: halosize
    integer, intent(in) :: nex
    integer, intent(in) :: ney
    integer, intent(in) :: nfp
    integer, intent(inout) :: vmapb2(  2*(nex + ney)*nfp*halosize )
    integer, intent(in) :: vmapb( 2*(nex + ney)*nfp )
 
    integer :: f
    integer :: iso, isso, i, ii
    integer :: ke, p, ph, pv
    !------------------------------------------------------------
  
    if ( halosize > nfp ) then
      log_info("MeshUtil2D_genPatchBoundaryMap_wide",*) "HaloSize should be <= Nfp. Check!"
      call prc_abort
    end if
 
    !--
    f = 1
    iso = 0; isso = 0
    do ke=1, nex
      do ph=1, halosize
        do p=1, nfp
          i  = iso  + p              + (ke-1)*nfp
          ii = isso + p + (ph-1)*nfp + (ke-1)*nfp*halosize
          vmapb2(ii) = vmapb(i) + (ph-1)*nfp
        end do
      end do
    end do
 
    !--
    f = 2
    iso  = iso  + nex * nfp
    isso = isso + nex * nfp*halosize
    do ke=1, ney
      do ph=1, halosize
        do p=1, nfp
          i  = iso  + p              + (ke-1)*nfp
          ii = isso + p + (ph-1)*nfp + (ke-1)*nfp*halosize
          vmapb2(ii) = vmapb(i) - (ph-1)
        end do
      end do
    end do
 
    !--
    f = 3
    iso  = iso  + ney * nfp
    isso = isso + ney * nfp*halosize
    do ke=1, nex
      do ph=1, halosize
        do p=1, nfp
          i  = iso  + p              + (ke-1)*nfp
          ii = isso + p + (ph-1)*nfp + (ke-1)*nfp*halosize
          vmapb2(ii) = vmapb(i) - (ph-1)*nfp
        end do
      end do
    end do
 
    !--
    f = 4
    iso  = iso  + nex * nfp
    isso = isso + nex * nfp*halosize
    do ke=1, ney
      do ph=1, halosize
        do p=1, nfp
          i  = iso  + p              + (ke-1)*nfp
          ii = isso + p + (ph-1)*nfp + (ke-1)*nfp*halosize
          vmapb2(ii) = vmapb(i) + (ph-1)
        end do
      end do
    end do
    return
  end subroutine meshutil2d_genpatchboundarymap_wide
 
!OCL SERIAL
  subroutine meshutil2d_buildglobalmap( &
    panelID_table, pi_table, pj_table,           &
    tileID_map, tileFaceID_map, tilePanelID_map, &
    Ntile, isPeriodicX, isPeriodicY,             &
    Ne_x, Ne_y )
 
    implicit none
 
    integer, intent(in) :: ntile
    integer, intent(out) :: panelid_table(ntile)
    integer, intent(out) :: pi_table(ntile)
    integer, intent(out) :: pj_table(ntile)
    integer, intent(out) :: tileid_map(4,ntile)
    integer, intent(out) :: tilefaceid_map(4,ntile)
    integer, intent(out) :: tilepanelid_map(4,ntile)
    logical, intent(in) :: isperiodicx
    logical, intent(in) :: isperiodicy
    integer, intent(in) :: ne_x
    integer, intent(in) :: ne_y
 
    integer :: ntileperpanel
    integer :: nvx, nvy
    integer, allocatable :: nodesid_2d(:,:)
    integer, allocatable :: etov(:,:)
    integer, allocatable :: etoe(:,:)
    integer, allocatable :: etof(:,:)
    integer :: i, j, f
    integer :: panelid
    integer :: tileid, tileid_r
    integer :: counter
 
    !-----------------------------------------------------------------------------
 
    ntileperpanel = ntile/1
    nvx = ne_x + 1
    nvy = ne_y + 1
    allocate( nodesid_2d(nvx, nvy) )
    allocate( etov(ntile,4), etoe(ntile,4), etof(ntile,4) )
 
    counter = 0
    do j = 1, nvy
    do i = 1, nvx
      counter = counter + 1
      nodesid_2d(i,j) = counter
    end do
    end do
 
    !----
 
    tileid = 0
    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
      etov(tileid,:) = (/ nodesid_2d(i,j  ), nodesid_2d(i+1,j  ),    &
        &                 nodesid_2d(i,j+1), nodesid_2d(i+1,j+1)/)
    end do
    end do
 
    call meshutil2d_genconnectivity( etoe, etof, &
      & etov, ntile, 4 )
    tileid_map(:,:) = transpose(etoe)
    tilefaceid_map(:,:) = transpose(etof)
 
    do tileid=1, ntile
    do f=1, 4
      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
          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
          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
          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)
          tilefaceid_map(3,tileid) = 1
        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 meshutil2d_buildglobalmap
 
end module scale_meshutil_2d