module ShapeFunctionClass
    use, intrinsic :: iso_fortran_env
    use MathClass
    use ArrayClass
    use IOClass
    implicit none
    
    type::ShapeFunction_
        
        real(real64),allocatable::Nmat(:)
        real(real64),allocatable::dNdgzi(:,:)
        real(real64),allocatable::dNdgzidgzi(:,:)
        real(real64),allocatable::gzi(:)
        real(real64),allocatable::GaussPoint(:,:)
        real(real64),allocatable::GaussIntegWei(:)
        real(real64),allocatable::Jmat(:,:),JmatInv(:,:)
        real(real64),allocatable::ElemCoord(:,:)
        real(real64),allocatable::ElemCoord_n(:,:)
        real(real64),allocatable::du(:,:)
        real(real64) :: detJ
        integer(int32) :: NumOfNode
        integer(int32) :: NumOfOrder
        integer(int32) :: NumOfDim=0
        integer(int32) :: NumOfGp=0
        integer(int32) :: GpID
        integer(int32) :: ierr
        integer(int32) :: currentGpID
        integer(int32) :: ElementID
        logical :: ReducedIntegration = .false.
        logical :: Empty=.true.
        character*70::ElemType
        character(len=60):: ErrorMsg
    contains
        procedure :: init => initShapeFunction
        procedure :: update => updateShapeFunction
        procedure :: SetType => SetShapeFuncType
        procedure :: GetAll  => GetAllShapeFunc
        procedure :: get  => GetAllShapeFunc
        procedure :: getOnlyNvec  => GetShapeFunction
        procedure :: Deallocate => DeallocateShapeFunction
        procedure :: getType => getShapeFuncType 
        procedure :: GetGaussPoint => GetGaussPoint
        procedure :: export => exportShapeFunction
        procedure :: remove => removeShapeFunction
        procedure :: save => saveShapeFunction
        procedure :: open => openShapeFunction
        procedure :: getNvec => getNvecShapeFunction
    end type ShapeFunction_

contains

! #####################################################
function getNvecShapeFunction(obj,x,y,z) result(nvec)
    class(Shapefunction_),intent(in) :: obj
    real(real64),optional,intent(in)::x,y,z
    real(real64),allocatable :: nvec(:)

    
end function
! #####################################################

