spec_grid.f90

Spectral Element Method for wave propagation and rupture dynamics.

  integer :: i,j,e,nel
  ! WARNING: dimension should be larger than the maximum expected
  ! number of neighbour elements 
  integer, dimension(10) :: etmp,jtmp,itmp

  nel = 0
  do e= 1,grid%nelem
  do j= 1,grid%ngll
  do i= 1,grid%ngll
    if (grid%ibool(i,j,e)==iglob) then
      nel = nel +1
      etmp(nel) = e
      itmp(nel) = i
      jtmp(nel) = j
    endif
  enddo
  enddo
  enddo

  if (associated(etab)) deallocate(etab)
  if (associated(itab)) deallocate(itab)
  if (associated(jtab)) deallocate(jtab)
  allocate(etab(nel),itab(nel),jtab(nel))
  etab = etmp(1:nel)
  itab = itmp(1:nel)
  jtab = jtmp(1:nel)

  end subroutine SE_node_belongs_to

!=====================================================================
!
  subroutine SE_find_point(coord,grid,e,xi,eta,new_coord)

  double precision, intent(in)  :: coord(NDIME)
  type(sem_grid_type), intent(in) :: grid
  integer, intent(out) :: e
  double precision, intent(out) :: xi,eta,new_coord(NDIME)

  integer :: iglob,k,istat
  integer, pointer, dimension(:) :: itab,jtab,etab

  call SE_find_nearest_node(coord,grid,iglob)
  call SE_node_belongs_to(iglob,etab,itab,jtab,grid)
  do k = 1,size(etab)
    xi=grid%xgll(itab(k))
    eta=grid%xgll(jtab(k))
    e = etab(k)
    call FE_find_point(coord,grid%fem,e,xi,eta,new_coord,istat)
    if (istat==0) exit
  enddo
  
  deallocate(itab,jtab,etab)
  if (istat>0) call IO_abort('SE_find_point: point not found')

  end subroutine SE_find_point

!=====================================================================
! get the list of bulk indices of the EDGE nodes of an ELEMENT, 
! assume counterclockwise numbering for the output list NODES

  subroutine SE_get_edge_nodes(grid,element,edge,nodes)

  type(sem_grid_type), intent(in) :: grid
  integer, intent(in) :: element,edge
  integer, intent(out) :: nodes(grid%ngll)

  select case (edge)
    case(edge_D); nodes = grid%ibool(:,1,element)
    case(edge_R); nodes = grid%ibool(grid%ngll,:,element)
    case(edge_U); nodes = grid%ibool(grid%ngll:1:-1,grid%ngll,element)
    case(edge_L); nodes = grid%ibool(1,grid%ngll:1:-1,element)
  end select

  end subroutine SE_get_edge_nodes


!=======================================================================
!
! Jacobian matrix  =  | dx/dxi   dx/deta |
!                     | dz/dxi   dz/deta |
!
  function SE_Jacobian_e(sem,e) result(jac)

  type(sem_grid_type), intent(in) :: sem
  integer, intent(in) :: e
  double precision :: jac(2,2,sem%ngll,sem%ngll)

  double precision, pointer :: coorg(:,:)
  integer :: i,j

  coorg => FE_GetElementCoord(sem%fem,e)
  do j = 1,sem%ngll
  do i = 1,sem%ngll
    jac(:,:,i,j) = matmul(coorg, sem%dshape(:,:,i,j) )
  enddo
  enddo
  deallocate(coorg)

  end function SE_Jacobian_e

!-----------------------------------------------------------------------

  function SE_Jacobian_eij(sem,e,i,j) result(jac)

  type(sem_grid_type), intent(in) :: sem
  integer, intent(in) :: e,i,j
  double precision :: jac(2,2)
  double precision, pointer :: coorg(:,:)

  coorg => FE_GetElementCoord(sem%fem,e)
  jac = matmul(coorg, sem%dshape(:,:,i,j) )
  deallocate(coorg)

  end function SE_Jacobian_eij


