scale_meshutil_1d.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Mesh / utility for 1D mesh
!!
!! @par Description
!!          A module useful for generating 1D mesh 
!!
!! @author Yuta Kawai, Team SCALE
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_meshutil_1d
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
 
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !  
  public :: meshutil1d_genlinedomain
  public :: meshutil1d_genconnectivity
  public :: meshutil1d_buildinteriormap
  public :: meshutil1d_genpatchboundarymap
  public :: meshutil1d_genpatchboundarymap_wide
  public :: meshutil1d_genperiodicboundarymap
  public :: meshutil1d_buildglobalmap
 
contains
!OCL SERIAL
  subroutine meshutil1d_genlinedomain( pos_v, EToV, & ! (out)
    ke_x, xmin, xmax,                               & ! (in)
    fx )                                              ! (in)
 
    implicit none
 
    integer, intent(in) :: ke_x
    real(rp), intent(in) :: xmin, xmax
    real(rp), intent(in), optional :: fx(ke_x+1)
 
    real(rp), intent(out) :: pos_v(ke_x+1,1)
    integer, intent(out) :: etov(ke_x,2)
 
    integer :: i
    integer :: n
    integer :: nvx
 
    !-----------------------------------------------------------------------------
 
    nvx = ke_x + 1
 
    do i=1, nvx
      if (present(fx) ) then
        pos_v(i,1) = fx(i)
      else
        pos_v(i,1) = (xmax - xmin)*dble(i - 1)/dble(nvx - 1) + xmin
      end if
    end do
 
    do i=1, ke_x
      etov(i,1:2) = (/ i, i+1 /)
    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 meshutil1d_genlinedomain
 