! #####################################################
subroutine openShapeFunction(obj,path,name)
    class(ShapeFunction_ ),intent(inout)::obj
    character(*),intent(in) :: path
    character(*),optional,intent(in) :: name
    type(IO_) :: f

    if(present(name) )then
        call system("mkdir -p "//trim(path)//"/"//trim(adjustl(name)))
        
        call f%open(trim(path)//"/"//trim(adjustl(name))//"/","ShapeFunction","prop" )
        
        call openArray(f%fh,obj%Nmat)
        call openArray(f%fh,obj%dNdgzi)
        call openArray(f%fh,obj%dNdgzidgzi)
        call openArray(f%fh,obj%gzi)
        call openArray(f%fh,obj%GaussPoint)
        call openArray(f%fh,obj%GaussIntegWei)
        call openArray(f%fh,obj%Jmat)
        call openArray(f%fh,obj%JmatInv)
        call openArray(f%fh,obj%ElemCoord)
        call openArray(f%fh,obj%ElemCoord_n)
        call openArray(f%fh,obj%du)
        read(f%fh,*)obj%detJ
        read(f%fh,*)obj%NumOfNode
        read(f%fh,*)obj%NumOfOrder
        read(f%fh,*)obj%NumOfDim
        read(f%fh,*)obj%NumOfGp
        read(f%fh,*)obj%GpID
        read(f%fh,*)obj%ierr
        read(f%fh,*)obj%currentGpID
        read(f%fh,*)obj%ReducedIntegration
        read(f%fh,*)obj%ElemType
        read(f%fh,*)obj%ErrorMsg
        
        call f%close()
    else

        call system("mkdir -p "//trim(path)//"/ShapeFunction")
        call f%open(trim(path)//"/","ShapeFunction/ShapeFunction","prop" )
        
        call openArray(f%fh,obj%Nmat)
        call openArray(f%fh,obj%dNdgzi)
        call openArray(f%fh,obj%dNdgzidgzi)
        call openArray(f%fh,obj%gzi)
        call openArray(f%fh,obj%GaussPoint)
        call openArray(f%fh,obj%GaussIntegWei)
        call openArray(f%fh,obj%Jmat)
        call openArray(f%fh,obj%JmatInv)
        call openArray(f%fh,obj%ElemCoord)
        call openArray(f%fh,obj%ElemCoord_n)
        call openArray(f%fh,obj%du)
        read(f%fh,*)obj%detJ
        read(f%fh,*)obj%NumOfNode
        read(f%fh,*)obj%NumOfOrder
        read(f%fh,*)obj%NumOfDim
        read(f%fh,*)obj%NumOfGp
        read(f%fh,*)obj%GpID
        read(f%fh,*)obj%ierr
        read(f%fh,*)obj%currentGpID
        read(f%fh,*)obj%ReducedIntegration
        read(f%fh,*)obj%ElemType
        read(f%fh,*)obj%ErrorMsg
        
        call f%close()
    endif

end subroutine
! #####################################################


! #####################################################
subroutine saveShapeFunction(obj,path,name)
    class(ShapeFunction_ ),intent(inout)::obj
    character(*),intent(in) :: path
    character(*),optional,intent(in) :: name
    type(IO_) :: f

    if(present(name) )then
        call system("mkdir -p "//trim(path)//"/"//trim(adjustl(name)))
        
        call f%open(trim(path)//"/"//trim(adjustl(name))//"/","ShapeFunction","prop" )
        
        call writeArray(f%fh,obj%Nmat)
        call writeArray(f%fh,obj%dNdgzi)
        call writeArray(f%fh,obj%dNdgzidgzi)
        call writeArray(f%fh,obj%gzi)
        call writeArray(f%fh,obj%GaussPoint)
        call writeArray(f%fh,obj%GaussIntegWei)
        call writeArray(f%fh,obj%Jmat)
        call writeArray(f%fh,obj%JmatInv)
        call writeArray(f%fh,obj%ElemCoord)
        call writeArray(f%fh,obj%ElemCoord_n)
        call writeArray(f%fh,obj%du)
        write(f%fh,*)obj%detJ
        write(f%fh,*)obj%NumOfNode
        write(f%fh,*)obj%NumOfOrder
        write(f%fh,*)obj%NumOfDim
        write(f%fh,*)obj%NumOfGp
        write(f%fh,*)obj%GpID
        write(f%fh,*)obj%ierr
        write(f%fh,*)obj%currentGpID
        write(f%fh,*)obj%ReducedIntegration
        write(f%fh,*)trim(obj%ElemType)
        write(f%fh,*)trim(obj%ErrorMsg)
        
        call f%close()
    else

        call system("mkdir -p "//trim(path)//"/ShapeFunction")
        call f%open(trim(path)//"/","ShapeFunction/ShapeFunction","prop" )
        
        call writeArray(f%fh,obj%Nmat)
        call writeArray(f%fh,obj%dNdgzi)
        call writeArray(f%fh,obj%dNdgzidgzi)
        call writeArray(f%fh,obj%gzi)
        call writeArray(f%fh,obj%GaussPoint)
        call writeArray(f%fh,obj%GaussIntegWei)
        call writeArray(f%fh,obj%Jmat)
        call writeArray(f%fh,obj%JmatInv)
        call writeArray(f%fh,obj%ElemCoord)
        call writeArray(f%fh,obj%ElemCoord_n)
        call writeArray(f%fh,obj%du)
        write(f%fh,*) obj%detJ
        write(f%fh,*) obj%NumOfNode
        write(f%fh,*) obj%NumOfOrder
        write(f%fh,*) obj%NumOfDim
        write(f%fh,*) obj%NumOfGp
        write(f%fh,*) obj%GpID
        write(f%fh,*) obj%ierr
        write(f%fh,*) obj%currentGpID
        write(f%fh,*) obj%ReducedIntegration
        write(f%fh,*) trim(obj%ElemType)
        write(f%fh,*) trim(obj%ErrorMsg)
        
        call f%close()
    endif

end subroutine
! #####################################################

subroutine removeShapeFunction(obj)
    class(ShapeFunction_),intent(inout)::obj

    if(allocated(obj%Nmat))then
        deallocate(obj%Nmat)
    endif
    if(allocated(obj%dNdgzi))then
        deallocate(obj%dNdgzi)
    endif
    if(allocated(obj%dNdgzidgzi))then
        deallocate(obj%dNdgzidgzi)
    endif
    if(allocated(obj%gzi))then
        deallocate(obj%gzi)
    endif
    if(allocated(obj%GaussPoint))then
        deallocate(obj%GaussPoint)
    endif
    if(allocated(obj%GaussIntegWei))then
        deallocate(obj%GaussIntegWei)
    endif
    if(allocated(obj%Jmat))then
        deallocate(obj%Jmat)
    endif
    if(allocated(obj%ElemCoord))then
        deallocate(obj%ElemCoord)
    endif
    if(allocated(obj%ElemCoord_n))then
        deallocate(obj%ElemCoord_n)
    endif
    if(allocated(obj%du))then
        deallocate(obj%du)
    endif
    obj%detJ=0.0d0
    obj%NumOfNode=0
    obj%NumOfOrder=0
    obj%NumOfDim=0
    obj%NumOfGp=0
    obj%GpID=0
    obj%ierr=0
    obj%currentGpID=0
    obj%ReducedIntegration = .false.
    obj%ElemType=" "
    obj%ErrorMsg=" "

end subroutine

!##################################################
subroutine initShapeFunction(obj,ElemType)
    class(ShapeFunction_),intent(inout) :: obj
    character(*),intent(in),optional :: ElemType

    obj%ElemType = ElemType

    call obj%SetType()

end subroutine
!##################################################


!##################################################
subroutine updateShapeFunction(obj,ElemType,NodCoord,ElemNod,ElemID,GpID)
    class(ShapeFunction_),intent(inout) :: obj
    character(*),optional,intent(in) :: ElemType
    integer(int32),intent(in) ::ElemNod(:,:),ElemID,GpID
    real(real64),intent(in) ::NodCoord(:,:)
    integer(int32) :: i,j

    if(present(ElemType) )then
        call obj%init(ElemType)
    endif
    call obj%GetAll(elem_id=i,nod_coord=NodCoord,elem_nod=ElemNod,OptionalGpID=j)

end subroutine
!##################################################


!##################################################
subroutine SetShapeFuncType(obj,NumOfDim,NumOfNodePerElem,ReducedIntegration)
    class(ShapeFunction_),intent(inout)::obj
    logical,optional,intent(in) :: ReducedIntegration
    character*70 ::  TrimedElemType
    integer(int32),optional,intent(in) :: NumOfDim,NumOfNodePerElem

    if(present(NumOfDim) )then
        if(present(NumOfNodePerElem) )then
            call obj%getType(NumOfDim,NumOfNodePerElem)
        endif
    endif
    
    if(present(ReducedIntegration) )then
        obj%ReducedIntegration=ReducedIntegration
    endif
    
    TrimedElemType=trim(obj%ElemType)
    if(trim(TrimedElemType)=="LinearRectangularGp4")then

        obj%NumOfNode   = 4
        obj%NumOfOrder  = 1
        obj%NumOfDim    = 2
        obj%NumOfGp     = 4
        

    elseif(trim(TrimedElemType)=="LinearHexahedralGp8")then
        obj%NumOfNode   = 8
        obj%NumOfOrder  = 1
        obj%NumOfDim    = 3
        obj%NumOfGp     = 8
        
    elseif(trim(TrimedElemType)=="Triangular")then
        obj%NumOfNode   = 3
        obj%NumOfOrder  = 1
        obj%NumOfDim    = 2
        obj%NumOfGp     = 1
    elseif(trim(TrimedElemType)=="LinearTetrahedral")then
        obj%NumOfNode   = 4
        obj%NumOfOrder  = 1
        obj%NumOfDim    = 3
        obj%NumOfGp     = 1
    else
        obj%ErrorMsg="ShapeFunctionClass.f90 :: SetShapeFuncType >> Element : "//TrimedElemType//"is not defined."
        print *, "ShapeFunctionClass.f90 :: SetShapeFuncType >> Element : ",TrimedElemType,"is not defined."
        return
    endif

    

end subroutine SetShapeFuncType
!##################################################

!##################################################
subroutine getShapeFuncType(obj,NumOfDim,NumOfNodePerElem)
    class(ShapeFunction_),intent(inout)::obj
    character*70 ::  TrimedElemType
    integer(int32),intent(in) :: NumOfDim,NumOfNodePerElem

    if(NumOfDim==2)then
        if(NumOfNodePerElem==4)then
            obj%ElemType="LinearRectangularGp4"
        else
            print *, "For 2D, NumOfNodePerElem = ",NumOfNodePerElem," is not set."
        endif
    elseif(NumOfDim==3)then
        if(NumOfNodePerElem==8)then
            obj%ElemType="LinearHexahedralGp8"
        else
            print *, "For 3D, NumOfNodePerElem = ",NumOfNodePerElem," is not set."
        endif
    else
        print *, "NumOfDim = ",NumOfDim," is not set."
    endif
    

end subroutine getShapeFuncType
!##################################################

!##################################################
subroutine GetAllShapeFunc(obj,elem_id,nod_coord,nod_coord_n,elem_nod,OptionalNumOfNode,OptionalNumOfOrder,&
    OptionalNumOfDim,OptionalNumOfGp,OptionalGpID,ReducedIntegration,NumOfDim,NumOfNodePerElem)
    class(ShapeFunction_),intent(inout)::obj
    integer(int32),optional,intent(in)::elem_nod(:,:),elem_id
    integer(int32),optional,intent(in) :: NumOfDim,NumOfNodePerElem
    real(real64),optional,intent(in)::nod_coord(:,:),nod_coord_n(:,:)
    integer(int32),optional,intent(in)::OptionalNumOfNode,OptionalNumOfOrder,&
    OptionalNumOfDim,OptionalNumOfGp,OptionalGpID
    logical,optional,intent(in) :: ReducedIntegration
    integer(int32) :: nd, ne

    nd = input(default=0, option=NumOfDim )
    ne = input(default=0, option=NumOfNodePerElem )
    if(present(nod_coord) )then
        nd = size(nod_coord,2)
    endif
    if(present(nod_coord) )then
        ne = size(elem_nod,2)
    endif

    if(obj%NumOfGp==0 .and. obj%NumOfDim==0)then
        
        if(nd==0 .or. ne==0)then
            print *, "ERROR :: GetAllShapeFunc >> please import NumOfDim and NumOfNodePerElem"
            stop
        endif
        call obj%SetType(NumOfDim=nd,NumOfNodePerElem=ne,ReducedIntegration=ReducedIntegration)
    endif
    if(present(ReducedIntegration) )then
        obj%ReducedIntegration=ReducedIntegration
    endif

    if(present(OptionalNumOfNode ) )then
        obj%NumOfNode=OptionalNumOfNode
    endif

    if(present(OptionalNumOfOrder ) )then
        obj%NumOfOrder=OptionalNumOfOrder
    endif

    if(present(OptionalNumOfDim ) )then
        obj%NumOfDim=OptionalNumOfDim
    endif

    if(present(OptionalNumOfGp ) )then
        obj%NumOfGp=OptionalNumOfGp
    endif
    
    if(present(OptionalGpID ) )then
        obj%GpID=OptionalGpID
    endif


    call GetGaussPoint(obj)
    call SetGaussPoint(obj)



    call GetShapeFunction(obj)
    call GetShapeFuncDer1(obj)
    !call GetShapeFuncDer2(obj)



    if(.not.present(nod_coord) .or. .not.present(elem_nod) )then
        
        obj%ErrorMsg="Mesh%NodCoord and Mesh%ElemNod is necessary to get Jmat"
        print *, obj%ErrorMsg
        
        return
    endif
    call GetElemCoord(obj,nod_coord,elem_nod,elem_id)
    if(present(nod_coord_n) )then
        call GetElemCoord_n(obj,nod_coord_n,elem_nod,elem_id)
        call Getdu(obj)
    endif



    call GetJmat(obj)
    obj%detJ = det_mat(obj%Jmat,size(obj%Jmat,1 ) )

end subroutine GetAllShapeFunc
!##################################################


!##################################################
subroutine DeallocateShapeFunction(obj)
    class(ShapeFunction_),intent(inout)::obj

    if( allocated(obj%Nmat            ) ) deallocate(obj%Nmat             )
    if( allocated(obj%dNdgzi          ) ) deallocate(obj%dNdgzi           )
    if( allocated(obj%dNdgzidgzi      ) ) deallocate(obj%dNdgzidgzi       )
    if( allocated(obj%gzi             ) ) deallocate(obj%gzi              )
    if( allocated(obj%GaussPoint      ) ) deallocate(obj%GaussPoint       )
    if( allocated(obj%GaussIntegWei   ) ) deallocate(obj%GaussIntegWei    )
    if( allocated(obj%Jmat            ) ) deallocate(obj%Jmat           )
    
    obj%ErrorMsg="All allocatable entities are deallocated"
end subroutine DeallocateShapeFunction
!##################################################

!--------------------------------------------------------
    subroutine GetGaussPoint(obj)

        class(ShapeFunction_),intent(inout)::obj

        !include "./GetGaussPoint.f90"

        if(.not.allocated(obj%GaussPoint) )then
            allocate(obj%GaussPoint(obj%NumOfDim,obj%NumOfGp ) )
        else
            if(size(obj%GaussPoint,1)/=obj%NumOfDim .or.&
            size(obj%GaussPoint,2)/=obj%NumOfGp)then
                deallocate(obj%GaussPoint)
                allocate(obj%GaussPoint(obj%NumOfDim,obj%NumOfGp ) )
            endif
        endif

        if(.not.allocated(obj%GaussIntegWei) )then
            allocate(obj%GaussIntegWei(obj%NumOfGp ) )
        else
            if(size(obj%GaussIntegWei)/=obj%NumOfGp)then
                deallocate(obj%GaussIntegWei)
                allocate(obj%GaussIntegWei(obj%NumOfGp ) )
            endif
        endif

        if(size(obj%GaussPoint,1)==2 .and. size(obj%GaussPoint,2)==4)then
            obj%GaussPoint(1,1) = -0.57735026918962576d0
            obj%GaussPoint(1,2) =  0.57735026918962576d0
            obj%GaussPoint(1,3) =  0.57735026918962576d0
            obj%GaussPoint(1,4) = -0.57735026918962576d0

            obj%GaussPoint(2,1) = -0.57735026918962576d0
            obj%GaussPoint(2,2) = -0.57735026918962576d0
            obj%GaussPoint(2,3) =  0.57735026918962576d0
            obj%GaussPoint(2,4) =  0.57735026918962576d0
            obj%GaussIntegWei(:)=1.0d0
            if(obj%ReducedIntegration .eqv. .true.)then
                obj%GaussPoint(:,:) =0.0d0
                obj%GaussIntegWei(:)=0.250d0
            endif
        elseif(size(obj%GaussPoint,1)==3 .and. size(obj%GaussPoint,2)==8)then
            obj%GaussPoint(1,1) = -0.57735026918962576d0
            obj%GaussPoint(1,2) =  0.57735026918962576d0
            obj%GaussPoint(1,3) =  0.57735026918962576d0
            obj%GaussPoint(1,4) = -0.57735026918962576d0
            obj%GaussPoint(1,5) = -0.57735026918962576d0
            obj%GaussPoint(1,6) =  0.57735026918962576d0
            obj%GaussPoint(1,7) =  0.57735026918962576d0
            obj%GaussPoint(1,8) = -0.57735026918962576d0

            obj%GaussPoint(2,1) = -0.57735026918962576d0
            obj%GaussPoint(2,2) = -0.57735026918962576d0
            obj%GaussPoint(2,3) =  0.57735026918962576d0
            obj%GaussPoint(2,4) =  0.57735026918962576d0
            obj%GaussPoint(2,5) = -0.57735026918962576d0
            obj%GaussPoint(2,6) = -0.57735026918962576d0
            obj%GaussPoint(2,7) =  0.57735026918962576d0
            obj%GaussPoint(2,8) =  0.57735026918962576d0

            obj%GaussPoint(3,1) = -0.57735026918962576d0
            obj%GaussPoint(3,2) = -0.57735026918962576d0
            obj%GaussPoint(3,3) = -0.57735026918962576d0
            obj%GaussPoint(3,4) = -0.57735026918962576d0
            obj%GaussPoint(3,5) =  0.57735026918962576d0
            obj%GaussPoint(3,6) =  0.57735026918962576d0
            obj%GaussPoint(3,7) =  0.57735026918962576d0
            obj%GaussPoint(3,8) =  0.57735026918962576d0
            obj%GaussIntegWei(:)=1.0d0
            
            if(obj%ReducedIntegration .eqv. .true.)then
                obj%GaussPoint(:,:) =0.0d0
                obj%GaussIntegWei(:)=0.1250d0
            endif

        elseif(size(obj%GaussPoint,2)==1 )then
            ! Triangular or Tetrahedral
            if(allocated(obj%GaussPoint))then
                deallocate(obj%GaussPoint)
            endif
            if(allocated(obj%GaussIntegWei))then
                deallocate(obj%GaussIntegWei)
            endif
            
            obj%GaussPoint(1,1) = 1.0d0
            obj%GaussPoint(2,1) = 1.0d0
            obj%GaussPoint(3,1) = 1.0d0
            obj%GaussIntegWei(:)= 1.0d0
        else
            print *, "ERROR :: ShapeFunctionClass/GetGaussPoint is not defined"
        endif


    end subroutine GetGaussPoint
!--------------------------------------------------------

!--------------------------------------------------------
    subroutine SetGaussPoint(obj)

        class(ShapeFunction_),intent(inout)::obj
        if(allocated(obj%gzi) ) then
            if(size(obj%gzi,1)/=obj%NumOfDim )then
            deallocate(obj%gzi)
                allocate(obj%gzi(obj%NumOfDim) )
            endif
        else
            allocate(obj%gzi(obj%NumOfDim) )
        endif
        
        
        if(obj%NumOfDim /= size(obj%GaussPoint,1) )then
            print *, "ERROR::SetGaussPoint",obj%NumOfDim, size(obj%GaussPoint,1)
            obj%ErrorMsg="ERROR::SetGaussPoint"
            obj%ierr=1
        else
            obj%gzi(:) = obj%GaussPoint(:,obj%GpID )
            
            obj%ErrorMsg="Succeed::SetGaussPoint"
            obj%ierr=0
        endif
        


    end subroutine SetGaussPoint
!--------------------------------------------------------


!--------------------------------------------------------
    subroutine GetShapeFunction(obj)
        class(ShapeFunction_),intent(inout)::obj

        
if(allocated(obj%Nmat) ) then
    if(size(obj%Nmat,1)/=obj%NumOfNode)then
        deallocate(obj%Nmat)
        allocate(obj%Nmat(obj%NumOfNode) )
    endif
else
    allocate(obj%Nmat(obj%NumOfNode) )
endif


if(obj%NumOfNode==1) then
    !#########################################################################
    !#######                                                           #######
    !#######                             +     (1)                     #######
    !#######                                                           #######
    !#######                                                           #######
    !#########################################################################
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction NumOfNode==1 "
    obj%ierr=1

elseif(obj%NumOfNode==2) then
    
    if(obj%NumOfDim==1) then
        
        !#########################################################################
        !#######                                                           #######
        !#######                +>----------------------->+                #######
        !#######            (1)                               (2)          #######
        !#######                                                           #######
        !#########################################################################
        
        obj%Nmat(1)=0.50d0*( 1.0d0-obj%gzi(1))
        obj%Nmat(2)=0.50d0*(-1.0d0+obj%gzi(1))
        
        obj%ErrorMsg="Succeed::GetShapeFunction "
        obj%ierr=0
    else
        print *, "ERROR::GetShapeFunction"
        obj%ErrorMsg="ERROR::GetShapeFunction NumOfNode==2 NumOfOrder/=1 "
        obj%ierr=1
    endif
elseif(obj%NumOfNode==3) then
    if(obj%NumOfDim==1) then
        !#########################################################################
        !#######                                                           #######
        !#######                +-----------+-------------+                #######
        !#######            (1)             (2)            (3)             #######
        !#######                                                           #######
        !#########################################################################
        !allocate(obj%Nmat(obj%NumOfNode,obj%NumOfDim) )
        !
        !obj%Nmat(1,1)=0.50d0*( 1.0d0-gzi(1))
        !obj%Nmat(1,2)=0.50d0*(-1.0d0+gzi(1))
        !obj%Nmat(1,3)=0.50d0*( 1.0d0-gzi(1))
        !
        !obj%ErrorMsg="Succeed::GetShapeFunction "
        !obj%ierr=0
        print *, "ERROR::GetShapeFunction"
        obj%ErrorMsg="ERROR::GetShapeFunction "
        obj%ierr=1
    elseif(obj%NumOfDim==2) then
        !#########################################################################
        !#######                                                           #######
        !#######          (1)   +-------------------------+                #######
        !#######                 \                       /  (3)            #######
        !#######                   \                   /                   #######
        !#######                     \               /                     #######
        !#######                       \           /                       #######
        !#######                         \       /                         #######
        !#######                           \   /                           #######        
        !#######                       (2)   +                             #######        
        !#########################################################################
        
        
        print *, "ERROR::GetShapeFunction"
        obj%ErrorMsg="ERROR::GetShapeFunction "
        obj%ierr=1
    else

        print *, "ERROR::GetShapeFunction"
        obj%ErrorMsg="ERROR::GetShapeFunction "
        obj%ierr=1
    endif
elseif(obj%NumOfNode==4) then
    if(obj%NumOfDim==1) then
        !#########################################################################
        !#######                                                           #######
        !#######                +-------+---------+-------+                #######
        !#######          (1)          (2)        (3)       (4)            #######
        !#######                                                           #######
        !#########################################################################
        !obj%ErrorMsg="Succeed::GetShapeFunction "
        !obj%ierr=0
        print *, "ERROR::GetShapeFunction"
        obj%ErrorMsg="ERROR::GetShapeFunction "
        obj%ierr=1
    elseif(obj%NumOfDim==2) then
        !#########################################################################
        !#######                                                           #######
        !#######           (1)  +-------------------------+  (4)           #######
        !#######                |                         |                #######
        !#######                |                         |                #######
        !#######                |                         |                #######
        !#######                |                         |                #######
        !#######                |                         |                #######
        !#######                |                         |                #######        
        !#######           (2)  +-------------------------+  (3)           #######        
        !#########################################################################
        
		obj%Nmat(1)=0.250d0*(1.0d0-obj%gzi(1))*(1.0d0-obj%gzi(2))
	    obj%Nmat(2)=0.250d0*(1.0d0+obj%gzi(1))*(1.0d0-obj%gzi(2))
		obj%Nmat(3)=0.250d0*(1.0d0+obj%gzi(1))*(1.0d0+obj%gzi(2))
		obj%Nmat(4)=0.250d0*(1.0d0-obj%gzi(1))*(1.0d0+obj%gzi(2))
	  
        obj%ErrorMsg="Succeed::GetShapeFunction "
        obj%ierr=0
    elseif(obj%NumOfDim==3) then
        !#########################################################################
        !#######                    (4)   +                                #######
        !#######                        /   \                              #######
        !#######                      /       \                            #######
        !#######                    /           \                          #######
        !#######                  /           ___+    (3)                  #######
        !#######                /  ___----'''     \                        #######
        !#######          (1)  +  -                 \                      #######
        !#######                '''-->---____       \                      #######        
        !#######                              -->--- +   (2)               #######        
        !#########################################################################
        
        print *, "ERROR::GetShapeFunction"
        obj%ErrorMsg="ERROR::GetShapeFunction "
        obj%ierr=1
    else
        print *, "ERROR::GetShapeFunction"
        obj%ErrorMsg="ERROR::GetShapeFunction "
        obj%ierr=1
    endif
    
    
elseif(obj%NumOfNode==5) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==6) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==7) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==8) then
    if(obj%NumOfDim==3) then
        !#########################################################################
        !#######           (8)   +-------------------------+ (7)           #######
        !#######                /|                        /|               #######
        !#######               / |                  (6)  / |               #######
        !#######          (5) +--|----------------------+  |               #######
        !#######              |  |                      |  |               #######
        !#######              |  +----------------------|--+ (3)           #######
        !#######              | / (4)                   | /                #######
        !#######              |/                        |/                 #######        
        !#######          (1) +-------------------------+ (2)              #######        
        !#########################################################################
        
        obj%Nmat(1)=0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 - obj%gzi(3))
	    obj%Nmat(2)=0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 - obj%gzi(3))
		obj%Nmat(3)=0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 - obj%gzi(3))
		obj%Nmat(4)=0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 - obj%gzi(3))
        obj%Nmat(5)=0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 + obj%gzi(3))
	    obj%Nmat(6)=0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 + obj%gzi(3))
		obj%Nmat(7)=0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 + obj%gzi(3))
		obj%Nmat(8)=0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 + obj%gzi(3))
	  

        obj%ErrorMsg="Succeed::GetShapeFunction "
        obj%ierr=0
    else
        print *, "ERROR::GetShapeFunction"
        obj%ErrorMsg="ERROR::GetShapeFunction "
        obj%ierr=1
    endif
    
