fem_grid.f90

Spectral Element Method for wave propagation and rupture dynamics.

  if (.not.associated(fem%VertexConn)) call FE_SetConnectivity(fem)
  k = fem%knods(n,e)
  ee => fem%VertexConn(k)%elem
  nn => fem%VertexConn(k)%node

end subroutine FE_GetVertexConn

!=======================================================================
! Get a table of the extreme nodes (node1,node2) 
! for each of the four edges of an element
! in counterclockwise order
!
  function FE_GetEdgeNodesTab(grid,e)

  integer :: FE_GetEdgeNodesTab(2,4)
  type(fem_grid_type), intent(in) :: grid
  integer :: e

  FE_GetEdgeNodesTab(1,:) = grid%knods(EdgeKnod1,e)
  FE_GetEdgeNodesTab(2,:) = grid%knods(EdgeKnod2,e)

  end function FE_GetEdgeNodesTab

!=======================================================================
!
! example:
!
!  double precision pointer :: ptr(:,:)
!  ptr => FE_GetElementCoord(grid,e)
!  ...
!  deallocate(ptr)
!
  function FE_GetElementCoord(grid,e) result(coorg)

  type(fem_grid_type), intent(in) :: grid
  integer, intent(in) :: e
  double precision, pointer :: coorg(:,:)

  allocate(coorg(NDIME,grid%ngnod))
  coorg = grid%coord(:,grid%knods(:,e))

  end function FE_GetElementCoord


!=======================================================================
!
  function FE_GetNodesPerElement(grid) result(ngnod)
  type(fem_grid_type), intent(in) :: grid
  integer :: ngnod
  ngnod = grid%ngnod
  end function FE_GetNodesPerElement

!=======================================================================
!
  integer function FE_GetNbElements(grid)
  type(fem_grid_type), intent(in) :: grid
  FE_GetNbElements = grid%nelem
  end function FE_GetNbElements

!=======================================================================
!
  integer function FE_GetNbNodes(grid)
  type(fem_grid_type), intent(in) :: grid
  FE_GetNbNodes = grid%npoin
  end function FE_GetNbNodes


!=======================================================================
!
! GET SHAPE FUNCTIONS
! s,t = local coordinates [-1:1]*[-1:1]
!
  function FE_GetShape(xi,eta,ngnod) result(shape)

  double precision, intent(in) :: xi,eta
  integer, intent(in) :: ngnod
  double precision :: shape(ngnod)

  if ( ngnod == 4 ) then
    shape = Q4_GetShape(xi,eta)
  else
    shape = Q9_GetShape(xi,eta)
  endif  
  
  end function FE_GetShape
  

!=======================================================================
!
! dshape(i,j) = derivatives of shape function #i (associated to node #i)
!               with respect to coordinate #j
!
  function FE_GetDerShape(xi,eta,ngnod) result(dshape)

  double precision, intent(in) :: xi,eta
  integer, intent(in) :: ngnod
  double precision :: dshape(ngnod,2)

  if ( ngnod == 4 ) then
    dshape = Q4_GetDerShape(xi,eta)
  else
    dshape = Q9_GetDerShape(xi,eta)
  endif  

  end function FE_GetDerShape


!=======================================================================
!
  function FE_GetNbBoundaries(grid) result(nb)

  type(fem_grid_type), intent(in) :: grid
  integer :: nb

  if (associated(grid%bnds)) then
    nb = size(grid%bnds)
  else
    nb = 0           
  endif

  end function FE_GetNbBoundaries

!=======================================================================
!
  subroutine FE_InquireBoundary(grid,i,nelem,tag,elems,edges)

  type(fem_grid_type), intent(in) :: grid
  integer, intent(in) :: i
  integer, intent(out), optional :: nelem,tag
  integer, pointer, optional :: elems(:), edges(:)

  type(bnd_grid_type), pointer :: b

  b => grid%bnds(i)
  if (present(nelem)) nelem = b%nelem
  if (present(tag))   tag   = b%tag
  if (present(elems)) elems => b%elem
  if (present(edges)) edges => b%edge
  
  end subroutine FE_InquireBoundary

