scale_mesh_base2d.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Mesh / Base 2D
!!
!! @par Description
!!      Base module to manage 2D meshes for element-based methods
!!
!! @author Yuta Kawai, Team SCALE
!<
#include "scaleFElib.h"
module scale_mesh_base2d
 
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
 
  use scale_localmesh_2d, only: &
    localmesh2d, localmesh2d_init, localmesh2d_final
 
  use scale_mesh_base, only: &
    meshbase, meshbase_init, meshbase_final
    
  use scale_element_base, only: elementbase2d
 
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  ! 
  !> Derived type to manage a computational mesh (base type for 2D domain)
  type, abstract, public, extends(meshbase) :: meshbase2d
    type(localmesh2d), allocatable :: lcmesh_list(:) !< Array of local meshes in each MPI process
    type(elementbase2d), pointer :: refelem2d        !< Pointer to an object with a reference element for 2D
  contains
    procedure(meshbase2d_generate), deferred :: generate 
    procedure :: getlocalmesh => meshbase2d_get_localmesh
  end type meshbase2d
 
  interface 
    subroutine meshbase2d_generate(this)
      import meshbase2d
      class(meshbase2d), intent(inout), target :: this
    end subroutine meshbase2d_generate
  end interface
 
  public :: meshbase2d_init, meshbase2d_final
  public :: meshbase2d_setgeometricinfo
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  integer, public :: meshbase2d_dimtype_num   = 4
  integer, public :: meshbase2d_dimtypeid_x   = 1
  integer, public :: meshbase2d_dimtypeid_y   = 2
  integer, public :: meshbase2d_dimtypeid_xy  = 3
  integer, public :: meshbase2d_dimtypeid_xyt = 4
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
 
contains
  !> Initialize an object to manage a 2D computational mesh
!OCL SERIAL
  subroutine meshbase2d_init(this, &
    refElem, NLocalMeshPerPrc,     &
    nprocs, myrank                 )
    
    implicit none
    class(meshbase2d), intent(inout) :: this
    class(elementbase2d), intent(in), target :: refelem
    integer, intent(in) :: nlocalmeshperprc
    integer, intent(in), optional :: nprocs
    integer, intent(in), optional :: myrank
 
    integer :: n
    !-----------------------------------------------------------------------------
    
    this%refElem2D => refelem
    call meshbase_init( this, &
      meshbase2d_dimtype_num, refelem, &
      nlocalmeshperprc, 4,             &
      nprocs                           )
 
    allocate( this%lcmesh_list(this%LOCAL_MESH_NUM) )
    do n=1, this%LOCAL_MESH_NUM
      call localmesh2d_init( this%lcmesh_list(n), n, refelem, myrank )
    end do
 
    call this%SetDimInfo( meshbase2d_dimtypeid_x, "x", "m", "X-coordinate" )
    call this%SetDimInfo( meshbase2d_dimtypeid_y, "y", "m", "Y-coordinate" )
    call this%SetDimInfo( meshbase2d_dimtypeid_xy, "xy", "m", "XY-coordinate" )
    call this%SetDimInfo( meshbase2d_dimtypeid_xyt, "xyt", "m", "XY-coordinate" )
 
    return
  end subroutine meshbase2d_init
 
  !> Finalize an object to manage a 2D computational mesh
!OCL SERIAL
  subroutine meshbase2d_final( this )
    implicit none
    class(meshbase2d), intent(inout) :: this
    integer :: n
    !-----------------------------------------------------------------------------
  
    if ( allocated ( this%lcmesh_list ) ) then 
      do n=1, this%LOCAL_MESH_NUM
        call localmesh2d_final( this%lcmesh_list(n), this%isGenerated )
      end do
  
      deallocate( this%lcmesh_list )
    end if
    
    call meshbase_final(this)
 
    return
  end subroutine meshbase2d_final
  
!OCL SERIAL
  subroutine meshbase2d_get_localmesh( this, id, ptr_lcmesh )
    use scale_localmesh_base, only: localmeshbase
    implicit none
    class(meshbase2d), target, intent(in) :: this
    integer, intent(in) :: id
    class(localmeshbase), pointer, intent(out) :: ptr_lcmesh
    !-------------------------------------------------------------
 
    ptr_lcmesh => this%lcmesh_list(id)
    return
  end subroutine meshbase2d_get_localmesh
 