elseif(obj%NumOfNode==9) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==10) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==11) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==12) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==13) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==14) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==15) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==16) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==17) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==18) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==19) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==20) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==21) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==22) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==23) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==24) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
else
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
endif

    end subroutine GetShapeFunction
!--------------------------------------------------------

!--------------------------------------------------------
    subroutine GetShapeFuncDer1(obj)
        class(ShapeFunction_),intent(inout)::obj

        

        if(allocated(obj%dNdgzi) ) then
            if(size(obj%dNdgzi,1)/=obj%NumOfDim .or. size(obj%dNdgzi,2)/=obj%NumOfNode)then
                deallocate(obj%dNdgzi)
                allocate(obj%dNdgzi(obj%NumOfDim,obj%NumOfNode) )
            endif
        else
            allocate(obj%dNdgzi(obj%NumOfDim,obj%NumOfNode) )
        endif
        
        
        
        if(obj%NumOfNode==1) then
            !#########################################################################
            !#######                                                           #######
            !#######                             +     (1)                     #######
            !#######                                                           #######
            !#######                                                           #######
            !#########################################################################
            print *, "ERROR::GetShapeFunction"
            obj%ErrorMsg="ERROR::GetShapeFunction NumOfNode==1 "
            obj%ierr=1
        
        elseif(obj%NumOfNode==2) then
            
            if(obj%NumOfDim==1) then
                
                !#########################################################################
                !#######                                                           #######
                !#######                +>----------------------->+                #######
                !#######            (1)                               (2)          #######
                !#######                                                           #######
                !#########################################################################
                
                obj%dNdgzi(1,1)= -0.50d0
                obj%dNdgzi(1,2)=  0.50d0
                
                obj%ErrorMsg="Succeed::GetShapeFunction "
                obj%ierr=0
            else
                print *, "ERROR::GetShapeFunction"
                obj%ErrorMsg="ERROR::GetShapeFunction NumOfNode==2 NumOfOrder/=1 "
                obj%ierr=1
            endif
        elseif(obj%NumOfNode==3) then
            if(obj%NumOfDim==1) then
                !#########################################################################
                !#######                                                           #######
                !#######                +-----------+-------------+                #######
                !#######            (1)             (2)            (3)             #######
                !#######                                                           #######
                !#########################################################################
                !allocate(obj%dNdgzi(obj%NumOfNode,obj%NumOfDim) )
                !
                !obj%dNdgzi(1,1)=0.50d0*( 1.0d0-gzi(1))
                !obj%dNdgzi(1,2)=0.50d0*(-1.0d0+gzi(1))
                !obj%dNdgzi(1,3)=0.50d0*( 1.0d0-gzi(1))
                !
                !obj%ErrorMsg="Succeed::GetShapeFunction "
                !obj%ierr=0
                print *, "ERROR::GetShapeFunction"
                obj%ErrorMsg="ERROR::GetShapeFunction "
                obj%ierr=1
            elseif(obj%NumOfDim==2) then
                !#########################################################################
                !#######                                                           #######
                !#######          (1)   +-------------------------+                #######
                !#######                 \                       /  (3)            #######
                !#######                   \                   /                   #######
                !#######                     \               /                     #######
                !#######                       \           /                       #######
                !#######                         \       /                         #######
                !#######                           \   /                           #######        
                !#######                       (2)   +                             #######        
                !#########################################################################
                
                
                print *, "ERROR::GetShapeFunction"
                obj%ErrorMsg="ERROR::GetShapeFunction "
                obj%ierr=1
            else
        
                print *, "ERROR::GetShapeFunction"
                obj%ErrorMsg="ERROR::GetShapeFunction "
                obj%ierr=1
            endif
        elseif(obj%NumOfNode==4) then
            if(obj%NumOfDim==1) then
                !#########################################################################
                !#######                                                           #######
                !#######                +-------+---------+-------+                #######
                !#######          (1)          (2)        (3)       (4)            #######
                !#######                                                           #######
                !#########################################################################
                !obj%ErrorMsg="Succeed::GetShapeFunction "
                !obj%ierr=0
                print *, "ERROR::GetShapeFunction"
                obj%ErrorMsg="ERROR::GetShapeFunction "
                obj%ierr=1
            elseif(obj%NumOfDim==2) then
                !#########################################################################
                !#######                                                           #######
                !#######           (4)  +-------------------------+  (3)           #######
                !#######                |                         |                #######
                !#######                |                         |                #######
                !#######                |                         |                #######
                !#######                |                         |                #######
                !#######                |                         |                #######
                !#######                |                         |                #######        
                !#######           (1)  +-------------------------+  (2)           #######        
                !#########################################################################
                
                obj%dNdgzi(1,1) = -0.250d0*(1.0d0-obj%gzi(2))
                obj%dNdgzi(1,2) =  0.250d0*(1.0d0-obj%gzi(2))
                obj%dNdgzi(1,3) =  0.250d0*(1.0d0+obj%gzi(2))
                obj%dNdgzi(1,4) = -0.250d0*(1.0d0+obj%gzi(2))
              
                obj%dNdgzi(2,1) = -0.250d0*(1.0d0-obj%gzi(1))
                obj%dNdgzi(2,2) = -0.250d0*(1.0d0+obj%gzi(1))
                obj%dNdgzi(2,3) =  0.250d0*(1.0d0+obj%gzi(1))
                obj%dNdgzi(2,4) =  0.250d0*(1.0d0-obj%gzi(1))
              
                obj%ErrorMsg="Succeed::GetShapeFunction "
                obj%ierr=0
            elseif(obj%NumOfDim==3) then
                !#########################################################################
                !#######                    (4)   +                                #######
                !#######                        /   \                              #######
                !#######                      /       \                            #######
                !#######                    /           \                          #######
                !#######                  /           ___+    (3)                  #######
                !#######                /  ___----'''     \                        #######
                !#######          (1)  +  -                 \                      #######
                !#######                '''-->---____       \                      #######        
                !#######                              -->--- +   (2)               #######        
                !#########################################################################
                
                print *, "ERROR::GetShapeFunction"
                obj%ErrorMsg="ERROR::GetShapeFunction "
                obj%ierr=1
            else
                print *, "ERROR::GetShapeFunction"
                obj%ErrorMsg="ERROR::GetShapeFunction "
                obj%ierr=1
            endif
            
            
        elseif(obj%NumOfNode==5) then
            print *, "ERROR::GetShapeFunction"
            obj%ErrorMsg="ERROR::GetShapeFunction "
            obj%ierr=1
        elseif(obj%NumOfNode==6) then
            print *, "ERROR::GetShapeFunction"
            obj%ErrorMsg="ERROR::GetShapeFunction "
            obj%ierr=1
        elseif(obj%NumOfNode==7) then
            print *, "ERROR::GetShapeFunction"
            obj%ErrorMsg="ERROR::GetShapeFunction "
            obj%ierr=1
        elseif(obj%NumOfNode==8) then
            if(obj%NumOfDim==3) then
                !#########################################################################
                !#######           (8)   +-------------------------+ (7)           #######
                !#######                /|                        /|               #######
                !#######               / |                  (6)  / |               #######
                !#######          (5) +--|----------------------+  |               #######
                !#######              |  |                      |  |               #######
                !#######              |  +----------------------|--+ (3)           #######
                !#######              | / (4)                   | /                #######
                !#######              |/                        |/                 #######        
                !#######          (1) +-------------------------+ (2)              #######        
                !#########################################################################
                
                obj%dNdgzi(1,1)= - 0.1250d0*(1.0d0 - obj%gzi(2))*(1.0d0 - obj%gzi(3))
                obj%dNdgzi(1,2)= + 0.1250d0*(1.0d0 - obj%gzi(2))*(1.0d0 - obj%gzi(3))
                obj%dNdgzi(1,3)= + 0.1250d0*(1.0d0 + obj%gzi(2))*(1.0d0 - obj%gzi(3))
                obj%dNdgzi(1,4)= - 0.1250d0*(1.0d0 + obj%gzi(2))*(1.0d0 - obj%gzi(3))
                obj%dNdgzi(1,5)= - 0.1250d0*(1.0d0 - obj%gzi(2))*(1.0d0 + obj%gzi(3))
                obj%dNdgzi(1,6)= + 0.1250d0*(1.0d0 - obj%gzi(2))*(1.0d0 + obj%gzi(3))
                obj%dNdgzi(1,7)= + 0.1250d0*(1.0d0 + obj%gzi(2))*(1.0d0 + obj%gzi(3))
                obj%dNdgzi(1,8)= - 0.1250d0*(1.0d0 + obj%gzi(2))*(1.0d0 + obj%gzi(3))
              
                obj%dNdgzi(2,1)= - 0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(3))
                obj%dNdgzi(2,2)= - 0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(3))
                obj%dNdgzi(2,3)= + 0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(3))
                obj%dNdgzi(2,4)= + 0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(3))
                obj%dNdgzi(2,5)= - 0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(3))
                obj%dNdgzi(2,6)= - 0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(3))
                obj%dNdgzi(2,7)= + 0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(3))
                obj%dNdgzi(2,8)= + 0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(3))
              
                obj%dNdgzi(3,1)= - 0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2))
                obj%dNdgzi(3,2)= - 0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2))
                obj%dNdgzi(3,3)= - 0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(2))
                obj%dNdgzi(3,4)= - 0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(2))
                obj%dNdgzi(3,5)= + 0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2))
                obj%dNdgzi(3,6)= + 0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2))
                obj%dNdgzi(3,7)= + 0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(2))
                obj%dNdgzi(3,8)= + 0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(2))
                
        
                obj%ErrorMsg="Succeed::GetShapeFunction "
                obj%ierr=0
            else
                print *, "ERROR::GetShapeFunction"
                obj%ErrorMsg="ERROR::GetShapeFunction "
                obj%ierr=1
            endif
            
        elseif(obj%NumOfNode==9) then
            print *, "ERROR::GetShapeFunction"
            obj%ErrorMsg="ERROR::GetShapeFunction "
            obj%ierr=1
        elseif(obj%NumOfNode==10) then
            print *, "ERROR::GetShapeFunction"
            obj%ErrorMsg="ERROR::GetShapeFunction "
            obj%ierr=1
        elseif(obj%NumOfNode==11) then
            print *, "ERROR::GetShapeFunction"
            obj%ErrorMsg="ERROR::GetShapeFunction "
            obj%ierr=1
        elseif(obj%NumOfNode==12) then
            print *, "ERROR::GetShapeFunction"
            obj%ErrorMsg="ERROR::GetShapeFunction "
            obj%ierr=1
        elseif(obj%NumOfNode==13) then
            print *, "ERROR::GetShapeFunction"
            obj%ErrorMsg="ERROR::GetShapeFunction "
            obj%ierr=1
        elseif(obj%NumOfNode==14) then
            print *, "ERROR::GetShapeFunction"
            obj%ErrorMsg="ERROR::GetShapeFunction "
            obj%ierr=1
        elseif(obj%NumOfNode==15) then
            print *, "ERROR::GetShapeFunction"
            obj%ErrorMsg="ERROR::GetShapeFunction "
            obj%ierr=1
        elseif(obj%NumOfNode==16) then
            print *, "ERROR::GetShapeFunction"
            obj%ErrorMsg="ERROR::GetShapeFunction "
            obj%ierr=1
        elseif(obj%NumOfNode==17) then
            print *, "ERROR::GetShapeFunction"
            obj%ErrorMsg="ERROR::GetShapeFunction "
            obj%ierr=1
        elseif(obj%NumOfNode==18) then
            print *, "ERROR::GetShapeFunction"
            obj%ErrorMsg="ERROR::GetShapeFunction "
            obj%ierr=1
        elseif(obj%NumOfNode==19) then
            print *, "ERROR::GetShapeFunction"
            obj%ErrorMsg="ERROR::GetShapeFunction "
            obj%ierr=1
        elseif(obj%NumOfNode==20) then
            print *, "ERROR::GetShapeFunction"
            obj%ErrorMsg="ERROR::GetShapeFunction "
            obj%ierr=1
        elseif(obj%NumOfNode==21) then
            print *, "ERROR::GetShapeFunction"
            obj%ErrorMsg="ERROR::GetShapeFunction "
            obj%ierr=1
        elseif(obj%NumOfNode==22) then
            print *, "ERROR::GetShapeFunction"
            obj%ErrorMsg="ERROR::GetShapeFunction "
            obj%ierr=1
        elseif(obj%NumOfNode==23) then
            print *, "ERROR::GetShapeFunction"
            obj%ErrorMsg="ERROR::GetShapeFunction "
            obj%ierr=1
        elseif(obj%NumOfNode==24) then
            print *, "ERROR::GetShapeFunction"
            obj%ErrorMsg="ERROR::GetShapeFunction "
            obj%ierr=1
        else
            print *, "ERROR::GetShapeFunction"
            obj%ErrorMsg="ERROR::GetShapeFunction "
            obj%ierr=1
        endif
    end subroutine GetShapeFuncDer1