!=====================================================================
!
  subroutine FE_GetElementSizes(dmax,dmin,e,grid)

    type (fem_grid_type), intent(in) :: grid
    double precision, intent(out) :: dmax,dmin 
    integer, intent(in) :: e

    double precision, pointer :: coorg(:,:)
    double precision :: x0,z0,x1,z1,x2,z2,x3,z3,hsq(4)

    coorg => FE_GetElementCoord(grid,e)
    x0 = coorg(1,1)
    z0 = coorg(2,1)
    x1 = coorg(1,2)
    z1 = coorg(2,2)
    x2 = coorg(1,3)
    z2 = coorg(2,3)
    x3 = coorg(1,4)
    z3 = coorg(2,4)
    deallocate(coorg)

    hsq(1) = (x1-x0)**2 + (z1-z0)**2 
    hsq(2) = (x2-x1)**2 + (z2-z1)**2 
    hsq(3) = (x3-x2)**2 + (z3-z2)**2 
    hsq(4) = (x3-x0)**2 + (z3-z0)**2

    dmax = sqrt( maxval(hsq) )
    dmin = sqrt( minval(hsq) )

  end subroutine FE_GetElementSizes

!=====================================================================
!
  function FE_GetElementArea(grid,e) result(area)

  type(fem_grid_type), intent(in) :: grid
  integer, intent(in) :: e

  double precision :: area,a,b,c,d,p,q
  double precision, pointer :: coorg(:,:)

  coorg => FE_GetElementCoord(grid,e)
  a = (coorg(1,2)-coorg(1,1))**2+(coorg(2,2)-coorg(2,1))**2
  b = (coorg(1,3)-coorg(1,2))**2+(coorg(2,3)-coorg(2,2))**2
  c = (coorg(1,4)-coorg(1,3))**2+(coorg(2,4)-coorg(2,3))**2
  d = (coorg(1,1)-coorg(1,4))**2+(coorg(2,1)-coorg(2,4))**2
  p = (coorg(1,1)-coorg(1,3))**2+(coorg(2,1)-coorg(2,3))**2
  q = (coorg(1,2)-coorg(1,4))**2+(coorg(2,2)-coorg(2,4))**2
  area = 0.25d0*sqrt( 4d0*p*q - (b+d-a-c)**2 )
  deallocate(coorg)

  end function FE_GetElementArea
  
!=====================================================================
!
  subroutine FE_find_point(coord,grid,e,xi,eta,new_coord,istatus)

  use utils, only : invert2
  use stdio, only : IO_abort

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

  integer, parameter :: ntrial=100

  double precision :: esize,tolf,tolx, x(NDIME),p(NDIME), &
                      fvec(NDIME),fjac(NDIME,NDIME)
  double precision , pointer :: coorg(:,:)
  double precision :: shap(grid%ngnod), dshap(grid%ngnod,NDIME)
  integer :: n

  coorg => FE_GetElementCoord(grid,e)
  esize = sqrt( FE_GetElementArea(grid,e)/PI )

  x(1) = xi
  x(2) = eta
  tolx=TINY_XABS
  tolf=TINY_XABS*esize
  p=0d0

  do n=1,ntrial

    shap = FE_getshape(x(1),x(2),grid%ngnod)
    fvec = matmul(coorg, shap) - coord

    dshap = FE_getdershape(x(1),x(2),grid%ngnod)
    fjac = matmul(coorg, dshap)

    if (sum(abs(fvec)) <= tolf) exit
    p=-fvec
    p = matmul(invert2(fjac),p)
    x=x+p
    if (sum(abs(p)) <= tolx) exit

  end do

  if (n>=ntrial) call IO_abort('FE_find_point: did not converge')
  if ( abs(xi)<1.d0+TINY_XABS .and. abs(eta)<1.d0+TINY_XABS ) then
    istatus = 0
  else
    istatus = 1
  endif
  xi=x(1)
  eta=x(2)
  shap = FE_getshape(xi,eta,grid%ngnod)
  new_coord = matmul(coorg, shap)
  deallocate(coorg)

  end subroutine FE_find_point

!=====================================================================
!
! reorder grid data by element permutation
  subroutine FE_reorder(grid,perm,perm_inv)

  type(fem_grid_type), intent(inout) :: grid
  integer :: perm(grid%nelem), perm_inv(grid%nelem)

  integer :: k

  grid%knods = grid%knods(:,perm)
  grid%tag = grid%tag(perm)
  do k=1,size(grid%bnds)
    grid%bnds(k)%elem=perm_inv( grid%bnds(k)%elem )
  enddo

  end subroutine FE_reorder

end module fem_grid