element.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
!       If you wish introduce new elements with new geometry or/and
!    new shape functions, you need do the following
!!
!!
!!-    Introduce new element ID corresponding to your element in this module.
!!-    Provide corresponding shape function and shape derivative in a new
!!    module and include this new module here.
!!-    Add the information of your new element into functions listed above.
!!-    If new quadrature method needed, do some modification in MODULE
!!    Quadrature.
!!-    If you introduce new surface geometry and do contact calculation also,
!!    you may need do some modification in MODULE mSurfElement.
module elementinfo
 
  use shape_line2n
  use shape_line3n
  use shape_tri3n
  use shape_tri6n
  use shape_quad4n
  use shape_quad8n
  use shape_quad9n
  use shape_hex8n
  use shape_hex20n
  use shape_tet4n
  use shape_tet10n
  use shape_prism6n
  use shape_prism15n
  implicit none
 
  integer, parameter, private :: kreal = kind(0.0d0)
 
  !-------------------------------------------
  !     Fllowing ID of element types
  !-------------------------------------------
  integer, parameter :: fe_unknown  = -1
 
  integer, parameter :: fe_line2n   = 111
  integer, parameter :: fe_line3n   = 112
  integer, parameter :: fe_tri3n    = 231
  integer, parameter :: fe_tri6n    = 232
  integer, parameter :: fe_tri6nc   = 2322
  integer, parameter :: fe_quad4n   = 241
  integer, parameter :: fe_quad8n   = 242
  integer, parameter :: fe_truss    = 301
  integer, parameter :: fe_tet4n    = 341
  integer, parameter :: fe_tet4n_pipi = 3414
  integer, parameter :: fe_tet10n   = 342
  integer, parameter :: fe_tet10nc  = 3422
  integer, parameter :: fe_prism6n  = 351
  integer, parameter :: fe_prism15n = 352
  integer, parameter :: fe_hex8n    = 361
  integer, parameter :: fe_hex20n   = 362
  integer, parameter :: fe_hex27n   = 363
 
  integer, parameter :: fe_if_line2n   = 511
 
  integer, parameter :: fe_beam2n   = 611
  integer, parameter :: fe_beam3n   = 612
  integer, parameter :: fe_beam341  = 641
 
  integer, parameter :: fe_tri6n_shell  = 732
  integer, parameter :: fe_dsg3_shell   = 733
  integer, parameter :: fe_mitc3_shell  = 731
  integer, parameter :: fe_mitc4_shell  = 741
  integer, parameter :: fe_mitc8_shell  = 742
  integer, parameter :: fe_mitc9_shell  = 743
 
  integer, parameter :: fe_mitc3_shell361  = 761
  integer, parameter :: fe_mitc4_shell361  = 781
 
  integer, parameter :: fe_nodesmooth_tet4n = 881
  integer, parameter :: fe_edgesmooth_tet4n = 891
 
  integer, parameter :: fe_tri3n_patch    = 1031
  integer, parameter :: fe_tri6n_patch    = 1032
  integer, parameter :: fe_quad4n_patch   = 1041
  integer, parameter :: fe_quad8n_patch   = 1042
  ! ---------------------------------------------
 
