fem_grid.f90

Spectral Element Method for wave propagation and rupture dynamics.

! SEM2DPACK version 2.2.11 -- A Spectral Element Method for 2D wave propagation and fracture dynamics,
!                             with emphasis on computational seismology and earthquake source dynamics.
! 
! Copyright (C) 2003-2007 Jean-Paul Ampuero
! All Rights Reserved
! 
! Jean-Paul Ampuero
! 
! ETH Zurich (Swiss Federal Institute of Technology)
! Institute of Geophysics
! Seismology and Geodynamics Group
! ETH H鰊ggerberg HPP O 13.1
! CH-8093 Z黵ich
! Switzerland
! 
! [e-mail]
! [phone] (office)
! [phone] (fax)
! 
! http://www.sg.geophys.ethz.ch/geodynamics/ampuero/
! 
! 
! This software is freely available for scientific research purposes. 
! If you use this software in writing scientific papers include proper 
! attributions to its author, Jean-Paul Ampuero.
! 
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License
! as published by the Free Software Foundation; either version 2
! of the License, or (at your option) any later version.
! 
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
! 
! You should have received a copy of the GNU General Public License
! along with this program; if not, write to the Free Software
! Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
! 
module fem_grid

!=======================================================================
!
!                       This module deals with quadrangle elements
!                       defined by 4 or 9 control nodes.
!                       The control nodes are defined as follows:
!
!                                       edge 3
!
!                               4 . . . . 7 . . . . 3
!                               .                   .
!                               .         t         .
!                               .                   .
!                   edge 4      8         9  s      6     edge 2
!                               .                   .
!                               .                   .
!                               .                   .
!                               1 . . . . 5 . . . . 2
!
!                                       edge 1
!
!                           Local coordinate system : s,t
!
!=======================================================================

  use elem_q4
  use elem_q9
  use bnd_grid
  use constants

  implicit none
  private

  type VertexConn_type
    integer, pointer :: elem(:)=>null(), node(:)=>null()
  end type VertexConn_type

!---- Finite element mesh (Q4 or Q9)
!     npoin = total number of control nodes
!     ngnod = number of control nodes per element (4 or 9)
!     nelem = total number of elements
!     coord = coordinates of the control nodes
!     knods = local to global numbering for the control nodes
!             (ignod,element) -> control node index
!     tag   = element tags
!     ref   = node tags
!     bnds  = list of domain boundaries
!     flat  = flag for cartesian grid (allows simplifications)
!
! vertex connectivity data: (in Q9 also includes internal nodes)
!     VertexConn(k)%elem(:) = elements sharing node k
!     VertexConn(k)%node(:) = local index of node k in each element
!
! edge connectivity data :
!     EdgeConn_elem(n,e) = element that shares edge#n of element#e
!     EdgeConn_edge(n,e) = edge# on matching element

  type fem_grid_type
    integer :: npoin=0, ngnod=0, nelem=0
    double precision, pointer :: coord(:,:) =>null()
    integer, pointer :: knods(:,:) =>null(), &
                        tag(:)     =>null(), &
                        ref(:)     =>null()
    type(bnd_grid_type), pointer :: bnds(:) =>null()
    logical :: flat = .false.
    integer, dimension(:,:), pointer :: EdgeConn_elem =>null(), &
                                        EdgeConn_edge =>null()
    type(VertexConn_type), pointer :: VertexConn(:) =>null()
  end type fem_grid_type

  !-- element edges tags Up, Down, Left, Right
  integer, parameter :: edge_D = 1
  integer, parameter :: edge_R = 2
  integer, parameter :: edge_U = 3
  integer, parameter :: edge_L = 4

  !-- control nodes defining the element edges
  integer, dimension(4), parameter :: EdgeKnod1 = (/  1, 2, 3, 4 /) &
                                     ,EdgeKnod2   = (/  2, 3, 4, 1 /)

  !-- domain side tags Up, Down, Left, Right
  integer, parameter :: side_D = 1
  integer, parameter :: side_R = 2
  integer, parameter :: side_U = 3
  integer, parameter :: side_L = 4

  public :: fem_grid_type &
          , FE_GetShape, FE_GetDerShape, FE_GetElementCoord &
          , FE_GetEdgeNodesTab, FE_GetNodesPerElement, FE_GetNbElements &
          , FE_GetNbBoundaries, FE_GetNbNodes, FE_InquireBoundary &
          , FE_SetConnectivity, FE_UnsetConnectivity, FE_GetEdgeConn, FE_GetVertexConn &
          , FE_find_point, FE_GetElementSizes, FE_reorder &
          , edge_D,edge_R,edge_U,edge_L &
          , side_D,side_R,side_U,side_L

contains

!=======================================================================
!