!OCL SERIAL
  subroutine meshutil1d_genconnectivity( EToE, EToF, & ! (out)
    etov, ne, nfaces )                                 ! (in)
    
    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,2)
 
    integer :: nodes(ne*nfaces,1)
    integer :: face_ids(ne*nfaces)
    integer :: spnodetonode(ne*nfaces,4)
    integer :: k
    integer :: f
    integer :: n
    integer :: n1, n2
    integer :: nnodes
    integer :: tmp
    integer :: spnodetonoderowtmp(4)
    integer :: nnodes_row
    integer :: matchl(2,4), matchr(2,4)
 
    real(rp) :: etoe_1d(ne*nfaces)
    real(rp) :: etof_1d(ne*nfaces)
    !-----------------------------------------------------------------------------
 
    nnodes = maxval( etov )
    nnodes_row = size(nodes,1)
 
    !---------
    do n=1, ne
       nodes(n     ,:) = etov(n,(/ 1 /))
       nodes(n+ne  ,:) = etov(n,(/ 2 /))
    end do
    nodes = nodes - 1
 
    !---------
 
    do n=1, ne
       etoe(n,:) = n
       etof(n,:) = (/ 1, 2 /)
    end do
 
    face_ids(:) = nodes(:,1) + 1
 
    do f=1, nfaces
    do k=1, ne
       n = k + (f-1)*ne
       spnodetonode(n,:) = (/ face_ids(n), n, etoe(k,f), etof(k,f) /)
       ! write(*,*) "face_id, n, EToE, EToF:", spNodeToNode(n,:)
    end do
    end do
 
    ! Sort row
    do n1=1, nnodes_row-1
    do n2=n1+1, nnodes_row
       if (spnodetonode(n1,1) > spnodetonode(n2,1)) then
          spnodetonoderowtmp(:) = spnodetonode(n1,:)
          spnodetonode(n1,:) = spnodetonode(n2,:)
          spnodetonode(n2,:) = spnodetonoderowtmp
       end if
    end do
    end do
 
    ! do n=1, Nnodes_row
    !    write(*,*) "(sorted) face_id, n, EToE, EToF:", spNodeToNode(n,:)
    ! end do
 
    !
    etoe_1d(:) = -1
    etof_1d(:) = -1
    do n=1, nnodes_row-1
       if ( spnodetonode(n,1) - spnodetonode(n+1,1) == 0 ) then
          matchl(:,:) = spnodetonode((/ n, n+1 /), :)
          matchr(:,:) = spnodetonode((/ n+1, n /), :)
 
          etoe_1d(matchl(:,2)) = matchr(:,3)
          etof_1d(matchl(:,2)) = matchr(:,4)
       end if
    end do
 
    do f=1, nfaces
    do k=1, ne
       n = k + (f-1)*ne
       if ( etoe_1d(n) /= -1 ) then
          etoe(k,f) = etoe_1d(n)
          etof(k,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 meshutil1d_genconnectivity
 
!OCL SERIAL
  subroutine meshutil1d_buildinteriormap( VMapM, VMapP, MapM, MapP,  & ! (out)
    pos_en, pos_ev, etoe, etof, etov, fmask, ne, np, nfp, nfaces, nv ) ! (in)
 
 
  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,1)
  real(rp), intent(in) :: pos_ev(nv,1)
  integer, intent(in) :: etoe(ne, nfaces)
  integer, intent(in) :: etof(ne, nfaces)
  integer, intent(in) :: etov(ne,2)
  integer, intent(in) :: fmask(nfp,2)
 
 
  integer :: k, k1, k2
  integer :: f, f1, f2
  integer :: p
  integer :: n
  integer :: idp, idm
  integer :: v1, v2
  integer :: nodeids(np,ne)
  real(rp) :: refd2
  real(rp) :: x1(nfp,nfp), x2(nfp,nfp)
  real(rp) :: dist(nfp,nfp)
  real(rp) :: x(np*ne)
  !-----------------------------------------------------------------------------
 
  do k=1, ne
  do p=1, np
     n = p + (k-1)*np
     nodeids(p,k) = n
     x(n) = pos_en(p,k,1)
  end do
  end do
 
  do k=1, ne
  do f=1, nfaces
  do p=1, nfp
     n = p + (f-1)*nfp + (k-1)*nfp*nfaces
     mapm(p,f,k) = n
     mapp(p,f,k) = n
     vmapm(p,f,k) = nodeids(fmask(p,f),k)
  end do
  end do
  end do
 
  vmapp(:,:,:) = -1
  do k1=1, ne
  do f1=1, nfaces
     k2 = etoe(k1,f1); f2 = etof(k1,f1)
 
     v1 = etov(k1,f1); v2 = etov(k1,1+mod(f1,nfaces))
     refd2 = (pos_ev(v1,1) - pos_ev(v2,1))**2
 
     x1(:,:) = spread( x(vmapm(:,f1,k1)), 2, nfp )
     x2(:,:) = spread( x(vmapm(:,f2,k2)), 1, nfp )
  
     dist(:,:) = (x1 - x2)**2
     do idp=1, nfp
     do idm=1, nfp
        if (dist(idm,idp)/refd2 < 1.0e-14_rp) then
           vmapp(idm,f1,k1) = vmapm(idp,f2,k2)
           mapp(idm,f1,k1) = idp + (f2-1)*nfp + (k2-1)*nfp*nfaces
        end if
     end do
     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=1,Ne
  !     write(*,*) "k=", k, "---"
  !     write(*,*) " VMapM:", VMapM(:,:,k)
  !     write(*,*) " VMapP:", VMapP(:,:,k)
  !  end do
 
!    write(*,*) "mapB:", mesh%mapB(:)
!    write(*,*) "vmapB:", mesh%vmapB(:)
 
    return
  end subroutine meshutil1d_buildinteriormap
 