!OCL SERIAL
  subroutine meshbase2d_setgeometricinfo( lcmesh, coord_conv, calc_normal )
    use scale_mesh_base, only: meshbase_setgeometricinfo
    implicit none
    type(localmesh2d), intent(inout) :: lcmesh
    interface
      subroutine coord_conv( x, y, xr, xs, yr, ys, &
        vx, vy, elem )
        import elementbase2d
        import rp
        type(elementbase2d), intent(in) :: elem
        real(rp), intent(out) :: x(elem%np), y(elem%np)
        real(rp), intent(out) :: xr(elem%np), xs(elem%np), yr(elem%np), ys(elem%np)
        real(rp), intent(in) :: vx(elem%nv), vy(elem%nv)        
      end subroutine coord_conv
      subroutine calc_normal( normal_fn, &
        Escale_f, fid, elem )
        import elementbase2d
        import rp
        type(elementbase2d), intent(in) :: elem
        real(rp), intent(out) :: normal_fn(elem%nfptot,2)
        integer, intent(in) :: fid(elem%nfp,elem%nfaces)
        real(rp), intent(in) :: escale_f(elem%nfptot,2,2)
      end subroutine calc_normal      
    end interface
 
    class(elementbase2d), pointer :: refelem
    integer :: ke
    integer :: f
    integer :: i, j
    integer :: d
    integer :: node_ids(lcmesh%refelem%nv)
    real(rp) :: vx(lcmesh%refelem%nv), vy(lcmesh%refelem%nv)
    real(rp) :: xr(lcmesh%refelem%np), xs(lcmesh%refelem%np)
    real(rp) :: yr(lcmesh%refelem%np), ys(lcmesh%refelem%np)
    integer :: fmask(lcmesh%refelem%nfptot)
    integer :: fid(lcmesh%refelem2d%nfp,lcmesh%refelem2d%nfaces)
    real(rp) :: escale_f(lcmesh%refelem%nfptot,2,2)
 
    !-----------------------------------------------------------------------------
 
    refelem => lcmesh%refElem2D
 
    call meshbase_setgeometricinfo(lcmesh, 2)
 
    allocate( lcmesh%G_ij(refelem%Np,lcmesh%Ne,2,2) )
    allocate( lcmesh%GIJ (refelem%Np,lcmesh%Ne,2,2) )
    allocate( lcmesh%lon(refelem%Np,lcmesh%Ne) )
    allocate( lcmesh%lat(refelem%Np,lcmesh%Ne) )
    !$acc enter data create(lcmesh%G_ij, lcmesh%GIJ, lcmesh%lon, lcmesh%lat)
 
    fmask(:) = reshape(refelem%Fmask, shape(fmask))
    do f=1, refelem%Nfaces
    do i=1, refelem%Nfp
       fid(i,f) = i + (f-1)*refelem%Nfp
    end do
    end do
 
    do ke=1, lcmesh%Ne
      node_ids(:) = lcmesh%EToV(ke,:)
      vx(:) = lcmesh%pos_ev(node_ids(:),1)
      vy(:) = lcmesh%pos_ev(node_ids(:),2)
      call coord_conv( &
        lcmesh%pos_en(:,ke,1), lcmesh%pos_en(:,ke,2), xr, xs, yr, ys, & ! (out)
        vx, vy, refelem )                                             ! (in)
      
      lcmesh%J(:,ke) = - xs*yr + xr*ys
 
      lcmesh%Escale(:,ke,1,1) =   ys/lcmesh%J(:,ke)
      lcmesh%Escale(:,ke,1,2) = - xs/lcmesh%J(:,ke)
      lcmesh%Escale(:,ke,2,1) = - yr/lcmesh%J(:,ke)
      lcmesh%Escale(:,ke,2,2) =   xr/lcmesh%J(:,ke)
      
      !* Face
      
      !
      !mesh%fx(:,n) = mesh%x(fmask(:),n)
      !mesh%fy(:,n) = mesh%y(fmask(:),n)
 
      ! Calculate normal vectors
      do j=1, 2
      do i=1, 2
        escale_f(:,i,j) = lcmesh%Escale(fmask(:),ke,i,j)
      end do
      end do
      
      call calc_normal( lcmesh%normal_fn(:,ke,:), & ! (out)
        escale_f, fid, refelem )                   ! (in)
      
      lcmesh%sJ(:,ke) = sqrt( lcmesh%normal_fn(:,ke,1)**2 + lcmesh%normal_fn(:,ke,2)**2 )
      do d=1, 2
        lcmesh%normal_fn(:,ke,d) = lcmesh%normal_fn(:,ke,d)/lcmesh%sJ(:,ke)
      end do
      lcmesh%sJ(:,ke) = lcmesh%sJ(:,ke)*lcmesh%J(fmask(:),ke)
 
      lcmesh%Fscale(:,ke) = lcmesh%sJ(:,ke)/lcmesh%J(fmask(:),ke)       
    end do
    !$acc update device(lcmesh%pos_en, lcmesh%normal_fn, lcmesh%sJ, lcmesh%J, lcmesh%Escale, lcmesh%Fscale)
 
    !$omp parallel
    !$omp workshare
    lcmesh%Gsqrt (:,:)     = 1.0_rp
    lcmesh%GIJ   (:,:,1,1) = 1.0_rp
    lcmesh%GIJ   (:,:,2,1) = 0.0_rp
    lcmesh%GIJ   (:,:,1,2) = 0.0_rp
    lcmesh%GIJ   (:,:,2,2) = 1.0_rp
    lcmesh%G_ij  (:,:,1,1) = 1.0_rp
    lcmesh%G_ij  (:,:,2,1) = 0.0_rp
    lcmesh%G_ij  (:,:,1,2) = 0.0_rp
    lcmesh%G_ij  (:,:,2,2) = 1.0_rp
    !$omp end workshare
    !$omp end parallel
    !$acc update device(lcmesh%Gsqrt, lcmesh%GIJ, lcmesh%G_ij)
 
    return
  end subroutine meshbase2d_setgeometricinfo
  
end module scale_mesh_base2d