scale_meshutil_1d Module Reference

module FElib / Mesh / utility for 1D mesh More...

Functions/Subroutines
subroutine, public	meshutil1d_genlinedomain (pos_v, etov, ke_x, xmin, xmax, fx)
subroutine, public	meshutil1d_genconnectivity (etoe, etof, etov, ne, nfaces)
subroutine, public	meshutil1d_buildinteriormap (vmapm, vmapp, mapm, mapp, pos_en, pos_ev, etoe, etof, etov, fmask, ne, np, nfp, nfaces, nv)
subroutine, public	meshutil1d_genpatchboundarymap (vmapb, mapb, vmapp, pos_en, xmin, xmax, fmask, ne, np, nfp, nfaces, nv)
subroutine, public	meshutil1d_genperiodicboundarymap (etoe, etof, vmapp, xperiod, x, fx, fmask, ne, np, nfp, nfaces, nv)
subroutine, public	meshutil1d_buildglobalmap (panelid_table, pi_table, tileid_map, tilefaceid_map, tilepanelid_map, ntile)

Detailed Description

module FElib / Mesh / utility for 1D mesh

Description

A module useful for generating 1D mesh

Author

Yuta Kawai, Team SCALE

Function/Subroutine Documentation

meshutil1d_genlinedomain()

subroutine, public scale_meshutil_1d::meshutil1d_genlinedomain	(	real(rp), dimension(ke_x+1,1), intent(out)	pos_v,
		integer, dimension(ke_x,2), intent(out)	etov,
		integer, intent(in)	ke_x,
		real(rp), intent(in)	xmin,
		real(rp), intent(in)	xmax,
		real(rp), dimension(ke_x+1), intent(in), optional	fx )

Definition at line 35 of file scale_meshutil_1d.F90.

 
    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

Referenced by scale_mesh_base1d::meshbase1d_setuplocaldom().

meshutil1d_genconnectivity()

subroutine, public scale_meshutil_1d::meshutil1d_genconnectivity	(	integer, dimension(ne, nfaces), intent(out)	etoe,
		integer, dimension(ne, nfaces), intent(out)	etof,
		integer, dimension(ne,2), intent(in)	etov,
		integer, intent(in)	ne,
		integer, intent(in)	nfaces )

Definition at line 90 of file scale_meshutil_1d.F90.

    
    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

Referenced by scale_mesh_base1d::meshbase1d_setuplocaldom(), and meshutil1d_buildglobalmap().

meshutil1d_buildinteriormap()

subroutine, public scale_meshutil_1d::meshutil1d_buildinteriormap	(	integer, dimension(nfp,nfaces,ne), intent(out)	vmapm,
		integer, dimension(nfp,nfaces,ne), intent(out)	vmapp,
		integer, dimension(nfp,nfaces,ne), intent(out)	mapm,
		integer, dimension(nfp,nfaces,ne), intent(out)	mapp,
		real(rp), dimension(np,ne,1), intent(in)	pos_en,
		real(rp), dimension(nv,1), intent(in)	pos_ev,
		integer, dimension(ne, nfaces), intent(in)	etoe,
		integer, dimension(ne, nfaces), intent(in)	etof,
		integer, dimension(ne,2), intent(in)	etov,
		integer, dimension(nfp,2), intent(in)	fmask,
		integer, intent(in)	ne,
		integer, intent(in)	np,
		integer, intent(in)	nfp,
		integer, intent(in)	nfaces,
		integer, intent(in)	nv )

Definition at line 198 of file scale_meshutil_1d.F90.

 
 
  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

Referenced by scale_mesh_base1d::meshbase1d_setuplocaldom().

meshutil1d_genpatchboundarymap()

subroutine, public scale_meshutil_1d::meshutil1d_genpatchboundarymap	(	integer, dimension(:), intent(inout), allocatable	vmapb,
		integer, dimension(:), intent(inout), allocatable	mapb,
		integer, dimension(nfp,nfaces,ne), intent(inout)	vmapp,
		real(rp), dimension(np,ne,1), intent(in)	pos_en,
		real(rp), intent(in)	xmin,
		real(rp), intent(in)	xmax,
		integer, dimension(nfp,2), intent(in)	fmask,
		integer, intent(in)	ne,
		integer, intent(in)	np,
		integer, intent(in)	nfp,
		integer, intent(in)	nfaces,
		integer, intent(in)	nv )

Definition at line 315 of file scale_meshutil_1d.F90.

 
    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

Referenced by scale_mesh_base1d::meshbase1d_setuplocaldom().

meshutil1d_genperiodicboundarymap()

subroutine, public scale_meshutil_1d::meshutil1d_genperiodicboundarymap	(	integer, dimension(ne, nfaces), intent(inout)	etoe,
		integer, dimension(ne, nfaces), intent(inout)	etof,
		integer, dimension(nfp,nfaces,ne), intent(inout)	vmapp,
		real(rp), intent(in)	xperiod,
		real(rp), dimension(np,ne), intent(in)	x,
		real(rp), dimension(nfp*ne), intent(in)	fx,
		integer, dimension(nfp,4), intent(in)	fmask,
		integer, intent(in)	ne,
		integer, intent(in)	np,
		integer, intent(in)	nfp,
		integer, intent(in)	nfaces,
		integer, intent(in)	nv )

Definition at line 416 of file scale_meshutil_1d.F90.

 
    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

meshutil1d_buildglobalmap()

subroutine, public scale_meshutil_1d::meshutil1d_buildglobalmap	(	integer, dimension(ntile), intent(out)	panelid_table,
		integer, dimension(ntile), intent(out)	pi_table,
		integer, dimension(2,ntile), intent(out)	tileid_map,
		integer, dimension(2,ntile), intent(out)	tilefaceid_map,
		integer, dimension(2,ntile), intent(out)	tilepanelid_map,
		integer, intent(in)	ntile )

Definition at line 520 of file scale_meshutil_1d.F90.

    
    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

References meshutil1d_genconnectivity().

Referenced by scale_mesh_base1d::meshbase1d_assigndomid().