!OCL SERIAL
  subroutine meshutil1d_genpatchboundarymap(  VMapB, MapB, VMapP, & ! (inout)
    pos_en, xmin, xmax, fmask, ne, np, nfp, nfaces, nv            ) ! (in)
 
    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,1)
    real(rp), intent(in) :: xmin, xmax
    integer, intent(in) :: fmask(nfp,2)
 
    integer :: k
    integer :: b
    integer :: f
    integer :: i, j
    real(rp) :: x
 
    real(rp), parameter :: node_tol = 1.0e-12_rp
 
    integer :: elemids(nfp*ne,2)
    integer :: faceids(nfp*ne,2)
    integer :: counterb(2)
    integer :: mapb_counter
    real(rp) :: rdomx
    !-----------------------------------------------------------------------------
 
    counterb(:) = 0d0
    rdomx = 1.0_rp/(xmax - xmin)
 
    do k=1, ne
    do f=1, nfaces
      x = sum(pos_en(fmask(:,f),k,1)/dble(nfp))
 
      call eval_domain_boundary(1, x, xmin, k, f, rdomx)
      call eval_domain_boundary(2, x, xmax, k, f, rdomx)
    end do
    end do
 
 
    allocate( mapb(sum(counterb*nfp))    )
    allocate( vmapb(size(mapb)) )
 
    mapb_counter = 1
    do b = 1, 2
      ! write(*,*) "LocalMesh boundary ID:", b
      ! write(*,*) counterB(b)
      ! write(*,*) ordInfo(1:counterB(1),b)
      ! write(*,*) elemIds(1:counterB(1),b)
      ! write(*,*) faceIds(1:counterB(1),b)
 
      do i=1, counterb(b)
        k = elemids(i,b); f = faceids(i,b)
        do j=1, nfp
          vmapp(j,f,k) = np*ne + mapb_counter
          vmapb(mapb_counter) = fmask(j,f) + (k-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, k_, f_, normalized_fac)
      implicit none
      integer, intent(in) :: domb_id
      real(rp), intent(in) :: r
      real(rp), intent(in) :: rbc
      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
        elemids(counterb(domb_id),domb_id) = k_
        faceids(counterb(domb_id),domb_id) = f_
      end if
 
      return
    end subroutine eval_domain_boundary
  end subroutine meshutil1d_genpatchboundarymap
 
!OCL SERIAL
  subroutine meshutil1d_genpatchboundarymap_wide( VMapB2, &
    VMapB, HaloSize,                                      &
    NeX, Np )
    use scale_prc, only: prc_abort
    implicit none
    integer, intent(in) :: halosize
    integer, intent(in) :: nex
    integer, intent(in) :: np
    integer, intent(inout) :: vmapb2(2*halosize)
    integer, intent(in) :: vmapb(2)
 
    integer :: f
    integer :: iso, isso, i, ii
    integer :: ke, p, ph, pv
    !------------------------------------------------------------
  
    if ( halosize > np ) then
      log_info("MeshUtil1D_genPatchBoundaryMap_wide",*) "HaloSize should be <= Np. Check!"
      call prc_abort
    end if
 
    !--
    f = 1
    iso = 0; isso = 0
    do ph=1, halosize
      i  = iso  + 1
      ii = isso + ph
      vmapb2(ii) = vmapb(i) + (ph-1)
    end do
 
    !--
    f = 2
    iso  = iso  + 1
    isso = isso + halosize
    do ph=1, halosize
      i  = iso  + 1
      ii = isso + ph
      vmapb2(ii) = vmapb(i) - (ph-1)
    end do
    return
  end subroutine meshutil1d_genpatchboundarymap_wide
 
!OCL SERIAL
  subroutine meshutil1d_genperiodicboundarymap( EToE, EToF, VMapP, & ! (inout)
    xperiod, x, fx, fmask, ne, np, nfp, nfaces, nv                 ) ! (in)
 
    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) :: etoe(ne, nfaces)
    integer, intent(inout) :: etof(ne, nfaces)  
    integer, intent(inout) :: vmapp(nfp,nfaces,ne)    
    real(rp), intent(in) :: xperiod
    real(rp), intent(in) :: x(np,ne)
    real(rp), intent(in) :: fx(nfp*ne)
    integer, intent(in) :: fmask(nfp,4)
 
    integer :: k, k1, k2
    integer :: f, f1, f2
    integer :: p
    integer :: n
 
    real(rp) :: cx1, cx2
    real(rp) :: dx_pbc
    integer :: idp, idm
    integer :: vidl(nfp), vidr(nfp)
    integer :: fidl(nfp), fidr(nfp)
 
    real(rp) :: refd2
    real(rp) :: x1(nfp,nfp), x2(nfp,nfp)
    real(rp) :: dist(nfp,nfp)
 
    real(rp), parameter :: node_tol = 1.0e-12_rp
    !-----------------------------------------------------------------------------
 
 
    do k1=1, ne
    do f1=1, nfaces
     cx1 = sum(x(fmask(:,f1),k1)/dble(nfp))
     
     k2 = etoe(k1,f1); f2 = etof(k1,f1)
     if (k2==k1) then
       do k2=1, ne
         if (k1 /= k2) then
           do f2=1, nfaces
             if (etoe(k2,f2) == k2) then
               cx2 = sum(x(fmask(:,f2),k2)/dble(nfp))
               
               dx_pbc = sqrt( (abs(cx1-cx2)-xperiod)**2 )
               
               if (dx_pbc < node_tol) then
 
                 etoe(k1,f1) = k2; etoe(k2,f2) = k1
                 etof(k1,f1) = f2; etof(k2,f2) = f1
 
                 !mesh%BCType(f1,k1) = BCTYPE_PERIODIC
                 !mesh%BCType(f2,k2) = BCTYPE_PERIODIC
 
                 do p=1, nfp
                   vidl(p) = fmask(p,f1) + (k1-1)*np
                   vidr(p) = fmask(p,f2) + (k2-1)*np
                   fidl(p) = p + (f1-1)*nfp + (k1-1)*nfp*nfaces
                   fidr(p) = p + (f2-1)*nfp + (k2-1)*nfp*nfaces
                 end do
 
                 x1(:,:) = spread( fx(fidl(:)), 2, nfp )
                 x2(:,:) = spread( fx(fidr(:)), 1, nfp )
                 
                 dist(:,:) = (abs(x1-x2) -xperiod)**2
 
                 do idp=1,nfp
                 do idm=1,nfp
                     if (dist(idm,idp) < node_tol) then
                       vmapp(idm,f1,k1) = vidr(idp)
                       vmapp(idp,f2,k2) = vidl(idm)
                     end if
                 end do
                 end do
 
               end if
             end if
           end do
         end if
       end do
     end if
 
    end do
    end do
 
   !-------
 