!--------------------------------------------------------

!--------------------------------------------------------
    subroutine GetShapeFuncDer2(obj)
        class(ShapeFunction_),intent(inout)::obj

        
        include "./GetShapeFuncDer2.f90"

    end subroutine GetShapeFuncDer2
!--------------------------------------------------------



!--------------------------------------------------------
subroutine GetElemCoord(obj,nod_coord,elem_nod,elem_id)
    class(ShapeFunction_),intent(inout)::obj
    integer(int32),intent(in)::elem_nod(:,:),elem_id
    real(real64),intent(in)::nod_coord(:,:)
    integer(int32)::i,j,k,n,m

    
    n=size(elem_nod,2)
    m=size(nod_coord,2)
    
    if(allocated(obj%ElemCoord) )then
        if(n/=size(obj%ElemCoord,1) .or. m/=size(obj%ElemCoord,2)  )then
            deallocate(obj%ElemCoord)
            allocate(obj%ElemCoord(n,m)  )
        endif
    else
        allocate(obj%ElemCoord(n,m)  )
    endif
    do j=1,n
        obj%ElemCoord(j,1:m)=nod_coord(elem_nod(elem_id,j),1:m )
        !print *, obj%ElemCoord(j,1:m)
    enddo
    
    return


end subroutine GetElemCoord
!--------------------------------------------------------