subroutine FE_SetConnectivity(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,LI_Get_Len

  type(fem_grid_type), intent(inout) :: grid

  type Vertex_Data_Type
    integer :: elem=0,node=0
  end type Vertex_Data_Type

  type Vertex_Type
    type(Link_Type) :: Link
    type(Vertex_Data_Type) :: Data
  end type Vertex_Type

  type Vertex_Ptr_Type
    type(Vertex_Type), pointer :: P
  end type Vertex_Ptr_Type

  type(List_Type)     :: Vertex_List(grid%npoin)
  type(Link_Ptr_Type) :: Link, Link2
  type(Vertex_Ptr_Type) :: Vertex, Vertex2

  integer :: i,j,k,e,n,ii,jj,ee,nn,k1,k2,n1,n2

 ! vertex data
  do k = 1,grid%npoin
    call LI_Init_List(Vertex_List(k))
  enddo
  do e = 1,grid%nelem
  do n = 1,grid%ngnod
    allocate( Vertex%P ) ! reallocation
    Vertex%P%Data%elem = e
    Vertex%P%Data%node = n
    Link = transfer(Vertex,Link)
    call LI_Add_To_Head(Link,Vertex_List(grid%knods(n,e)))
  enddo
  enddo

  ! convert the vertex linked list into arrays, for easier handling
  allocate(grid%VertexConn(grid%npoin))
  do k = 1,grid%npoin
    n = LI_Get_Len(Vertex_List(k))
    if (n==0) cycle  ! center point in Q9 element
    allocate(grid%VertexConn(k)%elem(n))
    allocate(grid%VertexConn(k)%node(n))
    n=0
    Link = LI_Remove_Head( Vertex_List(k) )
    do while ( LI_Associated(Link) )
      n = n+1
      Vertex = transfer(Link,Vertex)
      grid%VertexConn(k)%elem(n) = Vertex%P%Data%elem
      grid%VertexConn(k)%node(n) = Vertex%P%Data%node
      deallocate( Vertex%P )
      Link = LI_Remove_Head( Vertex_List(k) )
    enddo
  enddo

 ! edge data :
  allocate(grid%EdgeConn_elem(4,grid%nelem)) 
  allocate(grid%EdgeConn_edge(4,grid%nelem)) 
  grid%EdgeConn_elem = 0
  grid%EdgeConn_edge = 0
  do e = 1,grid%nelem
  do n = 1,4
    if ( grid%EdgeConn_elem(n,e)>0 ) cycle ! skip if already processed
    k1 = grid%knods(EdgeKnod1(n),e) 
    k2 = grid%knods(EdgeKnod2(n),e) 
   ! search another element (ee) shared by the two vertices of the current edge 
   ! ( Vertex_List(k1) and Vertex_List(k2) )
    nn = 0
    do n1=1,size(grid%VertexConn(k1)%elem)
      ! the fist vertex also belongs to element ee
      ee = grid%VertexConn(k1)%elem(n1)
      if (ee == e) cycle
      ! check if the second vertex also belongs to element ee
      ! if so: we found the matching edge
      do n2=1,size(grid%VertexConn(k2)%elem) 
        if (grid%VertexConn(k2)%elem(n2)==ee) then
         ! nn = edge# = first_control_node# 
         ! + because of the counterclockwise node numbering convention
         ! the second node of the current edge  
         ! is matched by the first node of the opposite edge:
          nn = grid%VertexConn(k2)%node(n2)
          grid%EdgeConn_elem(n,e) = ee
          grid%EdgeConn_edge(n,e) = nn
          grid%EdgeConn_elem(nn,ee) = e
          grid%EdgeConn_edge(nn,ee) = n
          exit
        endif
      enddo
      if (nn>0) exit
    enddo
  enddo
  enddo
   

end subroutine FE_SetConnectivity

!-----------------------------------------------------------------------
subroutine FE_UnsetConnectivity(grid)

  type(fem_grid_type), intent(inout) :: grid
  integer :: k

  deallocate(grid%EdgeConn_elem)
  deallocate(grid%EdgeConn_edge)
  do k=1,grid%npoin
    deallocate(grid%VertexConn(k)%elem)
    deallocate(grid%VertexConn(k)%node)
  enddo
  deallocate(grid%VertexConn)

end subroutine FE_UnsetConnectivity

!-----------------------------------------------------------------------
subroutine FE_GetEdgeConn(grid,e,n,ee,nn)

  type(fem_grid_type), intent(inout) :: grid
  integer, intent(in) :: e,n
  integer, intent(out) :: ee,nn

  if (.not.associated(grid%EdgeConn_elem)) call FE_SetConnectivity(grid)
  ee=grid%EdgeConn_elem(n,e)
  nn=grid%EdgeConn_edge(n,e)

end subroutine FE_GetEdgeConn

!-----------------------------------------------------------------------
subroutine FE_GetVertexConn(fem,e,n,ee,nn)

  type(fem_grid_type), intent(inout) :: fem
  integer, intent(in) :: e,n
  integer, pointer :: ee(:),nn(:)

  integer :: k