!    write(*,*) "Build MapInfo (Periodic BC): "
!    do k=1,mesh%Ne
!       write(*,*) "k=", k, "---"
!       write(*,*) " VMapM:", mesh%VMapM(:,:,k)
!       write(*,*) " VMapP:", mesh%VMapP(:,:,k)
!       write(*,*) "BCType:", mesh%BCType(:,k)
!    end do
 
    return
  end subroutine meshutil1d_genperiodicboundarymap
 
!OCL SERIAL
  subroutine meshutil1d_buildglobalmap( &
    panelID_table, pi_table,                     & ! (out)
    tileid_map, tilefaceid_map, tilepanelid_map, & ! (out)
    ntile )                                        ! (in)
    
    use scale_prc, only: prc_ismaster
    implicit none
 
    integer, intent(in) :: ntile
    integer, intent(out) :: panelid_table(ntile)
    integer, intent(out) :: pi_table(ntile)
    integer, intent(out) :: tileid_map(2,ntile)
    integer, intent(out) :: tilefaceid_map(2,ntile)
    integer, intent(out) :: tilepanelid_map(2,ntile)
 
    integer :: ntileperpanel
    integer :: ne, nv
    integer, allocatable :: nodesid_1d(:)
    integer, allocatable :: etov(:,:)
    integer, allocatable :: etoe(:,:)
    integer, allocatable :: etof(:,:)
    integer :: i, f
    integer :: panelid
    integer :: tileid, tileid_r
    integer :: counter
    real(dp) :: del
 
    !-----------------------------------------------------------------------------
 
    ntileperpanel = ntile
    ne = ntileperpanel
    nv = ne + 1
    allocate( nodesid_1d(nv) )
    allocate( etov(ntile,2), etoe(ntile,2), etof(ntile,2) )
 
    nodesid_1d(:) = -1
    counter = 0
    do i = 2, nv-1
      counter = counter + 1
      nodesid_1d(i) = counter
    end do
 
    !----
 
    tileid = 0
 
    do i = 1, ne
      tileid = tileid + 1
      panelid_table(tileid) = 1
      pi_table(tileid) = i
      etov(tileid,:) = (/ nodesid_1d(i), nodesid_1d(i+1) /)
    end do
 
    call meshutil1d_genconnectivity( etoe, etof, &
      & etov, ntile, 2 )
    tileid_map(:,:) = transpose(etoe)
    tilefaceid_map(:,:) = transpose(etof)
 
    do tileid=1, ntile
    do f=1, 2
      tileid_r = tileid_map(f,tileid)
      tilepanelid_map(f,tileid) = panelid_table(tileid_r)
    end do
    end do
 
    !--
 
    ! if (PRC_isMaster) then
    !
    !   write(*,*) "TotTile", tileID, 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
 
    return
  end subroutine meshutil1d_buildglobalmap
 
end module scale_meshutil_1d