!--------------------------------------------------------
subroutine GetElemCoord_n(obj,nod_coord_n,elem_nod,elem_id)
    class(ShapeFunction_),intent(inout)::obj
    integer(int32),intent(in)::elem_nod(:,:),elem_id
    real(real64),intent(in)::nod_coord_n(:,:)
    integer(int32)::i,j,k,n,m

    
    n=size(elem_nod,2)
    m=size(nod_coord_n,2)
    
    if(allocated(obj%ElemCoord_n) )then
        if(n/=size(obj%ElemCoord_n,1) .or. m/=size(obj%ElemCoord_n,2)  )then
            deallocate(obj%ElemCoord_n)
            allocate(obj%ElemCoord_n(n,m)  )
        endif
    else
        allocate(obj%ElemCoord_n(n,m)  )
    endif
    do j=1,n
        obj%ElemCoord_n(j,1:m)=nod_coord_n(elem_nod(elem_id,j),1:m )
        !print *, obj%ElemCoord_n(j,1:m)
    enddo
    
    return


end subroutine 
!--------------------------------------------------------



!--------------------------------------------------------
subroutine getdu(obj)
    class(ShapeFunction_),intent(inout)::obj
    integer(int32)::i,j,k,n,m

    n=size(obj%ElemCoord,1)
    m=size(obj%ElemCoord,2)
    
    if(allocated(obj%du) )then
        if(n/=size(obj%du,1) .or. m/=size(obj%du,2)  )then
            deallocate(obj%du)
            allocate(obj%du(n,m)  )
        endif
    else
        allocate(obj%du(n,m)  )
    endif
    do j=1,n
        obj%du(j,1:m)=obj%ElemCoord(j,1:m)-obj%ElemCoord_n(j,1:m)
    enddo
    
    return