contains
 
  !************************************
  !    Following geometric information
  !************************************
  integer(kind=kind(2)) function getspacedimension( etype )
    integer, intent(in) :: etype
 
    select case( etype)
      case (fe_line2n, fe_line3n, fe_if_line2n)
        getspacedimension = 1
      case (fe_tri3n, fe_tri6n, fe_tri6nc, fe_quad4n, fe_quad8n)
        getspacedimension = 2
      case default
        getspacedimension = 3
    end select
  end function
 
  integer(kind=kind(2)) function getnumberofnodes( etype )
    integer, intent(in) :: etype
 
    select case (etype)
      case (fe_line2n, fe_beam2n, fe_if_line2n)
        getnumberofnodes = 2
      case (fe_line3n, fe_beam3n)
        getnumberofnodes = 3
      case (fe_tri3n, fe_mitc3_shell, fe_mitc3_shell361, fe_tri3n_patch )
        getnumberofnodes = 3
      case ( fe_tri6n, fe_tri6nc, fe_tri6n_shell, fe_tri6n_patch )
        getnumberofnodes = 6
      case ( fe_quad4n, fe_mitc4_shell, fe_mitc4_shell361, fe_quad4n_patch )
        getnumberofnodes = 4
      case ( fe_quad8n, fe_mitc8_shell, fe_quad8n_patch)
        getnumberofnodes = 8
      case ( fe_mitc9_shell )
        getnumberofnodes = 9
      case ( fe_tet4n, fe_tet4n_pipi, fe_beam341 )
        getnumberofnodes = 4
      case ( fe_tet10n, fe_tet10nc )
        getnumberofnodes = 10
      case ( fe_prism6n )
        getnumberofnodes = 6
      case ( fe_prism15n )
        getnumberofnodes = 15
      case ( fe_hex8n )
        getnumberofnodes = 8
      case ( fe_hex20n )
        getnumberofnodes = 20
      case default
        getnumberofnodes = -1
        ! error message
    end select
  end function
 
  integer(kind=kind(2)) function getnumberofsubface( etype )
    integer, intent(in) :: etype
 
    select case (etype)
      case (fe_line2n, fe_line3n, fe_if_line2n)
        getnumberofsubface = 2
      case (fe_tri3n, fe_mitc3_shell, fe_tri6n, fe_tri6nc, fe_tri6n_shell, fe_mitc3_shell361  )
        getnumberofsubface = 3
      case ( fe_quad4n, fe_mitc4_shell, fe_quad8n, fe_mitc8_shell, fe_mitc9_shell, fe_mitc4_shell361  )
        getnumberofsubface = 4
      case ( fe_tet4n, fe_tet4n_pipi, fe_tet10n, fe_tet10nc, fe_beam341 )
        getnumberofsubface = 4
      case ( fe_prism6n, fe_prism15n )
        getnumberofsubface = 5
      case ( fe_hex8n, fe_hex20n)
        getnumberofsubface = 6
      case ( fe_tri3n_patch, fe_tri6n_patch, fe_quad4n_patch, fe_quad8n_patch )
        getnumberofsubface = 1
      case default
        getnumberofsubface = -1
        ! error message
    end select
  end function
 
  subroutine getsubface( intype, innumber, outtype, nodes )
    integer, intent(in)    :: intype
    integer, intent(in)    :: innumber
    integer, intent(out)   :: outtype
    integer, intent(out)   :: nodes(:)
 
    if( innumber>getnumberofsubface( intype ) ) stop "Error in getting subface"
    select case ( intype )
      case (fe_tet4n, fe_tet4n_pipi, fe_beam341)
        outtype = fe_tri3n
        select case ( innumber )
          case (1)
            nodes(1)=1; nodes(2)=2; nodes(3)=3
          case (2)
            nodes(1)=4; nodes(2)=2; nodes(3)=1
          case (3)
            nodes(1)=4; nodes(2)=3; nodes(3)=2
          case (4)
            nodes(1)=4; nodes(2)=1; nodes(3)=3
        end select
      case (fe_tet10n)
        outtype = fe_tri6n
        select case ( innumber )
          case (1)
            nodes(1)=1; nodes(2)=2; nodes(3)=3
            nodes(4)=5; nodes(5)=6; nodes(6)=7
          case (2)
            nodes(1)=4; nodes(2)=2; nodes(3)=1
            nodes(4)=9; nodes(5)=5; nodes(6)=8
          case (3)
            nodes(1)=4;  nodes(2)=3; nodes(3)=2
            nodes(4)=10; nodes(5)=6; nodes(6)=9
          case (4)
            nodes(1)=4; nodes(2)=1; nodes(3)=3
            nodes(4)=8; nodes(5)=7; nodes(6)=10
        end select
      case (fe_tet10nc)
        outtype = fe_tri6nc
        select case ( innumber )
          case (1)
            nodes(1)=1; nodes(2)=2; nodes(3)=3
            nodes(4)=5; nodes(5)=6; nodes(6)=7
          case (2)
            nodes(1)=4; nodes(2)=2; nodes(3)=1
            nodes(4)=9; nodes(5)=5; nodes(6)=8
          case (3)
            nodes(1)=4;  nodes(2)=3; nodes(3)=2
            nodes(4)=10; nodes(5)=6; nodes(6)=9
          case (4)
            nodes(1)=4; nodes(2)=1; nodes(3)=3
            nodes(4)=8; nodes(5)=7; nodes(6)=10
        end select
      case ( fe_hex8n )
        outtype = fe_quad4n
        select case ( innumber )
          case (1)
            nodes(1)=1; nodes(2)=2; nodes(3)=3; nodes(4)=4
          case (2)
            nodes(1)=8; nodes(2)=7; nodes(3)=6; nodes(4)=5
          case (3)
            nodes(1)=5; nodes(2)=6; nodes(3)=2; nodes(4)=1
          case (4)
            nodes(1)=6; nodes(2)=7; nodes(3)=3; nodes(4)=2
          case (5)
            nodes(1)=7; nodes(2)=8; nodes(3)=4; nodes(4)=3
          case (6)
            nodes(1)=8; nodes(2)=5; nodes(3)=1; nodes(4)=4
          case default
            ! error
        end select
      case (fe_hex20n)
        outtype = fe_quad8n
        select case ( innumber )
          case (1)
            nodes(1)=1;  nodes(2)=2;  nodes(3)=3;  nodes(4)=4
            nodes(5)=9;  nodes(6)=10; nodes(7)=11; nodes(8)=12
          case (2)
            nodes(1)=8;  nodes(2)=7;  nodes(3)=6;  nodes(4)=5
            nodes(5)=15; nodes(6)=14; nodes(7)=13; nodes(8)=16
          case (3)
            nodes(1)=5;  nodes(2)=6;  nodes(3)=2;  nodes(4)=1
            nodes(5)=13; nodes(6)=18; nodes(7)=9;  nodes(8)=17
          case (4)
            nodes(1)=6;  nodes(2)=7;  nodes(3)=3;  nodes(4)=2
            nodes(5)=14; nodes(6)=19; nodes(7)=10; nodes(8)=18
          case (5)
            nodes(1)=7;  nodes(2)=8;  nodes(3)=4;  nodes(4)=3
            nodes(5)=15; nodes(6)=20; nodes(7)=11; nodes(8)=19
          case (6)
            nodes(1)=8;  nodes(2)=5;  nodes(3)=1;  nodes(4)=4
            nodes(5)=16; nodes(6)=17; nodes(7)=12; nodes(8)=20
          case default
            ! error
        end select
      case (fe_prism6n)
        select case ( innumber )
          case (1)
            outtype = fe_tri3n
            nodes(1)=1; nodes(2)=2; nodes(3)=3
          case (2)
            outtype = fe_tri3n
            nodes(1)=6; nodes(2)=5; nodes(3)=4
          case (3)
            outtype = fe_quad4n
            nodes(1)=4; nodes(2)=5; nodes(3)=2; nodes(4)=1
          case (4)
            outtype = fe_quad4n
            nodes(1)=5; nodes(2)=6; nodes(3)=3; nodes(4)=2
          case (5)
            outtype = fe_quad4n
            nodes(1)=6; nodes(2)=4; nodes(3)=1; nodes(4)=3
        end select
      case (fe_prism15n)
        select case ( innumber )
          case (1)
            outtype = fe_tri6n
            nodes(1)=1; nodes(2)=2; nodes(3)=3
            nodes(4)=7; nodes(5)=8; nodes(6)=9
          case (2)
            outtype = fe_tri6n
            nodes(1)=6;  nodes(2)=5;  nodes(3)=4
            nodes(4)=11; nodes(5)=10; nodes(6)=12
          case (3)
            outtype = fe_quad8n
            nodes(1)=4;  nodes(2)=5;  nodes(3)=2;  nodes(4)=1
            nodes(5)=10; nodes(6)=14; nodes(7)=7;  nodes(8)=13
          case (4)
            outtype = fe_quad8n
            nodes(1)=5;  nodes(2)=6;  nodes(3)=3;  nodes(4)=2
            nodes(5)=11; nodes(6)=15; nodes(7)=8;  nodes(8)=14
          case (5)
            outtype = fe_quad8n
            nodes(1)=6;  nodes(2)=4;  nodes(3)=1;  nodes(4)=3
            nodes(5)=12; nodes(6)=13; nodes(7)=9;  nodes(8)=15
        end select
      case ( fe_tri3n, fe_mitc3_shell )
        outtype = fe_line2n
        select case (innumber )
          case (1)
            nodes(1) = 1; nodes(2)=2
          case (2)
            nodes(1) = 2; nodes(2)=3
          case (3)
            nodes(1) = 3; nodes(2)=1
        end select
      case ( fe_tri6n, fe_tri6nc, fe_tri6n_shell )
        outtype = fe_line3n
        select case (innumber )
          case (1)
            nodes(1) = 1; nodes(2)=2;  nodes(3)=4
          case (2)
            nodes(1) = 2; nodes(2)=3;  nodes(3)=5
          case (3)
            nodes(1) = 3; nodes(2)=1;  nodes(3)=6
        end select
      case ( fe_quad4n, fe_mitc4_shell )
        outtype = fe_line2n
        select case (innumber )
          case (1)
            nodes(1) = 1; nodes(2)=2
          case (2)
            nodes(1) = 2; nodes(2)=3
          case (3)
            nodes(1) = 3; nodes(2)=4
          case (4)
            nodes(1) = 4; nodes(2)=1
        end select
      case ( fe_quad8n, fe_mitc8_shell, fe_mitc9_shell )
        outtype = fe_line3n
        select case (innumber )
          case (1)
            nodes(1) = 1; nodes(2)=2;  nodes(3)=5
          case (2)
            nodes(1) = 2; nodes(2)=3;  nodes(3)=6
          case (3)
            nodes(1) = 3; nodes(2)=4;  nodes(3)=7
          case (4)
            nodes(1) = 4; nodes(2)=1;  nodes(3)=8
        end select
      case (fe_mitc3_shell361)
        select case ( innumber )
          case (1)
            outtype = fe_tri3n
            nodes(1)=1; nodes(2)=2; nodes(3)=3
          case (2)
            outtype = fe_tri3n
            nodes(1)=6; nodes(2)=5; nodes(3)=4
          case (3)
            outtype = fe_quad4n
            nodes(1)=4; nodes(2)=5; nodes(3)=2; nodes(4)=1
          case (4)
            outtype = fe_quad4n
            nodes(1)=5; nodes(2)=6; nodes(3)=3; nodes(4)=2
          case (5)
            outtype = fe_quad4n
            nodes(1)=6; nodes(2)=4; nodes(3)=1; nodes(4)=3
        end select
      case ( fe_mitc4_shell361 )
        outtype = fe_quad4n
        select case ( innumber )
          case (1)
            nodes(1)=1; nodes(2)=2; nodes(3)=3; nodes(4)=4
          case (2)
            nodes(1)=8; nodes(2)=7; nodes(3)=6; nodes(4)=5
          case (3)
            nodes(1)=5; nodes(2)=6; nodes(3)=2; nodes(4)=1
          case (4)
            nodes(1)=6; nodes(2)=7; nodes(3)=3; nodes(4)=2
          case (5)
            nodes(1)=7; nodes(2)=8; nodes(3)=4; nodes(4)=3
          case (6)
            nodes(1)=8; nodes(2)=5; nodes(3)=1; nodes(4)=4
          case default
            ! error
        end select
      case ( fe_tri3n_patch )
        outtype = fe_tri3n
        select case ( innumber )
          case (1)
            nodes(1)=1; nodes(2)=2; nodes(3)=3
          case default
            !error
        end select
      case ( fe_tri6n_patch )
        outtype = fe_tri6n
        select case ( innumber )
          case (1)
            nodes(1)=1; nodes(2)=2; nodes(3)=3
            nodes(4)=4; nodes(5)=5; nodes(6)=6
          case default
            !error
        end select
      case ( fe_quad4n_patch )
        outtype = fe_quad4n
        select case ( innumber )
          case (1)
            nodes(1)=1; nodes(2)=2; nodes(3)=3; nodes(4)=4
          case default
            !error
        end select
      case ( fe_quad8n_patch )
        outtype = fe_quad8n
        select case ( innumber )
          case (1)
            nodes(1)=1; nodes(2)=2; nodes(3)=3; nodes(4)=4
            nodes(5)=5; nodes(6)=6; nodes(7)=7; nodes(8)=8
          case default
            !error
        end select
      case default
        outtype = fe_unknown
        stop "element type not defined-sbs"
        ! error message
    end select
  end subroutine
 
  integer function numofquadpoints( fetype )
    integer, intent(in) :: fetype
    select case (fetype)
      case (fe_line2n, fe_tri3n, fe_tet4n, fe_beam2n , fe_beam341, fe_truss, fe_if_line2n )
        numofquadpoints = 1
      case ( fe_tri6n )
        numofquadpoints = 3
      case ( fe_tri6nc )
        numofquadpoints = 4
      case (fe_line3n )
        numofquadpoints = 2
      case ( fe_quad4n, fe_mitc4_shell, fe_mitc4_shell361 )
        numofquadpoints = 4
      case ( fe_quad8n, fe_mitc9_shell )
        numofquadpoints = 9
      case ( fe_hex8n )
        numofquadpoints = 8
      case ( fe_hex20n, fe_mitc8_shell )
        numofquadpoints = 27
      case ( fe_prism6n )
        numofquadpoints = 2
      case ( fe_mitc3_shell, fe_mitc3_shell361 )
        numofquadpoints = 3
      case ( fe_prism15n, fe_tri6n_shell )
        numofquadpoints = 9
      case ( fe_tet10n, fe_tet4n_pipi )
        numofquadpoints = 4
      case ( fe_tet10nc )
        numofquadpoints = 12
      case ( fe_nodesmooth_tet4n, fe_edgesmooth_tet4n )
        numofquadpoints = 1
      case default
        numofquadpoints = -1
        ! error message
        stop "element type not defined-np"
    end select
  end function
 
  integer function numofshellthicknessquadpoints( etype )
    integer, intent(in) :: etype
 
    select case (etype)
      case (fe_mitc3_shell, fe_mitc4_shell, fe_mitc3_shell361, fe_mitc4_shell361)
        numofshellthicknessquadpoints = 2
      case (fe_mitc9_shell)
        numofshellthicknessquadpoints = 3
      case default
        numofshellthicknessquadpoints = -1
    end select
  end function numofshellthicknessquadpoints
 
  subroutine getshellthicknessquadpoint( etype, np, zeta )
    use quadrature, only: gauss1d2, gauss1d3
    integer, intent(in) :: etype, np
    real(kind=kreal), intent(out) :: zeta
 
    if (np < 1 .or. np > numofshellthicknessquadpoints(etype)) then
      stop "invalid shell thickness quadrature point"
    endif
 
    select case (etype)
      case (fe_mitc3_shell, fe_mitc4_shell, fe_mitc3_shell361, fe_mitc4_shell361)
        zeta = gauss1d2(1, np)
      case (fe_mitc9_shell)
        zeta = gauss1d3(1, np)
    end select
  end subroutine getshellthicknessquadpoint
 
  real(kind=kreal) function getshellthicknessweight( etype, np )
    use quadrature, only: weight1d2, weight1d3
    integer, intent(in) :: etype, np
 
    if (np < 1 .or. np > numofshellthicknessquadpoints(etype)) then
      stop "invalid shell thickness quadrature point"
    endif
 
    select case (etype)
      case (fe_mitc3_shell, fe_mitc4_shell, fe_mitc3_shell361, fe_mitc4_shell361)
        getshellthicknessweight = weight1d2(np)
      case (fe_mitc9_shell)
        getshellthicknessweight = weight1d3(np)
    end select
  end function getshellthicknessweight
 
  subroutine getquadpoint( fetype, np, pos )
    use quadrature
    integer, intent(in)           :: fetype
    integer, intent(in)           :: np
    real(kind=kreal), intent(out) :: pos(:)    
 
    if( np<1 .or. np>numofquadpoints(fetype) ) then
      ! error
    endif
 
    select case (fetype)
      case (fe_tri3n)
        pos(1:2)=gauss2d4(:,np)
      case ( fe_tri6n, fe_mitc3_shell )
        pos(1:2)=gauss2d5(:,np)
      case (fe_tri6nc )
        pos(1:2)=gauss2d6(:,np)
      case ( fe_quad4n, fe_mitc4_shell )
        pos(1:2)=gauss2d2(:,np)
      case ( fe_quad8n, fe_mitc9_shell )
        pos(1:2)=gauss2d3(:,np)
      case ( fe_hex8n, fe_mitc4_shell361 )
        pos(1:3)=gauss3d2(:,np)
      case ( fe_hex20n, fe_mitc8_shell )
        pos(1:3)=gauss3d3(:,np)
      case ( fe_prism6n, fe_mitc3_shell361 )
        pos(1:3)=gauss3d7(:,np)
      case ( fe_prism15n, fe_tri6n_shell )
        pos(1:3)=gauss3d8(:,np)
      case ( fe_tet4n, fe_beam341 )
        pos(1:3)=gauss3d4(:,np)
      case ( fe_tet10n, fe_tet4n_pipi )
        pos(1:3)=gauss3d5(:,np)
      case ( fe_tet10nc )
        pos(1:3)=np
      case ( fe_line2n, fe_if_line2n )
        pos(1:1)=gauss1d1(:,np)
      case ( fe_line3n )
        pos(1:1)=gauss1d2(:,np)
      case default
        ! error message
        stop "element type not defined-qp"
    end select
  end subroutine
 
  real(kind=kreal) function getweight( fetype, np )
    use quadrature
    integer, intent(in)           :: fetype
    integer, intent(in)           :: np
    if( np<1 .or. np>numofquadpoints(fetype) ) then
      ! error
    endif
 
    select case (fetype)
      case (fe_tri3n)
        getweight = weight2d4(1)
      case ( fe_tri6n, fe_mitc3_shell )
        getweight = weight2d5(np)
      case ( fe_quad4n, fe_mitc4_shell )
        getweight = weight2d2(np)
      case ( fe_quad8n, fe_mitc9_shell )
        getweight = weight2d3(np)
      case ( fe_hex8n, fe_mitc4_shell361 )
        getweight = weight3d2(np)
      case ( fe_hex20n)
        getweight = weight3d3(np)
      case ( fe_prism6n, fe_mitc3_shell361 )
        getweight = weight3d7(np)
      case ( fe_prism15n )
        getweight = weight3d8(np)
      case ( fe_tet4n, fe_beam341 )
        getweight = weight3d4(1)
      case ( fe_tet10n, fe_tet4n_pipi )
        getweight = weight3d5(np)
      case ( fe_line2n, fe_if_line2n )
        getweight = weight1d1(1)
      case ( fe_line3n )
        getweight = weight1d2(np)
      case default
        getweight = 0.d0
        ! error message
    end select
  end function
 
  !************************************
  !    Following shape function information
  !************************************
  subroutine getshapederiv( fetype, localcoord, shapederiv )
    integer, intent(in)           :: fetype
    real(kind=kreal), intent(in)  :: localcoord(:)      
    real(kind=kreal), intent(out) :: shapederiv(:,:)    
 
    select case (fetype)
      case ( fe_tri3n, fe_mitc3_shell )
        !error check
        call shapederiv_tri3n(shapederiv(1:3,1:2))
      case (fe_tri6n)
        !error check
        call shapederiv_tri6n(localcoord,shapederiv(1:6,1:2) )
      case ( fe_quad4n, fe_mitc4_shell )
        !error check
        call shapederiv_quad4n(localcoord,shapederiv(1:4,1:2))
      case (fe_quad8n)
        !error check
        call shapederiv_quad8n(localcoord,shapederiv(1:8,1:2))
      case ( fe_mitc9_shell )
        !error check
        call shapederiv_quad9n(localcoord,shapederiv(1:9,1:2))
      case (fe_hex8n, fe_mitc4_shell361)
        ! error check
        call shapederiv_hex8n(localcoord,shapederiv(1:8,1:3))
      case (fe_hex20n)
        ! error check
        call shapederiv_hex20n(localcoord, shapederiv(1:20,1:3))
      case (fe_prism6n, fe_mitc3_shell361)
        call shapederiv_prism6n(localcoord,shapederiv(1:6,1:3))
      case (fe_prism15n)
        call shapederiv_prism15n(localcoord,shapederiv(1:15,1:3))
      case (fe_tet4n, fe_tet4n_pipi, fe_beam341)
        ! error check
        call shapederiv_tet4n(shapederiv(1:4,1:3))
      case (fe_tet10n)
        ! error check
        call shapederiv_tet10n(localcoord,shapederiv(1:10,1:3))
      case default
        ! error message
        stop "Element type not defined-sde"
    end select
  end subroutine
 
  subroutine getshape2ndderiv( fetype, localcoord, shapederiv )
    integer, intent(in)           :: fetype
    real(kind=kreal), intent(in)  :: localcoord(:)      
    real(kind=kreal), intent(out) :: shapederiv(:,:,:)  
 
    select case (fetype)
      case ( fe_tri3n, fe_mitc3_shell )
        !error check
        call shape2ndderiv_tri3n(shapederiv(1:3,1:2,1:2))
      case (fe_tri6n)
        !error check
        call shape2ndderiv_tri6n(shapederiv(1:6,1:2,1:2))
      case ( fe_quad4n, fe_mitc4_shell )
        !error check
        call shape2ndderiv_quad4n(shapederiv(1:4,1:2,1:2))
      case (fe_quad8n)
        !error check
        call shape2ndderiv_quad8n(localcoord,shapederiv(1:8,1:2,1:2))
      case default
        ! error message
        stop "Cannot calculate second derivatives of shape function"
    end select
  end subroutine
 
  subroutine getshapefunc( fetype, localcoord, func )
    integer, intent(in)           :: fetype
    real(kind=kreal), intent(in)  :: localcoord(:)     
    real(kind=kreal), intent(out) :: func(:)           
 
    select case (fetype)
      case ( fe_tri3n, fe_mitc3_shell )
        !error check
        call shapefunc_tri3n(localcoord,func(1:3))
      case (fe_tri6n)
        !error check
        call shapefunc_tri6n(localcoord,func(1:6))
      case ( fe_quad4n, fe_mitc4_shell )
        !error check
        call shapefunc_quad4n(localcoord,func(1:4))
      case (fe_quad8n)
        !error check
        call shapefunc_quad8n(localcoord,func(1:8))
      case (fe_hex8n, fe_mitc4_shell361)
        ! error check
        call shapefunc_hex8n(localcoord,func(1:8))
      case ( fe_mitc9_shell )
        !error check
        call shapefunc_quad9n(localcoord,func(1:9))
      case (fe_hex20n)
        ! error check
        call shapefunc_hex20n(localcoord,func(1:20))
      case (fe_prism6n, fe_mitc3_shell361)
        call shapefunc_prism6n(localcoord,func(1:6))
      case (fe_prism15n)
        call shapefunc_prism15n(localcoord,func(1:15))
      case (fe_tet4n, fe_tet4n_pipi, fe_beam341)
        ! error check
        call shapefunc_tet4n(localcoord,func(1:4))
      case (fe_tet10n)
        ! error check
        call shapefunc_tet10n(localcoord,func(1:10))
      case (fe_line2n, fe_if_line2n)
        !error check
        call shapefunc_line2n(localcoord,func(1:2))
      case (fe_line3n)
        !error check
        call shapefunc_line3n(localcoord,func(1:3))
      case default
        stop "Element type not defined-sf"
        ! error message
    end select
  end subroutine
 
 
  ! (Gaku Hashimoto, The University of Tokyo, 2012/11/15) <
  !####################################################################
  subroutine getnodalnaturalcoord(fetype, nncoord)
    !####################################################################
 
    integer, intent(in)             :: fetype
    real(kind = kreal), intent(out) :: nncoord(:, :)
 
    !--------------------------------------------------------------------
 
    select case( fetype )
      case( fe_tri3n, fe_mitc3_shell, fe_mitc3_shell361 )
 
        !error check
        call nodalnaturalcoord_tri3n( nncoord(1:3, 1:2) )
 
      case( fe_quad4n, fe_mitc4_shell, fe_mitc4_shell361 )
 
        !error check
        call nodalnaturalcoord_quad4n( nncoord(1:4, 1:2) )
 
      case( fe_mitc9_shell )
 
        !error check
        call nodalnaturalcoord_quad9n( nncoord(1:9, 1:2) )
 
      case default
 
        ! error message
        stop "Element type not defined-sde"
 
    end select
 
    !--------------------------------------------------------------------
 
    return
 
    !####################################################################
  end subroutine getnodalnaturalcoord
  !####################################################################
  ! > (Gaku Hashimoto, The University of Tokyo, 2012/11/15)
 
 
  subroutine getglobalderiv( fetype, nn, localcoord, elecoord, det, gderiv )
    integer, intent(in)           :: fetype
    integer, intent(in)           :: nn
    real(kind=kreal), intent(in)  :: localcoord(:)   
    real(kind=kreal), intent(in)  :: elecoord(:,:)   
    real(kind=kreal), intent(out) :: det             
    real(kind=kreal), intent(out) :: gderiv(:,:)     
 
    real(kind=kreal) :: dum, xj(3,3), xji(3,3), deriv(nn,3)
    integer          :: nspace
 
    nspace = getspacedimension( fetype )
    call getshapederiv( fetype, localcoord(:), deriv(1:nn,:) )
 
    if( nspace==2 ) then
      xj(1:2,1:2)=matmul( elecoord(1:2,1:nn), deriv(1:nn,1:2) )
      det=xj(1,1)*xj(2,2)-xj(2,1)*xj(1,2)
      if( det==0.d0 ) stop "Math error in GetGlobalDeriv! Determinant==0.0"
      dum=1.d0/det
      xji(1,1)= xj(2,2)*dum
      xji(1,2)=-xj(1,2)*dum
      xji(2,1)=-xj(2,1)*dum
      xji(2,2)= xj(1,1)*dum
    else
      !  JACOBI MATRIX
      xj(1:3,1:3)= matmul( elecoord(1:3,1:nn), deriv(1:nn,1:3) )
      !DETERMINANT OF JACOBIAN
      det=xj(1,1)*xj(2,2)*xj(3,3)                                             &
        +xj(2,1)*xj(3,2)*xj(1,3)                                             &
        +xj(3,1)*xj(1,2)*xj(2,3)                                             &
        -xj(3,1)*xj(2,2)*xj(1,3)                                             &
        -xj(2,1)*xj(1,2)*xj(3,3)                                             &
        -xj(1,1)*xj(3,2)*xj(2,3)
      if( det==0.d0 ) stop "Math error in GetGlobalDeriv! Determinant==0.0"
      ! INVERSION OF JACOBIAN
      dum=1.d0/det
      xji(1,1)=dum*( xj(2,2)*xj(3,3)-xj(3,2)*xj(2,3) )
      xji(1,2)=dum*(-xj(1,2)*xj(3,3)+xj(3,2)*xj(1,3) )
      xji(1,3)=dum*( xj(1,2)*xj(2,3)-xj(2,2)*xj(1,3) )
      xji(2,1)=dum*(-xj(2,1)*xj(3,3)+xj(3,1)*xj(2,3) )
      xji(2,2)=dum*( xj(1,1)*xj(3,3)-xj(3,1)*xj(1,3) )
      xji(2,3)=dum*(-xj(1,1)*xj(2,3)+xj(2,1)*xj(1,3) )
      xji(3,1)=dum*( xj(2,1)*xj(3,2)-xj(3,1)*xj(2,2) )
      xji(3,2)=dum*(-xj(1,1)*xj(3,2)+xj(3,1)*xj(1,2) )
      xji(3,3)=dum*( xj(1,1)*xj(2,2)-xj(2,1)*xj(1,2) )
    endif
 
    gderiv(1:nn,1:nspace)=matmul( deriv(1:nn,1:nspace), xji(1:nspace,1:nspace) )
  end subroutine
 
  real(kind=kreal) function getdeterminant( fetype, nn, localcoord, elecoord )
    integer, intent(in)           :: fetype
    integer, intent(in)           :: nn
    real(kind=kreal), intent(in)  :: localcoord(:)   
    real(kind=kreal), intent(in)  :: elecoord(:,:)   
 
    real(kind=kreal) :: xj(3,3), deriv(nn,3)
    integer          :: nspace
 
    nspace = getspacedimension( fetype )
    call getshapederiv( fetype, localcoord(:), deriv(1:nn,:) )
 
    if( nspace==2 ) then
      xj(1:2,1:2)=matmul( elecoord(1:2,1:nn), deriv(1:nn,1:2) )
      getdeterminant=xj(1,1)*xj(2,2)-xj(2,1)*xj(1,2)
    else
      xj(1:3,1:3)= matmul( elecoord(1:3,1:nn), deriv(1:nn,1:3) )
      getdeterminant=xj(1,1)*xj(2,2)*xj(3,3)                                             &
        +xj(2,1)*xj(3,2)*xj(1,3)                                             &
        +xj(3,1)*xj(1,2)*xj(2,3)                                             &
        -xj(3,1)*xj(2,2)*xj(1,3)                                             &
        -xj(2,1)*xj(1,2)*xj(3,3)                                             &
        -xj(1,1)*xj(3,2)*xj(2,3)
    endif
 
  end function
 
  subroutine getjacobian( fetype, nn, localcoord, elecoord, det, jacobian, inverse )
    integer, intent(in)           :: fetype
    integer, intent(in)           :: nn
    real(kind=kreal), intent(in)  :: localcoord(:)   
    real(kind=kreal), intent(in)  :: elecoord(:,:)   
    real(kind=kreal), intent(out) :: det             
    real(kind=kreal), intent(out) :: jacobian(:,:)   
    real(kind=kreal), intent(out) :: inverse(:,:)    
 
    real(kind=kreal) :: dum, deriv(nn,3)
    integer          :: nspace
 
    nspace = getspacedimension( fetype )
    call getshapederiv( fetype, localcoord(:), deriv(1:nn,:) )
 
    if( nspace==2 ) then
      jacobian(1:2,1:2)=matmul( elecoord(1:2,1:nn), deriv(1:nn,1:2) )
      det=jacobian(1,1)*jacobian(2,2)-jacobian(2,1)*jacobian(1,2)
      if( det==0.d0 ) stop "Math error in getJacobain! Determinant==0.0"
      dum=1.0/det
      inverse(1,1)= jacobian(2,2)*dum
      inverse(1,2)=-jacobian(1,2)*dum
      inverse(2,1)=-jacobian(2,1)*dum
      inverse(2,2)= jacobian(1,1)*dum
    else
      !  JACOBI MATRIX
      jacobian(1:3,1:3)= matmul( elecoord(1:3,1:nn), deriv(1:nn,1:3) )
      !DETERMINANT OF JACOBIAN
      det=jacobian(1,1)*jacobian(2,2)*jacobian(3,3)                                         &
        +jacobian(2,1)*jacobian(3,2)*jacobian(1,3)                                         &
        +jacobian(3,1)*jacobian(1,2)*jacobian(2,3)                                         &
        -jacobian(3,1)*jacobian(2,2)*jacobian(1,3)                                         &
        -jacobian(2,1)*jacobian(1,2)*jacobian(3,3)                                         &
        -jacobian(1,1)*jacobian(3,2)*jacobian(2,3)
      if( det==0.d0 ) stop "Math error in getJacobain! Determinant==0.0"
      ! INVERSION OF JACOBIAN
      dum=1.d0/det
      inverse(1,1)=dum*( jacobian(2,2)*jacobian(3,3)-jacobian(3,2)*jacobian(2,3) )
      inverse(1,2)=dum*(-jacobian(1,2)*jacobian(3,3)+jacobian(3,2)*jacobian(1,3) )
      inverse(1,3)=dum*( jacobian(1,2)*jacobian(2,3)-jacobian(2,2)*jacobian(1,3) )
      inverse(2,1)=dum*(-jacobian(2,1)*jacobian(3,3)+jacobian(3,1)*jacobian(2,3) )
      inverse(2,2)=dum*( jacobian(1,1)*jacobian(3,3)-jacobian(3,1)*jacobian(1,3) )
      inverse(2,3)=dum*(-jacobian(1,1)*jacobian(2,3)+jacobian(2,1)*jacobian(1,3) )
      inverse(3,1)=dum*( jacobian(2,1)*jacobian(3,2)-jacobian(3,1)*jacobian(2,2) )
      inverse(3,2)=dum*(-jacobian(1,1)*jacobian(3,2)+jacobian(3,1)*jacobian(1,2) )
      inverse(3,3)=dum*( jacobian(1,1)*jacobian(2,2)-jacobian(2,1)*jacobian(1,2) )
    endif
  end subroutine
 
  function surfacenormal( fetype, nn, localcoord, elecoord ) result( normal )
    integer, intent(in)           :: fetype
    integer, intent(in)           :: nn
    real(kind=kreal), intent(in)  :: localcoord(2)     
    real(kind=kreal), intent(in)  :: elecoord(3,nn)    
    real(kind=kreal) :: normal(3)
    real(kind=kreal) :: deriv(nn,2), gderiv(3,2)
 
    select case (fetype)
      case (fe_tri3n)
        !error check
        call shapederiv_tri3n(deriv(1:3,1:2))
      case (fe_tri6n)
        !error check
        call shapederiv_tri6n(localcoord,deriv(1:6,1:2))
      case (fe_quad4n)
        !error check
        call shapederiv_quad4n(localcoord,deriv(1:4,1:2))
      case (fe_quad8n)
        !error check
        call shapederiv_quad8n(localcoord,deriv(1:8,1:2))
      case default
        ! error message
        normal =0.d0
        return
    end select
 
    gderiv = matmul( elecoord, deriv )
    normal(1) = gderiv(2,1)*gderiv(3,2) - gderiv(3,1)*gderiv(2,2)
    normal(2) = gderiv(3,1)*gderiv(1,2) - gderiv(1,1)*gderiv(3,2)
    normal(3) = gderiv(1,1)*gderiv(2,2) - gderiv(2,1)*gderiv(1,2)
    !  normal = normal/dsqrt(dot_product(normal, normal))
  end function
 
  function edgenormal( fetype, nn, localcoord, elecoord ) result( normal )
    integer, intent(in)           :: fetype
    integer, intent(in)           :: nn
    real(kind=kreal), intent(in)  :: localcoord(1)     
    real(kind=kreal), intent(in)  :: elecoord(2,nn)    
    real(kind=kreal) :: normal(2)
    real(kind=kreal) :: deriv(nn,1), gderiv(2,1)
 
    select case (fetype)
      case (fe_line2n, fe_if_line2n)
        !error check
        call shapederiv_line2n(deriv(1:nn,:))
      case (fe_line3n)
        !error check
        call shapederiv_line3n(localcoord,deriv(1:nn,:))
      case default
        ! error message
        normal =0.d0
        return
    end select
 
    gderiv = matmul( elecoord, deriv )
    normal(1) = -gderiv(2,1)
    normal(2) = gderiv(1,1)
    !  normal = normal/dsqrt(dot_product(normal, normal))
  end function
 
  subroutine tangentbase( fetype, nn, localcoord, elecoord, tangent )
    integer, intent(in)           :: fetype
    integer, intent(in)           :: nn
    real(kind=kreal), intent(in)  :: localcoord(2)     
    real(kind=kreal), intent(in)  :: elecoord(3,nn)    
    real(kind=kreal), intent(out) :: tangent(3,2)      
    real(kind=kreal) :: deriv(nn,2)
 
    select case (fetype)
      case (fe_tri3n)
        !error check
        call shapederiv_tri3n(deriv(1:3,1:2))
      case (fe_tri6n)
        !error check
        call shapederiv_tri6n(localcoord,deriv(1:6,1:2))
      case (fe_tri6nc)
        !error check
        call shapederiv_tri6n(localcoord,deriv(1:6,1:2))
      case (fe_quad4n)
        !error check
        call shapederiv_quad4n(localcoord,deriv(1:4,1:2))
      case (fe_quad8n)
        !error check
        call shapederiv_quad8n(localcoord,deriv(1:8,1:2))
      case default
        ! error message
        tangent =0.d0
        return
    end select
 
    tangent = matmul( elecoord, deriv )
  end subroutine tangentbase
 
  subroutine curvature( fetype, nn, localcoord, elecoord, l2ndderiv, normal, curv )
    integer, intent(in)           :: fetype
    integer, intent(in)           :: nn
    real(kind=kreal), intent(in)  :: localcoord(2)     
    real(kind=kreal), intent(in)  :: elecoord(3,nn)    
    real(kind=kreal), intent(out) :: l2ndderiv(3,2,2)  
    real(kind=kreal), intent(in), optional  :: normal(3)     
    real(kind=kreal), intent(out), optional :: curv(2,2)     
    real(kind=kreal) :: deriv2(nn,2,2)
 
    select case (fetype)
      case (fe_tri3n)
        !error check
        call shape2ndderiv_tri3n(deriv2(1:3,1:2,1:2))
      case (fe_tri6n)
        !error check
        call shape2ndderiv_tri6n(deriv2(1:6,1:2,1:2))
      case (fe_tri6nc)
        !error check
        call shape2ndderiv_tri6n(deriv2(1:6,1:2,1:2))
        ! deriv2=0.d0
      case (fe_quad4n)
        !error check
        call shape2ndderiv_quad4n(deriv2(1:4,1:2,1:2))
      case (fe_quad8n)
        !error check
        call shape2ndderiv_quad8n(localcoord,deriv2(1:8,1:2,1:2))
      case default
        ! error message
        stop "Cannot calculate second derivatives of shape function"
    end select
 
    l2ndderiv(1:3,1,1) = matmul( elecoord(1:3,1:nn), deriv2(1:nn,1,1) )
    l2ndderiv(1:3,1,2) = matmul( elecoord(1:3,1:nn), deriv2(1:nn,1,2) )
    l2ndderiv(1:3,2,1) = matmul( elecoord(1:3,1:nn), deriv2(1:nn,2,1) )
    l2ndderiv(1:3,2,2) = matmul( elecoord(1:3,1:nn), deriv2(1:nn,2,2) )
    if( present(curv) ) then
      curv(1,1) = dot_product( l2ndderiv(:,1,1), normal(:) )
      curv(1,2) = dot_product( l2ndderiv(:,1,2), normal(:) )
      curv(2,1) = dot_product( l2ndderiv(:,2,1), normal(:) )
      curv(2,2) = dot_product( l2ndderiv(:,2,2), normal(:) )
    endif
  end subroutine curvature
 
  subroutine getelementcenter( fetype, localcoord )
    integer, intent(in)           :: fetype
    real(kind=kreal), intent(out) :: localcoord(:)     
 
    select case (fetype)
      case (fe_tri3n, fe_tri6n, fe_tri6nc)
        localcoord(1:2) = 1.d0/3.d0
      case (fe_quad4n, fe_quad8n)
        localcoord(1:2) = 0.d0
      case (fe_tet4n, fe_tet10n, fe_tet10nc)
        localcoord(1:3) = 1.d0/4.d0
      case (fe_prism6n, fe_prism15n)
        localcoord(1) = 1.d0/3.d0
        localcoord(2) = 1.d0/3.d0
        localcoord(3) = 0.d0
      case (fe_hex8n, fe_hex20n, fe_hex27n)
        localcoord(1:3) = 0.d0
      case default
        localcoord(:) = 0.d0
    end select
  end subroutine getelementcenter
 
  integer function isinsideelement( fetype, localcoord, clearance )
    integer, intent(in)              :: fetype
    real(kind=kreal), intent(inout)  :: localcoord(2)   
    real(kind=kreal), optional       :: clearance       
    real(kind=kreal) :: clr, coord3
 
    clr = 1.d-6
    if( present(clearance) ) clr = clearance
    if( dabs(localcoord(1))<clr ) localcoord(1)=0.d0
    if( dabs(localcoord(2))<clr ) localcoord(2)=0.d0
    if( dabs(dabs(localcoord(1))-1.d0)<clr )    &
      localcoord(1)=sign(1.d0,localcoord(1))
    if( dabs(dabs(localcoord(2))-1.d0)<clr )    &
      localcoord(2)=sign(1.d0,localcoord(2))
    isinsideelement = -1
    select case (fetype)
      case (fe_tri3n, fe_tri6n, fe_tri6nc)
        !error check
        coord3 = 1.d0-(localcoord(1)+localcoord(2))
        if( dabs(coord3)<clr ) coord3=0.d0
        if( localcoord(1)>=0.d0 .and. localcoord(1)<=1.d0 .and.   &
            localcoord(2)>=0.d0 .and. localcoord(2)<=1.d0 .and.   &
            coord3>=0.d0 .and. coord3<=1.d0 )  then
          isinsideelement = 0
          if( localcoord(1)==1.d0 ) then
            isinsideelement = 1
          elseif( localcoord(2)==1.d0 ) then
            isinsideelement = 2
          elseif( coord3==1.d0 ) then
            isinsideelement = 3
          elseif( coord3==0.d0 ) then
            isinsideelement = 12
          elseif( localcoord(1)==0.d0 ) then
            isinsideelement = 23
          elseif( localcoord(2)==0.d0 ) then
            isinsideelement = 31
          endif
        endif
      case (fe_quad4n, fe_quad8n)
        !error check
        if( all(dabs(localcoord)<=1.d0) ) then
          isinsideelement = 0
          if( localcoord(1)==-1.d0 .and. localcoord(2)==-1.d0 ) then
            isinsideelement = 1
          elseif( localcoord(1)==1.d0 .and. localcoord(2)==-1.d0 ) then
            isinsideelement = 2
          elseif( localcoord(1)==1.d0 .and. localcoord(2)==1.d0 ) then
            isinsideelement = 3
          elseif( localcoord(1)==-1.d0 .and. localcoord(2)==1.d0 ) then
            isinsideelement = 4
          elseif( localcoord(2)==-1.d0 ) then
            isinsideelement = 12
          elseif( localcoord(1)==1.d0 ) then
            isinsideelement = 23
          elseif( localcoord(2)==1.d0 ) then
            isinsideelement = 34
          elseif( localcoord(1)==-1.d0 ) then
            isinsideelement = 41
          endif
        endif
    end select
  end function isinsideelement
 
  integer function isinside3delement( fetype, localcoord, clearance )
    integer, intent(in)              :: fetype
    real(kind=kreal), intent(inout)  :: localcoord(3)   
    real(kind=kreal), optional       :: clearance       
    real(kind=kreal) :: clr, coord4
 
    integer :: idof
 
    clr = 1.d-6
    if( present(clearance) ) clr = clearance
    do idof=1,3
      if( dabs(localcoord(idof))<clr ) localcoord(idof)=0.d0
      if( dabs(dabs(localcoord(idof))-1.d0)<clr )    &
        &  localcoord(idof)=sign(1.d0,localcoord(idof))
    enddo
 
    isinside3delement = -1
    select case (fetype)
      case (fe_tet4n, fe_tet4n_pipi, fe_tet10n, fe_tet10nc)
        !error check
        coord4 = 1.d0-(localcoord(1)+localcoord(2)+localcoord(3))
        if( dabs(coord4)<clr ) coord4=0.d0
        isinside3delement = 0
        do idof=1,3
          if( localcoord(idof) < 0.d0 .or. localcoord(idof) > 1.d0 ) isinside3delement = -1
        enddo
        if( coord4 < 0.d0 .or. coord4 > 1.d0 ) isinside3delement = -1
      case (fe_prism6n, fe_prism15n)
        !error check
        coord4 = 1.d0-(localcoord(1)+localcoord(2))
        isinside3delement = 0
        do idof=1,2
          if( localcoord(idof) < 0.d0 .or. localcoord(idof) > 1.d0 ) isinside3delement = -1
        enddo
        if( localcoord(3) < -1.d0 .or. localcoord(3) > 1.d0 ) isinside3delement = -1
        if( coord4 < 0.d0 .or. coord4 > 1.d0 ) isinside3delement = -1
      case (fe_hex8n, fe_hex20n, fe_hex27n)
        if( all(dabs(localcoord)<=1.d0) ) isinside3delement = 0
      end select
  end function
 
  subroutine getvertexcoord( fetype, cnode, localcoord )
    integer, intent(in)            :: fetype
    integer, intent(in)            :: cnode
    real(kind=kreal), intent(out)  :: localcoord(2)   
 
    select case (fetype)
      case (fe_tri3n, fe_tri6n, fe_tri6nc)
        if( cnode==1 ) then
          localcoord(1) =1.d0
          localcoord(2) =0.d0
        elseif( cnode==2 ) then
          localcoord(1) =0.d0
          localcoord(2) =1.d0
        else
          localcoord(1) =0.d0
          localcoord(2) =0.d0
        endif
      case (fe_quad4n, fe_quad8n)
        if( cnode==1 ) then
          localcoord(1) =-1.d0
          localcoord(2) =-1.d0
        elseif( cnode==2 ) then
          localcoord(1) =1.d0
          localcoord(2) =-1.d0
        elseif( cnode==3 ) then
          localcoord(1) =1.d0
          localcoord(2) =1.d0
        else
          localcoord(1) =-1.d0
          localcoord(2) =1.d0
        endif
    end select
  end subroutine
 
  subroutine extrapolatevalue( lpos, fetype, nnode, pvalue, ndvalue )
    real(kind=kreal), intent(in)  :: lpos(:)        
    integer, intent(in)           :: fetype
    integer, intent(in)           :: nnode
    real(kind=kreal), intent(in)  :: pvalue(:)      
    real(kind=kreal), intent(out) :: ndvalue(:,:)   
 
    integer         :: i
    real(kind=kreal) :: shapefunc(nnode)
    call getshapefunc( fetype, lpos, shapefunc )
    do i=1,nnode
      ndvalue(i,:) = shapefunc(i)*pvalue(:)
    enddo
  end subroutine
 
  subroutine interapolatevalue( lpos, fetype, nnode, pvalue, ndvalue )
    real(kind=kreal), intent(in)  :: lpos(:)        
    integer, intent(in)           :: fetype
    integer, intent(in)           :: nnode
    real(kind=kreal), intent(out) :: pvalue(:)      
    real(kind=kreal), intent(in)  :: ndvalue(:,:)   
 
    integer         :: i
    real(kind=kreal) :: shapefunc(nnode)
    call getshapefunc( fetype, lpos, shapefunc )
    pvalue(:) = 0
    do i=1,nnode
      pvalue(:) = pvalue(:)+ shapefunc(i)*ndvalue(i,:)
    enddo
  end subroutine
 
  subroutine gauss2node( fetype, gaussv, nodev )
    integer, intent(in)           :: fetype
    real(kind=kreal), intent(in)  :: gaussv(:,:)  
    real(kind=kreal), intent(out) :: nodev(:,:)   
 
    integer :: i, ngauss, nnode
    real(kind=kreal) :: localcoord(3), func(100)
    ngauss = numofquadpoints( fetype )
    nnode = getnumberofnodes( fetype )
    ! error checking
    select case (fetype)
      case (fe_tri3n)
        !error check
        do i=1,nnode
          nodev(i,:) = gaussv(1,:)
        enddo
      case (fe_tri6n)
        !error check
        !   func(1:6) = ShapeFunc_tri6n(localcoord)
      case (fe_quad4n)
        !error check
        !  nodev(:,:) = gaussv(1,:)
      case (fe_quad8n)
        !error check
        call shapefunc_quad8n(localcoord,func(1:8))
      case (fe_hex8n, fe_mitc4_shell361)
        ! error check
        call shapefunc_hex8n(localcoord,func(1:8))
      case (fe_hex20n)
        ! error check
        call shapefunc_hex20n(localcoord,func(1:20))
      case (fe_prism6n, fe_mitc3_shell361)
        do i=1,3
          nodev(i,:) = gaussv(1,:)
        enddo
        do i=1,3
          nodev(i+3,:) = gaussv(2,:)
        enddo
      case (fe_prism15n)
        call shapefunc_prism15n(localcoord,func(1:15))
      case (fe_tet4n, fe_tet4n_pipi, fe_beam341)
        ! error check
        do i=1,nnode
          nodev(i,:) = gaussv(1,:)
        enddo
      case (fe_tet10n)
        ! error check
        call shapefunc_tet10n(localcoord,func(1:10))
      case default
        stop "Element type not defined"
        ! error message
    end select
  end subroutine
 
  real(kind=kreal) function getreferencelength( fetype, nn, localcoord, elecoord )
    integer, intent(in)         :: fetype
    integer, intent(in)         :: nn
    real(kind=kreal),intent(in) :: localcoord(2)    
    real(kind=kreal),intent(in) :: elecoord(3,nn)   
    real(kind=kreal) :: detjxy, detjyz, detjxz, detj
    detjxy = getdeterminant( fetype, nn, localcoord, elecoord(1:2,1:nn) )
    detjyz = getdeterminant( fetype, nn, localcoord, elecoord(2:3,1:nn) )
    detjxz = getdeterminant( fetype, nn, localcoord, elecoord(1:3:2,1:nn) )
    detj = dsqrt( detjxy **2 + detjyz **2 + detjxz **2 )
    getreferencelength = dsqrt( detj )
  end function getreferencelength
 
 
end module