!=======================================================================
!
!-- Compute weights(i,j) = wgll(i)*wgll(j)*dvol(i,j) for one element
!
  function SE_VolumeWeights_e(sem,e) result(dvol)

  type(sem_grid_type), intent(in) :: sem
  integer, intent(in) :: e
  double precision :: dvol(sem%ngll,sem%ngll)

  double precision :: jac(2,2,sem%ngll,sem%ngll)

  jac = SE_Jacobian_e(sem,e)
  dvol = jac(1,1,:,:)*jac(2,2,:,:) - jac(1,2,:,:)*jac(2,1,:,:)
  dvol = dvol*sem%wgll2

  end function SE_VolumeWeights_e

!-----------------------------------------------------------------------

  function SE_VolumeWeights_eij(sem,e,i,j) result(dvol)

  type(sem_grid_type), intent(in) :: sem
  integer, intent(in) :: e,i,j
  double precision :: dvol

  double precision :: jac(2,2)

  jac = SE_Jacobian_eij(sem,e,i,j)
  dvol = ( jac(1,1)*jac(2,2) - jac(1,2)*jac(2,1) )*sem%wgll2(i,j)

  end function SE_VolumeWeights_eij


!=======================================================================
!
! Inverse Jacobian matrix  =  | dxi/dx    dxi/dz |
!                             | deta/dx  deta/dz |
!
  function SE_InverseJacobian_e(sem,e) result(jaci)

  use utils, only : invert2

  type(sem_grid_type), intent(in) :: sem
  integer, intent(in) :: e
  double precision :: jaci(2,2,sem%ngll,sem%ngll)

  double precision :: jac(2,2,sem%ngll,sem%ngll)
  integer :: i,j

  jac = SE_Jacobian_e(sem,e)
  do j = 1,sem%ngll
  do i = 1,sem%ngll
    jaci(:,:,i,j) = invert2(jac(:,:,i,j))
  enddo
  enddo

  end function SE_InverseJacobian_e

!----------------------------------------------------------------------

  function SE_InverseJacobian_eij(sem,e,i,j) result(jaci)

  use utils, only : invert2

  type(sem_grid_type), intent(in) :: sem
  integer, intent(in) :: e,i,j
  double precision :: jaci(2,2)

  double precision :: jac(2,2)

  jac = SE_Jacobian_eij(sem,e,i,j)
  jaci = invert2(jac)

  end function SE_InverseJacobian_eij


!=======================================================================
!
subroutine SE_BcTopoInit(grid)
  type(sem_grid_type), intent(inout) :: grid
  integer :: i
  grid%bounds => grid%fem%bnds
  do i=1,size(grid%bounds)
    call BC_set_bulk_node(grid%bounds(i),grid)
  enddo
end subroutine SE_BcTopoInit

!=======================================================================
!
subroutine BC_set_bulk_node(bc,grid)

  use generic_list, only : Link_Ptr_Type,Link_Type,List_Type &
     ,LI_Init_List,LI_Add_To_Head,LI_Get_Head &
     ,LI_Get_Next,LI_Associated,LI_Remove_Head

  use utils, only: drank

  type(bnd_grid_type), intent(inout) :: bc
  type(sem_grid_type), intent(in) :: grid
  
  type BC_Node_Data_Type
    integer :: bulk_node,bc_node
  end type BC_Node_Data_Type

  type BC_Node_Type
    TYPE(Link_Type) :: Link
    type(BC_Node_Data_Type), pointer :: Data
  end type BC_Node_Type

  type BC_Node_Ptr_Type
    type(BC_Node_Type), pointer :: P
  end type BC_Node_Ptr_Type

  integer :: bulk_node,bulk_node_vect(grid%ngll),i,n,kloc,bc_npoin,bc_inode
  logical :: node_at_corner,new_node
  type(List_Type)      :: BC_Nodes_List,BC_Corners_List
  TYPE(Link_Ptr_Type)  :: Link,Sublink
  TYPE(BC_Node_Ptr_Type)  :: BC_node,BC_corner

! For reordering:
  double precision, allocatable :: BcCoord(:,:)
  integer, allocatable :: Isort(:),Iback(:)
  double precision :: Lx,Lz

  bc%ngnod = grid%ngll
  allocate(bc%ibool(bc%ngnod,bc%nelem))

!-----------------------------------------------------------------------
!  Build the database as a linked list

  bc_npoin = 0
  call LI_Init_List(BC_Nodes_List)
  call LI_Init_List(BC_Corners_List)

  do n=1,bc%nelem

