scale_meshutil_2d Module Reference

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

Functions/Subroutines
subroutine, public	meshutil2d_genrectdomain (pos_v, etov, ke_x, xmin, xmax, ke_y, ymin, ymax)
subroutine, public	meshutil2d_genconnectivity (etoe, etof, etov, ne, nfaces)
subroutine, public	meshutil2d_buildinteriormap (vmapm, vmapp, mapm, mapp, pos_en, pos_ev, etoe, etof, etov, fmask, ne, np, nfp, nfaces, nv)
subroutine, public	meshutil2d_genpatchboundarymap (vmapb, mapb, vmapp, pos_en, xmin, xmax, ymin, ymax, fmask, ne, np, nfp, nfaces, nv)
subroutine, public	meshutil2d_buildglobalmap (panelid_table, pi_table, pj_table, tileid_map, tilefaceid_map, tilepanelid_map, ntile, isperiodicx, isperiodicy, ne_x, ne_y)

Detailed Description

module FElib / Mesh / utility for 2D mesh

Description

A module useful for generating 2D mesh

Author

Yuta Kawai, Team SCALE

Function/Subroutine Documentation

meshutil2d_genrectdomain()

subroutine, public scale_meshutil_2d::meshutil2d_genrectdomain	(	real(rp), dimension((ke_x+1)*(ke_y+1),2), intent(out)	pos_v,
		integer, dimension(ke_x*ke_y,4), intent(out)	etov,
		integer, intent(in)	ke_x,
		real(rp), intent(in)	xmin,
		real(rp), intent(in)	xmax,
		integer, intent(in)	ke_y,
		real(rp), intent(in)	ymin,
		real(rp), intent(in)	ymax )

Definition at line 34 of file scale_meshutil_2d.F90.

 
    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

Referenced by scale_mesh_rectdom2d::meshrectdom2d_setuplocaldom().

meshutil2d_genconnectivity()

subroutine, public scale_meshutil_2d::meshutil2d_genconnectivity	(	integer, dimension(ne, nfaces), intent(out)	etoe,
		integer, dimension(ne, nfaces), intent(out)	etof,
		integer, dimension(ne,4), intent(in)	etov,
		integer, intent(in)	ne,
		integer, intent(in)	nfaces )

Definition at line 94 of file scale_meshutil_2d.F90.

    
    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

Referenced by scale_mesh_rectdom2d::meshrectdom2d_setuplocaldom(), meshutil2d_buildglobalmap(), and scale_meshutil_cubedsphere2d::meshutilcubedsphere2d_buildglobalmap().

meshutil2d_buildinteriormap()

subroutine, public scale_meshutil_2d::meshutil2d_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,2), intent(in)	pos_en,
		real(rp), dimension(nv,2), intent(in)	pos_ev,
		integer, dimension(ne, nfaces), intent(in)	etoe,
		integer, dimension(ne, nfaces), intent(in)	etof,
		integer, dimension(ne, 4), intent(in)	etov,
		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 211 of file scale_meshutil_2d.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,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(:)
 

Referenced by scale_mesh_rectdom2d::meshrectdom2d_setuplocaldom().

meshutil2d_genpatchboundarymap()

subroutine, public scale_meshutil_2d::meshutil2d_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,2), intent(in)	pos_en,
		real(rp), intent(in)	xmin,
		real(rp), intent(in)	xmax,
		real(rp), intent(in)	ymin,
		real(rp), intent(in)	ymax,
		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 343 of file scale_meshutil_2d.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,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
     

Referenced by scale_mesh_rectdom2d::meshrectdom2d_setuplocaldom().

meshutil2d_buildglobalmap()

subroutine, public scale_meshutil_2d::meshutil2d_buildglobalmap	(	integer, dimension(ntile), intent(out)	panelid_table,
		integer, dimension(ntile), intent(out)	pi_table,
		integer, dimension(ntile), intent(out)	pj_table,
		integer, dimension(4,ntile), intent(out)	tileid_map,
		integer, dimension(4,ntile), intent(out)	tilefaceid_map,
		integer, dimension(4,ntile), intent(out)	tilepanelid_map,
		integer, intent(in)	ntile,
		logical, intent(in)	isperiodicx,
		logical, intent(in)	isperiodicy,
		integer, intent(in)	ne_x,
		integer, intent(in)	ne_y )

Definition at line 453 of file scale_meshutil_2d.F90.

 
    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

References meshutil2d_genconnectivity().

Referenced by scale_mesh_rectdom2d::meshrectdom2d_setuplocaldom().