end subroutine 
!--------------------------------------------------------


!--------------------------------------------------------
subroutine GetJmat(obj)
    class(ShapeFunction_),intent(inout)::obj
    integer(int32) n


    n=size(obj%ElemCoord,2)
    if(allocated(obj%Jmat) )then
        if(n/=size(obj%Jmat,1) .or. n/=size(obj%Jmat,2)  )then
            deallocate(obj%Jmat)
            allocate(obj%Jmat(n,n)  )
            allocate(obj%JmatInv(n,n) )
    
        endif
    else
        allocate(obj%Jmat(n,n)  )
        allocate(obj%JmatInv(n,n) )
    endif
    
    obj%Jmat(:,:)=matmul(obj%dNdgzi,obj%ElemCoord)
    
    call inverse_rank_2(obj%Jmat,obj%JmatInv)
    
end subroutine GetJmat
!--------------------------------------------------------

subroutine exportShapeFunction(obj,restart,path)
    class(ShapeFunction_),intent(inout) :: obj
    logical,optional,intent(in) :: restart
    character(*),intent(in) :: path
    type(IO_) :: f

    if(present(restart) )then
        call system("mkdir -p "//trim(path)//"/ShapeFunction")
        call f%open(trim(path)//"/ShapeFunction/","ShapeFunction",".res")
        write(f%fh,*) obj%Nmat(:)
        write(f%fh,*) obj%dNdgzi(:,:)
        write(f%fh,*) obj%dNdgzidgzi(:,:)
        write(f%fh,*) obj%gzi(:)
        write(f%fh,*) obj%gaussPoint(:,:)
        write(f%fh,*) obj%gaussIntegWei(:)
        write(f%fh,*) obj%Jmat(:,:)
        write(f%fh,*) obj%JmatInv(:,:)
        write(f%fh,*) obj%ElemCoord(:,:)
        write(f%fh,*) obj%ElemCoord_n(:,:)
        write(f%fh,*) obj%du(:,:)
        write(f%fh,*) obj%detJ
        write(f%fh,*) obj%NumOfNode
        write(f%fh,*) obj%NumOfOrder
        write(f%fh,*) obj%NumOfDim
        write(f%fh,*) obj%NumOfGp
        write(f%fh,*) obj%GpID
        write(f%fh,*) obj%ierr
        write(f%fh,*) obj%currentGpID
        write(f%fh,*) obj%ReducedIntegration
        write(f%fh,'(A)' ) trim(obj%elemType)
        write(f%fh,'(A)' ) trim(obj%ErrorMsg)
        call f%close()
    endif

end subroutine

end module ShapeFunctionClass