!   get edge nodes, counterclockwise by default
    call SE_get_edge_nodes(grid,bc%elem(n),bc%edge(n) &
                          ,bulk_node_vect)

    do kloc = 1,bc%ngnod

      bulk_node = bulk_node_vect(kloc)
      new_node = .true.

!  When the node is an element corner (vertex) 
!  check if it is already in the BC list.
!  Note that the search is done in a BC_corners sublist.

      node_at_corner = kloc==1 .or. kloc==bc%ngnod
      if (node_at_corner) then
        SubLink = LI_Get_Head(BC_Corners_List)
        do while(LI_Associated(SubLink))
          BC_corner = transfer(SubLink,BC_corner)
          if (BC_corner%P%Data%bulk_node == bulk_node) then
            new_node = .false.
            bc_inode = BC_corner%P%Data%bc_node
            exit
          endif
          Sublink = LI_Get_Next(Sublink)
        enddo
      endif

!  If it is a new node, add it to the list of BC nodes ...
      if (new_node) then 

        bc_npoin = bc_npoin +1
        bc_inode = bc_npoin
        allocate(BC_node%P); allocate(BC_node%P%data)
        BC_node%P%data%bulk_node = bulk_node
        BC_node%P%data%bc_node = bc_inode
        Link = transfer(BC_node,Link)
        call LI_Add_To_Head(Link,BC_Nodes_List)
         
!  ... and eventually to the sublist of BC corners.
        if (node_at_corner) then
          allocate(BC_corner%P)
          BC_corner%P%data => BC_node%P%data
          SubLink = transfer(BC_corner,SubLink)
          call LI_Add_To_Head(SubLink,BC_Corners_List)
        endif

      endif

     ! Set the [ (gll,bc_element) -> bc_node ] table
      bc%ibool(kloc,n) = bc_inode

    enddo
  enddo

!-----------------------------------------------------------------------
!  Translate BC database from linked list to array storage

  bc%npoin = bc_npoin ! = LI_Get_Len(BC_Nodes_List)
  allocate(bc%node(bc%npoin)) 
  do i=1,bc%npoin
    Link = LI_Remove_Head(BC_Nodes_List)
    BC_node = transfer(Link,BC_node)
    bc%node(BC_node%P%Data%bc_node) = BC_node%P%Data%bulk_node
    deallocate(BC_node%P)
  enddo

!-----------------------------------------------------------------------
!  Sort BC nodes by increasing coordinate
!  Use the coordinate with the largest range
  allocate(BcCoord(2,bc%npoin))
  allocate(Isort(bc%npoin))
  BcCoord = grid%coord(:,bc%node)
  Lx = maxval(BcCoord(1,:)) - minval(BcCoord(1,:))
  Lz = maxval(BcCoord(2,:)) - minval(BcCoord(2,:))
  if (Lx>Lz) then
    call drank(BcCoord(1,:),Isort)
  else
    call drank(BcCoord(2,:),Isort)
  endif
!    write(51,'(I5,1X,I3)') ( bc%node(n),Isort(n),n=1,bc%npoin) 
  bc%node = bc%node(Isort)
!    write(51,'(I5)') bc%node
  allocate(Iback(bc%npoin))
  Iback(Isort) = (/ (i, i=1,bc%npoin) /) 
  do n=1,bc%nelem
    bc%ibool(:,n) = Iback(bc%ibool(:,n))
  enddo

  deallocate(BcCoord,Isort,Iback)

end subroutine BC_set_bulk_node


!=======================================================================
!
  subroutine SE_inquire(grid,element,edge &
                       ,itab,jtab,dim_t &
                       ,size_min,size_max)

  type (sem_grid_type), intent(in) :: grid
  integer, optional, intent(in) :: element,edge
  double precision, optional, intent(out) :: size_min,size_max
  integer, dimension(grid%ngll), optional, intent(out) :: itab,jtab
  integer, optional, intent(out) :: dim_t

  integer :: k
  double precision :: dmin,dmax

  !-- give the GLL local numbering ITAB,JTAB for a given EDGE
  if ( present(edge) .and. present(itab) .and. present(jtab) ) then
    select case(edge)
      case(edge_D); itab = (/ (k, k=1,grid%ngll) /)   ; jtab = 1
      case(edge_R); itab = grid%ngll                  ; jtab = (/ (k, k=1,grid%ngll) /)
      case(edge_U); itab = (/ (k, k=grid%ngll,1,-1) /); jtab = grid%ngll
      case(edge_L); itab = 1                          ; jtab = (/ (k, k=grid%ngll,1,-1) /)
    end select
  endif

  !-- give local dimension (xi,eta) of tangent direction to an EDGE
  if ( present(dim_t) ) then
    select case(edge)
      case(edge_D,edge_U); dim_t = 1
      case(edge_R,edge_L); dim_t = 2
    end select
  endif

  if (present(size_min) .or. present(size_max)) then
    call FE_GetElementSizes(dmax,dmin,element,grid%fem)
    if (present(size_min)) size_min = dmin
    if (present(size_max)) size_max = dmax
  endif
    
  end subroutine SE_inquire

!=======================================================================
!
  subroutine BC_inquire(bounds,tag,bc_topo_ptr)

  type (bnd_grid_type), pointer :: bounds(:),bc_topo_ptr
  integer, intent(in) :: tag
  
  integer :: i
 
 ! if no boundary exists with the requested tag a null pointer is returned
  nullify(bc_topo_ptr) 
  do i=1,size(bounds)
    if ( bounds(i)%tag == tag ) then
      bc_topo_ptr => bounds(i)
      return
    endif
  enddo

  end subroutine BC_inquire

!=======================================================================
! Normal to a boundary, pointing out of the element
! Normals are assembled --> "average" normal between elements
  subroutine BC_get_normal_and_weights(bc,grid,NORM,W,periodic)

  type (bnd_grid_type), intent(inout) :: bc
  type (sem_grid_type), intent(in) :: grid
  double precision, intent(out) :: NORM(bc%npoin,2)
  double precision, intent(out) :: W(bc%npoin)
  logical, intent(in) :: periodic

  double precision :: DGlobDLoc(2,2),DxzDt(2),Tang(2),Jac1D,SignTang
  integer :: iGLLtab(bc%ngnod),jGLLtab(bc%ngnod)
  integer :: BcElmt,kGLL,BcNode,LocDimTanToEdge

  NORM = 0d0
  W = 0d0
  do BcElmt = 1,bc%nelem
    call SE_inquire(grid, edge=bc%edge(BcElmt) &
                   ,itab=iGLLtab, jtab=jGLLtab, dim_t=LocDimTanToEdge)
   ! SignTang enforces the counterclockwise convention for tangent vector
    if (bc%edge(BcElmt)==edge_U .OR. bc%edge(BcElmt)==edge_L) then
      SignTang = -1d0
    else
      SignTang = 1d0
    endif
    do kGLL = 1,bc%ngnod
      DGlobDLoc = SE_Jacobian(grid,bc%elem(BcElmt),iGLLtab(kGLL),jGLLtab(kGLL))
      DxzDt = DGlobDLoc(:,LocDimTanToEdge)
      Jac1D = sqrt( DxzDt(1)**2 + DxzDt(2)**2 )
      Tang = SignTang*DxzDt/Jac1D
      BcNode = bc%ibool(kGLL,BcElmt)
                            ! assembled, outwards normal
      NORM(BcNode,:) = NORM(BcNode,:) + (/ Tang(2),-Tang(1) /)
      W(BcNode) = W(BcNode) + grid%wgll(kGLL)*Jac1D
    enddo
  enddo
  if (periodic) then
    NORM(1,:) = NORM(1,:) + NORM(bc%npoin,:)
    NORM(bc%npoin,:) = NORM(1,:)
    W(1) = W(1) + W(bc%npoin)
    W(bc%npoin) = W(1)
  endif
  do BcElmt = 1,bc%nelem ! fix interelement normals:
    BcNode = bc%ibool(1,BcElmt)
    NORM(BcNode,:) = NORM(BcNode,:) / sqrt( NORM(BcNode,1)**2 + NORM(BcNode,2)**2 )
    BcNode = bc%ibool(bc%ngnod,BcElmt)
    NORM(BcNode,:) = NORM(BcNode,:) / sqrt( NORM(BcNode,1)**2 + NORM(BcNode,2)**2 )
  enddo
  
  end subroutine BC_get_normal_and_weights


end module spec_grid